demo fmri rois

demo_up:

demo_fmri_rois

%% Demo: fMRI Region-Of-Interest (ROI) analyses
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% It is based on the following work:
% * Connolly et al (2012), Representation of biological classes in the human
%   brain. Journal of Neuroscience, doi 10.1523/JNEUROSCI.5547-11.2012
%
% Six categories (monkey, lemur, mallard, warbler, ladybug, lunamoth)
% during ten runs in an fMRI study. Using the General Linear Model response
% were estimated for each category in each run, resulting in 6*10=60
% t-values.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config = cosmo_config();
study_path = fullfile(config.tutorial_data_path, 'ak6');
output_path = config.output_data_path;

readme_fn = fullfile(study_path, 'README');
cosmo_type(readme_fn);

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

%% Example 1: split-half correlation measure (Haxby 2001-style)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subject_id = 's01';
mask_label = 'ev_mask';
data_path = fullfile(study_path, subject_id); % data from subject s01

% Define data locations and load data from even and odd runs
mask_fn = fullfile(data_path, [mask_label '.nii']); % whole brain

data_odd_fn = fullfile(data_path, 'glm_T_stats_odd.nii');
ds_odd = cosmo_fmri_dataset(data_odd_fn, 'mask', mask_fn, ...
                            'targets', 1:6, 'chunks', 1);

data_even_fn = fullfile(data_path, 'glm_T_stats_even.nii');
ds_even = cosmo_fmri_dataset(data_even_fn, 'mask', mask_fn, ...
                             'targets', 1:6, 'chunks', 2);

% Combine even and odd runs
ds_odd_even = cosmo_stack({ds_odd, ds_even});

% remove constant features
ds_odd_even = cosmo_remove_useless_data(ds_odd_even);

% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_odd_even);

% compute correlations
ds_corr = cosmo_correlation_measure(ds_odd_even);

% show result
fprintf(['Average correlation difference between matching and '...
         'non-matching categories in %s for %s is %.3f\n'], ...
        mask_label, subject_id, ds_corr.samples);

%% Example 2: split-half correlation measure with group analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subject_ids = {'s01', 's02', 's03', 's04', 's05', 's06', 's07', 's08'};
nsubjects = numel(subject_ids);

mask_label = 'vt_mask';

ds_corrs = cell(nsubjects, 1); % allocate space for output
for subject_num = 1:nsubjects
    subject_id = subject_ids{subject_num};

    % Code from here is pretty much identical to that above >>>

    % set path for this subject
    data_path = fullfile(study_path, subject_id);

    % Define data locations and load data from even and odd runs
    mask_fn = fullfile(data_path, [mask_label '.nii']); % whole brain

    data_odd_fn = fullfile(data_path, 'glm_T_stats_odd.nii');
    ds_odd = cosmo_fmri_dataset(data_odd_fn, 'mask', mask_fn, ...
                                'targets', 1:6, 'chunks', 1);

    data_even_fn = fullfile(data_path, 'glm_T_stats_even.nii');
    ds_even = cosmo_fmri_dataset(data_even_fn, 'mask', mask_fn, ...
                                 'targets', 1:6, 'chunks', 2);

    % Combine even and odd runs
    ds_odd_even = cosmo_stack({ds_odd, ds_even});

    % remove constant features
    ds_odd_even = cosmo_remove_useless_data(ds_odd_even);

    ds_corr = cosmo_correlation_measure(ds_odd_even);

    % <<< identical up to here

    % set targets and chunks for the output, so that cosmo_stat can be used
    % below
    ds_corr.sa.targets = 1;
    ds_corr.sa.chunks = subject_num;

    ds_corrs{subject_num} = ds_corr;
end

% combine the data from all subjects
ds_all = cosmo_stack(ds_corrs);

%%
samples = ds_all.samples; % get the correlations for all subjects

% run one-sample t-test again zero

% Using cosmo_stats
ds_t = cosmo_stat(ds_all, 't');     % t-test against zero
ds_p = cosmo_stat(ds_all, 't', 'p'); % convert to p-value

