Cluster-based permutation tests on event related fields

The objective of this tutorial is to give an introduction to the statistical analysis of event-related EEG and MEG data (denoted as MEEG data in the following) by means of cluster-based permutation tests. The tutorial starts with a long background section that sketches the background of permutation tests. Subsequently it is shown how to use FieldTrip to perform cluster-based permutation tests on actual axial and planar event-related fields in a between-trials and in a within-subjects design.

In this tutorial we will continue working on the dataset described in the preprocessing tutorials. Below we will repeat code to select the trials and preprocess the data as described in the first tutorials (trigger based trial selection, artifact rejection, event related averaging and planar gradient). We assume that the preprocessing and averaging steps of the analysis are already clear for the reader.

This tutorial is not covering statistical test on time-frequency representations. If you are interested in that, you can read the Cluster-based permutation tests on time-frequency data tutorial. If you are interested how parametric statistical tests can be used with FieldTrip, you can read the Parametric and non-parametric statistics on event-related fields tutorial.

In this tutorial we will consider a between-trials experiment, in which we analyze the data of a single subject. The statistical analysis for this experiment we perform both on axial and planar ERFs. The steps we perform are as follows:

• Preprocessing and time-locked analysis with the ft_definetrial, ft_preprocessing and ft_timelockanalysis functions
• (Calculation of the planar gradient with the ft_megplanar and ft_combineplanar functions)
• Permutation test with the ft_timelockstatistics function
• Plotting the result with the ft_topoplotER function

Figure 1. Analysis protocol of a between-trials experiment with axial data

Figure 2. Analysis protocol of a between-trials experiment with planar data

Subsequently, we consider a within-subjects experiment, in which we compare differences between (planar gradient) ERFs over all subjects. The steps we perform are as follows:

• Preprocessing and time-locked analysis with ft_definetrial, ft_preprocessing and ft_timelockanalysis functions
• Calculation of the planar gradient with the ft_megplanar and ft_combineplanar functions
• Make grandaverage with the ft_timelockgrandaverage function
• Permutation test with the ft_timelockstatistics function
• Plotting the result with the ft_topoplotER function

Figure 3. Analysis protocol of a within-subjects experiment with planar data

In a between-trials experiment, we analyze the data of a single subject. By means of a statistical test, we want to answer the question whether there is a systematic difference in the MEG recorded on trials with a fully congruent and trials with a fully incongruent sentence ending. (This comparison is the MEG-version of the classical comparison that showed the N400-component in EEG-data).

We first extract the trials of the fully incongruent condition

% find the interesting segments of data
cfg = [];                                           % empty configuration
cfg.dataset                 = 'Subject01.ds';       % name of CTF dataset  
cfg.trialdef.eventtype      = 'backpanel trigger';
cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 2;
cfg.trialdef.eventvalue     = 3;                    % trigger value for fully incongruent (FIC)
cfg = ft_definetrial(cfg);            

% remove the trials that have artifacts from the trl
cfg.trl([15, 36, 39, 42, 43, 49, 50, 81, 82, 84],:) = []; 

% preprocess the data
cfg.channel    = {'MEG', '-MLP31', '-MLO12'};        % read all MEG channels except MLP31 and MLO12
cfg.demean     = 'yes';
cfg.baselinewindow  = [-0.2 0];
cfg.lpfilter   = 'yes';                              % apply lowpass filter
cfg.lpfreq     = 35;                                 % lowpass at 35 Hz.

dataFIC_LP = ft_preprocessing(cfg);                      

save dataFIC_LP dataFIC_LP

Then we also extract the trails of the fully congruent condition.

% find the interesting segments of data
cfg = [];                                           % empty configuration
cfg.dataset                 = 'Subject01.ds';       % name of CTF dataset  
cfg.trialdef.eventtype      = 'backpanel trigger';
cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 2;
cfg.trialdef.eventvalue     = 9;                    % trigger value for fully incongruent (FC)
cfg = ft_definetrial(cfg);            

% remove the trials that have artifacts from the trl
cfg.trl([2, 3, 4, 30, 39, 40, 41, 45, 46, 47, 51, 53, 59, 77, 85],:) = []; 

% preprocess the data
cfg.channel    = {'MEG', '-MLP31', '-MLO12'};       % read all MEG channels except MLP31 and MLO12
cfg.demean     = 'yes';
cfg.baselinewindow  = [-0.2 0];
cfg.lpfilter   = 'yes';                              % apply lowpass filter
cfg.lpfreq     = 35;                                 % lowpass at 35 Hz.

dataFC_LP = ft_preprocessing(cfg);                      

save dataFC_LP dataFC_LP

The format of the preprocessed data is inconvenient for statistical analyses. A more convenient format is produced by ft_timelockanalysis with the option cfg.keeptrials = 'yes'. The output of ft_timelockanalysis contains an .avg field with the average event-related field, and a .trial field with the individual trial data. The output is stored in timelockFIC and timelockFC for the fully incongruent and the fully congruent condition.

Ft_timelockanalysis requires preprocessed data (see above), which is available from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/cluster_permutation_timelock/dataFIC_LP.mat and dataFC_LP.mat.

load dataFIC_LP
load dataFC_LP

cfg = [];
cfg.keeptrials = 'yes';
timelockFIC = ft_timelockanalysis(cfg, dataFIC_LP);
timelockFC  = ft_timelockanalysis(cfg, dataFC_LP);

Cluster-level permutation tests for event-related fields are performed by the function ft_timelockstatistics. This function takes as its arguments a configuration (cfg) and two or more data structures. These data structures must be produced by ft_timelockanalysis or ft_timelockgrandaverage, which all operate on preprocessed data. The argument list of ft_timelockstatistics must contain one data structure for every experimental condition. For comparing the data structures timelockFIC and timelockFC, you must call ft_timelockstatistics as follows: [stat] = ft_timelockstatistics(cfg, timelockFIC, timelockFC);

Some fields of the configuration (cfg), such as channel and latency, are not specific for ft_timelockstatistics; their role is similar in other FieldTrip functions. We first concentrate on the fields that are specific for ft_timelockstatistics.

cfg = [];
cfg.method = 'montecarlo';       % use the Monte Carlo Method to calculate the significance probability
cfg.statistic = 'indepsamplesT'; % use the independent samples T-statistic as a measure to evaluate 
                                 % the effect at the sample level
cfg.correctm = 'cluster';
cfg.clusteralpha = 0.05;         % alpha level of the sample-specific test statistic that will be used for thresholding
cfg.clusterstatistic = 'maxsum'; % test statistic that will be evaluated under the permutation distribution. 
cfg.minnbchan = 2;               % minimum number of neighborhood channels that is required for a selected 
cfg.neighbours = neighbours;
                                 % sample to be included in the clustering algorithm (default=0).
cfg.tail = 0;                    % -1, 1 or 0 (default = 0); one-sided or two-sided test
cfg.clustertail = 0;
cfg.alpha = 0.025;               % alpha level of the permutation test
cfg.numrandomization = 100;      % number of draws from the permutation distribution

design = zeros(1,size(timelockFIC.trial,1) + size(timelockFC.trial,1));
design(1,1:size(timelockFIC.trial,1)) = 1;
design(1,(size(timelockFIC.trial,1)+1):(size(timelockFIC.trial,1) + size(timelockFC.trial,1)))= 2;

cfg.design = design;             % design matrix
cfg.ivar  = 1;                   % number or list with indices, independent variable(s)

We now describe these options one-by-one.

• With cfg.method = 'montecarlo' we choose the Monte Carlo method for calculating the significance probability. This significance probability is a Monte Carlo estimate of the p-value under the permutation distribution.
• With cfg.statistic = 'indepsamplesT', we choose the independent samples T-statistic to evaluate the effect (the difference between the fully congruent and the fully incongruent condition) at the sample level. In cfg.statistic, many other test statistics can be specified. Which test statistic is appropriate depends on your research question and your experimental design. For instance, in a within-UO design, one must use the dependent samples T-statistic ('depsamplesT'). And if you want to compare more than two experimental conditions, you should choose an F-statistic ('indepsamplesF' or 'depsamplesF').
• We use cfg.clusteralpha to choose the critical value that will be used for thresholding the sample-specific T-statistics. With cfg.clusteralpha = 0.05, every sample-specific T-statistic is compared with the critical value of the univariate T-test with a critical alpha-level of 0.05. (This statistical test would have been the most appropriate test if we had observed a single channel at a single time-point.) The value of cfg.clusteralpha does not affect the false alarm rate of the statistical test at the cluster-level. It is a rational threshold for deciding whether a sample should be considered a member of some large cluster of samples (which may or may not be significant at the cluster-level).
• We use cfg.clusterstatistic to choose the test statistic that will be evaluated under the permutation distribution. This is the actual test statistic and it must be distinguished from the sample-specific T-statistics that are used for thresholding. With cfg.clusterstatistic = 'maxsum', the actual test statistic is the maximum of the cluster-level statistics. A cluster-level statistic is equal to the sum of the sample-specific T-statistics that belong to this cluster. Taking the largest of these cluster-level statistics of the different clusters produces the actual test statistic.
• The value of cfg.minnbchan is a tuning parameter that determines the way the clusters are formed. More specifically, we use cfg.minnbchan to specify the minimum number of neighborhood channels that is required for a selected sample (i.e., a sample who's T-value exceeds the threshold) to be included in the clustering algorithm. With cfg.minnbchan = 0 (the default), it sometimes happens that two clusters are spatially connected via a narrow bridge of samples. Because they are connected, these two clusters are considered as a single cluster. If clusters are interpreted as reflecting spatially distinct sources, such a combined cluster does not make much sense. To suppress this type of combined clusters, one can choose to ignore all selected samples (on the basis of their T-values) if they have less than some minimum number of neighbors that were also selected. This minimum number is assigned to cfg.minnbchan. This number must be chosen independently of the data.
• cfg.neighbours is a structure that you need to have previously created using ft_neighbourselection.
• We use cfg.tail to choose between a one-sided and a two-sided statistical test. Choosing cfg.tail = 0 affects the calculations in three ways. First, the sample-specific T-values are thresholded from below as well as from above. This implies that both large negative and large positive T-statistics are selected for later clustering. Second, clustering is performed separately for thresholded positive and thresholded negative T-statistics. And third, the critical value for the cluster-level test statistic (determined by cfg.alpha; see further) is now two-sided: negative cluster-level statistics must be compared with the negative critical value, and positive cluster-level statistics must be compared with the positive critical value.
• We use cfg.alpha to control the false alarm rate of the permutation test (the probability of falsely rejecting the null hypothesis). The value of cfg.alpha determines the critical values with which we must compare the test statistic (i.e., the maximum and the minimum cluster-level statistic). Note that if you want to run a two-sided test, you have to split the critical alpha value by setting cfg.correcttail = 'alpha'; i.e. this sets cfg.alpha = 0.025, corresponding to a false alarm rate of 0.05 in a two-sided test. The field cfg.alpha is not crucial. This is because the output of ft_timelockstatistics (see further) contains a p-value for every cluster (calculated under the permutation distribution of the maximum/minimum cluster-level statistic). Instead of the critical values, we can also use these p-values to determine the significance of the clusters.
• We use cfg.numrandomization to control the number of draws from the permutation distribution. Remember that ft_timelockstatistics approximates the permutation distribution by means of a histogram. This histogram is a so-called Monte Carlo approximation of the permutation distribution. This Monte Carlo approximation is used to calculate the p-values and the critical values that are shown in the output of ft_timelockstatistics. In this tutorial, we use cfg.numrandomization = 100. This number is too small for actual applications. As a rule of thumb, use cfg.numrandomization = 500, and double this number if it turns out that the p-value differs from the critical alpha-level (0.05 or 0.01) by less than 0.02.
• We use cfg.design to store information about the UOs. The content of cfg.design must be a matrix. Consider the hypothetical case that your subject has performed 5 trials in the first condition and 4 trials in the second condition. Then the correct design matrix looks like this: cfg.design = [1 1 1 1 1 2 2 2 2].
• We use cfg.ivar to indicate the row of the design matrix that contains the independent variable. For a between-trials statistical analysis, the minimum design matrix contains only a single row, and cfg.ivar is superfluous. However, in other statistical analyses (e.g., those for a within-UO design) the design matrix must contain more than one row, and specifying cfg.ivar is essential.

We now briefly discuss the configuration fields that are not specific for ft_timelockstatistics:

cfg.channel = {'MEG'};          % cell-array with selected channel labels
cfg.latency = [0 1];            % time interval over which the experimental 
                                % conditions must be compared (in seconds)
cfg_neighb.method = 'distance';         
cfg.neighbours = ...            % specifies with which sensors other sensors
ft_prepare_neighbours(...       % can form clusters
cfg_neighb, dataFC_LP);

With these two options, we select the spatio-temporal dataset involving all MEG channels and the time interval between 0 and 1 second. The two experimental conditions will only be compared on this selection of the complete spatio-temporal dataset. Also, feel free to consult our frequenly asked questions about prepare_neighbours

One should be aware of the fact that the sensitivity of ft_timelockstatistics (i.e., the probability of detecting an effect) depends on the length of the time interval that is analyzed, as specified in cfg.latency. For instance, assume that the difference between the two experimental conditions extends over a short time interval only (e.g., between 0.3 and 0.4 sec.). If it is known in advance that this short time interval is the only interval where an effect is likely to occur, then one should limit the analysis to this time interval (i.e., choose cfg.latency = [0.3 0.4]). Choosing a time interval on the basis of prior information about the time course of the effect will increase the sensitivity of the statistical test. If there is no prior information, then one must compare the experimental conditions over the complete time interval. This is accomplished by choosing cfg.latency = 'all'. In this tutorial, we choose cfg.latency = [0 1] because EEG-studies have shown that the strongest effect of semantic incongruity is observed in the first second after stimulus presentation.

Now, run ft_timelockstatistics to compare timelockFIC and timelockFC using the configuration described above.

[stat] = ft_timelockstatistics(cfg, timelockFIC, timelockFC);

Save the output to disk:

save stat_ERF_axial_FICvsFC stat;

The output can also be obtained from stat_ERF_axial_FICvsFC.mat. If you need to reload the statistics output, use:

load stat_ERF_axial_FICvsFC

The output of ft_timelockstatistics has separate fields for positive and negative clusters. For the positive clusters, the output is given in the following pair of fields: stat.posclusters and stat.posclusterslabelmat. The field stat.posclusters is an array that provides the following information for every cluster:

• The field clusterstat contains the cluster-level statistic (the sum of the T-values in this cluster).
• The field prob contains the proportion of draws from the permutation distribution with a maximum cluster-level statistic that is larger than clusterstat.

The elements in the array stat.posclusters are sorted according to their p-value: the cluster with the smallest p-value comes first, followed by the cluster with the second-smallest, etc. Thus, if the k-th cluster has a p-value that is larger than the critical alpha-level (e.g., 0.025), then so does the (k+1)-th. Type stat.posclusters(k) on the Matlab command line to see the information for the k-th cluster.

The field stat.posclusterslabelmat is a spatiotemporal matrix. This matrix contains numbers that identify the clusters to which the (channel,time)-pairs (the samples) belong. For example, all (channel,time)-pairs that belong to the third cluster, are identified by the number 3. As will be shown in the following, this information can be used to visualize the topography of the clusters.

For the negative clusters, the output is given in the following pair of fields: stat.negclusters and stat.negclusterslabelmat. These fields contain the same type of information as stat.posclusters and stat.posclusterslabelmat, but now for the negative clusters.

By inspecting stat.posclusters and stat.negclusters, it can be seen that only the first positive and the first negative cluster have a p-value less than the critical alpha-level of 0.025. This critical alpha-level corresponds to a false alarm rate of 0.05 in a two-sided test. By typing stat.posclusters(1) on the Matlab command line, you should obtain the following:

stat.posclusters(1)
ans = 
         prob: 0
  clusterstat: 8.2721e+03

And by typing stat.negclusters(1), you should obtain the following:

stat.negclusters(1)
ans = 
         prob: 0
  clusterstat: -1.0251e+04

It is possible that the p-values in your output are a little bit different from 0. This is because ft_timelockstatistics calculated as a Monte Carlo approximation of the permutation p-values: the p-value for the k-th positive cluster is calculated as the proportion of random draws from the permutation distribution in which the maximum of the cluster-level statistics is more larger than stat.posclusters(k).clusterstat.

To plot the results of the permutation test, we use the plotting function ft_topoplotER. In doing so, we will plot a topography of the difference between the two experimental conditions (FC and FIC). Atop that, and for each timestep of interest, we'll highlight the sensors which are members of significant clusters. First, however, we must calculate the difference between conditions.

cfg = [];
cfg.keeptrials = 'no';  % now keep only the subject-wise average, not the single trials
avgFIC = ft_timelockanalysis(cfg, dataFIC_LP);
avgFC  = ft_timelockanalysis(cfg, dataFC_LP);

% Copy the entire timelockanalysis structure to preserve all
% the information it holds in addition to ERP averages.
raweffectFICvsFC     = avgFIC;
% Then take the difference of the averages.
raweffectFICvsFC.avg = avgFIC.avg - avgFC.avg;  

We then construct a boolean matrix indicating membership in the significant clusters. This matrix has size [Number_of_MEEG_channels × Number_of_temporal_samples], like stat.posclusterslabelmat. We'll make two such matrices: one for positive clusters (named pos), and one for negative (neg). All (channel,time)-pairs belonging to the significant clusters will be coded in the new boolean matrix as 1, and all those that don't will be coded as 0.

In plotting significant clusters, we must of course first determine which clusters are reliable.

% Make a vector of all p-values associated with the clusters from ft_timelockstatistics.
pos_cluster_pvals = [stat.posclusters(:).prob];
% Then, find which clusters are significant, outputting their indices as held in stat.posclusters
pos_signif_clust = find(pos_cluster_pvals < stat.cfg.alpha);
% (stat.cfg.alpha is the alpha level we specified earlier for cluster comparisons; In this case, 0.025)
% make a boolean matrix of which (channel,time)-pairs are part of a significant cluster
pos = ismember(stat.posclusterslabelmat, pos_signif_clust);

% and now for the negative clusters...
neg_cluster_pvals = [stat.negclusters(:).prob];
neg_signif_clust = find(neg_cluster_pvals < stat.cfg.alpha);
neg = ismember(stat.negclusterslabelmat, neg_signif_clust);

Alternatively, we can manually select which clusters we want to plot. If we only want to see the extext of the first (i.e. most significant) positive and negative clusters, for instance, we can do so as follows:

pos = stat.posclusterslabelmat == 1;	% or == 2, or 3, etc.
neg = stat.negclusterslabelmat == 1;

To plot a sequence of twenty topographic plots equally spaced between 0 and 1 second, we define the vector j of time steps. These time intervals correspond to the samples m in stat and in the variables pos and neg. m and j must, therefore, have the same length.

To be sure that your sample-based time windows align with your time windows in seconds, check the following:

timestep = 0.05;		% timestep between time windows for each subplot (in seconds)
sampling_rate = dataFC_LP.fsample;	% Data has a temporal resolution of 300 Hz
sample_count = length(stat.time);
					% number of temporal samples in the statistics object
j = [0:timestep:1];   % Temporal endpoints (in seconds) of the ERP average computed in each subplot
m = [1:timestep*sampling_rate:sample_count];  % temporal endpoints in MEEG samples

To plot the data use the following for-loop:

for k = 1:20;
     subplot(4,5,k);
     cfg = [];   
     cfg.xlim=[j(k) j(k+1)];   % time interval of the subplot
     cfg.zlim = [-2.5e-13 2.5e-13];
   % If a channel reaches this significance, then
   % the element of pos_int with an index equal to that channel
   % number will be set to 1 (otherwise 0).
   
   % Next, check which channels are significant over the
   % entire time interval of interest.
     pos_int = all(pos(:, m(k):m(k+1)), 2);
     neg_int = all(neg(:, m(k):m(k+1)), 2);

     cfg.highlight = 'on';
   % Get the index of each significant channel
     cfg.highlightchannel = find(pos_int | neg_int);
     cfg.comment = 'xlim';   
     cfg.commentpos = 'title';   
     cfg.layout = 'CTF151.lay';
     ft_topoplotER(cfg, raweffectFICvsFC);   
end

In this for-loop, cfg.xlim defines the time interval of each subplot. The variables pos_int and neg_int boolean vectors indicating which channels of pos and neg are significant in the time interval of interest. This is defined in cfg.highlight. The for-loop plots 20 subplots covering a time interval of 50 ms each. Running this for-loop creates the following figure:

Figure 4: Raw effect (FIC-FC) on the ERFs of subject 1, significant clusters are highlighted.

To perform the permutation test using synthetic planar gradient data, the data must first be converted using the functions ft_megplanar and ft_combineplanar. These functions were described in the tutorial on event-related fields. After running these functions, the statistical analysis using ft_timelockstatistics involves the same configuration options as for ordinary event-related averages (no synthetic planar gradients). There is only one additional step, which is needed to add the gradiometer structure to one of the planar gradient data sets.

cfg = [];
cfg.planarmethod = 'sincos';
timelockFIC_planar = ft_megplanar(cfg, timelockFIC);
timelockFC_planar  = ft_megplanar(cfg, timelockFC);
       
timelockFIC_planar_cmb = ft_combineplanar(cfg, timelockFIC_planar);
timelockFC_planar_cmb  = ft_combineplanar(cfg, timelockFC_planar);
  
timelockFIC_planar_cmb.grad = timelockFIC.grad;  % add the gradiometer structure
timelockFC_planar_cmb.grad  = timelockFC.grad;
     

Having calculated synthetic planar gradient data, one can use the same configuration parameters as used for the analysis of the original data.

cfg = [];
cfg.channel = {'MEG'};
cfg.latency = [0 1];
 
cfg.method = 'montecarlo';
cfg.statistic = 'indepsamplesT';
cfg.correctm = 'cluster';
cfg.clusteralpha = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan = 2;
cfg.tail = 0;
cfg.clustertail = 0;
cfg.alpha = 0.025;
cfg.numrandomization = 100;

design = zeros(1,size(timelockFIC_planar_cmb.trial,1) + size(timelockFC_planar_cmb.trial,1));
design(1,1:size(timelockFIC_planar_cmb.trial,1)) = 1;
design(1,(size(timelockFIC_planar_cmb.trial,1)+1):(size(timelockFIC_planar_cmb.trial,1) + size(timelockFC_planar_cmb.trial,1)))= 2;

cfg.design = design;
cfg.ivar = 1;

[stat] = ft_timelockstatistics(cfg, timelockFIC_planar_cmb, timelockFC_planar_cmb)

save stat_ERF_planar_FICvsFC stat

The output can also be obtained from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/cluster_permutation_timelock/stat_ERF_planar_FICvsFC.mat. If you need to reload the statistics output, use:

load stat_ERF_planar_FICvsFC

We now calculate the raw effect in the average with planar gradient data using the following configuration:

cfg = [];
cfg.keeptrials = 'no';   % now only the average, not the single trials
avgFIC_planar = ft_timelockanalysis(cfg, timelockFIC_planar);
avgFC_planar  = ft_timelockanalysis(cfg, timelockFC_planar);
cfg = [];
avgFIC_planar_cmb = ft_combineplanar(cfg, avgFIC_planar);
avgFC_planar_cmb  = ft_combineplanar(cfg, avgFC_planar);

% make a copy of the structure to contain the raw effect and all other information
% and subtract the averages from each other
raweffectFICvsFC     = avgFIC_planar_cmb;
raweffectFICvsFC.avg = avgFIC_planar_cmb.avg - avgFC_planar_cmb.avg;  

Using the following configuration for ft_topoplotER we can plot the raw effect and highlight the channels belonging to the significant clusters:

figure;  
timestep = 0.05;		%(in seconds)
sampling_rate = dataFC_LP.fsample;
sample_count = length(stat.time);
j = [0:timestep:1];   % Temporal endpoints (in seconds) of the ERP average computed in each subplot
m = [1:timestep*sampling_rate:sample_count];  % temporal endpoints in MEEG samples

pos_cluster_pvals = [stat.posclusters(:).prob];
pos_signif_clust = find(pos_cluster_pvals < stat.cfg.alpha);
pos = ismember(stat.posclusterslabelmat, pos_signif_clust);

% Remember to do the same for negative clusters if you want them!

for k = 1:20;
     subplot(4,5,k);   
     cfg = [];
     cfg.xlim =[j(k) j(k+1)];
     cfg.zlim = [-1.0e-13 1.0e-13];   
     pos_int = all(pos(:, m(k):m(k+1)), 2);
     cfg.highlight = 'on';
     cfg.highlightchannel = find(pos_int);
     cfg.comment = 'xlim';
     cfg.commentpos = 'title';
     cfg.layout = 'CTF151.lay';
     ft_topoplotER(cfg, raweffectFICvsFC);
end

Figure 5: Raw effect (FIC-FC) on the planar gradient ERFs of subject 1, the significant clusters are highlighted.

We now consider experiments involving multiple subjects that are each observed in multiple experimental conditions. Typically, every subject is observed in a large number of trials, each one belonging to one experimental condition. Usually, for every subject, averages are computed over all trials belonging to each of the experimental conditions. Thus, for every subject, the data are summarized in an array of condition-specific averages. The permutation test that is described in this section informs us about the following null hypothesis: the probability distribution of the condition-specific averages is independent of the experimental conditions.

We now describe how we can statistically test the difference between the event-related averages for fully incongruent (FIC) and the fully congruent (FC) sentence endings. For this analysis we use planar gradient data. For convenience we will not do the reading-in and preprocessing steps on all subjects. Instead we begin by loading the timelock structures containing the event-related averages (of the planar gradient data) of all ten subjects. The data is available from the FieldTrip FTP server (GA_ERF_orig.mat).

load GA_ERF_orig;

The event-related averages for the fully incongruent and the fully congruent sentence endings are stored in GA_FIC and GA_FC. These averages were calculated using the function ft_timelockgrandaverage. In the configuration for ft_timelockgrandaverage, the option cfg.keepindividual = ‘yes’ has to be chosen. (This has been done in the example data.) If this option is not used, then only the average ERF of all subjects will be calculated.

We now perform the permutation test using ft_timelockstatistics. The configuration settings for this analysis differ from the previous settings in several fields:

1. We have to select a different measure to evaluate the effect at sample level (in cfg.statistic)
2. The design matrix is different (i.c., it now contains two lines instead of one)
3. The so-called unit variable has to be defined.

The configuration looks as follows:

cfg = [];
cfg.channel = {'MEG'};
cfg.latency = [0 1];

cfg.method = 'montecarlo';
cfg.statistic = 'depsamplesT';
cfg.correctm = 'cluster';
cfg.clusteralpha = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan = 2;
cfg.neighbours = neighbours;
cfg.tail = 0;
cfg.clustertail = 0;
cfg.alpha = 0.025;
cfg.numrandomization = 500;

subj = 10;
design = zeros(2,2*subj);
for i = 1:subj
  design(1,i) = i;
end
for i = 1:subj
  design(1,subj+i) = i;
end
design(2,1:subj)        = 1;
design(2,subj+1:2*subj) = 2;

cfg.design = design;
cfg.uvar  = 1;
cfg.ivar  = 2;

We now describe the differences between this configuration and the configuration for a between-trials experiment.

• Instead of an independent samples T-statistic, we use the dependent samples T-statistic to evaluate the effect at the sample level (cfg.statistic = 'depsamplesT'). This is because we are dealing with a within-UO instead of a between-UO design.
• The design matrix in a within-UO design is different from the design matrix in a between-UO design. In the design matix for a within-UO design, you have to specify the unit variable. The unit variable specifies the units that have produced the different condition-specific data structures. For example, consider a hypothetical study with 4 subjects and 2 experimental conditions. The design matrix may then look like this: design = [1 2 3 4 1 2 3 4; 1 1 1 1 2 2 2 2 ]. The first row of this matrix is the unit variable: it specifies that the first subject produced the first and the fifth data structure, the second subject produced the second and the sixth data structure, etc. The second row of the design matrix is the independent variable.
• Because the design matrix contains both a unit variable and an independent variable, it has to be specified in the configuration which row contains which variable. This information is passed in the fields cfg.uvar (for the unit variable) and cfg.ivar (for the independent variable).

Now, use the configuration above to perform the following statistical analysis:

[stat] = ft_timelockstatistics(cfg, GA_FIC, GA_FC)

save stat_ERF_planar_FICvsFC_GA stat

The output can also be obtained from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/cluster_permutation_timelock/stat_ERF_planar_FICvsFC_GA.mat. If you need to reload the statistics output, use:

load stat_ERF_planar_FICvsFC_GA

From inspection of stat.posclusters and stat.negclusters, we observe that there is only one significant positive cluster and no significant negative cluster.

To plot the results:

GA_FICvsFC = GA_FIC;
GA_FICvsFC.avg = GA_FIC.avg - GA_FC.avg;

figure;  
timestep = 0.05;      %(in seconds)
sampling_rate = dataFC_LP.fsample;
sample_count = length(stat.time);
j = [0:timestep:1];   % Temporal endpoints (in seconds) of the ERP average computed in each subplot
m = [1:timestep*sampling_rate:sample_count];  % temporal endpoints in MEEG samples

pos_cluster_pvals = [stat.posclusters(:).prob];
pos_signif_clust = find(pos_cluster_pvals < stat.cfg.alpha);
pos = ismember(stat.posclusterslabelmat, pos_signif_clust);

for k = 1:20;
     subplot(4,5,k);   
     cfg = [];   
     cfg.xlim=[j(k) j(k+1)];   
     cfg.zlim = [-1.0e-13 1.0e-13];   
     pos_int = all(pos(:, m(k):m(k+1)), 2);
     cfg.highlight = 'on';
     cfg.highlightchannel = find(pos_int);       
     cfg.comment = 'xlim';   
     cfg.commentpos = 'title';   
     ft_topoplotER(cfg, GA_FICvsFC);
end  

Figure 6: Raw effect (FIC-FC) on the grand average planar gradient ERFs the significant cluster is highlighted

In this tutorial, it was shown how to do non-parametric statistics on axial and planar ERFs in a between-trials and in within-subjects design. It was also shown how to plot the results.

If you are interested in parametric tests in FieldTrip, you can read the Parametric and non-parametric statistics on event-related fields tutorial. If you are interested in how to do the same statistics on time-frequency representations, you can read the Cluster-based permutation tests on time-frequency data tutorial.

If you would like to read more about issues related to statistical analysis, you can read the following FAQ's as well:
Why should I use the cfg.correcttail option when using statistics_montecarlo?
What is the idea behind statistical inference at the second-level?

If you would like to read about neighborhood selection, you can read the following FAQ's:
How can I define my own neighbourhood templates or updating an already existing template?
How can I define neighbouring sensors?
How does ft_prepare_neighbours work?

Or you can look also at the following example scripts:
Apply clusterrandanalysis on TFRs of power that were computed with BESA
Source statistics
Stratify the distribution of one variable that differs in two conditions
Using the MATLAB statistics toolbox

This tutorial was last tested by Jan-Mathijs, with version 20110409 of FieldTrip, using Matlab 2009b on a 64-bit Linux platform.