Skip to content
Permalink
Browse files

EXC: exercise for time generalization multiple comparsion correction

  • Loading branch information...
nno committed Jul 21, 2019
1 parent 3c1c844 commit a9cfb10e12ad06f2ec910fcfae7e5f60b8d19a10
Showing with 238 additions and 4 deletions.
  1. +16 −4 doc/source/ex_multiple_comparisons.rst
  2. +222 −0 examples/run_meeg_time_generalization_mcc.m
@@ -15,16 +15,17 @@ A few approaches for whole-brain multiple-comparison methods have been proposed

- Bonferroni correction; typically too conservative.
- FDR: can lead to invalid inferences.
- Parametric cluster-baed methods (SPM, AFNI, FSL): can be too liberal.
- Parametric cluster-based methods (SPM, AFNI, FSL): can be too liberal.
- Fixed-treshold cluster-based approach: requires choosing an uncorrected (feature-wise) threshold.

In CoSMoMVPA we supply an implementation for Threshold-Free Cluster Enhancement (:cite:`SN09`). Optionally this can also be used with first-level null-data (:cite:`SCT13`), although that is not done in this exercise.


Exercises
+++++++++
Exercise with fMRI data
+++++++++++++++++++++++


Note: Since data in the present dataset is not in MNI space, we cannot do group analysis. This exercise therefore considers data from one subject using data in ten chunks corresponding to ten runs.
Note: Since data in the ak6 dataset is not in MNI space, we cannot do group analysis. This exercise therefore considers data from one subject using data in ten chunks corresponding to ten runs.
However, the appraoch can also be used (in exactly the same way) for group analysis, where each chunk corresponds to one subject.

In this exercise, load data from ten runs from the ``ak6`` datasets. Compute, for each chunk, the effect of stimulus presentation versus baseline. Then, using :ref:`cosmo_cluster_neighborhood` and :ref:`cosmo_montecarlo_cluster_stat`, compute a TFCE zscore map corrected for multiple comparisons.
@@ -35,4 +36,15 @@ Template: :ref:`run_multiple_comparison_correction_skl`

Check your answers here: :ref:`run_multiple_comparison_correction` / :pb:`multiple_comparison_correction`


Exercise with time-generalization data
++++++++++++++++++++++++++++++++++++++
In the MEG datasets coming with CoSMoMVPA there is only data from one participant. For this exercise we split the data in ten parts and treat each as coming from one pseudo-participant.

After splitting the data in ten parts, run the time generalization measure on each part, looking for pairs of time points from which the data shows category information. After the necessary transpositions run cosmo_montecarlo_cluster_stat to obtain a group map corrected for multiple comparisons.

Template: :ref:`run_meeg_time_generalization_mcc_skl`

Check your answers here: :ref:`run_meeg_time_generalization_mcc` / :pb:`run_meeg_time_generalization_mcc`

.. include:: links.txt
@@ -0,0 +1,222 @@
%% MEEG time generalization multiple comparison correction
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm where a participant saw
% images of six object categories.
%
% The code presented here can be adapted for other MEEG analyses, but
% there please note:
% * the current examples do not perform baseline corrections or signal
% normalizations, which may reduce discriminatory power.
%
% Note: running this code requires FieldTrip.
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #


%% get timelock data in CoSMoMVPA format

% set configuration
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'meg_obj6');

% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);

% reset citation list
cosmo_check_external('-tic');

% load preprocessed data
data_fn=fullfile(data_path,'meg_obj6_s00.mat');
data_tl=load(data_fn);

% convert to cosmomvpa struct and show the dataset
ds=cosmo_meeg_dataset(data_tl);
cosmo_disp(ds);

% set the targets (trial condition)
ds.sa.targets=ds.sa.trialinfo(:,1); % 6 categories

% set the chunks (independent measurements)
% all trials are here considered to be independent
nsamples=size(ds.samples,1);
ds.sa.chunks=(1:nsamples)';

% in addition give a label to each trial
index2label={'body','car','face','flower','insect','scene'};
ds.sa.labels=cellfun(@(x)index2label(x),num2cell(ds.sa.targets));

% just to check everything is ok
cosmo_check_dataset(ds);

%% Select subset of sensors and time points

% Select posterior gradiometers
sensor_posterior_planar={'MEG1632', 'MEG1642', 'MEG1732', 'MEG1842', ...
'MEG1912', 'MEG1922', 'MEG1942', 'MEG2232', ...
'MEG2312', 'MEG2322', 'MEG2342', 'MEG2432', ...
'MEG2442', 'MEG2512', 'MEG2532',...
'MEG1633', 'MEG1643', 'MEG1733', 'MEG1843', ...
'MEG1913', 'MEG1923', 'MEG1943', 'MEG2233', ...
'MEG2313', 'MEG2323', 'MEG2343', 'MEG2433', ...
'MEG2443', 'MEG2513', 'MEG2533'};

msk=cosmo_dim_match(ds,'chan',sensor_posterior_planar,...
'time',@(t)t>=0 & t<=.3);

ds_sel=cosmo_slice(ds,msk,2);
ds_sel=cosmo_dim_prune(ds_sel);

