Giter Site home page Giter Site logo

cvmanova's Introduction

MVPA by cross-validated MANOVA

This is an implementation for Matlab 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

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 errors. The functions use the generated SPM.mat file and the data files referenced therein, as well as the analysis brain mask image. Other files generated during estimation can be deleted.

Searchlight analysis

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 (see below). permute specifies whether permutation values should be computed and defaults to false.

The searchlight radius 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 0 leads to a searchlight size of 1 voxel, 1 to 7 voxels, 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 the searchlight radius. For a table of searchlight radii leading to different searchlight sizes, run slSize.

The result of the analysis are estimates of a measure of multivariate 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.

While D is the main measure of multivariate effect size of cross-validated MANOVA and the statistic whose values should be reported, for the purpose of statistical inference the standardized pattern distinctness may be better suited (see Eq. 17 in the paper). Statistical parametric maps of standardized D are written to images with filenames of the form

spmDs_C####_P####.nii

Additionally, an image of the numbers of voxels contained in each searchlight is written to VPSL.nii, and the analysis parameters are saved to cmsParameters.mat.

Region-of-interest analysis

The paper describes cross-validated MANOVA only for searchlight analyses. This package contains additional code for ROI analyses. The function

[D, p] = cvManovaRegion(dirName, regions, Cs, permute)

performs ROI-based analysis on a set of voxels, specified by the parameter regions in the form of a logical volume or an image filename, or a cell array to process several ROIs at once.

Since the result of ROI analysis is not an image but a scalar, it is directly returned by the function.

Note that the order of parameters of cvManovaRegion has changed with respect to version 2. Please adjust your code accordingly!

Contrasts

In cross-validated MANOVA, effects of interest are specified in the form of contrast vectors or matrices, in the same way as for univariate analysis. 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, but identical to SPM’s internal format.

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. This makes it easy to ignore regressors that may be present in only some sessions or subjects, e.g. modeling error trials, as long as they are trailing.

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]) is equivalent to

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]' };

describing the two main effects and the interaction.

Note that the resulting contrasts may have to be modified in the case of including HRF derivatives or using FIR models, or if a factor is nested in another one.

Remarks

– 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, even if they do not enter the contrast. Because sub-effects can be selected through the mechanism of contrasts, it is neither necessary nor advisable to use different GLMs as the basis of different MVPA analyses of the same data set.

– The fMRI model specification must include the modeling of temporal autocorrelations in order to correctly estimate the pattern distinctness. For this, the option ‘serial correlations’ in SPM has to be kept at the default value AR(1), or changed to the newer FAST.

– 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.

– For computational efficiency, the functions read 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. Peak memory usage is about twice the amount to load the data set, which is (number of in-mask voxels) × (number of scans) × 8 bytes.

cvManovaSearchlight contains a checkpointing mechanism. If the computation is interrupted for reasons other than an internal error and then restarted with the same parameters, it picks up at the last checkpoint. Intermediate results are stored in a file cmsCheckpoint####.mat in the same directory as SPM.mat, where #### is a hash encoding the parameters. The file is deleted after the computation finishes successfully.

Example and test script

The subdirectory cvManovaTest contains a script cvManovaTest that analyses the data of Haxby et al. (2001), both as a test of the implementation and as an example for how to use it.

Regularization

For very large searchlight sizes or for a large ROI, it is possible that there are so many voxels that adequate estimation of the error covariance matrix is no longer possible. The resulting numerical instability of the method may be remedied by using regularization of the matrix. This implementation contains the option to apply shrinkage towards the covariance matrix diagonal, by supplying the shrinkage parameter lambda as an additional parameter to cvManovaSearchlight or cvManovaRegion.

Note however that with regularization, D is no longer an unbiased estimator of the true pattern distinctness (unless it is zero). It is therefore recommended to avoid regularization, and rather reduce the number of voxels. The recommended searchlight radius for cross-validated MANOVA is 3, leading to a searchlight of 123 voxels. If regularization is used, the parameter should be kept very small, e.g. 0.001.

The implementation contains a hard-coded limit on the number of voxels within a searchlight or ROI regardless of regularization, of 90% of the available error degrees of freedom. That is already a rather large threshold, which one should normally not get close to.

Note that previous experimental code to estimate the optimal shrinkage parameter based on the method of Schäfer and Strimmer (2005) in version 2 has been removed, because it proved to be unreliable.

Negative pattern distinctness?

Even though pattern distinctness D is a ratio of explained multivariate variance, or a generalized squared distance, the values provided by cross-validated MANOVA may be negative. That raises the question how such values should be interpreted.

The simple answer is: Negative values do not have an interpretation per se, and they can never be significantly below zero, so there is no problem for reporting.

The longer answer is that the values produced by the algorithm are only estimates of the true pattern distincess, and these estimates randomly vary around the true value because of a finite amount of data. The true pattern distinctness can never be below zero. However, the estimator was designed to be unbiased (correct on average), and that implies that if the true value is zero or close to it, estimates vary around zero, and therefore about half of them have to be below zero.

