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])
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])
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]:
In [10]:
plt.hist(amp_bin_dev,alpha=0.5)
plt.hist(amp_ste_dev,alpha=0.5)
plt.legend(['binaural','stereo'])
Out[10]:
In [11]:
plt.boxplot([lat_bin_dev,lat_ste_dev])
Out[11]:
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')
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')
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')
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')
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()
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)
In [108]:
evoked = mne.combine_evoked([ev_ste_dev, -ev_ste_std],
weights='equal') # calculate difference wave
evoked.plot_joint(title="Stereo Deviant vs. Std")#, ts_args=time_unit,topomap_args=time_unit) # show difference wave
plt.show()
for clust,pval in zip(clusters_ste,clusters_pb_ste):
if pval < 0.005:
evoked.plot_image(mask=clust.T,mask_style='contour')
plt.show()
clustersize = np.sum(clust)
print(clustersize)
extraire les electrodes (max) et timepoints across tous les clusters. On se restreint à 300-400ms.
In [295]:
def extract_electrodes_times(clusters,clusters_pb,tmin_ind=500,tmax_ind=640,alpha=0.005):
ch_list_temp = []
time_list_temp = []
for clust,pval in zip(clusters,clusters_pb):
if pval < alpha:
for j,curline in enumerate(clust[tmin_ind:tmax_ind]):
for k,el in enumerate(curline):
if el:
ch_list_temp.append(evoked.ch_names[k])
time_list_temp.append(evoked.times[j+tmin_ind])
return np.unique(ch_list_temp),np.unique(time_list_temp)
In [296]:
channels_deviance_ste,times_deviance_ste=extract_electrodes_times(clusters_ste,clusters_pb_ste)
In [297]:
channels_deviance_bin,times_deviance_bin=extract_electrodes_times(clusters_bin,clusters_pb_bin)
In [300]:
print(channels_deviance_bin),print(times_deviance_bin)
Out[300]:
In [301]:
print(channels_deviance_ste),print(times_deviance_ste)
Out[301]:
In [314]:
np.intersect1d(times_deviance_bin,times_deviance_ste)
Out[314]:
In [316]:
times_union = np.union1d(times_deviance_bin,times_deviance_ste)
In [324]:
ch_union = np.unique(np.hstack([channels_deviance_bin,channels_deviance_ste]))
In [326]:
plt.plot(times_union,np.ones_like(times_union),'*')
Out[326]:
In [346]:
plt.plot(times_deviance_bin,np.ones_like(times_deviance_bin),'*')
Out[346]:
In [347]:
plt.plot(times_deviance_ste,np.ones_like(times_deviance_ste),'*')
Out[347]:
Isolons les trois périodes
In [343]:
t1 = times_union[:32]
print(t1)
In [344]:
t2 = times_union[32:70]
print(t2)
In [345]:
t3 = times_union[71:]
print(t3)
In [362]:
evoked_diff_ste = mne.combine_evoked([ev_ste_dev, -ev_ste_std],weights='equal').crop(t1[0],t1[-1])
evoked_diff_ste.plot_topomap(time_unit='ms')
plt.show()
evoked_diff_ste = mne.combine_evoked([ev_ste_dev, -ev_ste_std],weights='equal').crop(t2[0],t2[-1])
evoked_diff_ste.plot_topomap(time_unit='ms')
plt.show()
evoked_diff_ste = mne.combine_evoked([ev_ste_dev, -ev_ste_std],weights='equal').crop(t3[0],t3[-1])
evoked_diff_ste.plot_topomap(time_unit='ms')
plt.show()
In [363]:
evoked_diff_bin = mne.combine_evoked([ev_bin_dev, -ev_bin_std],weights='equal').crop(t1[0],t1[-1])
evoked_diff_bin.plot_topomap(time_unit='ms')
plt.show()
evoked_diff_bin = mne.combine_evoked([ev_bin_dev, -ev_bin_std],weights='equal').crop(t2[0],t2[-1])
evoked_diff_bin.plot_topomap(time_unit='ms')
plt.show()
evoked_diff_bin = mne.combine_evoked([ev_bin_dev, -ev_bin_std],weights='equal').crop(t3[0],t3[-1])
evoked_diff_bin.plot_topomap(time_unit='ms')
plt.show()
In [449]:
epochs_bin_dev.plot_sensors(show_names=True)
plt.show()
In [12]:
#Selecting channels
epochs_bin_dev_ch = epochs_bin_dev.pick_channels(ch_union)
epochs_bin_std_ch = epochs_bin_std.pick_channels(ch_union)
epochs_ste_dev_ch = epochs_ste_dev.pick_channels(ch_union)
epochs_ste_std_ch = epochs_ste_std.pick_channels(ch_union)
X_diff = [epochs_bin_dev_ch.get_data().transpose(0, 2, 1) - epochs_bin_std_ch.get_data().transpose(0, 2, 1),
epochs_ste_dev_ch.get_data().transpose(0, 2, 1) - epochs_ste_std_ch.get_data().transpose(0, 2, 1)]
t1_ind = evoked.time_as_index(t1)
t2_ind = evoked.time_as_index(t2)
t3_ind = evoked.time_as_index(t3)
X_diff_ste_bin = X_diff[1]-X_diff[0]
In [9]:
from scipy import stats
from mne.stats import bonferroni_correction,fdr_correction
def ttest_amplitude(X,times_ind,ch_names,times):
# Selecting time points and averaging over time
amps = X[:,times_ind,:].mean(axis=1)
T, pval = ttest_1samp(amps, 0)
alpha = 0.05
n_samples, n_tests= amps.shape
threshold_uncorrected = stats.t.ppf(1.0 - alpha, n_samples - 1)
reject_bonferroni, pval_bonferroni = bonferroni_correction(pval, alpha=alpha)
threshold_bonferroni = stats.t.ppf(1.0 - alpha / n_tests, n_samples - 1)
reject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method='indep')
mask_fdr = pval_fdr < 0.05
mask_bonf = pval_bonferroni < 0.05
print('FDR from %02f to %02f' % ((times[times_ind[0]]),times[times_ind[-1]]))
for i,curi in enumerate(mask_fdr):
if curi:
print("Channel %s, T = %0.2f, p = %0.3f " % (ch_names[i], T[i],pval_fdr[i]))
print('Bonferonni from %02f to %02f' % ((times[times_ind[0]]),times[times_ind[-1]]))
for i,curi in enumerate(mask_bonf):
if curi:
print("Channel %s, T = %0.2f, p = %0.3f " % (ch_names[i], T[i],pval_bonferroni[i]))
return T,pval,pval_fdr,pval_bonferroni
In [463]:
T,pval,pval_fdr,pval_bonferroni = ttest_amplitude(X_diff_ste_bin,t1_ind,epochs_bin_dev_ch.ch_names,times=epochs_bin_dev_ch.times)
In [464]:
T,pval,pval_fdr,pval_bonferroni = ttest_amplitude(X_diff_ste_bin,t2_ind,epochs_bin_dev_ch.ch_names,times=epochs_bin_dev_ch.times)
In [465]:
T,pval,pval_fdr,pval_bonferroni = ttest_amplitude(X_diff_ste_bin,t3_ind,epochs_bin_dev_ch.ch_names,times=epochs_bin_dev_ch.times)
In [466]:
from scipy import stats
from mne.stats import bonferroni_correction,fdr_correction
T, pval = ttest_1samp(X_diff_ste_bin, 0)
alpha = 0.05
n_samples, n_tests,_ = X_diff_ste_bin.shape
threshold_uncorrected = stats.t.ppf(1.0 - alpha, n_samples - 1)
reject_bonferroni, pval_bonferroni = bonferroni_correction(pval, alpha=alpha)
threshold_bonferroni = stats.t.ppf(1.0 - alpha / n_tests, n_samples - 1)
reject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method='indep')
#threshold_fdr = np.min(np.abs(T)[reject_fdr])
In [467]:
masking_mat = pval<0.05
Tbis = np.zeros_like(T)
Tbis[masking_mat] = T[masking_mat]
plt.matshow(Tbis.T,cmap=plt.cm.RdBu_r)
plt.colorbar()
plt.show()
plt.matshow(-np.log10(pval).T)
plt.colorbar()
Out[467]:
In [ ]: