In [1]:
import scipy.io as sio
from matplotlib import pyplot as plt 
%matplotlib inline
import numpy as np 
import os 

import mne
import numpy as np
import scipy.io as sio

In [2]:
# coding: utf-8

def _loadftfile(path):

    filecontents = sio.whosmat(path)
    
    strucname = filecontents[0][0]

    mat = sio.loadmat(path, squeeze_me=True, struct_as_record=False)
    matstruct = mat[strucname]
    return matstruct 


def _matstruc2mne(matstruct,ch_names=None):
    
    if ch_names is None:
        ch_names=list(matstruct.label)
        
    myinfo = mne.create_info(ch_names=ch_names,sfreq=1/(matstruct.time[1] - matstruct.time[0]),ch_types='eeg')
    ev_arr = mne.EvokedArray(matstruct.individual.mean(axis=0),myinfo,tmin=-0.2) ### Specific to this dataset, 200ms baseline
    ev_arr.set_montage(mne.channels.read_montage("EGI_256"))
    return ev_arr


def _matstruc2latency(matstruct,peak_tmin,peak_tmax,ch_names=None):
    
    if ch_names is None:
        ch_names=list(matstruct.label)
        
    myinfo = mne.create_info(ch_names=ch_names,sfreq=1/(matstruct.time[1] - matstruct.time[0]),ch_types='eeg')
    
    all_chpeaks = []
    all_lat = []
    all_amp = []
    for mat in matstruct.individual:
    
        ev_arr = mne.EvokedArray(mat,myinfo,tmin=-0.2) ### Specific to this dataset, 200ms baseline
        ev_arr.set_montage(mne.channels.read_montage("EGI_256"))
        chpeak,lat,amp = ev_arr.get_peak(tmin=peak_tmin,tmax=peak_tmax,return_amplitude=True,mode='neg')
        all_lat.append(lat)
        all_chpeaks.append(chpeak)
        all_amp.append(amp)
    return all_chpeaks,all_lat,all_amp


def _matstruc2mne_epochs(matstruct,ch_names=None):
    
    if ch_names is None:
        ch_names=list(matstruct.label)
        
    myinfo = mne.create_info(ch_names=ch_names,sfreq=1/(matstruct.time[1] - matstruct.time[0]),ch_types='eeg')
    ev_arr = mne.EpochsArray(matstruct.individual,myinfo,tmin=-0.2) ### Specific to this dataset, 200ms baseline
    ev_arr.set_montage(mne.channels.read_montage("EGI_256"))
    return ev_arr

In [3]:
tcrop = 0.7

In [5]:
matfile = '/home/nfarrugi/datasets/eeg4sounds/result-eeg4sounds/oddball/grav/grav_bin_dev.mat'

mat_bin_dev = _loadftfile(matfile)

matfile = '/home/nfarrugi/datasets/eeg4sounds/result-eeg4sounds/oddball/grav/grav_bin_std.mat'

mat_bin_std = _loadftfile(matfile)



matfile = '/home/nfarrugi/datasets/eeg4sounds/result-eeg4sounds/oddball/grav/grav_ste_dev.mat'

mat_ste_dev = _loadftfile(matfile)

matfile = '/home/nfarrugi/datasets/eeg4sounds/result-eeg4sounds/oddball/grav/grav_ste_std.mat'

mat_ste_std = _loadftfile(matfile)

In [6]:
ev_bin_dev = _matstruc2mne(mat_bin_dev).crop(tmax=tcrop)
ev_bin_std = _matstruc2mne(mat_bin_std).crop(tmax=tcrop)

ev_ste_dev = _matstruc2mne(mat_ste_dev).crop(tmax=tcrop)
ev_ste_std = _matstruc2mne(mat_ste_std).crop(tmax=tcrop)

mne.equalize_channels([ev_bin_dev,ev_ste_dev,ev_bin_std,ev_ste_std])


Identifying common channels ...
Dropped the following channels:
['E73', 'E144', 'E9', 'E211', 'E201', 'E64', 'E203', 'E49', 'E63', 'E245', 'E158', 'E75', 'E251', 'E168', 'E227', 'E169', 'E160', 'E150', 'E183', 'E253', 'E59', 'E179', 'E35', 'E192', 'E12', 'E238', 'E213', 'E177', 'E256', 'E125', 'E29', 'E138', 'E34', 'E47', 'E248', 'E53', 'E231']

