permutation testing: opening thoughts

This is the first of what will probably be a continuing series of musing, meandering posts about permutation testing for group fMRI MVPA.

I'm a big fan of permutation testing, especially for MVPA: our questions, analyses, sample sizes, distributions, etc. are so far from what was envisioned for parametric stats that using a "standard" parametric test "out of the box" gives me pause (or, in practice, lots of assumption-checking). I like that with permutation tests it's possible to design a test that gets at exactly the hypothesis we're asking, while keeping the relevant variance structures (e.g. how many runs, type of classifier, number of examples in each cross-validation fold) intact.

But, the exact way two different people design a permutation test for any particular question will probably be different. And this gets worse in group analyses.

I'll be writing a series of posts illustrating some of the common ways permutation tests are done, both single subject and group analyses. My intent is that these posts will highlight and explain how to carry out what I think are the most useful and common methods, but also share my perception of where some of the open questions lie. As always, please feel free to comment, especially pointing out where your practice/observation/experience/conclusions differ from mine.

more searchlights: the Princeton mvpa toolbox

MS Al-Rawi kindly sent me coordinates of the searchlights created by the Princeton MVPA toolbox, which I made into this image. As before, the center voxel is black, the one-voxel-radius searchlight red, the two-voxel-radius searchlight purple, and the three-voxel-radius searchlight blue.

The number of voxels in each is:
1-voxel radius: 6 surrounding voxels (7 total)
2-voxel radius: 32 surrounding voxels (33 total)
3-voxel radius: 122 surrounding (123 total)
4-voxel radius: 256 surrounding (257 total)

From a quick look at the code generating the coordinates, it appears to "draw" a sphere around the center voxel, then retrieve the voxels falling within the sphere. (anyone disagree with this characterization?)

This "draw a sphere first" strategy explains why I couldn't come up with this number of voxels for a three-voxel-radius searchlight if the shape follows additive rules (e.g. "add voxels that share a face with the next-size-smaller radius"). It's another example of how much it can vary when different people actually implement a description: to me, it was more natural to think of searchlights of different radii as growing iteratively, a rather different solution than the one used in the Princeton MVPA toolbox.

my searchlight, and some (final?) thoughts

Since I've been thinking about searchlight shapes, here's a diagram of the two-voxel radius searchlight I'm currently using in a few projects. As in the previous post, the center voxel is black, the one-voxel radius surrounds red, and the two-voxel radius surrounds purple.

This two-voxel radius searchlight has 13 + 21 + 24 + 21 + 13 = 92 voxels in the surround, which makes it larger than the three-voxel radius faces-must-touch searchlight (at 86 surrounding voxels).

So how should we grow searchlights, edge-face-corner, edge-face, faces-only? Beats me. I can imaging comparing the results of a few different shapes in a particular dataset, but I doubt there's a single shape that is always best. But the searchlight shape should be added to our list of how to describe a searchlight analysis.

PS - I use R code to generate lists of adjacent voxels for running searchlight analyses; contact me if you'd like a copy/details.

pyMVPA searchlight shapes

Michael Hanke listed the number of voxels for different searchlight radii in pyMVPA. He also wrote that pyMVPA considers neighboring voxels to be those that share two edges (i.e., a face).

Here's a diagram of following the faces-touching rule. The center voxel is black, the one-voxel radius surrounds red, the two-voxel radius surrounds purple, and the three-voxel radius surrounds blue.

This gives 1 + 9 + 21 + 24 + 21 + 9 + 1 = 86 surrounding voxels for the three-voxel radius.

But Michael wrote that a three-voxel radius searchlight in pyMVPA has 123 voxels (counting the center). If so, at left is how I could come by that count: 13 + 28 + 40 + 28 + 13 = 122.

This shape is not quite symmetrical: there are four voxels to the left, right, front, and back of the center voxel but only two above and below.

Anyone know if this is what pyMVPA does?

UPDATE: Neither of these sketches is correct, see pyMVPA searchlight: solved.

searchlight shapes

How exactly the searchlight was shaped is one of those details that is usually not mentioned in searchlight analysis papers but is not obvious, mostly since voxels are cubes.

Nikolaus Kriegeskorte (the father of searchlight analysis) uses the illustration on the left.

Expanded, I think this would look like the picture at left. There are 32 voxels surrounding the center; the center voxel in black and surrounding in blue.

But there are other ways of thinking of a searchlight. Here are three ways of defining a one-voxel radius searchlight.

Personally, I implemented the middle version (edges and faces, not corners, 18 voxels) as a one-voxel radius searchlight in my code. Larger radii are calculated iteratively (so for a two-voxel radius, combine the one-voxel radius searchlights defined by treating all 18 surrounding voxels present in the one-voxel radius searchlight as centers). I don't think one version of defining surrounds is necessarily better than others, but we do need to specify the voxels we are including.


update 11 September: I changed the caption of the searchlight shape images for clarity.

what MVPA detects

MVPA is not a universal solution for fMRI; it's always necessary to carefully think about what types of activity changes you want to detect.

