version of 2014-2-18

README.text

@@ -0,0 +1,101 @@
+% MVPA by cross-validated MANOVA
+% Carsten Allefeld
+% v1, 2013-12-20
+
+This text documents the Matlab implementation for SPM of the method introduced
+by Carsten Allefeld and John-Dylan Haynes, "Searchlight-based multi-voxel
+pattern analysis of fMRI by cross-validated MANOVA", NeuroImage, 89:345–357,
+2014.
+
+# prerequisites
+
+The cross-validated MANOVA is based on a ("first-level") multivariate General
+Linear Model. This model has to be specified and estimated in SPM before using
+these functions.
+
+Estimation of the model is necessary in order to access SPM's estimates of
+various fMRI data properties, especially the temporal correlation of the
+residuals. The functions use the generated `SPM.mat` file and the data files
+referenced therein, as well as the mask image (`mask.hdr` and `mask.img`). The
+ beta images generated during estimation can be deleted.
+
+The functions may work with other versions, but have been specifically
+developed and tested with SPM8 and Matlab 7.11–7.14.
+
+
+# interface
+
+The main interface is given by the function
+
+>       cvManovaSearchlight(dirName, slRadius, Cs, permute)
+
+which computes the cross-validated MANOVA on a searchlight. `dirName` is the
+name of the directory where the `SPM.mat` file is located, `slRadius` is the
+radius of the searchlight in voxels, and `Cs` is a cell array whose elements
+are contrast matrices. `permute` specifies whether permutation values should
+be computed and defaults to `false`.
+
+Simple ("t-like") contrasts are specified as a *column vector*, complex
+("F-like") contrasts as a matrix of several columns. Please note that this is
+the transpose of the format used in the SPM user interface.
+
+The rows of a contrast matrix correspond to the model regressors for each
+session *separately*, i.e. other than in SPM the contrast should not be
+explicitly replicated for several sessions. Instead, the program performs the
+replication internally, assuming that (at least the leading) regressors for
+each session model the same effects. If there are fewer rows in a contrast
+matrix than there are regressors for a session, the matrix is zero-padded.
+
+The searchlight radius $r$ is interpreted such that every voxel is included
+for which the distance from the center voxel is *smaller than or equal* to the
+radius. This means that $r = 0$ leads to a searchlight size of 1 voxel,
+$r = 1$ to 7 voxels, $r = 2$ to 33 voxels, and so on. This definition may
+differ from the one used in other implementations of MVPA algorithms and in
+publications. Note that it is possible to use fractional values for $r$.
+
+The result of the analysis are estimates of a multivariate measure of effect
+size, the pattern discriminability $D$, which is intended as a drop-in
+replacement for the conventional measure of classification accuracy.
+Statistical parametric maps of $D$ are written to images with filenames
+of the form
+
+>       spmD_C####_P####.nii
+
+enumerating all contrasts and permutations, in the same directory as the
+`SPM.mat` file. Additionally, an image of the numbers of voxels contained in
+each searchlight is written to `VPSL.nii`.
+
+To ease the specification of contrasts, the utility function `contrasts` can
+be used to generate contrast matrices for all main effects and interactions of
+a factorial design, in a form suitable for use with `cvManovaSearchlight`.
+
+For example, `Cs = contrasts([2 3])` results in
+
+>       Cs = { [ 1  1  1 -1 -1 -1]'
+>              [ 1 -1 -0  1 -1 -0 ; -0  1 -1 -0  1 -1]'
+>              [ 1 -1 -0 -1  1  0 ; -0  1 -1  0 -1  1]' };
+
+For further documentation, please refer to the help texts in the `m`-files.
+
+
+# remarks
+
+— The functions are optimized for the computation of several contrasts (and
+permutations) in one run. One call of `cvManovaSearchlight` with several
+contrasts will take substantially less time than several calls for each
+contrast separately.
+
+— The function reads the complete data set into memory. The analysis should
+therefore be run on a computer with a sufficient amount of main memory, and
+using other memory-intensive programs at the same time should be avoided.
+
+— The estimation of $D$ is based on the GLM residuals and therefore depends on
+a properly specified model. That means that all effects that are known to
+systematically occur should be included in the model. Because sub-effects
+can be selected through the mechanism of constrasts, it is neither necessary
+nor advisable to use different GLMs as the basis of different MVPA analyses.
+
+— Depending on the data set, it may be possible to perform the analysis for
+searchlight radii of up to 5 or 6. However, a large radius leads to very long
+computation times as well as decreased numerical precision. The recommended
+radius is 3, resulting in a searchlight of 123 voxels.

contrasts.m

@@ -0,0 +1,60 @@
+function [cMatrix, cName] = contrasts(fLevel, fName)
+
+% generate contrasts for main effects and interactions of a factorial design
+% 
+% [cMatrix, cName] = contrasts(fLevel)
+% [cMatrix, cName] = contrasts(fLevel, fName)
+%
+% fLevel:   number of levels for each factor (row vector)
+% fName:    names of factors (cell array)
+% cMatrix:  contrast matrices (cell array)
+% cName:    contrast names (cell array)
+%
+% the first factor is the one being enumerated slowest.
+% contrasts are not orthonormalized!
+%
+% Copyright (C) 2013 Carsten Allefeld
+%
+% This program is free software: you can redistribute it and/or modify it
+% under the terms of the GNU General Public License as published by the
+% Free Software Foundation, either version 3 of the License, or (at your
+% option) any later version. This program is distributed in the hope that
+% it will be useful, but without any warranty; without even the implied
+% warranty of merchantability or fitness for a particular purpose. See the
+% GNU General Public License <http://www.gnu.org/licenses/> for more details.
+
+% number of factors
+nf = size(fLevel, 2);
+
+if nargin < 2
+    % generate generic names for factors
+    fName = num2cell(char((1 : nf) + '@'));
+end
+
+% generate sorted list of contrast signatures
+nc = 2 ^ nf - 1;
+cs = bitget(ones(nf, 1) * (1 : nc), (1 : nf)' * ones(1, nc))';
+[~, ind] = sort(sum(cs, 2));
+cs = cs(ind, :);
+cs = logical(cs);
+
+% compute contrast elements
+e = cell(2, nf);
+for fi = 1 : nf
+    e{1, fi} = ones(fLevel(fi), 1);
+    e{2, fi} = -(diff(eye(fLevel(fi)))');
+end
+
+% compute contrasts
+cMatrix = cell(nf, 1);
+cName = cell(nf, 1);
+for ci = 1 : nc
+    % contrast matrix
+    cMatrix{ci} = 1;
+    for fi = 1 : nf
+        cMatrix{ci} = kron(cMatrix{ci}, e{cs(ci, fi) + 1, fi});
+    end
+    % contrast name
+    l = sprintf('×%s', fName{cs(ci, :)});
+    cName{ci} = l(2:end);
+end

cvManovaSearchlight.m

@@ -0,0 +1,119 @@
+function cvManovaSearchlight(dirName, slRadius, Cs, permute)
+
+% cross-validated MANOVA on searchlight
+%
+% cvManovaSearchlight(dirName, slRadius, Cs, permute = false)
+%
+% dirName:  directory where the SPM.mat file referring to an estimated
+%           model is located
+% slRadius: radius of the searchlight sphere, in voxels
+% Cs:       cell array of contrast matrices
+% permute:  whether to compute permutation values of the test statistic
+%
+% Output files are written to the same directory:
+% spmD_C####_P####.nii:
+%           images of the pattern discriminability D
+%           contrast and permutation are identified by numbers
+% VPSL.img  an image of the number of voxels within each searchlight
+% cms.mat   a record of the analysis parameters
+%
+% Copyright (C) 2013 Carsten Allefeld
+%
+% This program is free software: you can redistribute it and/or modify it
+% under the terms of the GNU General Public License as published by the
+% Free Software Foundation, either version 3 of the License, or (at your
+% option) any later version. This program is distributed in the hope that
+% it will be useful, but without any warranty; without even the implied
+% warranty of merchantability or fitness for a particular purpose. See the
+% GNU General Public License <http://www.gnu.org/licenses/> for more details.
+
+
+fprintf('\n\ncvManovaSearchlight\n\n')
+
+if nargin < 4
+    permute = false;
+end
+
+nContrasts = numel(Cs);
+outpattern = 'spmD_C%04d_P%04d.nii';
+
+if dirName(end) ~= '/', dirName = [dirName '/']; end
+
+% simplify saving images
+wd = cd;
+cd(dirName)
+
+% load data, design matrix etc.
+[Y, X, mask, misc] = loadDataSPM(dirName);
+
+% determine per-run design and data matrices
+nRuns = misc.m;
+Xrun = cell(nRuns, 1);
+Yrun = cell(nRuns, 1);
+% for each run
+for ri = 1 : nRuns
+    Yrun{ri} = Y(misc.sRow{ri}, :);
+    Xrun{ri} = X(misc.sRow{ri}, misc.sCol{ri});
+end
+clear Y X
+
+% check contrasts
+for ci = 1 : nContrasts
+    if size(Cs{ci}, 2) > rank(Cs{ci})
+        error('contrast %d is misspecified!', ci)
+    end
+    for ri = 1 : nRuns
+        if inestimability(Cs{ci}, Xrun{ri}) > 1e-6
+            error('contrast %d is not estimable in run %d!', ci, ri)
+        end
+    end
+end
+
+% precomputation
+fprintf('\nprecomputing GLM runwise\n')
+[XXs, betas, xis] = cvManova_precompute(Xrun, Yrun);
+clear Xrun Yrun
+
+% determine voxels per searchlight, and save as image
+fprintf('\ncomputing voxels per searchlight image\n')
+p = runSearchlight(slRadius, mask, @(vi)(size(vi, 1)));
+saveMRImage(reshape(p, size(mask)), 'VPSL.nii', misc.v2mm, ...
+    'voxels per searchlight');
+
+% error check
+pMax = max(p(:));
+if pMax > 0.9 * (misc.m - 1) * misc.fE
+    error('insufficient amount of data for searchlight size %d!', pMax)
+end
+
+% bias correction factor
+factor = ((misc.m - 1) * misc.fE - p - 1) / ((misc.m - 1) * misc.n);
+clear p
+
+% run searchlight
+fprintf('\ncomputing cross-validated MANOVA on searchlight\n')
+mDl = runSearchlight(slRadius, mask, @cvManova_compute, ...
+    XXs, betas, xis, Cs, permute);
+
+% separate contrast and permutation dimensions
+nContrasts = numel(Cs);
+nPerms = size(mDl, 2) / nContrasts;
+mDl = reshape(mDl, [], nContrasts, nPerms);
+
+% compute the unbiased estimate of the pattern discriminability D
+D = bsxfun(@times, factor, mDl);
+clear mDl
+
+% save results
+for ci = 1 : nContrasts
+    for pi = 1 : nPerms
+        fn = sprintf(outpattern, ci, pi);
+        saveMRImage(reshape(D(:, ci, pi), size(mask)), fn, misc.v2mm, ...
+            'pattern discriminability');
+    end
+end
+
+% analysis parameters
+save cms.mat slRadius Cs permute misc nPerms
+
+cd(wd)