In [6]:
epochs_bin_dev = _matstruc2mne_epochs(mat_bin_dev).crop(tmax=tcrop)
epochs_bin_std = _matstruc2mne_epochs(mat_bin_std).crop(tmax=tcrop)

epochs_ste_dev = _matstruc2mne_epochs(mat_ste_dev).crop(tmax=tcrop)
epochs_ste_std = _matstruc2mne_epochs(mat_ste_std).crop(tmax=tcrop)

mne.equalize_channels([epochs_bin_dev,epochs_bin_std,epochs_ste_dev,epochs_ste_std])


19 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
tmax is not in epochs time interval. tmax is set to epochs.tmax
19 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
tmax is not in epochs time interval. tmax is set to epochs.tmax
19 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
tmax is not in epochs time interval. tmax is set to epochs.tmax
<ipython-input-6-43cba28dbd19>:1: RuntimeWarning: tmax is not in epochs time interval. tmax is set to epochs.tmax
  epochs_bin_dev = _matstruc2mne_epochs(mat_bin_dev).crop(tmax=tcrop)
<ipython-input-6-43cba28dbd19>:2: RuntimeWarning: tmax is not in epochs time interval. tmax is set to epochs.tmax
  epochs_bin_std = _matstruc2mne_epochs(mat_bin_std).crop(tmax=tcrop)
<ipython-input-6-43cba28dbd19>:4: RuntimeWarning: tmax is not in epochs time interval. tmax is set to epochs.tmax
  epochs_ste_dev = _matstruc2mne_epochs(mat_ste_dev).crop(tmax=tcrop)
19 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
tmax is not in epochs time interval. tmax is set to epochs.tmax
Identifying common channels ...
Dropped the following channels:
['E201', 'E245', 'E231', 'E211', 'E144', 'E73', 'E47', 'E251', 'E34', 'E192', 'E59', 'E213', 'E160', 'E138', 'E248', 'E12', 'E256', 'E168', 'E169', 'E35', 'E150', 'E9', 'E64', 'E49', 'E53', 'E158', 'E29', 'E227', 'E63', 'E183', 'E177', 'E125', 'E75', 'E238', 'E253', 'E179', 'E203']
<ipython-input-6-43cba28dbd19>:5: RuntimeWarning: tmax is not in epochs time interval. tmax is set to epochs.tmax
  epochs_ste_std = _matstruc2mne_epochs(mat_ste_std).crop(tmax=tcrop)

Latency analysis


In [7]:
peak_tmin = 0.30
peak_tmax = 0.35
ch_lat_bin_dev,lat_bin_dev,amp_bin_dev = _matstruc2latency(mat_bin_dev,peak_tmin=peak_tmin,peak_tmax=peak_tmax)
#ev_bin_std = _matstruc2latency(mat_bin_std)

ch_lat_ste_dev,lat_ste_dev,amp_ste_dev = _matstruc2latency(mat_ste_dev,peak_tmin=peak_tmin,peak_tmax=peak_tmax)
#ev_ste_std = _matstruc2latency(mat_ste_std)

In [9]:
plt.hist(lat_bin_dev,alpha=0.5)
plt.hist(lat_ste_dev,alpha=0.5)
plt.legend(['binaural','stereo'])


Out[9]:
<matplotlib.legend.Legend at 0x7f65602a7198>

In [10]:
plt.hist(amp_bin_dev,alpha=0.5)
plt.hist(amp_ste_dev,alpha=0.5)
plt.legend(['binaural','stereo'])


Out[10]:
<matplotlib.legend.Legend at 0x7f65623cdf28>

In [11]:
plt.boxplot([lat_bin_dev,lat_ste_dev])


Out[11]:
{'whiskers': [<matplotlib.lines.Line2D at 0x7f65601a5a58>,
  <matplotlib.lines.Line2D at 0x7f65601a5ef0>,
  <matplotlib.lines.Line2D at 0x7f65601b7860>,
  <matplotlib.lines.Line2D at 0x7f65601b7c88>],
 'caps': [<matplotlib.lines.Line2D at 0x7f65601ad320>,
  <matplotlib.lines.Line2D at 0x7f65601ad748>,
  <matplotlib.lines.Line2D at 0x7f65601bf0f0>,
  <matplotlib.lines.Line2D at 0x7f65601bf518>],
 'boxes': [<matplotlib.lines.Line2D at 0x7f65601a5908>,
  <matplotlib.lines.Line2D at 0x7f65601b7438>],
 'medians': [<matplotlib.lines.Line2D at 0x7f65601adb70>,
  <matplotlib.lines.Line2D at 0x7f65601bf940>],
 'fliers': [<matplotlib.lines.Line2D at 0x7f65601adf98>,
  <matplotlib.lines.Line2D at 0x7f65601bfd68>],
 'means': []}

