%% 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 variable 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 <<< %%%%