%%
% subsample time dimension to speed up the analysis
subsample_time_factor=3; % take every 3rd time point

% take every subsample_time_factor-th
% hint: use ds.fa.time to find the desired features, then use cosmo_slice
% and cosmo_dim_prune.
% >@@>
msk = mod(ds_sel.fa.time,subsample_time_factor)==1;
ds_sel=cosmo_slice(ds_sel,msk,2);
ds_sel=cosmo_dim_prune(ds_sel);
% <@@<

% to illustrate group analysis, we use data from a single participant
% and divide it in ten parts. Each part represents a pseudo-participant.
n_pseudo_participants = 10;
ds_sel.sa.subject_id=cosmo_chunkize(ds_sel,n_pseudo_participants);
ds_cell=cosmo_split(ds_sel, 'subject_id');

%%
% apply the cosmo_dim_generalization_measure to the data from each
% pseudo-participant
group_cell=cell(n_pseudo_participants,1);

for k=1:n_pseudo_participants
ds_subj=ds_cell{k};
ds_subj.sa=rmfield(ds_subj.sa,'subject_id');

ds_subj=cosmo_balance_dataset(ds_subj);
ds_subj.sa.chunks=cosmo_chunkize(ds_subj,2);
ds_subj_tr=cosmo_dim_transpose(ds_subj,'time',1);

% use a custom measure that computes a one-way ANOVA F-value and
% then converts this to a z-score
measure=@(d,opt)cosmo_stat(d,'F','z');

ds_time_gen=cosmo_dim_generalization_measure(ds_subj_tr,...
'measure',@cosmo_correlation_measure,...
'dimension','time');


group_cell{k} = ds_time_gen;
end

%%
% show an element of group_cell. What is the size of .samples?
cosmo_disp(group_cell{1});

%%
% To do group analysis, the above format will not work.
% We want a dataset ds_group with size n_pseudo_participants x NF
% where NF is the number of features.

% allocate a cell group_cell_tr with the same size of group_cell
% >@@>
group_cell_tr=cell(size(group_cell));
% <@@<

for k=1:numel(group_cell)
% take data from the k-th participant and store
% in a varibale ds_time_gen
ds_time_gen=group_cell{k};

% change 'train_time' and 'test_time' from being sample dimensions
% to become feature dimensions.
% Hint: use cosmo_dim_transpose.
% >@@>
ds_time_gen_tr=cosmo_dim_transpose(ds_time_gen,...
{'train_time','test_time'},2);
% <@@<

% set chunks and targets for a one-sample t-test against zero,
% so that across participants: all targets have the same value, and
% all chunks have different values.
% >@@>
ds_time_gen_tr.sa.chunks=k;
ds_time_gen_tr.sa.targets=1;
% <@@<

% store ds_time_gen_tr as the k-th element in group_cell_tr
% >@@>
group_cell_tr{k}=ds_time_gen_tr;
% <@@<
end

% show an element of group_cell_tr. What is the size of .samples?
cosmo_disp(group_cell_tr{1});

%%
% stack the elements in group_cell_tr into a dataset ds_group
% >@@>
ds_group=cosmo_stack(group_cell_tr);
% <@@<

%%
% define a clustering neighborhood and store the result in a struct
% called nbrhood
% Hint: use cosmo_cluster_neighborhood
nbrhood=cosmo_cluster_neighborhood(ds_group);

%%

% run multiple comparison correction using cosmo_montecarlo_cluster_stat
% with 1000 iterations, for a t-test against h0_mean=0.
opt=struct();
opt.niter=1000;
opt.h0_mean=0;
ds_tfce=cosmo_montecarlo_cluster_stat(ds_group,nbrhood,opt);

%%
% >@@>
ds_group=cosmo_stack(group_cell_tr);
% <@@<

% extract the values from the ds_tfce, using cosmo_unflatten.
% store the array, and the dimension labels and values, into variables
% arr, dim_labels, and dim_values.
% Hint: use cosmo_unflatten.
% >@@>
[arr, dim_labels, dim_values]=cosmo_unflatten(ds_tfce);
% <@@<

% Reshape to arr to be 2-dimensional and store the result in arr_2d,
% then visualize the array
% >@@>
arr_2d=squeeze(arr);
% <@@<
clim=[-1,1]*max(abs(arr_2d(:)));
imagesc(arr_2d,clim);

% add axis labels
nticks=5;

ytick=round(linspace(1, numel(dim_values{1}), nticks));
ylabel(strrep(dim_labels{1},'_',' '));
set(gca,'Ytick',ytick,'YTickLabel',dim_values{1}(ytick));

xtick=round(linspace(1, numel(dim_values{2}), nticks));
xlabel(strrep(dim_labels{2},'_',' '));
set(gca,'Xtick',xtick,'XTickLabel',dim_values{2}(xtick));
colorbar();

% bonus: add markers indicating significance
% >@@>
z_min=1.96;
[i,j]=find(abs(arr_2d)>z_min);
hold on;
scatter(j,i,'o','k');
hold off;
% <@@<

0 comments on commit a9cfb10

Please sign in to comment.
You can’t perform that action at this time.