run meeg time generalization mcc

run_up:

run_meeg_time_generalization_mcc

%% 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 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.
    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;