diff --git a/EEG4_Step04_ExploratoryERPPermutation.m b/EEG4_Step04_ExploratoryERPPermutation.m new file mode 100644 index 0000000..c51f48d --- /dev/null +++ b/EEG4_Step04_ExploratoryERPPermutation.m @@ -0,0 +1,171 @@ +%% B_analysis +% Perform lmeEEG on simulated EEG data +clear all;clc +load('E:\Exp02\EEGCode\Exp02_chanloc_.mat'); +load('E:\Exp02\EEGCode\Exp02_time_.mat'); +load('E:\Exp02\EEGCode\Exp02_sEEG_1.mat'); +DesignM = readtable('E:\Exp02\EEGCode\Exp02_DesignM.csv'); +DesignM.Properties.VariableNames = {'SubjID', 'Markers'}; +DesignM.Properties.VariableNames{'Markers'} = 'Marker3'; +unique_SubjID = unique(DesignM.SubjID); +Table = readtable('SemCatEStork.csv'); +Table.Properties.VariableNames{'x___Target'} = 'Target'; +Table2 = readtable('FreqDistra.csv'); +Table = innerjoin(Table, Table2,'Keys', {'Distractor'}); + +repeatedTable = repmat(Table, 36, 1); +SubjList = [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,... + 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, ... + 45, 46]; +repeatedTable.SubjID = reshape(repelem(SubjList, 102), 3672, 1); + +mergedTable = innerjoin(DesignM, repeatedTable,'Keys', {'SubjID', 'Marker3'}); + + + +Order = readtable('Exp02_SeqSI.csv'); +repeatedOrder = innerjoin(mergedTable, Order, 'Keys', {'SubjID', 'Marker3'}); + + +rowsToRemove = repeatedOrder.SubjID == 12 | repeatedOrder.SubjID == 44 |... + repeatedOrder.ExpTrialList == 1 | repeatedOrder.ExpTrialList == 2 |... + repeatedOrder.ExpTrialList == 3 | repeatedOrder.ExpTrialList == 4 |... + repeatedOrder.ExpTrialList == 74 | repeatedOrder.ExpTrialList == 75 |... + repeatedOrder.ExpTrialList == 76 | repeatedOrder.ExpTrialList == 77; + +repeatedOrder(rowsToRemove, :) = []; +sEEG = sEEG(:, :, ~rowsToRemove); + +channelinfo = EEG.chanlocs; + +clear DesignM EEG mergedTable Order repeatedTable rowsToRemove ... + SubjList Table unique_SubjID +save('well_preparedData.mat'); +disp('Data is well-organized!!!!!!!!!!'); +%% STEP 1 +% Conduct mixed models on each channel/timepoint combination. +ID = nominal(repeatedOrder.SubjID); Item=nominal(repeatedOrder.Target); +%CON = categorical(repeatedOrder.Condition); +JSD = categorical(repeatedOrder.JSD); +Classifier = categorical(repeatedOrder.ClassifierCongruency); +Freq = log(repeatedOrder.Frequency); +Stroke = (repeatedOrder.NumbersofStorks - mean(repeatedOrder.NumbersofStorks))/std(repeatedOrder.NumbersofStorks); +CongruencySemanticCategories = categorical(repeatedOrder.CongruencySemanticCategories); +mEEG = nan(size(sEEG)); + +for ch = 1:size(sEEG,1) + for tpoint = 1:size(sEEG,2) + EEG = double(squeeze(sEEG(ch,tpoint,:))); + EEG = table(EEG, CongruencySemanticCategories,Stroke, Freq, JSD, Classifier, ID, Item); + m = fitlme(EEG,'EEG~CongruencySemanticCategories + Stroke + Freq + JSD + Classifier + +(1|ID)+(1|Item)'); + + mEEG(ch,tpoint,:) = fitted(m,'Conditional',0)+residuals(m); + Coeff = fixedEffects(m); + Coeff_1 = Coeff(2:4); X = designMatrix(m); X_1 = X(:, 2:4); + coeff_CongruencySemanticCategories_1 = Coeff(strcmp(m.CoefficientNames, 'CongruencySemanticCategories_1')); + Coeff_Stroke = Coeff(strcmp(m.CoefficientNames, 'Stroke')); + coeff_Freq = Coeff(strcmp(m.CoefficientNames, 'Freq')); + Coeff_Intercept = Coeff(strcmp(m.CoefficientNames, '(Intercept)')); + +% m = fitlme(EEG,'EEG~Freq+Stroke+Animacy+Classifier+Animacy:Classifier+(1|ID)+(1|Item)'); + mEEG(ch,tpoint,:) = fitted(m,'Conditional',0)+residuals(m) - X_1 * Coeff_1; + + end +end + +% Extract design matrix X +EEG = double(squeeze(sEEG(1,1,:))); +% EEG = table(EEG,Freq,Stroke,Animacy,Classifier,ID,Item); +EEG = table(EEG,Classifier,ID,Item); +m = fitlme(EEG,'EEG~Classifier+(1|ID)+(1|Item)'); +% m = fitlme(EEG,'EEG~Freq+Stroke+Animacy+Classifier+Animacy:Classifier+(1|ID)+(1|Item)'); + +X = designMatrix(m); +clear EEG EEGall +save('./Output/marginalResult_Classifier.mat'); +clear EEG repeatedOrder tpoint; +disp('the marginalization stage is done!!!!!!!!!!'); +%% Step 2 +% Perform mass univariate linear regressions on ???marginal??? EEG data. +%f_obs = nan(size(mEEG,1),size(mEEG,2)); + +t_obs = nan(size(mEEG,1),size(mEEG,2),size(X,2)); +beta_obs = nan(size(mEEG,1),size(mEEG,2),size(X,2)); +se_obs = nan(size(mEEG,1),size(mEEG,2),size(X,2)); +for ch = 1:size(mEEG,1) + parfor tpoint = 1:size(mEEG,2) + EEG = squeeze(mEEG(ch,tpoint,:)); + [t_obs(ch,tpoint,:), beta_obs(ch,tpoint,:), se_obs(ch,tpoint,:)]=lmeEEG_regress(EEG,X); + mdl_summary = anova(fitlm(X, EEG), 'summary'); + f_obs(ch,tpoint, :)= mdl_summary.F(2); + end +end +save('./Output/observed_stage_classifier.mat'); +disp('the observation stage is done!!!!!!!!!!'); +%% Step 3 permutations +% Perform permutation testing ... +nperms=999; +num_rows = size(X,1); +% rperms = randperm(num_rows); + +% [rperms] = lmeEEG_permutations4(nperms, Animacy, Classifier, ID, Item); +f_perms = nan(nperms,size(mEEG,1),size(mEEG,2)); +t_perms = nan(nperms,size(mEEG,1),size(mEEG,2),size(X,2)); +beta_perms = nan(nperms,size(mEEG,1),size(mEEG,2),size(X,2)); +se_perms = nan(nperms,size(mEEG,1),size(mEEG,2),size(X,2)); +perm_t=nan(nperms,size(mEEG,1),size(mEEG,2)); +for p =1:nperms + XX = X(randperm(num_rows),:); +% XX = X(rperms(:,p),:); + for ch = 1:size(mEEG,1) + parfor tpoint = 1:size(mEEG,2) + EEG = squeeze(mEEG(ch,tpoint,:)); + tic + [t_perms(p,ch,tpoint,:), beta_perms(p,ch,tpoint,:),se_perms(p,ch,tpoint,:)]=lmeEEG_regress(EEG,XX); + perm_t(p,ch,tpoint)=toc; + mdl2_summary = anova(fitlm(XX, EEG),'summary'); + f_perms(p,ch,tpoint) = mdl2_summary.F(2); + end + end +end + +save('./Output/permutate_stage_classifier.mat'); + +disp('the permutation stage is done!!!!!!!!!!'); +%% ... and apply TFCE. +X = designMatrix(m); +for i = 2:size(X,2) + if ndims(t_obs) == 3 + Results.(matlab.lang.makeValidName(m.CoefficientNames{i})) = lmeEEG_TFCE(squeeze(t_obs(:,:,i)),squeeze(t_perms(:,:,:,i)),channelinfo,[0.66 2]); + elseif ndims(t_obs) == 4 + Results.(matlab.lang.makeValidName(m.CoefficientNames{i})) = lmeEEG_TFCE(squeeze(t_obs(:,:,:,i)),squeeze(t_perms(:,:,:,:,i)),channelinfo,[0.66 2]); + end +end + +Results2 = lmeEEG_TFCE(squeeze(f_obs(:,:)),squeeze(f_perms(:,:,:)),channelinfo,[0.5 1]); +save('./Output/tfce_stage_classifierr.mat'); +%% +mT = Results2.Obs; +mT2 = mT; +mT(not(Results2.Mask))=0; +tick_labels = reshape({channelinfo.labels}, 32, 1); + + +figure, +imagesc(mT) +xlim([0 230]) +set(gca,'ytick',1:32,'FontSize',15,'FontName','Arial'); +set(gca,'TickLength',[0 0]); +set(gca,'XTick',linspace(1,230,10),'XTickLabel',-200:100:700,'FontSize',15,'FontName','Arial'); +yticklabels(tick_labels); +hc=colorbar; +xlabel("Time (ms)","FontWeight","bold","FontSize",30, "FontName","Arial"); +ylabel(hc,'f-value','FontWeight','bold','FontSize',30,'FontName','Arial'); + +cmap2=cmap; cmap2(129,:)=[.8 .8 .8]; +set(gca,'clim',[-20 20],'colormap',cmap2) + +set(gca,'color','none') + +a = get(hc,'YTickLabel') +set(hc, 'colormap', cmap,'YTickLabel',a,'FontSize',8,'FontName','Arial'); diff --git a/EEG_Step01_PreproEEG_Batch.m b/EEG_Step01_PreproEEG_Batch.m new file mode 100644 index 0000000..94b7459 --- /dev/null +++ b/EEG_Step01_PreproEEG_Batch.m @@ -0,0 +1,23 @@ +clear all; clc; close all +SavePath1 = 'E:\Exp02\DATA\EEG_clean01\';%place of the data +FilePath='E:\Exp02\DATA\EEG\'; +mkdir(SavePath1); +sub_name = data_load(FilePath, 'bdf'); +%% STEP01 PREPROCESSING: get the dataset into the memory and run ica +for sub = 1:length(sub_name)% a loop for all data\ + filename = [FilePath,sub_name{sub}]; + EEG = pop_biosig(filename); + EEG = pop_chanedit(EEG, 'lookup','E:\\Exp02\\EEGCode\\standard_1020.elc'); +% EEG = pop_reref( EEG, [37 38] ); + EEG = pop_eegfiltnew(EEG, 0.1,30,333792,0,[],0); + EEG = pop_eegfiltnew(EEG, 48,52,31690,1,[],0); + EEG = pop_resample( EEG, 256); + EEG = pop_epoch( EEG, { '11' '12' '13' '14' '15' '16' '17' '18' '19' '20' '21' '22' '23' '24' '25' '26' '27' '28' '29' '30' '31' '32' '33' '34' '35' '36' '37' '38' '39' '40' '41' '42' '43' '44' '45' '46' '47' '48' '49' '50' '51' '52' '53' '54' '55' '56' '57' '58' '59' '60' '61' '62' '63' '64' '65' '66' '67' '68' '69' '70' '71' '72' '73' '74' '75' '76' '77' '78' '79' '80' '81' '82' '83' '84' '85' '86' '87' '88' '89' '90' '91' '92' '93' '94' '95' '96' '97' '98' '99' '100' '101' '102' '103' '104' '105' '106' '107' '108' '109' '110' '111' '112' }, [-0.2 0.7], 'newname', 'BDF file epochs', 'epochinfo', 'yes'); + + EEG = pop_rmbase( EEG, [-200 0]); + EEG = pop_select( EEG,'nochannel',{'EXG1' 'EXG2' 'EXG3' 'EXG4' 'EXG7' 'EXG8' 'GSR1' 'GSR2' 'Erg1' 'Erg2' 'Resp' 'Plet' 'Temp'}); + +% EEG = pop_runica(EEG, 'extended',1,'interupt','on'); + EEG = pop_saveset( EEG, 'filename',sub_name{sub},'filepath',SavePath1); + fprintf('*****the subject %d has been done',sub); +end diff --git a/EEG_Step02_PreproEEG.m b/EEG_Step02_PreproEEG.m new file mode 100644 index 0000000..73858cd --- /dev/null +++ b/EEG_Step02_PreproEEG.m @@ -0,0 +1,14 @@ +clear all; clc; close all; +FilePath = 'E:\Exp02\DATA\EEG_clean02\';%place of the data +SavePath1='E:\Exp02\DATA\EEG_clean03\'; +mkdir(SavePath1); +sub_name = data_load(FilePath, 'set'); +for sub = 1:length(sub_name)% a loop for all data\ + filename = [FilePath,sub_name{sub}]; + EEG = pop_loadset(filename); + EEG = pop_reref( EEG, [33 34] ); + EEG = pop_select( EEG,'nochannel',{'EXG5' 'EXG6'}); + EEG = pop_runica(EEG, 'extended',1,'interupt','on'); + EEG = pop_saveset( EEG, 'filename',sub_name{sub},'filepath',SavePath1); + fprintf('*****the subject %d has been done',sub); +end diff --git a/EEG_Step03_PreproRemovthreshold.m b/EEG_Step03_PreproRemovthreshold.m new file mode 100644 index 0000000..0b37963 --- /dev/null +++ b/EEG_Step03_PreproRemovthreshold.m @@ -0,0 +1,55 @@ +clear all; clc; close all; +FilePath = 'E:\Exp02\DATA\EEG_clean04\';%place of the data +SavePath1='E:\Exp02\DATA\EEG_clean05\'; +mkdir(SavePath1); +AllDATA = []; +DesignM = []; +sub_name = data_load(FilePath, 'set'); +n = 1; +for sub = 1:length(sub_name)% a loop for all data\ + filename = [FilePath,sub_name{sub}]; + EEG = pop_loadset(filename); + EEG = pop_rmbase( EEG, [-200 0]); + EEG = pop_eegthresh(EEG,1,[1:32] ,-100,100,-0.19922,0.7,0,1); + EEG = pop_epoch( EEG, { '11' '12' '13' '14' '15' '16' '17' '18' '19' '20' '21' '22' '23' '24' '25' '26' '27' '28' '29' '30' '31' '32' '33' '34' '35' '36' '37' '38' '39' '40' '41' '42' '43' '44' '45' '46' '47' '48' '49' '50' '51' '52' '53' '54' '55' '56' '57' '58' '59' '60' '61' '62' '63' '64' '65' '66' '67' '68' '69' '70' '71' '72' '73' '74' '75' '76' '77' '78' '79' '80' '81' '82' '83' '84' '85' '86' '87' '88' '89' '90' '91' '92' '93' '94' '95' '96' '97' '98' '99' '100' '101' '102' '103' '104' '105' '106' '107' '108' '109' '110' '111' '112' }, [-0.2 0.7], 'newname', 'BDF file epochs', 'epochinfo', 'yes'); + EEG = pop_saveset( EEG, 'filename',sub_name{sub},'filepath',SavePath1); + Markers = []; + for i = 1:length(EEG.event) + mar = EEG.event(i).type; + Markers = [Markers; mar]; + end + Markers(find(Markers>178)) = []; + Markers = unique(Markers); + + ElecName =[]; + for i = 1:length(EEG.chanlocs) + ele = EEG.chanlocs(i).labels; + ElecName = [ElecName; string(ele)]; + end + splitedCell = strsplit(sub_name{sub}, '.'); + splitedCell2 = strsplit(splitedCell{1}, '_'); + subjID = str2num(splitedCell2{2}); + Design = [repmat(subjID, length(Markers), 1), Markers]; + DesignM = [DesignM; Design]; + j = n + length(Markers)-1; + sEEG(:, :, n:j) = EEG.data; + n = j+1; + for trl = 1:size(EEG.data, 3) + Frame = [repmat(subjID, length(ElecName), 1), repmat(Markers(trl),length(ElecName), 1), ElecName, EEG.data(:, :, trl)]; + AllDATA =[AllDATA; Frame]; + end + + + fprintf('*****the subject %d has been done/n',subjID); +end +cell2csv('Exp02_EEG.csv', AllDATA, ','); +save('Exp02_EEG.mat', 'AllDATA', '-v7.3'); +time = EEG.times; +csvwrite('Exp02_EEGtime.csv', time) + +csvwrite('Exp02_DesignM.csv', DesignM) + +save('Exp_02DesignM_.mat', 'DesignM') +save('Exp_02sEEG_1.mat', 'sEEG') +save('Exp02_time_.mat', 'time') +save('Exp02_chanloc_.mat', 'EEG') diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..6bf17e3 --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2025 Yufang Wang + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. \ No newline at end of file diff --git a/data_load.m b/data_load.m new file mode 100644 index 0000000..6926978 --- /dev/null +++ b/data_load.m @@ -0,0 +1,38 @@ +function [ rest_eyes ] = data_load( filename,key_words ) +% 自动读取filename路径下的与key_words有关的文件 +% -- input -- +% filename - 想要读取的文件的路径 +% key_words - 想要读取的文件关键字 +% -- output -- +% 输出所要得到的文件的名字 + +rest = dir(filename); +num = 1; +num1 = 1; +while num + if(rest(num1).bytes == 0) + rest(num1) = []; + else + num = 0; + end +end +for i = 1:length(rest) + rest_name{i} = rest(i).name; +end +if ~isempty(key_words) + for i = 1:length(rest_name) + ind = strfind(rest_name{i},key_words); + if (ind ~= 0) + kk(i) = i; + end + end + kk = kk(find(kk~=0)); + rest_eyes = cell(1,length(kk)); + for i = 1:length(kk) + rest_eyes{i} = rest_name{kk(i)}; + end +else + rest_eyes = rest_name; +end + +