This has an exact analogue in the case of cross-validated classification accuracy. The true accuracy can never be below chance level, but estimated accuracies can be.

In some cases, the estimated value of pattern distinctness strongly suggests that the true value is below zero, too. This is most likely the result of a violation of the assumption underlying cross-validation, that the different parts of the data (sessions) are generated in exactly the same way. This assumption may not hold if there are unmodelled confounds in the data, or problems with the design itself. Strongly negative estimated values of pattern distinctness therefore suggest that you should re-check your design matrix, or the experimental design itself. Again, the same problem will most likely occur with cross-validated classification accuracy computed from the same data.


Feel free to contact me with questions and comments. Bug reports and feature requests can also be submitted via the GitHub issue tracker.

This software was developed with SPM8 and SPM12 under Matlab 7.11–8.5 (R2010b–R2015a), but later versions should work, too. It is copyrighted © 2013–2018 by Carsten Allefeld and released under the terms of the GNU General Public License, version 3 or later.

This file is part of v3 of cvmanova, see
https://github.com/allefeld/cvmanova/releases

cvmanova's People

Contributors

allefeld avatar remi-gau avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar

cvmanova's Issues

searchlight size error

Hi Carsten

After loading a dataset
100 of 365 volumes loaded
200 of 365 volumes loaded
300 of 365 volumes loaded
365 of 365 volumes loaded
whitening
high-pass-filtering
df: 365 - 11 - 10 = 344 [SPM: trRV = 344 erdf = 344]

I get the error : "Error using cvManovaSearchlight (line 67)
data insufficient for searchlight of size 257!"

I've tried different searchlight sizes without success.

---> if pMax > fEMin * 0.9 % ensures decent numerical precision

seems that the images doesn't meet this criteria....?

thanks and best, mike

Normalization

Hi,

We were wondering to include normalization of the data (Pereira 2009, Example preprocessing section, https://www.sciencedirect.com/science/article/pii/S1053811908012263) right after data loading.

I think the line 36 in cvManovaRegion.m,
[Ys, Xs, ~, misc] = loadDataSPM(dirName, regions);
is the data loading part

I was wondering if I add following lines? If it does not make sense, can you suggest better way?

Yss{1}=transpose(zscore(transpose(Ys{1})));
Yss{2}=transpose(zscore(transpose(Ys{2})));
Yss{3}=transpose(zscore(transpose(Ys{3})));
Ys=Yss; clearvars Yss;

Thank you.

cvManovaCore contrast error

Hi,

I have 16 conditions and I set up 1 main effect and 1 interaction effect contrasts

% set up contrasts
Cs = {};
% 1) main effect of threat
Cs{1} = [ -1 1 -1 1 0 0 0 0 0 0 0 0 0 0 0 0]';
% 2) interaction
Cs{2} = [ -1 1 1 -1 0 0 0 0 0 0 0 0 0 0 0 0]';

and I ran cvManovaSearchlight,

cvManovaSearchlight(dir, 3, Cs)

but I had error message as follows,

_computing cross-validated MANOVA on searchlight
Error using cvManovaCore (line 82)
contrast 1 is not estimable in session 1!

Error in runSearchlight (line 73)
res = fun(nan(0, 1), varargin{:});

