demo surface tfce

demo_up:

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, ...
                                       'metric', 'dijkstra');

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