%% 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');