fprintf(['correlation difference in %s at group level: '...
         '%.3f +/- %.3f, %s=%.3f, p=%.5f (using cosmo_stat)\n'], ...
        mask_label, mean(samples), std(samples), ...
        ds_t.sa.stats{1}, ds_t.samples, ds_p.samples);

% Using matlab's stat toolbox (if present)
if cosmo_check_external('@stats', false)
    [h, p, ci, stats] = ttest(samples);
    fprintf(['Correlation difference in %s at group level: '...
             '%.3f +/- %.3f, t_%d=%.3f, p=%.5f (using matlab stats '...
             'toolbox)\n'], ...
            mask_label, mean(samples), std(samples), stats.df, stats.tstat, p);
else
    fprintf('Matlab stats toolbox not found\n');
end

%% Example 3: comparison of four classifiers in two regions of interest
config = cosmo_config();
data_path = fullfile(config.tutorial_data_path, 'ak6', 's01');
data_fn = fullfile(data_path, 'glm_T_stats_perrun.nii');

% Define classifiers and mask labels
classifiers = {@cosmo_classify_nn, ...
               @cosmo_classify_naive_bayes, ...
               @cosmo_classify_lda};

% Use svm classifiers, if present
svm_name2func = struct();
svm_name2func.matlabsvm = @cosmo_classify_matlabsvm;
svm_name2func.libsvm = @cosmo_classify_libsvm;
svm_name2func.svm = @cosmo_classify_svm;
svm_names = fieldnames(svm_name2func);
for k = 1:numel(svm_names)
    svm_name = svm_names{k};
    if cosmo_check_external(svm_name, false)
        classifiers{end + 1} = svm_name2func.(svm_name);
    else
        warning('Classifier %s skipped because not available', svm_name);
    end
end

mask_labels = {'vt_mask', 'ev_mask'};

%
nclassifiers = numel(classifiers);
nmasks = numel(mask_labels);

labels = {'monkey', 'lemur', 'mallard', 'warbler', 'ladybug', 'lunamoth'};
nlabels = numel(labels);

% little helper function to replace underscores by spaces
underscore2space = @(x) strrep(x, '_', ' ');

for j = 1:nmasks
    mask_label = mask_labels{j};
    mask_fn = fullfile(data_path, [mask_label '.nii']);
    ds = cosmo_fmri_dataset(data_fn, 'mask', mask_fn, ...
                            'targets', repmat(1:6, 1, 10), ...
                            'chunks', floor(((1:60) - 1) / 6) + 1);

    % remove constant features
    ds = cosmo_remove_useless_data(ds);

    % print dataset
    fprintf('Dataset input:\n');
    cosmo_disp(ds);

    % Define partitions
    partitions = cosmo_nfold_partitioner(ds);

    % print dataset
    fprintf('Partitions:\n');
    cosmo_disp(partitions);

    % show result for each classifier
    for k = 1:nclassifiers
        classifier = classifiers{k};

        % get predictions for each fold
        [pred, accuracy] = cosmo_crossvalidate(ds, classifier, partitions);

        % get confusion matrix for each fold
        confusion_matrix_folds = cosmo_confusion_matrix(ds.sa.targets, pred);

        % sum confusion for each ground-truth target and prediction,
        % resulting in an nclasses x nclasses matrix
        confusion_matrix = sum(confusion_matrix_folds, 3);
        figure;
        imagesc(confusion_matrix, [0 10]);
        cfy_label = underscore2space(func2str(classifier));
        title_ = sprintf('%s using %s: accuracy=%.3f', ...
                         underscore2space(mask_label), cfy_label, accuracy);
        title(title_);
        set(gca, 'XTick', 1:nlabels, 'XTickLabel', labels);
        set(gca, 'YTick', 1:nlabels, 'YTickLabel', labels);
        ylabel('target');
        xlabel('predicted');
    end
end

%% Show citation information
cosmo_check_external('-cite');