MEG:ANOVA dans Fieldtrip

From WikiMEG
Jump to: navigation, search

Analyse non paramétrique de données MEG à deux facteurs

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.

Ces scripts présentés lors du Club MEG du Mardi 26 Janvier 2016 seront prochainement téléchargeables avec des données exemples.

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';
% ins_temporalClusterPermuTest_between.m
function result = ins_temporalClusterPermuTest_between(data, def, cfg, structTemplate)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
%
% Nonparametric MEG data analysis of effects between two groups
% on two levels. 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_between(
%                               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
groupDesc = tabulate(def.list.group);
nGroupLevel1 = groupDesc{1,2};
nGroupLevel2 = groupDesc{2,2};
groupLevels = unique(def.list.group);
 
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';    % independant samples
% cfg.parameter = string          (default = 'trial' or 'avg')
 
% design
design = zeros(1,nGroupLevel1*2 + nGroupLevel2*2);
design(1,1:nGroupLevel1*2) = 1;
design(1,(nGroupLevel1*2 + 1):(nGroupLevel1*2 + nGroupLevel2*2))= 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]);
    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 (each subject has 2 cond)
							      % 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
        % set1.GM: group level 1, set2.GM: group level 2
        set1.GM = ft_timelockgrandaverage(cfgGM,component{1,1,1:ind_component(1,1)},...
            component{1,2,1:ind_component(1,2)});
        set2.GM = ft_timelockgrandaverage(cfgGM,component{2,1,1:ind_component(2,1)},...
            component{2,2,1:ind_component(2,2)});
 
        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_group1_group2 = ft_timelockstatistics(cfg,set1.GM,set2.GM);
        disp([roi,' ',other,' ',group,' ',subject,' mask: ',num2str(sum(stat_group1_group2.mask))]);
 
        % figure
        subplot(1,2,nplot);
        nplot = nplot + 1;
 
        set1.legend = groupLevels{1};
        set2.legend = groupLevels{2};
        stat_group1_group2.title = [roi,' ',other];
        stat_group1_group2.ylabel = def.value;
        res = ins_temporalClusterPermuTest_graph(set1, set2, stat_group1_group2, def, cfg);
    end
    filename = ['./output/f_between_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig'];
    disp([' -> saving ' filename]);
    saveas(gcf,filename);
end
 
%%
result = 'END of script ins_temporalClusterPermuTest_between.m';
% ins_temporalClusterPermuTest_within_simple.m
function result = ins_temporalClusterPermuTest_within_simple(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) for each group. 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_simple(
%                               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)
% Nov/2015 (this version)
% http://ins.univ-amu.fr
 
%% dimensions and grand average parameters
groupLevels = unique(def.list.group);
 
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 : in the loop on criteria because it may be specific for each group
 
%% 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 74 1419 740]);
    nplot = 1;
 
    for group_ind = 1:length(groupLevels)
        group = char(groupLevels(group_ind));
        groupSubjects = def.list.subject(strcmp(def.list.group,group));
 
        % design
        nSubjects = length(groupSubjects);
        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)
        % end design
 
        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(groupSubjects)
                subject = char(groupSubjects(subject_ind));
 
                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 = [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))]);
 
            % figure
            subplot(2,2,nplot);
            nplot = nplot + 1;
 
            set1.legend = def.list.cond{1};
            set2.legend = def.list.cond{2};
            stat_cond1_cond2.title = [roi,' ',group,' ',other];
            stat_cond1_cond2.ylabel = def.value;
            res = ins_temporalClusterPermuTest_graph(set1, set2, stat_cond1_cond2, def, cfg);
        end
    end
    filename = ['./output/f_within_simple_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig'];
    disp([' -> saving ' filename]);
    saveas(gcf,filename);
end
 
