Friday, April 10, 2015

Automating SPM Contrasts

Manually typing in contrasts in SPM is a grueling process that can have a wide array of unpleasant side effects, including diplopia, lumbago, carpal tunnel syndrome, psychosis, violent auditory and visual hallucinations, hives, and dry mouth. These symptoms are only compounded by the number of regressors in your model, and the number of subjects in your study.

Fortunately, there is a simply way to automate all of this - provided that each subject has the same number of runs, and that the regressors in each run are structured the same way. If they are, though, the following approach will work.

First, open up SPM and click on the TASKS button in the upper right corner of the Graphics window. The button is marked "TASKS" in capital letters, because they really, really want you to use this thing, and mitigate all of the damage and harm in your life caused by doing things manually. You then select the Stats menu, then Contrast Manager. The options from there are straightforward, similar to what you would do when opening up the Results section from the GUI and typing in contrasts manually.

When specifying the contrast vector, take note of how many runs there are per subject. This is because we want to take the average parameter estimate for each regressor we are considering; one can imagine a scenario where one of the regressors occurs in every run, but the other regressor only happens in a subset of runs, and this more or less puts them on equal footing. In addition, comparing the average parameter or contrast estimate across subjects is easier to interpret.

Once you have the settings to your satisfaction, save it out as a .mat file - for example, 'RunContrasts.mat'. This can then be loaded from the command line:

load('RunContrasts')

Which will put a structure called "jobs" in your workspace, which contains all of the code needed to run a first-level contrast. The only part of it we need to change when looping over subjects is the spmmat field, which can be done with code like the following:

subjList=[207 208]; %And so on, including however many subjects you want

for subj=subjList

    jobs{1}.stats{1}.con.spmmat =     {['/data/hammer/space4/MultiOutcome2/fmri/' num2str(subj) '/RESULTS/model_multiSess/SPM.mat']} %This could be modified so that the path is a variable reflecting where you put your SPM.mat file
    spm_jobman('run', jobs)

end

This is demonstrated in the following pair of videos; the first, showing the general setup, and the second showing the execution from the command line.





