Demonstrations - full listings¶
demo_fmri_correlation_searchlight¶
%% Demo: fMRI searchlights with split-half correlations, classifier, and representational similarity analysis
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% - 'ak6' is based on the following work (please cite if you use it):
% 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();
ak6_study_path=fullfile(config.tutorial_data_path,'ak6');
% show readme information
readme_fn=fullfile(ak6_study_path,'README');
cosmo_type(readme_fn);
% reset citation list
cosmo_check_external('-tic');
% set result directory
output_path=config.output_data_path;
%% Example: split-half correlation measure (Haxby 2001-style)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This example uses the 'ak6' dataset
% define data filenames & load data from even and odd runs
data_path=fullfile(ak6_study_path,'s01'); % data from subject s01
mask_fn=fullfile(data_path, 'brain_mask.nii'); % whole brain mask
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});
% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_odd_even);
% Use cosmo_correlation_measure.
% This measure returns, by default, a split-half correlation measure
% based on the difference of mean correlations for matching and
% non-matching conditions (a la Haxby 2001).
measure=@cosmo_correlation_measure;
% define spherical neighborhood with radius of 3 voxels
radius=3; % voxels
nbrhood=cosmo_spherical_neighborhood(ds_odd_even,'radius',radius);
% Run the searchlight with a 3 voxel radius
corr_results=cosmo_searchlight(ds_odd_even,nbrhood,measure);
% print output
fprintf('Dataset output:\n');
cosmo_disp(corr_results);
% Plot the output
cosmo_plot_slices(corr_results);
% Define output location
output_fn=fullfile(output_path,'corr_searchlight.nii');
% Store results to disc
cosmo_map2fmri(corr_results, output_fn);
% Show citation information
cosmo_check_external('-cite');
demo_fmri_distatis¶
%% Demo: DISTATIS
%
% 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.
%
% This example shows the application of DISTATIS, which tries to find an
% optimal 'compromise' dissimilarity matrix across a set of observations
% (participants)
%
% Reference:
% - Abdi, H., Valentin, D., O?Toole, A. J., & Edelman, B. (2005).
% DISTATIS: The analysis of multiple distance matrices. In
% Proceedings of the IEEE Computer Society: International conference
% on computer vision and pattern recognition, San Diego, CA, USA
% (pp. 42?47).
%
% # 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');
%% Preprocessing for DISTATIS: RSM analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subject_ids={'s01','s02','s03','s04','s05','s06','s07','s08'};
nsubjects=numel(subject_ids);
mask_label='vt_mask';
ds_rsms=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']); % vt mask
% Use odd runs only
data_fn=fullfile(data_path,'glm_T_stats_odd.nii');
ds=cosmo_fmri_dataset(data_fn,'mask',mask_fn,...
'targets',1:6,'chunks',1);
ds_rsm=cosmo_dissimilarity_matrix_measure(ds);
% set chunks (one chunk per subject)
ds_rsm.sa.chunks=subject_num*ones(size(ds_rsm.samples,1),1);
ds_rsms{subject_num}=ds_rsm;
end
% combine data from all subjects
all_ds=cosmo_stack(ds_rsms);
%% Run DISTATIS
distatis=cosmo_distatis(all_ds);
%% show comprimise distance matrix
[compromise_matrix,dim_labels,values]=cosmo_unflatten(distatis,1);
labels={'monkey', 'lemur', 'mallard', 'warbler', 'ladybug', 'lunamoth'};
n_labels=numel(labels);
figure();
imagesc(compromise_matrix)
title('DSM');
set(gca,'YTick',1:n_labels,'YTickLabel',labels);
set(gca,'XTick',1:n_labels,'XTickLabel',labels);
ylabel(dim_labels{1});
xlabel(dim_labels{2});
colorbar
% skip if stats toolbox is not present
if cosmo_check_external('@stats',false)
figure();
hclus = linkage(compromise_matrix);
dendrogram(hclus,'labels',labels,'orientation','left');
title('dendrogram');
figure();
F = cmdscale(squareform(compromise_matrix));
text(F(:,1), F(:,2), labels);
title('2D MDS plot');
mx = max(abs(F(:)));
xlim([-mx mx]); ylim([-mx mx]);
end
%% show citation information
cosmo_check_external('-cite');
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');
demo_fmri_searchlight_lda¶
%% Demo: fMRI searchlights LDA classifier
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
% A participant made finger pressed with the index and middle finger of
% the right hand during 4 runs in an fMRI study. Each run was divided in
% 4 blocks with presses of each finger and analyzed with the GLM,
% resulting in 2*4*4=32 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();
digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);
output_path=config.output_data_path;
% reset citation list
cosmo_check_external('-tic');
%% LDA classifier searchlight analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This analysis identified brain regions where the categories can be
% distinguished using an odd-even partitioning scheme and a Linear
% Discriminant Analysis (LDA) classifier.
data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig.HEAD');
mask_fn=fullfile(data_path,'brain_mask+orig.HEAD');
targets=repmat(1:2,1,16)'; % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/8)+1; % run labels: 1 1 1 1 1 1 1 1 2 2 ... 4 4
ds_per_run = cosmo_fmri_dataset(data_fn, 'mask', mask_fn,...
'targets',targets,'chunks',chunks);
% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_per_run);
% Use the cosmo_cross_validation_measure and set its parameters
% (classifier and partitions) in a measure_args struct.
measure = @cosmo_crossvalidation_measure;
measure_args = struct();
% Define which classifier to use, using a function handle.
% Alternatives are @cosmo_classify_{svm,matlabsvm,libsvm,nn,naive_bayes}
measure_args.classifier = @cosmo_classify_lda;
% Set partition scheme. odd_even is fast; for publication-quality analysis
% nfold_partitioner is recommended.
% Alternatives are:
% - cosmo_nfold_partitioner (take-one-chunk-out crossvalidation)
% - cosmo_nchoosek_partitioner (take-K-chunks-out " ").
measure_args.partitions = cosmo_oddeven_partitioner(ds_per_run);
% print measure and arguments
fprintf('Searchlight measure:\n');
cosmo_disp(measure);
fprintf('Searchlight measure arguments:\n');
cosmo_disp(measure_args);
% Define a neighborhood with approximately 100 voxels in each searchlight.
nvoxels_per_searchlight=100;
nbrhood=cosmo_spherical_neighborhood(ds_per_run,...
'count',nvoxels_per_searchlight);
% Run the searchlight
lda_results = cosmo_searchlight(ds_per_run,nbrhood,measure,measure_args);
% print output dataset
fprintf('Dataset output:\n');
cosmo_disp(lda_results);
% Plot the output
cosmo_plot_slices(lda_results);
% Define output location
output_fn=fullfile(output_path,'lda_searchlight+orig');
% Store results to disc
cosmo_map2fmri(lda_results, output_fn);
% Show citation information
cosmo_check_external('-cite');
demo_fmri_searchlight_naive_bayes¶
%% Demo: fMRI searchlight with naive bayes classifier
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
% A participant made finger pressed with the index and middle finger of
% the right hand during 4 runs in an fMRI study. Each run was divided in
% 4 blocks with presses of each finger and analyzed with the GLM,
% resulting in 2*4*4=32 t-values
%
% This example uses the cosmo_naive_bayes_classifier_searchlight, which is
% a fast alternative to using the regular searchlight with a
% crossvalidation measure and a classifier
%
% # 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();
digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);
output_path=config.output_data_path;
% reset citation list
cosmo_check_external('-tic');
%% LDA classifier searchlight analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This analysis identified brain regions where the categories can be
% distinguished using an odd-even partitioning scheme and a Naive Bayes
% classifier.
data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig.HEAD');
mask_fn=fullfile(data_path,'brain_mask+orig.HEAD');
targets=repmat(1:2,1,16)'; % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/8)+1; % run labels: 1 1 1 1 1 1 1 1 2 2 ... 4 4
ds_per_run = cosmo_fmri_dataset(data_fn, 'mask', mask_fn,...
'targets',targets,'chunks',chunks);
% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_per_run);
% set parameters for naive bayes searchlight (partitions) in a
% measure_args struct.
measure_args = struct();
% Set partition scheme. odd_even is fast; for publication-quality analysis
% nfold_partitioner is recommended.
% Alternatives are:
% - cosmo_nfold_partitioner (take-one-chunk-out crossvalidation)
% - cosmo_nchoosek_partitioner (take-K-chunks-out " ").
measure_args.partitions = cosmo_oddeven_partitioner(ds_per_run);
% print measure and arguments
fprintf('Searchlight measure arguments:\n');
cosmo_disp(measure_args);
% Define a neighborhood with approximately 100 voxels in each searchlight.
nvoxels_per_searchlight=100;
nbrhood=cosmo_spherical_neighborhood(ds_per_run,...
'count',nvoxels_per_searchlight);
%% Run the searchlight
%
nb_results = cosmo_naive_bayes_classifier_searchlight(ds_per_run,...
nbrhood,measure_args);
%% Show results
% print output dataset
fprintf('Dataset output:\n');
cosmo_disp(nb_results);
% Plot the output
cosmo_plot_slices(nb_results);
% Define output location
output_fn=fullfile(output_path,'naive_bayes_searchlight+orig');
% Store results to disc
cosmo_map2fmri(nb_results, output_fn);
% Show citation information
cosmo_check_external('-cite');
demo_fmri_searchlight_rsm¶
%% Demo: fMRI searchlights with representational similarity analysis
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% - 'ak6' is based on the following work (please cite if you use it):
% 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.
%
% The example shows a searchlight analysis matching local neural similarity
% patterns to three different target similarity matrices
%
% # 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();
ak6_study_path=fullfile(config.tutorial_data_path,'ak6');
% show readme information
readme_fn=fullfile(ak6_study_path,'README');
cosmo_type(readme_fn);
output_path=config.output_data_path;
% reset citation list
cosmo_check_external('-tic');
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This example uses the 'ak6' dataset
% In this example only one sample (response estimate) per condition (class)
% per feature (voxel) is used. Here this is done using t-stats from odd
% runs. One could also use output from a GLM based on an entire
% scanning session experiment.
%
% define data filenames & load data from even and odd runs
data_path=fullfile(ak6_study_path,'s01'); % data from subject s01
mask_fn=fullfile(data_path, 'brain_mask.nii'); % whole brain mask
data_fn=fullfile(data_path,'glm_T_stats_odd.nii');
ds=cosmo_fmri_dataset(data_fn,'mask',mask_fn,...
'targets',1:6,'chunks',1);
%% Set animal species & class
ds.sa.labels={'monkey','lemur','mallard','warbler','ladybug','lunamoth'}';
ds.sa.animal_class=[1 1 2 2 3 3]';
% simple sanity check to ensure all attributes are set properly
cosmo_check_dataset(ds);
% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds);
%% Define feature neighorhoods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For the searchlight, define neighborhood for each feature (voxel).
nvoxels_per_searchlight=100;
% The neighborhood defined here is used three times (one for each target
% similarity matrix), so it is not recomputed for every searchlight call.
fprintf('Defining neighborhood for each feature\n');
nbrhood=cosmo_spherical_neighborhood(ds,'count',nvoxels_per_searchlight);
% print neighborhood
fprintf('Searchlight neighborhood definition:\n');
cosmo_disp(nbrhood);
%% Simple RSM searchlight
nsamples=size(ds.samples,1);
target_dsm=zeros(nsamples);
% define 'simple' target structure where primates (monkey, lemur),
% birds (mallard, warbler) and insects (ladybug, lunamoth) are the same
% (distance=0), all other pairs are equally dissimilar (distance=1).
% Given the ds.sa.targets assignment, pairs (1,2), (3,4) and (5,6) have
% distance 0, all others distance 1.
animal_class=ds.sa.animal_class;
for row=1:nsamples
for col=1:nsamples
same_animal_class=animal_class(row)==animal_class(col);
if same_animal_class
target_dsm(row,col)=0;
else
target_dsm(row,col)=1;
end
end
end
fprintf('Using the following target dsm\n');
disp(target_dsm);
imagesc(target_dsm)
set(gca,'XTick',1:nsamples,'XTickLabel',ds.sa.labels,...
'YTick',1:nsamples,'YTickLabel',ds.sa.labels)
% set measure
measure=@cosmo_target_dsm_corr_measure;
measure_args=struct();
measure_args.target_dsm=target_dsm;
% print measure and arguments
fprintf('Searchlight measure:\n');
cosmo_disp(measure);
fprintf('Searchlight measure arguments:\n');
cosmo_disp(measure_args);
% run searchlight
ds_rsm_binary=cosmo_searchlight(ds,nbrhood,measure,measure_args);
% Note: when these results are used for group analysis across multiple
% participants, it may be good to Fisher-transform the correlation values,
% so that they are more normally distributed. This can be done by:
%
% ds_rsm_binary.samples=atanh(ds_rsm_binary.samples);
% show results
cosmo_plot_slices(ds_rsm_binary);
% store results
output_fn=fullfile(output_path,'rsm_binary.nii');
cosmo_map2fmri(ds_rsm_binary,output_fn);
%% Using another RSM
% This example is very similar to the previous example.
% - This example uses a different target representational similarity
% matrix. The code below allows for identifying regions that show a
% linear dissimilarity across animal class, with primates of distance 1
% from birdsand distance 2 from insects, and insects distance 1 from
% birds. Graphically:
%
% +------------------+------------------+
% primates birds insects
% {monkey,lemur} {mallard,warbler} {ladybug,lunamoth}
%
% - It uses a Spearman rather than Pearson correlation measure to match
% the neural similarity to the target similarity measure
animal_class=ds.sa.animal_class;
% compute absolute difference between each pair of samples
target_dsm=abs(bsxfun(@minus,animal_class,animal_class'));
fprintf('Using the following target dsm\n');
disp(target_dsm);
imagesc(target_dsm)
set(gca,'XTick',1:nsamples,'XTickLabel',ds.sa.labels,...
'YTick',1:nsamples,'YTickLabel',ds.sa.labels)
% set measure
measure=@cosmo_target_dsm_corr_measure;
measure_args=struct();
measure_args.target_dsm=target_dsm;
% 'Spearman' requires matlab with stats toolbox; if not present use
% 'Pearson'
if cosmo_check_external('@stats',false)
measure_args.type='Spearman';
else
measure_args.type='Pearson';
fprintf('Matlab stats toolbox not present, using %s correlation\n',...
measure_args.type)
end
% run searchlight
ds_rsm_linear=cosmo_searchlight(ds,nbrhood,measure,measure_args);
% Note: when these results are used for group analysis across multiple
% participants, it may be good to Fisher-transform the correlation values,
% so that they are more normally distributed. This can be done by:
%
% ds_rsm_linear.samples=atanh(ds_rsm_linear.samples);
% show results
cosmo_plot_slices(ds_rsm_linear);
% store results
output_fn=fullfile(output_path,'rsm_linear.nii');
cosmo_map2fmri(ds_rsm_linear,output_fn);
%% Using a behavioural RSM
% This example is very similar to the one above, but now the target
% similarity structure is based on behavioural similarity ratings
% load behavioural similarity matrix from disc
behav_model_fn=fullfile(ak6_study_path,'models','behav_sim.mat');
behav_model=importdata(behav_model_fn);
target_dsm=behav_model;
fprintf('Using the following target dsm\n');
disp(target_dsm);
imagesc(target_dsm)
set(gca,'XTick',1:nsamples,'XTickLabel',ds.sa.labels,...
'YTick',1:nsamples,'YTickLabel',ds.sa.labels)
% set measure
measure=@cosmo_target_dsm_corr_measure;
measure_args=struct();
measure_args.target_dsm=target_dsm;
% run searchlight
ds_rsm_behav=cosmo_searchlight(ds,nbrhood,measure,measure_args);
% Note: when these results are used for group analysis across multiple
% participants, it may be good to Fisher-transform the correlation values,
% so that they are more normally distributed. This can be done by:
%
% ds_rsm_behav.samples=atanh(ds_rsm_behav.samples);
% show results
cosmo_plot_slices(ds_rsm_behav);
% store results
output_fn=fullfile(output_path,'rsm_behav.nii');
cosmo_map2fmri(ds_rsm_behav,output_fn);
% Show citation information
cosmo_check_external('-cite');
demo_meeg_dataset_operations¶
%% MEEG dataset operations
%
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
% * assignment of targets (labels of conditions) is based here on
% stimulation periods versus pre-stimulation periods. In typical
% analyses the targets should be based on different trial conditions, for
% example as set a FieldTrip .trialinfo field.
% * assignment of chunks (parts of the data that are assumed to be
% independent) is done randomly. In typical analyses they can be based on
% different runs, with one chunk value per run.
% * the time window used for analyses is rather small. This means that in
% particular for time-freq analysis a lot of data is missing, especially
% for early and late timepoints in the lower frequency bands. For typical
% analyses it may be preferred to use a wider time window.
% * 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_20hz');
% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);
% reset citation list
cosmo_check_external('-tic');
% load data
data_fn=fullfile(data_path,'subj102_B01_20Hz_timelock.mat');
data_tl=load(data_fn);
% convert to cosmomvpa struct
ds_tl=cosmo_meeg_dataset(data_tl);
% set the target (trial condition)
ds_tl.sa.targets=ds_tl.sa.trialinfo(:,1); % 1=pre, 2=post
% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tl.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tl.sa.targets));
% just to check everything is ok
cosmo_check_dataset(ds_tl);
%% Print some information about dataset
% print the size. Note that the number of features is the number of time
% points times the number of channels.
[ns,nf]=size(ds_tl.samples);
fprintf('There are %d samples and %d features\n', ns, nf);
fdim_labels=ds_tl.a.fdim.labels;
fdim_values=ds_tl.a.fdim.values;
ndim=numel(fdim_labels);
for dim=1:ndim
dim_label=fdim_labels{dim};
nvalues=numel(fdim_values{dim});
fprintf('Dimension %d (%s) has %d values\n', dim, dim_label, nvalues);
end
%% Simple data manipulations
% select five trials from post-stim period
post_mask=cosmo_match(ds_tl.sa.labels,'post');
ds_tl_post=cosmo_slice(ds_tl,post_mask);
select_trials=2:2:18;
% the last argument indicates to slice along the first dimension.
% The value 1 is the default and can be omitted
ds_trial=cosmo_slice(ds_tl_post, select_trials, 1);
% select a single channel (multiple channels are also possible but
% makes plotting the samples directly a bit more tricky)
chan={'MEG1042'};
chan_msk=cosmo_dim_match(ds_trial,'chan', chan);
% slice again, this time along the feature dimension (dimension 2),
% to get a new dataset with just one channel
ds_trial_chan=cosmo_slice(ds_trial, chan_msk, 2);
% plot five individual trials
figure();
subplot(2,1,1);
plot(data_tl.time, ds_trial_chan.samples');
title(sprintf('%d trials of sensor %s', numel(select_trials), chan{1}));
% use original full data to select data in just one channel
chan_msk=cosmo_dim_match(ds_tl,'chan', chan);
ds_chan=cosmo_slice(ds_tl, chan_msk, 2);
ds_mean=cosmo_fx(ds_chan,@(x)mean(x,1),'targets');
% plot signal time course averaged over all samples
% note: the 'pre' data has been shifted in time to match the 'post' data
subplot(2,1,2);
[time_dim, time_index]=cosmo_dim_find(ds_mean,'time');
assert(time_dim==2); % it is a feature dimension
plot(ds_mean.a.fdim.values{time_index},ds_mean.samples);
title(sprintf('average of sensor %s', chan{1}));
legend(ds_mean.sa.labels);
% Show citation information
cosmo_check_external('-cite');
demo_meeg_timefreq_searchlight¶
%% MEEG time-frequency searchlight
%
% This example shows MVPA analyses performed on MEEG data, using a
% searchlight across the time, frequency and channel dimensions
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
% * assignment of targets (labels of conditions) is based here on
% stimulation periods versus pre-stimulation periods. In typical
% analyses the targets should be based on different trial conditions, for
% example as set a FieldTrip .trialinfo field.
% * assignment of chunks (parts of the data that are assumed to be
% independent) is based on a trial-by-trial basis. For cross-validation,
% the number of chunks is reduced to two to speed up the analysis.
% * the time window used for analyses is rather small. This means that in
% particular for time-freq analysis a lot of data is missing, especially
% for early and late timepoints in the lower frequency bands. For typical
% analyses it may be preferred to use a wider time window.
% * 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_20hz');
% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);
% reset citation list
cosmo_check_external('-tic');
% load data
data_fn=fullfile(data_path,'subj102_B01_20Hz_timefreq.mat');
data_tf=load(data_fn);
% convert to cosmomvpa struct
ds_tf=cosmo_meeg_dataset(data_tf);
% set the target (trial condition)
ds_tf.sa.targets=ds_tf.sa.trialinfo(:,1); % 1=pre, 2=post
% set the chunks (independent measurements)
% in this dataset, the first half of the samples (in order)
% are the post-trials;
% the second half the pre-trials
ds_tf.sa.chunks=[(1:145) (1:145)]';
% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tf.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tf.sa.targets));
% just to check everything is ok
cosmo_check_dataset(ds_tf);
% get rid of features with at least one NaN value across samples
fa_nan_mask=sum(isnan(ds_tf.samples),1)>0;
fprintf('%d / %d features have NaN\n', ...
sum(fa_nan_mask), numel(fa_nan_mask));
ds_tf=cosmo_slice(ds_tf, ~fa_nan_mask, 2);
%% set MVPA parameters
fprintf('The input has feature dimensions %s\n', ...
cosmo_strjoin(ds_tf.a.fdim.labels,', '));
% set chunks
% again for speed just two chunks
% (targets were already set above)
nchunks=2;
ds_tf.sa.chunks=cosmo_chunkize(ds_tf, nchunks);
% define neighborhood parameters for each dimension
% channel neighborhood uses meg_combined_from_planar, which means that the
% input are planar channels but the output has combined-planar channels.
% to use the magnetometers, use 'meg_axial'
chan_type='meg_combined_from_planar';
chan_count=10; % use 10 channel locations (relative to the combined
% planar channels)
% as we use meg_combined_from_planar there are
% 20 channels in each searchlight because
% gradiometers are paired
time_radius=2; % 2*2+1=5 time bines
freq_radius=4; % 4*2+1=9 freq bins
% define the neighborhood for each dimensions
chan_nbrhood=cosmo_meeg_chan_neighborhood(ds_tf, 'count', chan_count, ...
'chantype', chan_type);
freq_nbrhood=cosmo_interval_neighborhood(ds_tf,'freq',...
'radius',freq_radius);
time_nbrhood=cosmo_interval_neighborhood(ds_tf,'time',...
'radius',time_radius);
% cross neighborhoods for chan-time-freq searchlight
nbrhood=cosmo_cross_neighborhood(ds_tf,{chan_nbrhood,...
freq_nbrhood,...
time_nbrhood});
% print some info
nbrhood_nfeatures=cellfun(@numel,nbrhood.neighbors);
fprintf('Features have on average %.1f +/- %.1f neighbors\n', ...
mean(nbrhood_nfeatures), std(nbrhood_nfeatures));
% only keep features with at least 10 neighbors
% (some have zero neighbors - in particular, those with low frequencies
% early or late in time)
center_ids=find(nbrhood_nfeatures>10);
% for illustration purposes use the split-half measure because it is
% relatively fast - but clasifiers can also be used
measure=@cosmo_correlation_measure;
% split-half, as there are just two chunks
% (when using a classifier, do not use 'half' but the number of chunks to
% leave out for testing, e.g. 1).
measure_args=struct();
measure_args.partitions=cosmo_nchoosek_partitioner(ds_tf,'half');
%% run searchlight
sl_tf_ds=cosmo_searchlight(ds_tf,nbrhood,measure,measure_args,...
'center_ids',center_ids);
%% visualize results
% deduce layout from output
layout=cosmo_meeg_find_layout(sl_tf_ds);
fprintf('The output uses layout %s\n', layout.name);
% map to FT struct for visualization
sl_tf_ft=cosmo_map2meeg(sl_tf_ds);
% show figure
figure()
cfg = [];
if cosmo_wtf('is_octave')
% GNU Octave does not show data when labels are shown
cfg.interactive='no';
cfg.showlabels='no';
else
% Matlab supports interactive viewing and labels
cfg.interactive = 'yes';
cfg.showlabels = 'yes';
end
cfg.zlim='maxabs';
cfg.layout = layout;
ft_multiplotTFR(cfg, sl_tf_ft);
% Show citation information
cosmo_check_external('-cite');
demo_meeg_timelock_searchlight¶
%% MEEG time-lock searchlight
%
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
% Using a time-channel neighborhood, a searchlight map is computed
% indicating in time and space (channel) the pre-stimulation period can be
% distinguished from the peri/post stimulation period.
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
% * assignment of targets (labels of conditions) is based here on
% stimulation periods versus pre-stimulation periods. In typical
% analyses the targets should be based on different trial conditions, for
% example as set a FieldTrip .trialinfo field.
% * assignment of chunks (parts of the data that are assumed to be
% independent) is based on a trial-by-trial basis. For cross-validation,
% the number of chunks is reduced to two to speed up the analysis.
% * 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_20hz');
% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);
% reset citation list
cosmo_check_external('-tic');
% load data
data_fn=fullfile(data_path,'subj102_B01_20Hz_timelock.mat');
data_tl=load(data_fn);
% convert to cosmomvpa struct
ds_tl=cosmo_meeg_dataset(data_tl);
% set the target (trial condition)
ds_tl.sa.targets=ds_tl.sa.trialinfo(:,1); % 1=pre, 2=post
% set the chunks (independent measurements)
% in this dataset, the first half of the samples (in order)
% are the post-trials;
% the second half the pre-trials
ds_tl.sa.chunks=[(1:145) (1:145)]';
% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tl.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tl.sa.targets));
% just to check everything is ok
cosmo_check_dataset(ds_tl);
% get rid of features with at least one NaN value across samples
fa_nan_mask=sum(isnan(ds_tl.samples),1)>0;
fprintf('%d / %d features have NaN\n', ...
sum(fa_nan_mask), numel(fa_nan_mask));
ds_tl=cosmo_slice(ds_tl, ~fa_nan_mask, 2);
%% set MVPA parameters
fprintf('The input has feature dimensions %s\n', ...
cosmo_strjoin(ds_tl.a.fdim.labels,', '));
% set chunks
% again for speed just two chunks
% (targets were already set above)
nchunks=2;
ds_tl.sa.chunks=cosmo_chunkize(ds_tl, nchunks);
% define neighborhood parameters for each dimension
% channel neighborhood uses meg_combined_from_planar, which means that the
% input are planar channels but the output has combined-planar channels.
% to use the magnetometers, use 'meg_axial'
chan_type='meg_combined_from_planar';
chan_count=10; % use 10 channel locations (relative to the combined
% planar channels)
% as we use meg_combined_from_planar there are
% 20 channels in each searchlight because
% gradiometers are paired
time_radius=2; % 2*2+1=5 time bines
% define the neighborhood for each dimensions
chan_nbrhood=cosmo_meeg_chan_neighborhood(ds_tl, 'count', chan_count, ...
'chantype', chan_type);
time_nbrhood=cosmo_interval_neighborhood(ds_tl,'time',...
'radius',time_radius);
% cross neighborhoods for chan-time searchlight
nbrhood=cosmo_cross_neighborhood(ds_tl,{chan_nbrhood,...
time_nbrhood});
% print some info
nbrhood_nfeatures=cellfun(@numel,nbrhood.neighbors);
fprintf('Features have on average %.1f +/- %.1f neighbors\n', ...
mean(nbrhood_nfeatures), std(nbrhood_nfeatures));
% only keep features with at least 10 neighbors
center_ids=find(nbrhood_nfeatures>10);
% for illustration purposes use the split-half measure because it is
% relatively fast - but clasifiers can also be used
measure=@cosmo_correlation_measure;
% split-half, as there are just two chunks
% (when using a classifier, do not use 'half' but the number of chunks to
% leave out for testing, e.g. 1).
measure_args=struct();
measure_args.partitions=cosmo_nchoosek_partitioner(ds_tl,'half');
%% run searchlight
sl_tl_ds=cosmo_searchlight(ds_tl,nbrhood,measure,measure_args,...
'center_ids',center_ids);
%% visualize timeseries results
% deduce layout from output
layout=cosmo_meeg_find_layout(sl_tl_ds);
fprintf('The output uses layout %s\n', layout.name);
% map to FT struct for visualization
sl_tl_ft=cosmo_map2meeg(sl_tl_ds);
figure();
cfg = [];
cfg.interactive = 'yes';
cfg.zlim=[-1 1];
cfg.layout = layout;
% show figure with plots for each sensor
ft_multiplotER(cfg, sl_tl_ft);
%% visualize topology results
% show figure with topology for 0 to 600ms after stimulus onset in bins of
% 100 ms
figure();
cfg.xlim=-0.1:0.1:0.5;
ft_topoplotER(cfg, sl_tl_ft);
%% Show citation information
cosmo_check_external('-cite');
demo_meeg_timeseries_classification¶
%% MEEG timeseries classification
%
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
%
% * assignment of targets (labels of conditions) is based here on
% stimulation periods versus pre-stimulation periods. In typical
% analyses the targets should be based on different trial conditions, for
% example as set a FieldTrip .trialinfo field.
% * assignment of chunks (parts of the data that are assumed to be
% independent) is based on a trial-by-trial basis. For cross-validation,
% the number of chunks is reduced to four to speed up the analysis.
% * the time window used for analyses is rather small. This means that in
% particular for time-freq analysis a lot of data is missing, especially
% for early and late timepoints in the lower frequency bands. For typical
% analyses it may be preferred to use a wider time window.
% * 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_20hz');
% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);
% reset citation list
cosmo_check_external('-tic');
% load data
data_fn=fullfile(data_path,'subj102_B01_20Hz_timelock.mat');
data_tl=load(data_fn);
% convert to cosmomvpa struct
ds_tl=cosmo_meeg_dataset(data_tl);
% set the target (trial condition)
ds_tl.sa.targets=ds_tl.sa.trialinfo(:,1); % 1=pre, 2=post
% set the chunks (independent measurements)
% in this dataset, the first half of the samples (in order)
% are the post-trials;
% the second half the pre-trials
ds_tl.sa.chunks=[(1:145) (1:145)]';
% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tl.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tl.sa.targets));
% just to check everything is ok
cosmo_check_dataset(ds_tl);
%% Prepare MVPA
% reset chunks: use four chunks
nchunks=4;
ds_tl.sa.chunks=cosmo_chunkize(ds_tl,nchunks);
% do a take-one-fold out cross validation.
% except when using a splithalf correlation measure it is important that
% the partitions are *balanced*, i.e. each target (or class) is presented
% equally often in each chunk
partitions=cosmo_nchoosek_partitioner(ds_tl,1);
partitions=cosmo_balance_partitions(partitions, ds_tl);
npartitions=numel(partitions);
fprintf('There are %d partitions\n', numel(partitions.train_indices));
fprintf('# train samples:%s\n', sprintf(' %d', cellfun(@numel, ...
partitions.train_indices)));
fprintf('# test samples:%s\n', sprintf(' %d', cellfun(@numel, ...
partitions.test_indices)));
%% Run time-series searchlight on magneto- and gradio-meters seperately
% try two different classification approaches:
% 1) without averaging the samples in the train set
% 2) by averaging 5 samples at the time in the train set, and re-using
% every sample 3 times.
% (Note: As of July 2015, there is no clear indication in the literature
% which approach is 'better'. These two approaches are used here to
% illustrate how they can be used with a searchlight).
average_train_args_cell={{},... % { 1
{'average_train_count',5,... % { 2
'average_train_resamplings',3}}; % { 2
n_average_train_args=numel(average_train_args_cell);
% compute and plot accuracies for magnetometers and gradiometers separately
chantypes={'meg_axial','meg_planar'};
% in the time searchlight analysis, select the time-point itself, the two
% timepoints after it, and the two timepoints before it
time_radius=2;
nchantypes=numel(chantypes);
ds_chantypes=cosmo_meeg_chantype(ds_tl);
plot_counter=0;
for j=1:n_average_train_args
% define the measure and its argument.
% here a simple naive baysian classifier is used.
% Alternative are @cosmo_classify_{svm,nn,lda}.
measure=@cosmo_crossvalidation_measure;
measure_args=struct();
measure_args.classifier=@cosmo_classify_naive_bayes;
measure_args.partitions=partitions;
% add the options to average samples to the measure arguments.
% (if no averaging is desired, this step can be left out.)
average_train_args=average_train_args_cell{j};
measure_args=cosmo_structjoin(measure_args, average_train_args);
for k=1:nchantypes
parent_type=chantypes{k};
% find feature indices of channels matching the parent_type
chantype_idxs=find(cosmo_match(ds_chantypes,parent_type));
% define mask with channels matching those feature indices
chan_msk=cosmo_match(ds_tl.fa.chan,chantype_idxs);
% slice the dataset to select only the channels matching the channel
% types
ds_tl_sel=cosmo_dim_slice(ds_tl, chan_msk, 2);
ds_tl_sel=cosmo_dim_prune(ds_tl_sel); % remove non-indexed channels
% define neighborhood over time; for each time point the time
% point itself is included, as well as the two time points before
% and the two time points after it
nbrhood=cosmo_interval_neighborhood(ds_tl_sel,'time',...
'radius',time_radius);
% run the searchlight using the measure, measure arguments, and
% neighborhood defined above.
% Note that while the input has both 'chan' and 'time' as feature
% dimensions, the output only has 'time' as the feature dimension
sl_map=cosmo_searchlight(ds_tl_sel,nbrhood,measure,measure_args);
fprintf('The output has feature dimensions: %s\n', ...
cosmo_strjoin(sl_map.a.fdim.labels,', '));
plot_counter=plot_counter+1;
subplot(n_average_train_args,nchantypes,plot_counter);
time_values=sl_map.a.fdim.values{1}; % first dim (channels got nuked)
plot(time_values,sl_map.samples);
ylim([.4 .8])
xlim([min(time_values),max(time_values)]);
ylabel('classification accuracy (chance=.5)');
xlabel('time');
if isempty(average_train_args)
postfix=' no averaging';
else
postfix=' with averaging';
end
descr=sprintf('%s - %s', strrep(parent_type,'_',' '), postfix);
title(descr);
end
end
% Show citation information
cosmo_check_external('-cite');
demo_meeg_timeseries_generalization¶
%% MEEG time-by-time transfer classification
%
% This example shows MVPA generalization across time. The input is a
% time-locked dataset (chan x time), the output is either:
%
% * (train_time x test_time) indicating for each combination of time points
% in the training and test set how similar the patterns are
% * (train_time x test_time x chan), which uses a searchlight and
% indicates for each combination of time points in the training and
% test set how similar the patterns are in the neighborhood of each
% channel.
%
% Results can be visualized in FieldTrip.
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
%
% * assignment of targets (labels of conditions) is based here on
% stimulation periods versus pre-stimulation periods. In typical
% analyses the targets should be based on different trial conditions, for
% example as set a FieldTrip .trialinfo field.
% * the time window used for analyses is rather small. This means that in
% particular for time-freq analysis a lot of data is missing, especially
% for early and late timepoints in the lower frequency bands. For typical
% analyses it may be preferred to use a wider time window.
% * 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. #
% set configuration
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'meg_20hz');
% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);
% reset citation list
cosmo_check_external('-tic');
% load data
data_fn=fullfile(data_path,'subj102_B01_20Hz_timelock.mat');
data_tl=load(data_fn);
% convert to cosmomvpa struct
ds_tl=cosmo_meeg_dataset(data_tl);
% set the target (trial condition)
ds_tl.sa.targets=ds_tl.sa.trialinfo(:,1); % 1=pre, 2=peri/post
% set the chunks (independent measurements)
% in this dataset, the first half of the samples (in order)
% are the post-trials;
% the second half the pre-trials
ds_tl.sa.chunks=[(1:145) (1:145)]';
% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tl.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tl.sa.targets));
% just to check everything is ok
cosmo_check_dataset(ds_tl);
%% reduce the number of time points
% to allow this example to run relatively fast, only use the 'even'
% timepoints and only use timepoints between 0 and 300 ms relative to
% stimulus onset.
% for a publication-quality analysis this step could be omitted.
msk=cosmo_dim_match(ds_tl,'time',@(x) 0<x &x<.3) & ...
mod(ds_tl.fa.time,2)==0;
ds=cosmo_slice(ds_tl,msk,2);
ds=cosmo_dim_prune(ds);
%% Set up chunks
% use half of the data for training and half for testing. Here both chunks
% come from the same session, but it is also possible to use different
% sessions for training and testing.
%
% For example, the train data can contain responses to auditory stimuli
% with the words 'face' and 'house', while the test session measures
% responses to visual stimuli. In that case the chunks should be assigned
% 'manually' (without cosmo_chunkize)
nchunks=2; % two chunks are required for this analysis
ds.sa.chunks=cosmo_chunkize(ds,nchunks);
%% time-by-time generalization on magneto- and gradio-meters seperately
% compute and plot accuracies for magnetometers and gradiometers separately
chan_types={'meg_axial','meg_planar'};
nchan_types=numel(chan_types);
% set arguments for the cosmo_dim_generalization_measure
measure_args=struct();
% the cosmo_dim_generalization_measure requires that another
% measure (here: the crossvalidation measure) is specified. The
% specified measure is applied for each combination of time points
measure_args.measure=@cosmo_crossvalidation_measure;
% When used ordinary, the cosmo_crossvalidation_measure itself
% requires two arguments:
% - classifier (here: LDA)
% - partitions
% However, because the cosmo_dim_generalization_measure defines
% the partitions itself, they are not set here.
measure_args.classifier=@cosmo_classify_lda;
% define the dimension over which generalization takes place
measure_args.dimension='time';
% define the radius for the time dimension. Here not just a single
% time-point is used, but also the time-point before it and the time-point
% after it.
measure_args.radius=1;
% try both channel types (axial and planar)
for k=1:nchan_types
use_chan_type=chan_types{k};
% get layout for mag or planar labels, and select only those channels
ds_chan_types=cosmo_meeg_chantype(ds);
ds_chan_idxs=find(cosmo_match(ds_chan_types,use_chan_type));
feature_msk=cosmo_match(ds.fa.chan,ds_chan_idxs);
ds_sel=cosmo_slice(ds, feature_msk, 2);
% make 'time' a sample dimension
% (this necessary for cosmo_dim_generalization_measure)
ds_time=cosmo_dim_transpose(ds_sel,'time',1);
fprintf('The input for channel type %s is:\n', use_chan_type)
cosmo_disp(ds_time);
% run transfer across time with the searchlight neighborhood
%cdt_ds=cosmo_cartesian_dim_transfer(ds_time,'time',measure,...
% 'args',measure_args);
cdt_ds=cosmo_dim_generalization_measure(ds_time,measure_args);
fprintf('The output is:\n')
cosmo_disp(cdt_ds);
% unflatten the data to get train_time x test_time matrix
[data, labels, values]=cosmo_unflatten(cdt_ds,1);
% show the results
figure()
imagesc(data, [.3 .7]);
title(sprintf('classification accuracy for %s', use_chan_type));
colorbar();
nticks=5;
ytick=round(linspace(1, numel(values{1}), nticks));
ylabel(strrep(labels{1},'_',' '));
set(gca,'Ytick',ytick,'YTickLabel',values{1}(ytick));
xtick=round(linspace(1, numel(values{2}), nticks));
xlabel(strrep(labels{2},'_',' '));
set(gca,'Xtick',xtick,'XTickLabel',values{2}(xtick));
colorbar();
end
%% time-by-time generalization using a searchlight with correlation measure
% how many channel locations in each searchlight
nchannel_locations_searchlight=10;
% make 'time' a sample dimension (necessary for cartesian_dim_transfer)
ds_time=cosmo_dim_transpose(ds_sel,'time',1);
% use planar channels as input; output is
% like combine-planar
use_chan_type='meg_combined_from_planar';
% define searchlight neighborhood
nbrhood=cosmo_meeg_chan_neighborhood(ds_time,...
'count',nchannel_locations_searchlight,...
'chantype',use_chan_type);
fprintf('The input is:\n')
cosmo_disp(ds_time);
fprintf('The neighborhood is:\n');
cosmo_disp(nbrhood);
% set the measure to be the dim generalization measure
measure=@cosmo_dim_generalization_measure;
% define the arguments for the measure. One of them is called 'measure',
% because that measure is applied to combinations of samples in the
% train and test set
measure_args=struct();
measure_args.measure=@cosmo_correlation_measure;
measure_args.radius=1;
measure_args.dimension='time';
% run transfer across time with the searchlight neighborhood
cdt_ds=cosmo_searchlight(ds_time,nbrhood,measure,measure_args);
fprintf('The output is:\n')
cosmo_disp(cdt_ds)
% move {train,test}_time from being sample dimensions to feature
% dimensions, so they can be mapped to a fieldtrip struct
cdt_tf_ds=cosmo_dim_transpose(cdt_ds,{'train_time','test_time'},2);
fprintf('The output after transposing is:\n')
cosmo_disp(cdt_ds)
% trick fieldtrip into thinking this is a time-freq-chan dataset, by
% renaming train_time and test_time to freq and time, respectively
cdt_tf_ds=cosmo_dim_rename(cdt_tf_ds,'train_time','freq');
cdt_tf_ds=cosmo_dim_rename(cdt_tf_ds,'test_time','time');
%% Show correlation time-by-time generalization searchlight results
% determine layout
layout=cosmo_meeg_find_layout(cdt_tf_ds);
% convert to fieldtrip format
ft=cosmo_map2meeg(cdt_tf_ds);
% show train_time x test_time x chan figure
%
% * train_time is on the vertical ('freq') axis
% * test_time on the horizontal ('time') axis)
figure();
cfg = [];
cfg.interactive = 'yes';
cfg.showlabels = 'yes';
cfg.zlim=[-1 1];
cfg.layout = layout;
ft_multiplotTFR(cfg, ft);
% show train_time x test_time figure for selected channels
figure();
cfg = [];
cfg.interactive = 'yes';
cfg.showlabels = 'yes';
cfg.zlim=[-1 1];
cfg.layout = layout;
cfg.channel={'MEG0322+0323', 'MEG0332+0333', 'MEG0412+0413', ...
'MEG0422+0423', 'MEG0632+0633', 'MEG0642+0643'};
ft_singleplotTFR(cfg,ft);
% show channel topology for timewindow over train_time x test_time
figure();
cfg = [];
cfg.interactive = 'yes';
cfg.xlim=[.2 .25];
cfg.ylim=[.2 .25];
cfg.zlim=[-1 1];
cfg.layout = layout;
ft_topoplotTFR(cfg, ft);
%% show citation information
cosmo_check_external('-cite');
demo_surface_searchlight_lda¶
%% Demo: fMRI surface-based searchlights with LDA classifier
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
% A participant made finger pressed with the index and middle finger of
% the right hand during 4 runs in an fMRI study. Each run was divided in
% 4 blocks with presses of each finger and analyzed with the GLM,
% resulting in 2*4*4=32 t-values
%
% The example shows four possible searchlight analyses, covering typical
% use cases:
% 1) single or twin surfaces
% + Caret and BrainVoyager use a single surface; a parameter 'offsets'
% is used to define which voxels are considered to be in the "grey
% matter" (but this may be not so precise)
% + FreeSurfer uses twin surfaces (pial and white), and voxels in
% between or on them are considered to be in the grey matter
% 2) lower resolution output map
% + in the canonical surface-based searchlight, each node on the input
% surface(s) is assigned a measure value (accuracy, in this example)
% + it is also possible to have output in a lower resolution version
% than the input surfaces; this reduces both the execution time
% (a Good Thing) and spatial precision (a Bad Thing). Two approaches
% are illustrated to use a lower resolution surface for output:
% 1) from MapIcosahedron, with a lower value for the number of
% divisions of the triangles
% 2) using a surface subsampling approach, implemented by
% surfing_subsample_surface
%
% In all cases a searchlight is run with a 100 voxel searchlight, using a
% disc for which the metric radius varies from node to node. For a fixed
% metric radius of the disc, use a positive value for 'radius' below.
% Distances are measured across the cortical surface using a geodesic
% distance metric.
%
% This example requires the surfing toolbox, github.com/nno/surfing
%
% This example may take quite some time to run. For faster execution, set
% ld=16 (instead of ld=64) below
%
% If you use this code for a publication, please cite:
% Oosterhof, N.N., Wiestler, T, Downing, P.E., & Diedrichsen, J. (2011)
% A comparison of volume-based and surface-based information mapping.
% Neuroimage. DOI:10.1016/j.neuroimage.2010.04.270
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
%% Check externals
cosmo_check_external('surfing');
%% 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();
digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);
output_path=config.output_data_path;
% reset citation list
cosmo_check_external('-tic');
% resolution parameter for input surfaces
% 64 is for high-quality results; use 16 for fast execution
ld=64;
% Twin surface case (FS)
pial_fn=fullfile(digit_study_path,...
sprintf('ico%d_mh.pial_al.asc', ld));
white_fn=fullfile(digit_study_path,...
sprintf('ico%d_mh.smoothwm_al.asc', ld));
% Single surface case (Caret/BV)
intermediate_fn=fullfile(digit_study_path,...
sprintf('ico%d_mh.intermediate_al.asc', ld));
% Used for visualization
inflated_fn=fullfile(digit_study_path,...
sprintf('ico%d_mh.inflated_alCoMmedial.asc', ld));
%%
% Set parameters
% Searchlight radius: select 100 features in each searchlight
% (to use a fixed radius of 8mm, use:
% cosmo_surficial_neighborhood(...,'radius',8)
% in the code below)
feature_count=100;
% Single surface case: select voxels that are 3 mm or closer to the surface
% on the white-matter side, up to voxels that are 2 mm from the surface on
% the pial matter side
single_surf_offsets=[-2 3];
% Single surface case: number of iterations to downsample surface
lowres_output_onesurf_niter=10;
% Twin surface case: number of linear divisions from MapIcosahedron
lowres_output_twosurf_icold=16;
lowres_intermediate_fn=fullfile(digit_study_path,...
sprintf('ico%d_mh.intermediate_al.asc',...
lowres_output_twosurf_icold));
% Use the cosmo_cross_validation_measure and set its parameters
% (classifier and partitions) in a measure_args struct.
measure = @cosmo_crossvalidation_measure;
measure_args = struct();
% Define which classifier to use, using a function handle.
% Alternatives are @cosmo_classify_{svm,nn,naive_bayes}
measure_args.classifier = @cosmo_classify_lda;
%% Load functional data
data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig');
targets=repmat(1:2,1,16)'; % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/8)+1; % run labels: 1 1 1 1 1 1 1 1 2 2 ... 4 4
ds = cosmo_fmri_dataset(data_fn,'targets',targets,'chunks',chunks);
% remove zero elements
zero_msk=all(ds.samples==0,1);
ds = cosmo_slice(ds, ~zero_msk, 2);
fprintf('Dataset has %d samples and %d features\n', size(ds.samples));
% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds);
%% Set partition scheme. odd_even is fast; for publication-quality analysis
% nfold_partitioner is recommended.
% Alternatives are:
% - cosmo_nfold_partitioner (take-one-chunk-out crossvalidation)
% - cosmo_nchoosek_partitioner (take-K-chunks-out " ").
measure_args.partitions = cosmo_oddeven_partitioner(ds);
% print measure and arguments
fprintf('Searchlight measure:\n');
cosmo_disp(measure);
fprintf('Searchlight measure arguments:\n');
cosmo_disp(measure_args);
%% Read inflated surface
[v_inf,f_inf]=surfing_read(inflated_fn);
fprintf('The inflated surface has %d vertices, %d faces\n',...
size(v_inf,1), size(f_inf,1))
%% Run four types of searchlights
for one_surf=[true,false]
if one_surf
desc='1surf';
else
desc='2surfs';
end
for lowres_output=[false,true]
if lowres_output
desc=sprintf('%s_lowres', desc);
end
fprintf('\n\n *** Starting analysis with %s *** \n\n\n', desc)
% define searchlight surface paramters for each type of analysis
if one_surf && lowres_output
% single surface (Caret/BV) with lower-res output
surf_def={intermediate_fn,single_surf_offsets,...
lowres_output_onesurf_niter};
elseif one_surf && ~lowres_output
% single surface (Caret/BV) with original-res output
surf_def={intermediate_fn,single_surf_offsets};
elseif ~one_surf && lowres_output
% single surface (FS) with lower-res output
surf_def={white_fn,pial_fn,lowres_intermediate_fn};
elseif ~one_surf && ~lowres_output
% single surface (FS) with original-res output
surf_def={white_fn,pial_fn};
else
assert(false); % should never get here
end
% Define the feature neighborhood for each node on the surface
% - nbrhood has the neighborhood information
% - vo and fo are vertices and faces of the output surface
% - out2in is the mapping from output to input surface
fprintf('Defining neighborhood with %s\n', desc);
[nbrhood,vo,fo,out2in]=cosmo_surficial_neighborhood(ds,surf_def,...
'count',feature_count);
% print neighborhood
fprintf('Searchlight neighborhood definition:\n');
cosmo_disp(nbrhood);
fprintf('The output surface has %d vertices, %d nodes\n',...
size(vo,1), size(fo,1));
% Run the searchlight
lda_results = cosmo_searchlight(ds,nbrhood,measure,measure_args);
% print searchlight output
fprintf('Dataset output:\n');
cosmo_disp(lda_results);
% Apply the node mapping from the surifical neighborhood
% to the high-res inflated surface.
% (This example shows how such a mapping can be applied to new
% surfaces)
if lowres_output
v_inf_out=v_inf(out2in,:);
f_inf_out=fo;
else
v_inf_out=v_inf;
f_inf_out=f_inf;
end
% visualize the surfaces, if the afni matlab toolbox is present
if cosmo_check_external('afni',false)
nvertices=size(v_inf_out,1);
opt=struct();
for show_edge=[false, true]
opt.ShowEdge=show_edge;
if show_edge
t='with edges';
else
t='without edges';
end
header=strrep([desc ' ' t],'_',' ');
DispIVSurf(vo,fo,1:nvertices,lda_results.samples',0,opt);
title(sprintf('Original %s', header));
DispIVSurf(v_inf_out,f_inf_out,1:nvertices,...
lda_results.samples',0,opt);
title(sprintf('Inflated %s', header));
end
else
fprintf('skip surface display; no afni matlab toolbox\n');
end
if lowres_output && one_surf
% in this example only this case a new surface was generated.
% To aid visualization using external tools, store it to disc.
% The surface is stored in ASCII, GIFTI and BV SRF
% formats, if the required externals are present
surf_output_fn=fullfile(output_path,['inflated_' desc]);
% AFNI/SUMA ASC
surfing_write([surf_output_fn '.asc'],v_inf_out,f_inf_out);
% BV SRF
if cosmo_check_external('neuroelf',false)
surfing_write([surf_output_fn '.srf'],v_inf_out,f_inf_out);
end
% GIFTI
if cosmo_check_external('gifti',false)
surfing_write([surf_output_fn '.gii'],v_inf_out,f_inf_out);
end
end
% store searchlight results
data_output_fn=fullfile(output_path,['lda_' desc]);
if cosmo_check_external('afni',false)
cosmo_map2surface(lda_results, [data_output_fn '.niml.dset']);
end
if cosmo_check_external('neuroelf',false)
cosmo_map2surface(lda_results, [data_output_fn '.smp']);
end
% store voxel counts (how often each voxel is in a neighborhood)
% take a random sample (the first one) from the input dataset
% and count how often each voxel was selected.
% If everything works, then voxels in the grey matter have high
% voxel counts but voxels outside it low or zero counts.
% Thus, this can be used as a sanity check that can be visualized
% easily.
vox_count_ds=cosmo_slice(ds,1);
vox_count_ds.samples(:)=0;
ncenters=numel(nbrhood.neighbors);
for k=1:ncenters
idxs=nbrhood.neighbors{k}; % feature indices in neigborhood
vox_count_ds.samples(idxs)=vox_count_ds.samples(idxs)+1;
end
vox_count_output_fn=fullfile(output_path,['vox_count_' desc]);
% store voxel count results
cosmo_map2fmri(vox_count_ds, [vox_count_output_fn '.nii']);
if cosmo_check_external('afni',false)
cosmo_map2fmri(vox_count_ds, [vox_count_output_fn '+orig']);
end
end
end
% Show citation information
cosmo_check_external('-cite');
demo_surface_tfce¶
%% Demo: Threshold-Free Cluster Enhancement (TFCE) on surface dataset
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
% A participant made finger pressed with the index and middle finger of
% the right hand during 4 runs in an fMRI study. Each run was divided in
% 4 blocks with presses of each finger and analyzed with the GLM,
% resulting in 2*4*4=32 t-values
%
% This example illustrates the use of Threshold-Free Cluster Enhancement
% with a permutation test to correct for multiple comparisons.
%
% TFCE reference: Stephen M. Smith, Thomas E. Nichols, Threshold-free
% cluster enhancement: Addressing problems of smoothing, threshold
% dependence and localisation in cluster inference, NeuroImage, Volume 44,
% Issue 1, 1 January 2009, Pages 83-98.
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
%% Check externals
cosmo_check_external({'surfing','afni'});
%% 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();
digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);
output_path=config.output_data_path;
% resolution parameter for input surfaces
% 64 is for high-quality results; use 16 for fast execution
ld=16;
% reset citation list
cosmo_check_external('-tic');
% load single surface
intermediate_fn=fullfile(digit_study_path,...
sprintf('ico%d_mh.intermediate_al.asc', ld));
[vertices,faces]=surfing_read(intermediate_fn);
%% Load functional data
data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig');
targets=repmat(1:2,1,16)'; % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/4)+1; % half-run: 1 1 1 1 2 2 2 2 3 3 ... 8 8
vol_ds = cosmo_fmri_dataset(data_fn,'targets',targets,'chunks',chunks);
%% Map univariate response data to surface
% this measure averages the data near each node to get a surface dataset
radius=0;
surf_band_range=[-2 2]; % get voxel data within 2mm from surface
surf_def={vertices,faces,[-2 2]};
nbrhood=cosmo_surficial_neighborhood(vol_ds,surf_def,'radius',radius);
measure=@(x,opt) cosmo_structjoin('samples',mean(x.samples,2),'sa',x.sa);
surf_ds=cosmo_searchlight(vol_ds,nbrhood,measure);
fprintf('Univariate surface data:\n');
cosmo_disp(surf_ds);
%% Average data in each chunk
% for this example only consider the samples in the first condition
% (targets==1), and average the samples in each chunk
%
% for group analysis: set chunks to (1:nsubj)', assuming each sample is
% data from a single participant
surf_ds=cosmo_slice(surf_ds,surf_ds.sa.targets==1);
surf_ds=cosmo_average_samples(surf_ds);
fn_surf_ds=fullfile(output_path, 'digit_target1.niml.dset');
% save to disc
cosmo_map2surface(surf_ds, fn_surf_ds);
fprintf('Input data saved to %s\n', fn_surf_ds);
%% Run Threshold-Free Cluster Enhancement (TFCE)
% All data is prepared; surf_ds has 8 samples and 5124 nodes. We want to
% see if there are clusters that show a significant difference from zero in
% their response. Thus, .sa.targets is set to all ones (the same
% condition), whereas .sa.chunks is set to (1:8)', indicating that all
% samples are assumed to be independent.
%
% (While this is a within-subject analysis, exactly the same logic can be
% applied to a group-level analysis)
% define neighborhood for each feature
% (cosmo_cluster_neighborhood can be used also for meeg or volumetric
% fmri datasets)
cluster_nbrhood=cosmo_cluster_neighborhood(surf_ds,...
'vertices',vertices,'faces',faces);
fprintf('Cluster neighborhood:\n');
cosmo_disp(cluster_nbrhood);
opt=struct();
% number of null iterations. for publication-quality, use >=1000;
% 10000 is even better
opt.niter=250;
% in this case we run a one-sample test against a mean of 0, and it is
% necessary to specify the mean under the null hypothesis
% (when testing classification accuracies, h0_mean should be set to chance
% level, assuming a balanced crossvalidation scheme was used)
opt.h0_mean=0;
% this example uses the data itself (with resampling) to obtain cluster
% statistcs under the null hypothesis. This is (in this case) somewhat
% conservative due to how the resampling is performed.
% Alternatively, and for better estimates (at the cost of computational
% cost), one can generate a set of (say, 50) datasets using permuted data
% e.g. using cosmo_randomize_targets), put them in a cell and provide
% them as the null argument.
opt.null=[];
fprintf('Running multiple-comparison correction with these options:\n');
cosmo_disp(opt);
% Run TFCE-based cluster correction for multiple comparisons.
% The output has z-scores for each node indicating the probablity to find
% the same, or higher, TFCE value under the null hypothesis
tfce_ds=cosmo_montecarlo_cluster_stat(surf_ds,cluster_nbrhood,opt);
%% Show results
fprintf('TFCE z-score dataset\n');
cosmo_disp(tfce_ds);
nfeatures=size(tfce_ds.samples,2);
percentiles=(1:nfeatures)/nfeatures*100;
plot(percentiles,sort(tfce_ds.samples))
title('sorted TFCE z-scores');
xlabel('feature percentile');
ylabel('z-score');
nvertices=size(vertices,1);
disp_opt=struct();
disp_opt.DataRange=[-2 2];
DispIVSurf(vertices,faces,1:nvertices,tfce_ds.samples',0,disp_opt);
% store results
fn_tfce_ds=fullfile(output_path, 'digit_target1_tfce.niml.dset');
cosmo_map2surface(tfce_ds, fn_tfce_ds);
surf_fn=fullfile(output_path, 'digit_intermediate.asc');
surfing_write(surf_fn,vertices,faces);
% show citation information
cosmo_check_external('-cite');