In [6]:
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 [7]:
# 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 [8]:
tcrop = 0.7
In [9]:
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 [10]:
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 [117]:
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 [13]:
X_bin = [epochs_bin_dev.get_data().transpose(0, 2, 1),
epochs_bin_std.get_data().transpose(0, 2, 1)]
X_ste = [epochs_ste_dev.get_data().transpose(0, 2, 1),
epochs_ste_std.get_data().transpose(0, 2, 1)]
Analyse de clusters
In [14]:
nperm = 1000
T_obs_bin,clusters_bin,clusters_pb_bin,H0_bin = mne.stats.spatio_temporal_cluster_test(X_bin,threshold=None,n_permutations=nperm,out_type='mask')
T_obs_ste,clusters_ste,clusters_pb_ste,H0_ste = mne.stats.spatio_temporal_cluster_test(X_ste,threshold=None,n_permutations=nperm,out_type='mask')
On récupère les channels trouvés grace a l'analyse de clusters
In [17]:
def extract_electrodes_times(clusters,clusters_pb,tmin_ind=500,tmax_ind=640,alpha=0.005,evoked = ev_bin_dev):
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)
channels_deviance_ste,times_deviance_ste=extract_electrodes_times(clusters_ste,clusters_pb_ste)
channels_deviance_bin,times_deviance_bin=extract_electrodes_times(clusters_bin,clusters_pb_bin)
print(channels_deviance_bin),print(times_deviance_bin)
print(channels_deviance_ste),print(times_deviance_ste)
times_union = np.union1d(times_deviance_bin,times_deviance_ste)
ch_union = np.unique(np.hstack([channels_deviance_bin,channels_deviance_ste]))
In [18]:
print(ch_union)
In [98]:
#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)]
X_diff_ste_bin = X_diff[1]-X_diff[0]
In [20]:
epochs_bin_dev_ch.plot_sensors(show_names=True)
plt.show()
In [146]:
roi = ['E117','E116','E108','E109','E151','E139','E141','E152','E110','E131','E143','E154','E142','E153','E140','E127','E118']
roi_frontal = ['E224','E223','E2','E4','E5','E6','E13','E14','E15','E20','E21','E27','E28','E30','E36','E40','E41']
In [147]:
len(roi_frontal),len(roi)
Out[147]:
In [70]:
from scipy.stats import ttest_1samp
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_bonferronia
In [102]:
def ttest_amplitude_roi(X,times_ind,ch_names_roi,times):
print(X.shape)
# Selecting time points and averaging over time
amps = X[:,times_ind,:].mean(axis=1)
# averaging over channels
amps = amps.mean(axis=1)
T, pval = ttest_1samp(amps, 0)
alpha = 0.05
n_samples, _, n_tests= X.shape
print('Uncorrected from %02f to %02f' % ((times[times_ind[0]]),times[times_ind[-1]]))
print("T = %0.2f, p = %0.3f " % (T,pval))
return T,pval,pval_fdr,pval_bonferroni
Tests de 280 a 440, par fenetres de 20 ms avec chevauchement de 10 ms
In [36]:
toi = np.arange(0.28,0.44,0.001)
toi_index = ev_bin_dev.time_as_index(toi)
In [38]:
wsize = 20
wstep = 10
In [44]:
toi
Out[44]:
Printing and preparing all time windows
In [75]:
all_toi_indexes = []
for i in range(14):
print(toi[10*i],toi[10*i + 20])
cur_toi_ind = range(10*i+1,(10*i+21))
all_toi_indexes.append(ev_bin_dev.time_as_index(toi[cur_toi_ind]))
print(toi[10*14],toi[10*14 + 19])
cur_toi_ind = range(10*14+1,(10*14+19))
all_toi_indexes.append(ev_bin_dev.time_as_index(toi[cur_toi_ind]))
Tests on each time window
In [76]:
for cur_timewindow in all_toi_indexes:
T,pval,pval_fdr,pval_bonferroni = ttest_amplitude(X_diff_ste_bin,cur_timewindow,epochs_bin_dev_ch.ch_names,times=epochs_bin_dev_ch.times)
On a channel subset (ROI) - average over channels
Parietal roi
In [161]:
#Selecting channels
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])
epochs_bin_dev_ch = epochs_bin_dev.pick_channels(roi)
epochs_bin_std_ch = epochs_bin_std.pick_channels(roi)
epochs_ste_dev_ch = epochs_ste_dev.pick_channels(roi)
epochs_ste_std_ch = epochs_ste_std.pick_channels(roi)
X_diff_roi = [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)]
X_diff_ste_bin_roi = X_diff_roi[1]-X_diff_roi[0]
for cur_timewindow in all_toi_indexes:
T,pval,pval_fdr,pval_bonferroni = ttest_amplitude_roi(X_diff_ste_bin_roi,cur_timewindow,roi,times=epochs_bin_dev_ch.times)
In [162]:
grav_bin_dev = epochs_bin_dev_ch.average()
grav_bin_std = epochs_bin_std_ch.average()
grav_ste_dev = epochs_ste_dev_ch.average()
grav_ste_std = epochs_ste_std_ch.average()
evoked_bin = mne.combine_evoked([grav_bin_dev, -grav_bin_std],
weights='equal')
evoked_ste = mne.combine_evoked([grav_ste_dev, -grav_ste_std],
weights='equal')
mne.viz.plot_compare_evokeds([grav_bin_std,grav_bin_dev,grav_ste_std,grav_ste_dev],picks=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])
plt.show()
mne.viz.plot_compare_evokeds([evoked_bin,evoked_ste],picks=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])
plt.show()
Frontal roi
In [163]:
#Selecting channels
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])
epochs_bin_dev_ch = epochs_bin_dev.pick_channels(roi_frontal)
epochs_bin_std_ch = epochs_bin_std.pick_channels(roi_frontal)
epochs_ste_dev_ch = epochs_ste_dev.pick_channels(roi_frontal)
epochs_ste_std_ch = epochs_ste_std.pick_channels(roi_frontal)
X_diff_roi = [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)]
X_diff_ste_bin_roi = X_diff_roi[1]-X_diff_roi[0]
for cur_timewindow in all_toi_indexes:
T,pval,pval_fdr,pval_bonferroni = ttest_amplitude_roi(X_diff_ste_bin_roi,cur_timewindow,roi,times=epochs_bin_dev_ch.times)
In [164]:
grav_bin_dev = epochs_bin_dev_ch.average()
grav_bin_std = epochs_bin_std_ch.average()
grav_ste_dev = epochs_ste_dev_ch.average()
grav_ste_std = epochs_ste_std_ch.average()
evoked_bin = mne.combine_evoked([grav_bin_dev, -grav_bin_std],
weights='equal')
evoked_ste = mne.combine_evoked([grav_ste_dev, -grav_ste_std],
weights='equal')
mne.viz.plot_compare_evokeds([grav_bin_std,grav_bin_dev,grav_ste_std,grav_ste_dev],picks=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])
plt.show()
mne.viz.plot_compare_evokeds([evoked_bin,evoked_ste],picks=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])
plt.show()
In [131]:
mne.viz.plot_compare_evokeds?
In [ ]:
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]:
a faire :