MEG:ANOVA dans Fieldtrip
From WikiMEG
Contents
Analyse non paramétrique de données MEG à deux facteurs (groupes et conditions)
Analyse non paramétrique de données MEG avec deux facteurs (groupes et % conditions) à deux niveaux chacun. Test par permutation basé sur des clusters temporels sur les données regroupées en régions d'intérêt.
Scripts présentés lors du Club MEG du Mardi 26 Janvier 2016.
Le script principal : ins_temporalClusterPermuTest_WithinBetween.m
% ins_temporalClusterPermuTest_WithinBetween.m function result = ins_temporalClusterPermuTest_WithinBetween(data, def, cfg, structTemplate) % http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%% % % Nonparametric MEG data analysis with two factors (groups and conditions) % on two levels each. Permutation test based on temporal clusters on the % data grouped into regions of interest. % [Analyse non paramétrique de données MEG avec deux facteurs (groupes et % conditions) à deux niveaux chacun. Test par permutation basé sur des % clusters temporels sur les données regroupées en régions d'intérêt.] % % DEPENDANCES: % This function uses the FieldTrip MATLAB software toolbox % (http://www.fieldtriptoolbox.org/) % % USAGE: % result = ins_temporalClusterPermuTest_WithinBetween( % data, % def, % cfg, % structTemplate) % % INPUTS: % data = array <1xN struct> % where N = person x roi x cond x other criterium (side) % example of struct : person = 'sujet01' % group = 'control' % condition = 'condition_A' % roi = 'roi_3' % side = 'left' % time = <1x193 double> % avgabs = <1x193 double> % def = structure arrays with correspondance between waited % field names and data field names, and list of levels % for each categorical field. Waited field names are : % subject, group, cond, roi, other (other selection criterium), % time, value. % cfg = configuration parameters (cf. example). % structTemplate = template to reconstruct the data structure waited by % the FieldTrip function ft_timelockgrandaverage. % (cf. http://www.fieldtriptoolbox.org/reference/ft_timelockgrandaverage) % % OUTPUTS: % result = %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % REFERENCES: % * Cluster-based permutation tests on event related fields. % http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock % ________________________________ % Bernard Giusiano & Sophie Chen % INSERM UMR 1106 Institut de Neurosciences des Systèmes % Oct/2015 (first version) % Oct/2015 (this version) % http://ins.univ-amu.fr %% global parameters cfg.method = 'montecarlo' ; % use the Monte Carlo Method to calculate the significance probability cfg.correctm = 'cluster'; % Apply multiple-comparison correction % Interaction between groups and conditions if def.stats{1} disp('*** 1 - Interaction between groups and conditions ***'); resultat = ins_temporalClusterPermuTest_interaction(data, def, cfg, structTemplate); disp(['=== ' resultat ' ===']); end % Main effect of conditions (within) if def.stats{2} disp('*** 2 - Main effect of conditions (within) ***'); resultat = ins_temporalClusterPermuTest_within(data, def, cfg, structTemplate); disp(['=== ' resultat ' ===']); end % Main effect of groups (between) if def.stats{3} disp('*** 3 - Main effect of groups (between) ***'); resultat = ins_temporalClusterPermuTest_between(data, def, cfg, structTemplate); disp(['=== ' resultat ' ===']); end % Simple effect of conditions by group if def.stats{4} disp('*** 4 - Simple effect of conditions by group ***'); resultat = ins_temporalClusterPermuTest_within_simple(data, def, cfg, structTemplate); disp(['=== ' resultat ' ===']); end % Simple effect of groups by condition if def.stats{5} disp('*** 5 - Simple effect of groups by condition ***'); resultat = ins_temporalClusterPermuTest_between_simple(data, def, cfg, structTemplate); disp(['=== ' resultat ' ===']); end result = 'END of script ins_temporalClusterPermuTest_WithinBetween.m';
Exemple d'appel du script ci-dessus : RunningExample_TCPTWB.m
% RunningExample_TCPTWB.m % Example of utilization of the function % ins_temporalClusterPermuTest_WithinBetween.m % Bernard Giusiano & Sophie Chen - oct 2015 %% directory of scripts and example data : cd('D:\_Recherche\DEV_temporalClusterPermuTest_WithinBetween\scripts ins_TCPTWB'); % cd('C:\Users\dr00245\Dropbox\_RECHERCHE\INS_MEG_MCI3\clustersTemporels'); load('structurePorteuseDataMEG.mat'); % -> variable 'structurePorteuseDataMEG' (-> paramètre structTemplate) % data loading load('mydata.mat'); % makink directory for output currentFolder = pwd; outputFolder = sprintf('%s\\output', currentFolder); %%% maj dans les fonctions %%%% if ~exist(outputFolder, 'dir') mkdir(outputFolder); end %% definition of the experiment mydef.subject = 'person'; mydef.list.subject = {'01','02','03','04','05','06','07','11',... '12','13','14','15','16','17','18','19'}; mydef.group = 'group'; mydef.list.group = {'subject','subject','subject','subject','subject','subject','subject','subject',... 'control','control','control','control','control','control','control','control'}; mydef.cond = 'condition'; mydef.list.cond = {'condition_A','condition_B'}; mydef.roi = 'roi'; % mydef.list.roi = {'roi_1','roi_2','roi_3','roi_4','roi_5','roi_6','roi_7'}; mydef.list.roi = {'roi_1','roi_5','roi_6','roi_7'}; % select 4 roi only mydef.other = 'side'; % mydef.list.other = {'L','R'}; mydef.list.other = {'L'}; % select left side only mydef.time = 'time'; mydef.value = 'avgabs'; mydef.stats = {... % select statistical calculations to be made (true) or not (false): true,... % Interaction between groups and conditions false,... % Main effect of conditions (within) true,... % Main effect of groups (between) true,... % Simple effect of conditions by group false,... % Simple effect of groups by condition }; %% global parameters % for other parameters values cf. http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock cfg = []; cfg.latency = [-0.3 0.9]; % time interval over which the experimental % conditions must be compared (in seconds) or 'all' (default = 'all') cfg.numrandomization = 10000; % Number of draws from the permutation distribution cfg.clusteralpha = 0.05; % alpha level of the sample-specific statistic that will be used for thresholding % the values used to construct clusters cfg.tail = 0; % -1, 1 or 0 (default = 0); one-sided or two-sided test cfg.clusterstatistic = 'wcm'; % How to combine the single samples that belong to a cluster % -> test statistic that will be evaluated under the permutation distribution % cfg.clusterstatistic = 'maxsum', 'maxsize', 'wcm' (default = 'maxsum') % option 'wcm' refers to 'weighted cluster mass', % a statistic that combines cluster size and % intensity; see Hayasaka & Nichols (2004) NeuroImage % for details % When cfg.clusterstatistic = 'wcm': cfg.wcm_weight = 0.5; % cf. Hayasaka & Nichols (2004) NeuroImage (page 62) cfg.clustertail = 0; % -1, 1 or 0 (default = 0 : two-tail test) cfg.correcttail = 'alpha'; % correct p-values or alpha-values when doing a two-sided test % cf. http://www.fieldtriptoolbox.org/faq/why_should_i_use_the_cfg.correcttail_option_when_using_statistics_montecarlo cfg.alpha = 0.05; % 0.05/10; % alpha level of the permutation test / 10 tests % ATTENTION: % alpha SHOULD BE CORRECTED based on the number of tests % example of correction for 7 roi, 2 groups, 2 conditions, 2 sides: % - interactions: 7 roi * 2 sides = 14 tests % - main effects: 2 factors (groups and conditions) * 7 roi * 2 sides = 28 tests % - simple effects: 2 groups * 2 cond * 7 roi * 2 sides = 56 tests % => 14 + 28 + 56 = 98 tests => Bonferroni's correction: 0.05/98 = 0.0005 % => cfg.alpha = 0.0005; for a global alpha threshold really equal to 0.05 % Other parameters: % cfg.clusterthreshold = 'nonparametric_common'; % method for single-sample threshold cfg.avgoverchan = 'yes'; % average over channels => neighbours no needed cfg.neighbours = {}; % cf. http://mailman.science.ru.nl/pipermail/fieldtrip/2013-July/006853.html %% trace %%%%%%%%%%%%%%%% % traces are saved during execution of script ins_temporalClusterPermuTest_graph.m log_filename = ['./output/logProbaMin_' datestr(now, 'yyyymmdd_HHMM') '.txt']; cfg.fid = fopen(log_filename, 'wt'); %%%%%%%%%%%%%%%%%%%%%%%%% %% resultat = ins_temporalClusterPermuTest_WithinBetween(mydata, mydef, cfg, structurePorteuseDataMEG); disp(['### ' resultat ' ###']); disp('*** END of script RunningExample_temporalClusterPermuTest_WithinBetween.m ***'); %% trace %%%%%%%%%%%%%%%% fclose(cfg.fid); %%%%%%%%%%%%%%%%%%%%%%%%%
Les autres scripts : fonctions appelées par le script principal
% ins_temporalClusterPermuTest_interaction.m function result = ins_temporalClusterPermuTest_interaction(data, def, cfg, structTemplate) % http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%% % % Nonparametric MEG data analysis of interaction between two factors (groups and conditions) % on two levels each. Permutation test based on temporal clusters on the % data grouped into regions of interest. % % DEPENDANCES: % This function uses the FieldTrip MATLAB software toolbox % (http://www.fieldtriptoolbox.org/) % % USAGE: % result = ins_temporalClusterPermuTest_interaction( % data, % def, % cfg, % structTemplate) % % INPUTS: % cf. function ins_temporalClusterPermuTest_WithinBetween.m % % OUTPUTS: % result = %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % REFERENCES: % * Cluster-based permutation tests on event related fields. % http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock % * Eric Maris. How to test an interaction effect using cluster-based permutation tests? % http://www.fieldtriptoolbox.org/faq/how_can_i_test_an_interaction_effect_using_cluster-based_permutation_tests % ________________________________ % Bernard Giusiano & Sophie Chen % INSERM UMR 1106 Institut de Neurosciences des Systèmes % Sept/2015 (first version) % Oct/2015 (this version) % http://ins.univ-amu.fr %% dimensions and grand average parameters groupDesc = tabulate(def.list.group); nGroupLevel1 = groupDesc{1,2}; nGroupLevel2 = groupDesc{2,2}; sizeSamples = size(data(1,1).(def.time)); nSamples = sizeSamples(2); cfgGM = []; cfgGM.keepindividual = 'yes'; cfgGM.latency = cfg.latency; %% parameters specific to this function cfg.statistic = 'indepsamplesT'; % cfg.parameter = string (default = 'trial' or 'avg') % design design = zeros(1,nGroupLevel1 + nGroupLevel2); design(1,1:nGroupLevel1) = 1; design(1,(nGroupLevel1 + 1):(nGroupLevel1 + nGroupLevel2))= 2; cfg.design = design; % design matrix cfg.ivar = 1; % number or list with indices indicating the independent variable(s) %% loop on criteria for roi_ind = 1:length(def.list.roi) roi = char(def.list.roi(roi_ind)); fig = figure; set(fig,'Units','pixels'); set(fig,'Position',[10 415 1419 399]); % AxesTitleFont = ['arial']; nplot = 1; for other_ind = 1:length(def.list.other) other = char(def.list.other(other_ind)); component = cell(2,2,1); % 2 groups and 2 conditions % 3° dim -> length(def.list.subject) but initialized to 1 ind_component = zeros(2,2); for subject_ind = 1:length(def.list.subject) subject = char(def.list.subject(subject_ind)); group = char(def.list.group(subject_ind)); group_ind = find(ismember(groupDesc(:,1),group)); for cond_ind = 1:length(def.list.cond) cond = char(def.list.cond(cond_ind)); disp([roi,' ',cond,' ',other,' ',group,' ',subject]); object_ind = find(strcmp({data.(def.roi)},roi) & strcmp({data.(def.cond)},cond)... & strcmp({data.(def.other)},other) & strcmp({data.(def.subject)},subject)); objectstat = [group cond subject]; eval([objectstat ' = structTemplate;']); eval([objectstat '.trial{1,1} = data(' num2str(object_ind) ').(def.value);']); eval([objectstat '.time{1,1} = data(' num2str(object_ind) ').(def.time);']); eval([objectstat '.fsample = nSamples;']); eval([objectstat '.sampleinfo = [1 nSamples];']); eval([objectstat '.trialinfo = 1;']); ind_component(group_ind,cond_ind) = ind_component(group_ind,cond_ind) + 1; component{group_ind,cond_ind,ind_component(group_ind,cond_ind)} = eval(objectstat); end end % grand average by group and condition % GA11: group level 1 - cond level 1, GA12: group level 1 - cond level 2 % GA21: group level 2 - cond level 1, GA22: group level 2 - cond level 2 GA11 = ft_timelockgrandaverage(cfgGM,component{1,1,1:ind_component(1,1)}); GA12 = ft_timelockgrandaverage(cfgGM,component{1,2,1:ind_component(1,2)}); GA21 = ft_timelockgrandaverage(cfgGM,component{2,1,1:ind_component(2,1)}); GA22 = ft_timelockgrandaverage(cfgGM,component{2,2,1:ind_component(2,2)}); set1.GM = GA11; set1.GM.individual = GA11.individual - GA12.individual; set2.GM = GA11; set2.GM.individual = GA21.individual - GA22.individual; set1.mean = squeeze(mean(set1.GM.individual,1)); set1.stdmean = squeeze(std(set1.GM.individual,1)/sqrt(nGroupLevel1)); set2.mean = squeeze(mean(set2.GM.individual,1)); set2.stdmean = squeeze(std(set2.GM.individual,1)/sqrt(nGroupLevel2)); % test calculus stat_dif1_dif2 = ft_timelockstatistics(cfg,set1.GM,set2.GM); disp([roi,' ',cond,' ',other,' mask: ',num2str(sum(stat_dif1_dif2.mask))]); % figures subplot(1,2,nplot); nplot = nplot + 1; set1.legend = [def.list.cond{1} '-' def.list.cond{2} ' ' groupDesc{1,1}]; set2.legend = [def.list.cond{1} '-' def.list.cond{2} ' ' groupDesc{2,1}]; stat_dif1_dif2.title = [roi,' ',other]; stat_dif1_dif2.ylabel = ['difference between averages (' def.value ')']; res = ins_temporalClusterPermuTest_graph(set1, set2, stat_dif1_dif2, def, cfg); end filename = ['./output/f_interaction_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig']; disp([' -> saving ' filename]); saveas(gcf,filename); end %% result = 'END of script ins_temporalClusterPermuTest_interaction.m';
% ins_temporalClusterPermuTest_within.m function result = ins_temporalClusterPermuTest_within(data, def, cfg, structTemplate) % http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%% % % Nonparametric MEG data analysis of effects within two conditions % on two levels (paired samples). Permutation test based on temporal clusters on the % data grouped into regions of interest. % % DEPENDANCES: % This function uses the FieldTrip MATLAB software toolbox % (http://www.fieldtriptoolbox.org/) % % USAGE: % result = ins_temporalClusterPermuTest_within( % data, % def, % cfg, % structTemplate) % % INPUTS: % cf. function ins_temporalClusterPermuTest_WithinBetween.m % % OUTPUTS: % result = %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % REFERENCES: % * Cluster-based permutation tests on event related fields. % http://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock % ________________________________ % Bernard Giusiano & Sophie Chen % INSERM UMR 1106 Institut de Neurosciences des Systèmes % Sept/2015 (first version) % Oct/2015 (this version) % http://ins.univ-amu.fr %% dimensions and grand average parameters sizeSamples = size(data(1,1).(def.time)); nSamples = sizeSamples(2); cfgGM = []; cfgGM.keepindividual = 'yes'; cfgGM.latency = cfg.latency; %% parameters specific to this function cfg.statistic = 'depsamplesT'; % paired samples % cfg.parameter = string (default = 'trial' or 'avg') % design nSubjects = length(def.list.subject); design = zeros(2,2*nSubjects); for i = 1:nSubjects design(1,i) = i; end for i = 1:nSubjects design(1,nSubjects+i) = i; end design(2,1:nSubjects) = 1; design(2,nSubjects+1:2*nSubjects) = 2; cfg.design = design; % design matrix cfg.uvar = 1; % indice of the unit variable cfg.ivar = 2; % number or list with indices indicating the independent variable(s) %% loop on criteria for roi_ind = 1:length(def.list.roi) roi = char(def.list.roi(roi_ind)); fig = figure; set(fig,'Units','pixels'); set(fig,'Position',[10 415 1419 399]); nplot = 1; for other_ind = 1:length(def.list.other) other = char(def.list.other(other_ind)); % components of grand average component = cell(2,1); % 2 conditions. 2° dim -> length(def.list.subject) but initialized to 1 ind_component = zeros(2); for subject_ind = 1:length(def.list.subject) subject = char(def.list.subject(subject_ind)); for cond_ind = 1:length(def.list.cond) cond = char(def.list.cond(cond_ind)); disp([roi,' ',cond,' ',other,' ',subject]); object_ind = find(strcmp({data.(def.roi)},roi) & strcmp({data.(def.cond)},cond)... & strcmp({data.(def.other)},other) & strcmp({data.(def.subject)},subject)); objectstat = [cond subject]; eval([objectstat ' = structTemplate;']); eval([objectstat '.trial{1,1} = data(' num2str(object_ind) ').(def.value);']); eval([objectstat '.time{1,1} = data(' num2str(object_ind) ').(def.time);']); eval([objectstat '.fsample = nSamples;']); eval([objectstat '.sampleinfo = [1 nSamples];']); eval([objectstat '.trialinfo = 1;']); ind_component(cond_ind) = ind_component(cond_ind) + 1; component{cond_ind,ind_component(cond_ind)} = eval(objectstat); end end % grand average by condition % set1.GM: cond level 1, set2.GM: cond level 2 set1.GM = ft_timelockgrandaverage(cfgGM,component{1,1:ind_component(1)}); set2.GM = ft_timelockgrandaverage(cfgGM,component{2,1:ind_component(2)}); set1.mean = squeeze(mean(set1.GM.individual,1)); set1.stdmean = squeeze(std(set1.GM.individual,1)/sqrt(nSubjects)); set2.mean = squeeze(mean(set2.GM.individual,1)); set2.stdmean = squeeze(std(set2.GM.individual,1)/sqrt(nSubjects)); % test calculus stat_cond1_cond2 = ft_timelockstatistics(cfg,set1.GM,set2.GM); disp([roi,' ',cond,' ',other,' ',subject,' mask: ',num2str(sum(stat_cond1_cond2.mask))]); % figures subplot(1,2,nplot); nplot = nplot + 1; set1.legend = def.list.cond{1}; set2.legend = def.list.cond{2}; stat_cond1_cond2.title = [roi,' ',other]; stat_cond1_cond2.ylabel = def.value; res = ins_temporalClusterPermuTest_graph(set1, set2, stat_cond1_cond2, def, cfg); end filename = ['./output/f_within_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig']; disp([' -> saving ' filename]); saveas(gcf,filename); end %% result = 'END of script ins_temporalClusterPermuTest_within.m';