21 comments:

  1. You don't even have to create your .mat with a gui if in the end you're using a Matlab script. For example here is what i used for a two-sample t-test between patients and control groups:
    matlabbatch = [];
    matlabbatch{1}.spm.stats.con.spmmat = cellstr(strcat(save_folder,'/SPM.mat'));
    matlabbatch{1}.spm.stats.con.consess{1}.tcon.name = 'Control > Patients';
    matlabbatch{1}.spm.stats.con.consess{1}.tcon.convec = [1 -1];
    matlabbatch{1}.spm.stats.con.consess{1}.tcon.sessrep = 'none';
    matlabbatch{1}.spm.stats.con.consess{2}.tcon.name = 'Patients > Control';
    matlabbatch{1}.spm.stats.con.consess{2}.tcon.convec = [-1 1];
    matlabbatch{1}.spm.stats.con.consess{2}.tcon.sessrep = 'none';
    matlabbatch{1}.spm.stats.con.delete = 1;
    spm_jobman('run',matlabbatch);

    Of course this is valid for subject level contrasts. Makes it easier when you have a little thing to change as well.

    ReplyDelete
    Replies
    1. Hey Rodolphe,

      That's true, and very useful; I show it with the GUI since for newcomers that seems to be a more intuitive way to set up the structure. But either way would work.

      -Andy

      Delete
  2. Dear Andrew Jahn,

    First of all, thank you for all your efforts.


    I was wondering about the weight contrasts in first level...

    I'v seen that in several papers the rule is perform a [1 -1] (condition 1 > condition) and then use this contrasts in a second level analysis.

    But, my question relay if I can perform contrasts with a weight of [1 0] to identify voxels whose activation increased in one condition and then use this img contrasts in second level.

    Is this correct?


    Thank you in advance!

    Cheers!

    ReplyDelete
    Replies
    1. Hi Richard,

      Yes, that's right; if you simply have a [1 0] contrast vector, then that will give you the equivalent of a simple effect, i.e. only looking at a single condition (or cell, if your design is factorial). However, contrasts between conditions are more highly recommended since if the conditions are controlled for confounds and differ only on a particular quality of interest, than the resulting contrast isolates that quality. This is also referred to as the assumption of pure insertion, which isn't perfect, but is often a better approach than just looking at simple effects.


      Best,

      -Andy

      Delete
  3. Hey Andrew, thanks for all the posts, it's been very informative for getting used to interacting with SPM through the command line. I was wondering though if you could elaborate a bit more on setting the contrasts to .2 or -.2 versus 1 or -1. I understand you are summing them up across the number of runs so, you end up with 1 or -1 versus -5 or 5. But are your results going to be much different if you ran the contrasts at -1 and 1 over 5 runs than -.2 or .2 over 5 runs? I realize I could answer my own question with a bit of time but wanted to know more about why you do it that way. Is that common in the field?

    Thanks a bunch!

    ReplyDelete
    Replies
    1. Hey Tony,

      I think the reason is that if you have different numbers of runs for different subjects, it can be biased. I'm not 100% sure on this, but the convention is to have the weights sum to 1 and -1, and I believe that is why.

      -Andy

      Delete
  4. Hi Andrew,

    Thanks for your wonderful efforts in order to help us getting better at this complicated thing that is neuroimaging.

    I am relatively comfortable with writing T-contrasts, but I am trying for the very first time to write an F contrast and I can't find my answer anywhere.

    I ran a longitudinal VBM analysis using the VBM8 manual protocol. I have two groups of individuals that have been scanned twice over 2 years. I chose the flexible factorial design with factors 1: Subject, 2: Group, 3: Time, and I also entered ICV as a covariate.

    In the contrast manager, I entered a T contrast for the Main Group Effect (0 1), a T contrast for the Main Time Effect (0 0 1), but I don't know how to write the Group X Time interaction. Do you know how to write this?

    Also, am I supposed to include my covariate in my contrasts?

    Best regards,

    Nick

    ReplyDelete
    Replies
    1. Hey Nick,

      I do what I can! For these kinds of interactions, one of the best guides I know of is written by Glascher and Gitelman; if you include subject as a factor and then have a 2x2 design, your interaction term would look something like -0.5 0.5 0.5 -0.5 (nonlinear, parabola-shaped in this case). Check out the Design 1 graphic on page 8; I think that should give you what you want: http://www.sbirc.ed.ac.uk/cyril/download/Contrast_Weighting_Glascher_Gitelman_2008.pdf

      If you have a single number covariate, such as age (or ICV, in your case), then yes, you can add it as a covariate at the 2nd level. If you wanted to test whether any of the betas correlate significantly with ICV, you would just put a 1 in that column, and leave the rest as zeros.

      Best,

      -Andy

      Delete
  5. Hi Andy,

    How would the weights be affected if the size of the regressors within each run changes? I.e. A subject sees 4 images fitting Regressor A, but 5 images fitting Regressor B in run 1, and then sees 6 images for Reg A and 5 images for Reg B in run 2.
    Would A > B be set up as [0.5 -0.5 0.5 -0.5] or would it have to account for the variance within each run?

    Thanks

    ReplyDelete
    Replies
    1. Hi there,

      You can still use that contrast vector; it will take the average of the parameter estimate across the runs, and any variability not explained by the parameter will go into the error term.

      Best,

      -Andy

      Delete
  6. Hi, Andy. I've just started doing cognitive neuroscience research, and my advisor asked me to find out what contrast analysis we can do on fMRI data to discover trends (using fMRI contrasts). Any help?

    ReplyDelete
    Replies
    1. Hey there,

      Are you talking about trends as in linear and quadratic trends? It depends on how many regressors you have, but assuming that you have 4, a linear trend analysis would look like this:

      [-1.5 -0.5 0.5 1.5]

      Note that this is the result of the command (1,2,3,4) - mean(1,2,3,4), a formula you can apply to any number of regressors.

      For more details, see this PDF by Jeanette Mumford: http://mumford.bol.ucla.edu/lin_cont_illus.pdf


      Best,

      -Andy

      Delete
  7. Hi Andy,
    I'm a fMRI beginner.
    I read some information from the Internet and my own lab from prior students about the contrast.
    However, there're some different at the first level contrast manager.
    In most information from the Internet, they set contrasts equal to conditions, it means, if I have 2 condition, my contrast will be [condition 1, condition 2].
    However, in my lab's manual, they have conditions and time derivatives, it means, if I have 2 condition, my contrast will be [condition 1, condition 1 time derivatives, condition 2, condition 2 time derivatives].
    I rarely found information (because I haven't read enough papers) about the time derivatives.
    Thanks....

    ReplyDelete
    Replies
    1. Hey Miss Life,

      A time derivative is a regressor that captures any variability in the onset of the HRF. (I can't visualize what a time derivative would look like, so I don't understand it that well.) If you think there could be a lot of variability in the onset of the HRF in whatever region you're looking at, then including the time derivative could give you a better model fit.

      However, a paper from Della-Maggiore (2002) showed that including the time derivative decreases power, suggesting that just using the canonical HRF may be good enough to explain most BOLD responses. That paper is pretty old though, and I haven't kept up with what has happened since then.

      I would try doing a couple of models both with and without the time derivative, and see what your contrast maps look like. If there's no significant difference, I would omit the time derivatives, since I like to simplify things.

      -Andy

      Delete
  8. Hi Andy,

    How would the process go when you have different number of conditions in each run and they differ from subject to subject? Thanks for your help

    ReplyDelete
  9. Hi Andy,

    Love your blog/videos -- very helpful!

    I am running into a 'Maximum recursion limit' error while trying to add group-level contrasts to a whole brain analysis. I have 3 groups in the analysis, and a hypothesis that two of the groups are more similar than the other. I've created a t-contrast [2 -1 -1], and when I try run it I always get the error: Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N) to change the limit. Be aware that exceeding
    your available stack space can crash MATLAB and/or your computer. Error in spm_create_image.

    Is there a different way that I should be specifying the contrast? I have also tried [1 -.5 -.5], which didn't work either.

    Many thanks!

    ReplyDelete
    Replies
    1. Hey there,

      I haven't seen that error before, are you getting it by using the code in the blog post, or through using the contrast manager? Try it with the contrast manager first, and if that works, then we'll see about debugging the code.


      Best,

      -Andy

      Delete
    2. Hi Andy,

      I got the error using the contrast manager (I haven't tried your code). If you think that the contrasts I entered look appropriate, then maybe there is something going on with my version of SPM/Matlab.

      I am able to run contrasts between two group (e.g. [1 0 -1]), but not 3.

      Any thoughts are helpful, and I don't want to take up too much of your time!

      Delete
    3. Hey Ann,

      Try using an F-contrast; e.g.:

      [2 -1 -1
      -1 2 -1
      -1 -1 2]

      Which will test all three contrasts simultaneously. I'm not sure if this will work or not, but if it does, the resulting map won't reveal which one of those contrasts is driving the effect; you'll have to extract the parameter estimates from an ROI and plot them.


      Best,

      -Andy

      Delete
  10. Hi Andy,

    Thanks a lot for all your content, it is really helpful and well explained. I was wondering, is it possible to automate the process when you have a choice based model? I am talking of models where there is the possibility that all conditions (and therefore, all the contrasts) don't exist for some participants that never selected a particular option. Is it possible to create a script that does this? I've managed to do it with many "if" at the end, to create blocks of consessions for each case, but I was wondering if there is a "cleaner way to do it". Thanks for your help!

    Best,
    Andrea

    ReplyDelete