For example, suppose these are timecourses for two people completing a task with two types of trials. It's clear that in both people there is a very strong effect of the trials: for person #1 the BOLD goes up during trial type A and down during trial type B; the reverse is true for person #2.

"Standard" MVPA (e.g. linear svm) will detect both of these patterns equally well: there is a consistent difference between trial types A and B in both people. In addition, the difference in direction is usually not reflected in the analysis: often only each subject's accuracy is taken to the second level.

This can be a feature, or a bug, depending on the hypotheses: If you want to identify regions with consistent differences in activation in each person, regardless of what those differences are, it's a feature. If you want to identify regions with a particular sort of difference in activation, it can be a bug.

My suggestion is usually that if what you want to detect is a difference in the amount of BOLD (e.g. "there will be more BOLD in this brain region during this condition compared to that condition") then it's probably best to look at some sort of mass-univariate/GLM/SPM analysis. But if you want to detect consistent differences regardless of the amount or direction of BOLD change (e.g. "the BOLD in this brain region will be different in my conditions"), then MVPA is more suitable.

Note also that a linear svm is perfectly happy to detect areas in which adjacent voxels have opposite changes in BOLD - the two timecourses above can be within the same ROI yet be detected quite well as an informative area. As before, this can be a feature or a bug. So, again, if you want to detect consistent regional differences in the overall amount of BOLD, you probably don't want to use "standard" MVPA.

Nipype: comments from Yaroslav

Yaroslav Halchenko was kind enough to provide some detailed feedback on my previous post, which he allowed me to excerpt here.

"a decent understanding of both nipype and the programs it's calling (from fsl or whatever)."

Yaroslav: yes and no. As with any software -- decent understanding is desired but not strictly necessary. nipype comes with "pre-cooked" workflows for some standard analyses so you would just need to tune the script/parameters up to input your data/design.

"as most useful to a skilled programmer wanting to run several versions of an analysis" 

Yaroslav: That is one of the cases. I could argue that the most useful is an efficient utilization of computing resources. i.e running many computations in parallel (e.g. per each functional run/ subject etc). 

also -- noone runs their analysis ones with all parameters decided a priory and not 'tuning them up' -- that is also where nipype comes useful since it would be smart to recompute only what needs to be recomputed. Thus possibly helping to avoid human errors.

 as for a "skilled programmer" -- well... somewhat -- scripting nipype may be could be easier but altogether it is not that bad actually. Just as with any scripting/hacking with practice comes (greater) efficiency. 

Jo: I saw mentions of using it in parallel. Our computer support people had a devil of a time (and basically failed) with python on our shared systems in the past. But we're also having problems with matlab/spm scripts tying up way too much of the systems and not getting properly managed. So that might be worth a closer look for us.

"It's often necessary to check initial output before moving to the next step in the pipeline; these errors would need to be caught somehow (go back and check intermediate files?)" 

Yaroslav: you could code your "pipeline" incrementally indeed and check all the intermediate results... or just run it once and then check them and rerun pipeline if you would need to change any of the parameters. I do not see it being much different from what you would get with other software -- once again boiling down to how particular/anal you want to be about checking your data/results. Also you could code/assemble pieces which would provide you with reports of the specific processing steps so you could later on glance over them rapidly.

"Nipype will need to be updated any time one of the underlying programs is updated." 

Yaroslav: again -- yes and no -- nipype (interfaces) would need to be updated only if corresponding programs' command line interface gets changed in some non-backward compatible fashion. Most of the software suites avoid such evil actions; but yes -- nipype would need to be (manually) updated to accommodate new cmdline options, and at times remove some if they were removed from the original software.

Jo: And I would think the defaults. It'd be troublesome if, say, a default Nipype motion-correction routine used a different motion-correction algorithm specification than the GUI version.

Yaroslav: well -- they are trying to match the defaults and as long as it doesn't change and you do not mix GUI/non-gui -- how do you know that GUI's is different/better.

Jo: I don't know that the GUI options are better. But practically speaking, many people (at least with spm) use the default options for most choices, so I was thinking we'd need to confirm that the options are the same running it through the GUI or Nipype, just for consistency.

"I've had terrible troubles before with left/right image flipping and shifting centerpoints when changing from one software package to another (and that's not even thinking about changing from one atlas to another). Trying to make packages that have different standard spaces play nicely together strikes me as very non-trivial..."

Yaroslav:  RIGHT ON -- that is indeed one of the problems I have heard about (mixing AFNI and FSL IIRC) and to say the truth I am not sure how far nipype could help assuring correct orientation... if that is possible at all -- it would be a good question to nipype folks I guess. BUT many tools do play nicely together without doing all kinds of crazy flips, so, although it is a valid concern, it is manageable (run pipeline and verify correct orientation by inspecting results through the steps) and might be limited to only a limited set of abusers.

"... Second, would there be a benefit to using Nipype if the pipeline only includes one program (e.g. why Nipype instead of spm8 batching)?"

Yaroslav:   well -- not everyone is using spm8 and those tools might lack back processing/dependencies tracking etc... and even with spm8 -- you might like eventually to add few additional steps provided by other toolkits and would be ready to do so with ease (given you also check for possible flip).