%% 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.
%%%% >>> Your code here <<< %%%%
% 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
%%%% >>> Your code here <<< %%%%
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.
%%%% >>> Your code here <<< %%%%
% 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.
%%%% >>> Your code here <<< %%%%
% store ds_time_gen_tr as the k-th element in group_cell_tr
%%%% >>> Your code here <<< %%%%
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
%%%% >>> Your code here <<< %%%%
%%
% 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);
%%
%%%% >>> Your code here <<< %%%%
% 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.
%%%% >>> Your code here <<< %%%%
% Reshape to arr to be 2-dimensional and store the result in arr_2d,
% then visualize the array
%%%% >>> Your code here <<< %%%%
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
%%%% >>> Your code here <<< %%%%