Error in cvManovaSearchlight (line 69)
[D, p] = runSearchlight(['cmsCheckpoint' uid
'.mat'], slRadius, mask, ..._

Can you help me to handle with this? Thanks.

cvManovaCore

integrating as much stuff as possible into cvManova_compute

matlabbatch error

When running cvManovaTest, it echos this errors. Thanks very much.

cvManovaTest
analyzing data of subj1 from Haxby et al. (2001)
37 cvManovaTest_model
estimating model


08-Aug-2019 18:32:14 - Running job #1

08-Aug-2019 18:32:14 - Running 'Model estimation'

SPM12: spm_spm (v7120) 18:32:16 - 08/08/2019

08-Aug-2019 18:32:16 - Failed 'Model estimation'
错误使用 spm_spm (line 309)
Data have not been specified.
In file "/Volumes/VMware Shared Folders/E/vm_mac_matlab/spm12/spm_spm.m" (v7120), function "spm_spm" at line 309.
In file "/Volumes/VMware Shared Folders/E/vm_mac_matlab/spm12/config/spm_run_fmri_est.m" (v7354), function "spm_run_fmri_est" at line 36.

The following modules did not run:
Failed: Model estimation

错误使用 MATLABbatch system
Job execution failed. The full log of this run can be found in MATLAB command
window, starting with the lines (look for the line showing the exact #job as
displayed in this error message)

Running job #1

in 4d-files, volume indicator suffix ",1" is ignored

removed by fullfile?
no, probably ignored by

    % read data directly (faster)
    y = V.private.dat(:, :, :, V.n(1));
    y = reshape(y, [], V.dim(3));
    y = bsxfun(@plus, bsxfun(@times, y, V.pinfo(1, :)), V.pinfo(2, :));
    Y(i, :) = y(mask);

Minor bug in fletcher16.m

Dear Dr. Allefeld,
When the ‘date’ in 'spm.mat' contains chinese characters(it is due to the system's language setting), Funcation 'Checksum' will prompt the UID is ‘invalid data’.
I removed the call to the ‘date’ by modifying the code in cvManovaSearchlight.m(line 56,57):

    uid = gencode({‘SPM.mat’ , ...
    slRadius, Cs, permute, lambda});

after that it runs smoothly, maybe there is a better way to solve the bug.
Thanks sincerely for your work.

best
Liu

Pattern and activation

Hi,

This is not the coding issue but an open question. I was wondering if cvMANOVA can capture both 'activation' and 'neural pattern information' of interest. I ask this question because we want to compare standard searchlight and cvMANOVA, and we found standard searchlight was sensitive to typical activation region and others, while cvMANOVA was not sensitive to typical activation region. Is there any way how the method handles differences in activation magnitude?

Thank you.

error loading data in SPM.mat

Hi there

I'm receiving an error while loading the SPM.mat file using cvManovaSearchlight.m:

loading data
via X:\SPM.mat
153926 in-mask voxels
no region mask
reading images
Undefined function 'diff' for input arguments of type 'cell'.

Error in loadDataSPM (line 105)
pattern(~all(diff(SPM.xY.P) == 0)) = '?';

Error in cvManovaSearchlight (line 47)
[Ys, Xs, mask, misc] = loadDataSPM(dirName);

thanks for help!
cheers, mike

Issue in contrasts with parametric regressors

Dear Dr. Allefeld,
I was trying to run the cvMANOVA on my dataset and I encountered an issue.
My task is a 2x2 factorial design (factors: category (2) and familiarity (2)), with 9 runs (sessions) in total.
For each session, I have the following regressors:

• 4 parametric regressors modelling each factorial combination of category and familiarity. For each of these 4 aggressors,

we also added a parametric modulator.
• 2 regressors modelling auditory cues and subject’s response
• 6 confound regressors modelling motion correction parameters.

In addition, in some participants we have an additional regressor modelling errors (or missed responses), that we inserted at the very end, so that it can be zero-padded by your function as reported in the readme.txt file.
The issue we are encountering is relative to setting up the contrasts for the cvMANOVA. We are interested for the moment in computing category and familiarity main effects, both for “regular” regressors and for their parametric counterpart, to understand where in the brain neural response is modulated by our parameter.
We set up the contrasts as follows:

%Main effect
Cs{1} = [1 0 -1 0 1 0 -1 0 0 0];

%Parametric main effect
Cs{2} = [0 1 0 -1 0 1 0 -1 0 0];

In the second contrast, we aim to replicate the first but this time using the parametric regressors factorXparameterXbf, as reported in SPM.
When we run the cvMANOVA, we encounter the following error:

Error using cvManovaCore (line 82)
contrast 2 is not estimable in session X!

This happens only with the parametric contrast, with the “regular” main effect the computation starts just fine. We tried on multiple subjects, the same error message appears and the session number that creates the error varies systematically.
Following another question we found in this thread (#27) and thinking the error might be due to the additional regressor present in some participants, we tried and run the cvMANOVA only in participants which did not commit any error (so with an equal number of regressors in each run), but the same error message appeared.
Here I attach an SPM.mat file of a participant without the additional error regressor (https://we.tl/t-F2tSL0yLdv).
Thanks for your help.
Best,

Flavio

cvManovaRegion

Hi,

I have another issue with cvManovaRegion. The issue is that this works well with some ROI mask, but does not work with some other ROI masks. In my case, I created three ROI region masks, (cell array of 3D volumes of 1 and 0). ROI 1 and 2 have some overlappings and ROI 3 is union of ROIs 1 and 2. The problem is that cvManovaRegion works only for ROI 2, but and ROI 1 and 3, and the error message is that,

data insufficient for the 1259 voxels of region 1!
data insufficient for the 1526 voxels of region 3!

Can you please help? Thanks.

minor bugs in cvManovaTest

Hi Allefeld,

Great work and very much appreciated for making the toolbox available!

I cloned the repo today and went through the example with the cvManovaTest scripts, where I encountered 3 minor problems:

  1. websave was introduced in MATLAB R2014b, so I needed to download the dataset manually since I am using MATLAB R2014a.
    see line 8 in cvManoveTest_getdata.m

  2. sprintf runs in trouble in Windows, since the escape character is the same as the file separator.
    see line 10 of cvManovaTest_preprocess.m
    see line 19 of cvManovaTest_model.m

  3. cvManovaTest_preprocess.m uses the same variables (e.g. nRuns) as cvManoveTest_getdata.m. In case the design has been created earlier in another session (exist(fnDesign, 'file') evaluates true in cvManovaTest_getdata.m), then these variables are missing during running cvManovaTest_preprocess.m.

After fixing these minor issues, the code runs smoothly.

Thanks and best,
Agoston

profile memory usage

Really hard to do. It looks like peak memory usage is about double the amount for the loaded data.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.