Cluster tests


In [7]:
T_obs_bin,clusters_bin,clusters_pb_bin,H0_bin = mne.stats.spatio_temporal_cluster_test(X_bin,threshold=None,n_permutations=1000,out_type='mask')


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-24b28e85018e> in <module>()
----> 1 T_obs_bin,clusters_bin,clusters_pb_bin,H0_bin = mne.stats.spatio_temporal_cluster_test(X_bin,threshold=None,n_permutations=1000,out_type='mask')

NameError: name 'X_bin' is not defined

In [91]:
T_obs_ste,clusters_ste,clusters_pb_ste,H0_ste = mne.stats.spatio_temporal_cluster_test(X_ste,threshold=None,n_permutations=1000,out_type='mask')


Using a threshold of 21.093024
stat_fun(H1): min=0.000000 max=31.585274
Running initial clustering
Found 18 clusters
Permuting 999 times...
[....................................... ] 99.30%  \   Computing cluster p-values
Done.

In [92]:
T_obs_dev,clusters_dev,clusters_pb_dev,H0_dev = mne.stats.spatio_temporal_cluster_test(X_dev,threshold=12,n_permutations=1000,out_type='mask')


stat_fun(H1): min=0.000000 max=14.638737
Running initial clustering
Found 3 clusters
Permuting 999 times...
[....................................... ] 99.30%  \   Computing cluster p-values
Done.

In [93]:
T_obs_stan,clusters_stan,clusters_pb_stan,H0_stan = mne.stats.spatio_temporal_cluster_test(X_stan,threshold=None,n_permutations=1000,out_type='mask')


Using a threshold of 21.093024
stat_fun(H1): min=0.000000 max=9.579573
Running initial clustering
Found 0 clusters
No clusters found, returning empty H0, clusters, and cluster_pv
<ipython-input-93-ddeda34da736>:1: RuntimeWarning: No clusters found, returning empty H0, clusters, and cluster_pv
  T_obs_stan,clusters_stan,clusters_pb_stan,H0_stan = mne.stats.spatio_temporal_cluster_test(X_stan,threshold=None,n_permutations=1000,out_type='mask')

In [109]:
evoked = mne.combine_evoked([ev_bin_dev, -ev_ste_dev],
                            weights='equal')  # calculate difference wave
time_unit = dict(time_unit="ms")
evoked.plot_joint(title="Deviants Binaural vs. Stereo")#, ts_args=time_unit,topomap_args=time_unit)  # show difference wave
plt.show()

for clust,pval in zip(clusters_dev,clusters_pb_dev):
    if pval < 0.1:
        evoked.plot_image(mask=clust.T,mask_style='contour')
        plt.show()


time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-109-7b40f5617e24>:4: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_joint(title="Deviants Binaural vs. Stereo")#, ts_args=time_unit,topomap_args=time_unit)  # show difference wave
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-109-7b40f5617e24>:4: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_joint(title="Deviants Binaural vs. Stereo")#, ts_args=time_unit,topomap_args=time_unit)  # show difference wave

In [106]:
evoked = mne.combine_evoked([ev_bin_dev, -ev_bin_std],
                            weights='equal')  # calculate difference wave
time_unit = dict(time_unit="ms")
evoked.plot_joint(title="Binaural Deviant vs. Std")#, ts_args=time_unit,topomap_args=time_unit)  # show difference wave
plt.show()

for clust,pval in zip(clusters_bin,clusters_pb_bin):
    if pval < 0.005:
        evoked.plot_image(mask=clust.T,mask_style='contour')
        plt.show()
        clustersize = np.sum(clust)
        print(clustersize)


time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:4: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_joint(title="Binaural Deviant vs. Std")#, ts_args=time_unit,topomap_args=time_unit)  # show difference wave
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:4: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_joint(title="Binaural Deviant vs. Std")#, ts_args=time_unit,topomap_args=time_unit)  # show difference wave
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
51
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
33
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
137
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
56
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
74
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
187
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
84
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
45
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
82
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
29
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')
333
time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
<ipython-input-106-f3698ca74301>:9: DeprecationWarning: time_unit defaults to "ms" in 0.16 but will change to "s" in 0.17, set it explicitly to avoid this warning
  evoked.plot_image(mask=clust.T,mask_style='contour')