%%
result = 'END of script ins_temporalClusterPermuTest_within_simple.m';
% ins_temporalClusterPermuTest_between_simple.m
function result = ins_temporalClusterPermuTest_between_simple(data, def, cfg, structTemplate)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
%
% Nonparametric MEG data analysis of effects between two groups
% on two levels for each condition. 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_between_simple(
%                               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)
% Nov/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};
groupLevels = unique(def.list.group);
 
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';  % independant samples
% 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 74 1419 740]);
    nplot = 1;
 
    for cond_ind = 1:length(def.list.cond)
        cond = char(def.list.cond(cond_ind));
 
        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));
 
                group = char(def.list.group(subject_ind));
                group_ind = find(ismember(groupDesc(:,1),group));
 
                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 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) = ind_component(group_ind) + 1;
                component{group_ind,ind_component(group_ind)} = eval(objectstat);
            end
 
            % grand average by group
            % set1.GM: group level 1, set2.GM: group 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(nGroupLevel1));
            set2.mean = squeeze(mean(set2.GM.individual,1));
            set2.stdmean = squeeze(std(set2.GM.individual,1)/sqrt(nGroupLevel2));
 
            % test calculus
            stat_group1_group2 = ft_timelockstatistics(cfg,set1.GM,set2.GM);
            disp([roi,' ',other,' ',group,' ',subject,' mask: ',num2str(sum(stat_group1_group2.mask))]);
 
            % figure
            subplot(2,2,nplot);
            nplot = nplot + 1;
 
            set1.legend = groupLevels{1};
            set2.legend = groupLevels{2};
            stat_group1_group2.title = [roi,' ',cond,' ',other];
            stat_group1_group2.ylabel = def.value;
            res = ins_temporalClusterPermuTest_graph(set1, set2, stat_group1_group2, def, cfg);
        end
    end
    filename = ['./output/f_between_simple_' roi '_' datestr(now, 'yyyymmdd_HHMM') '.fig'];
    disp([' -> saving ' filename]);
    saveas(gcf,filename);
end
 
%%
result = 'END of script ins_temporalClusterPermuTest_between_simple.m';

ins_temporalClusterPermuTest_graph.m

% ins_temporalClusterPermuTest_graph.m
function result = ins_temporalClusterPermuTest_graph(set1, set2, stat_set1_set2, def, cfg)
% http://meg.univ-amu.fr/wiki/Main_Page %%%%%%%%%%%%%
%
% Specific function to draw graphics of functions
% ins_temporalClusterPermuTest_XXX
%
% DEPENDANCES:
% This function uses the function 'jbfill' to shade area between two curves
% (http://www.mathworks.com/matlabcentral/fileexchange/13188-shade-area-between-two-curves/content/jbfill.m)
%
% USAGE:
% result = ins_temporalClusterPermuTest_graph(
%                               set1,
%                               set2,
%                               stat_set1_set2,
%                               def,
%                               cfg)
%
% INPUTS:
% cf. function ins_temporalClusterPermuTest_XXX.m
%
% OUTPUTS:
% result = true or false if error
% ________________________________
% Bernard Giusiano & Sophie Chen
% INSERM UMR 1106 Institut de Neurosciences des Systèmes
% Sept/2015 (first version)
% Nov/2015 (this version)
% http://ins.univ-amu.fr
%
result = false;
 
%% trace %%%%%%%%%%%%%%%%
p = min(stat_set1_set2.prob);
st=dbstack;
fprintf(cfg.fid, '%s;%s;%6.4f\n', st(2,1).name, stat_set1_set2.title, p);
%%%%%%%%%%%%%%%%%%%%%%%%%
 
plot(set1.GM.time,set1.mean,'b','linewidth',2);
hold on;
plot(set2.GM.time,set2.mean,'r','linewidth',2);
 
plot(set1.GM.time,stat_set1_set2.stat,'g','linewidth',0.5);
line(cfg.latency,[0 0],'Color',[0.1 0.1 0.1],'linestyle','--');  % xlim
 
y1 = stat_set1_set2.mask;
y2(~y1) = nan ;
y2(y1) = 0.0 ;
plot(set1.GM.time,y2,'Color',[1.0,0.4,0.0],'linewidth',10); % clusters significatifs en orange
 
legend(set1.legend,set2.legend,'stat diff','Location','NorthWest');
 
x = jbfill(set1.GM.time,(set1.mean+1.96*set1.stdmean).',(set1.mean-1.96*set1.stdmean).','b','b',0,0.05);
x = jbfill(set2.GM.time,(set2.mean+1.96*set2.stdmean).',(set2.mean-1.96*set2.stdmean).','r','r',0,0.05);
 
title(stat_set1_set2.title,'FontWeight','bold','FontSize',14)
xlabel(def.time);ylabel(stat_set1_set2.ylabel);
hold off;
 
%%
result = true;