Allefeld 2014: Searchlight-based multi-voxel pattern analysis of fMRI by cross-validated MANOVA

A recent paper by Carsten Allefeld and John-Dylan Haynes, "Searchlight-based multi-voxel pattern analysis of fMRI by cross-validated MANOVA (see full citation below), caught my eye. The paper advocates using a MANOVA-related statistic for searchlight analysis instead of classification-based statistics (like linear SVM accuracy). Carsten implemented the full procedure for SPM8 and Matlab; the code is available on his website.

In this post I'm going to describe the statistic proposed in the paper, leaving the discussion of when (sorts of hypotheses, dataset structures) a MANOVA-type statistic might be most suitable for a (possible) later post. There's quite a bit more in the paper (and to the method) than what's summarized here!

MANOVA-related statistics have been used/proposed for searchlight analysis before, including Kriegeskorte's original paper and implementations in BrainVoyager and pyMVPA. From what I can tell (please let me know otherwise), this previous MANOVA-searchlights fit the MANOVA on the entire dataset at once: all examples/timepoints, no cross-validation. Allefeld and Haynes propose doing MANOVA-type searchlights a bit differently: “cross-validated MANOVA” and “standardized pattern distinctness”.

Most of the paper's equations review multivariate statistics and the MANOVA; the “cross-validated MANOVA” and “standardized pattern distinctness” proposed in the paper are in equations 14 to 17:

• Equation 14 is the equation for the Hotelling-Lawley Trace statistic, which Allefeld refers to as D ("pattern distinctness").
• Equation 15 shows how Allefeld and Haynes propose to calculate the statistic in a "cross-validated" way. Partitioning on the runs, they obtain a D for each partition by calculating the residual sum-of-squares matrix (E) and the first part of the H equation from the not-left-out-runs, but the second part of the H equation from the left-out run.
• Equation 16 averages the D from each "cross-validation" fold, then multiplies the average by a correction factor calculated from the number of runs, voxels, and timepoints.
• Finally, equation 17 is the equation for “standardized pattern distinctness”: dividing the value from equation 16 by the square root of the number of voxels in the searchlight.

To understand the method a bit better I coded up a two-class version in R, using the same toy dataset as my searchlight demo. Note that this is a minimal example to show how the "cross-validation" works, not necessarily what would be useful for an actual analysis, and not showing all parts of the procedure.

The key part from the demo code is below. The dataset a matrix, with the voxels (for a single searchlight) in the columns and the examples (volumes) in the rows. There are two classes, "a" and "b". For simplicity, the "left-out run" is called "test" and the others "train", though this is not training and testing as meant in machine learning. train.key is a vector giving the class labels for each row of the training dataset.

For Hotelling's T2 we first calculate the "pooled" sample covariance matrix in the usual way, but using the training data only:
S123 <- ((length(which(train.key == "a"))-1) * var(train.data[which(train.key == "a"),]) + (length(which(train.key == "b"))-1) * var(train.data[which(train.key == "b"),])) / (length(which(train.key == "a")) + length(which(train.key == "b")) - 2);

To make the key equation more readable we store the total number of "a" and "b" examples:
a.count <- length(which(train.key == "a")) + length(which(test.key == "a"));
b.count <- length(which(train.key == "b")) + length(which(test.key == "b"));

and the across-examples mean vectors:
a.test.mean <- apply(test.data[which(test.key == "a"),], 2, mean);
b.test.mean <- apply(test.data[which(test.key == "b"),], 2, mean) ;  
a.train.mean <- apply(train.data[which(train.key == "a"),], 2, mean); 
b.train.mean <- apply(train.data[which(train.key == "b"),], 2, mean);

now we can calculate the Hotelling's T2 (D) for this "cross-validation fold" (note that solve(S123) returns the inverse of matrix S123):
((a.count*b.count)/(a.count+b.count)) * (t(a.train.mean-b.train.mean) %*% solve(S123) %*% (a.test.mean-b.test.mean));

The key is that, paralleling equation 15, the covariance matrix is computed from the training data, multiplied on the left by the mean difference vector from the training data, then on the right by the mean difference vector from the testing data.

Should we think of this way of splitting the Hotelling-Lawley Trace calculation as cross-validation? It is certainly similar: a statistic is computed on data subsets, then combined over the subsets. It feels different to me though, partly because the statistic is calculated from the "training" and "testing" sets together, and partly because I'm not used to thinking in terms of covariance matrices. I'd like to explore how the statistic behaves with different cross-validation schemes (e.g. partitioning on participants or groups of runs), and how it compares to non-cross-validated MANOVA. It'd also be interesting to compare the statistic's performance to those that don't model covariance, such as Gaussian Naive Bayes.

Interesting stuff; I hope this post helps you understand the procedure, and to keep us all thinking about the statistics we choose for our analyses.


Allefeld C, & Haynes JD (2014). Searchlight-based multi-voxel pattern analysis of fMRI by cross-validated MANOVA. NeuroImage, 89, 345-57 PMID: 24296330

demo: R code to perform a searchlight analysis

Here is an R script for performing a single-subject searchlight analysis with R. You are welcome to adapt the code for your own use so long as I am cited; the code is a demo that performs an entire single-subject searchlight analysis on the sample dataset, but will need to be edited for other datasets. This is an R script, not a package, because it is not "plug-and-play", though I hope it is straightforward and clear.

The demo takes a 4d nifti as input and creates a 3d nifti output image in which the value in each voxel is the classification accuracy of the voxel's searchlight. The code creates iterative-shell shaped searchlights, but can be altered for other shapes.

The code performs two preliminary steps before doing the searchlight analysis. These just need to be done once for each dataset (e.g. if multiple subjects have been spatially normalized, all use the same lookup and 3d files). The first creates a 3d.rds file, which assigns an integer to each brain voxel and serves to go between each voxel's 3d coordinates and the list of searchlight surrounds, which is stored in a lookup table. Each row of the lookup table corresponds to the brain voxel integer labels stored in 3d.rds; the columns list the voxel-name integers making up the searchlight for that row's voxel. The searchlight analysis proceeds by going through each voxel in turn, looking up its searchlight voxels in the lookup table, pulling the values for those voxels from the 4d input image, then writing the classification accuracy into the voxel's place in the 3d output image.

I hope that this code and/or logic is useful for people setting up their own searchlight analyses, whether in R or another language. I believe that most searchlight analysis implementations use a similar logic, though some others are probably more optimized for speed! If you give it a try, please let me know how it goes.

I created the demo dataset from a little bit of a real dataset, labeling volumes to have strong information in motor and visual areas. This is what my final searchlight accuracy map looked like, plotted on the ch2bet template in MRIcroN and scaled for accuracies of  0.5 to 1.

UPDATE 13 May 2015

The searchlight demo code above has code for creating searchlights shaped as iterative shells. It is now perhaps most common to make "spherical" searchlights, so this updated searchlight demo has code to create searchlights of the same shape as used by pyMVPA, BrainVoyager, and the Princeton toolbox. I also cleaned up a few comments, and changed the lookup table and 3d image to plain text and NIfTI, instead of .rds, for access in multiple programs, not just R.

This shows the accuracy maps produced from this little demo dataset using 2-voxel radius searchlights, shaped as iterative shells (right) and "spheres" (left). The areas with the highest accuracies are the same in both, as we'd hoped, but the accuracies tend to be a bit higher on the right (shells), likely because 2-voxel radius spherical searchlights are smaller (have fewer voxels) than 2-voxel radius iterative shell searchlights. (Note that this doesn't mean one shape is better than the other, they're simply different shapes, and it's probably best for consistency that we all use the same shape, which means spheres.)

nice doughnuts: pretty searchlight quirks

More than a year ago I posted an explanation about how strange ring-type edge effects can show up in searchlight analysis information maps, but this morning I saw the most vivid example of the phenomenon I've ever run across in real data.

These images show slices from two leave-one-subject-out searchlight analyses (linear SVM). The subjects were viewing all sorts of movie stimuli, and for a positive control I classified whether the people were viewing stimuli showing people moving, or whether the stimuli showed still people. This analysis should work very well, reporting massive information (high classification accuracy, in this case) in motor, somatosensory, and visual areas.

At right are a few slices of the information map that results when the examples are normalized within each searchlight (i.e. to a mean of 0, standard deviation of 1). The highest-classifying searchlights are in the proper brain areas (and with properly high accuracies), but in strange, hollow, shapes.


And here's what it looks like when the images are not normalized within each searchlight. The same areas have accurate searchlights, but now they make solid "blobs", with the most highly accurate searchlights (yellow) in the centers - where normalized searchlights had chance accuracy.

Read that previous post for a full explanation. This is the sort of data we'd expect to have the strongest edge/center effects: focal areas with high information, areas large enough to allow searchlights to fit inside the middle of the informative areas. It's also likely that the activations in these areas are all in the same direction (higher activity during the moving stimuli than the non-moving stimuli), creating consistent differences in the amount of BOLD throughout the areas. Normalizing each example within each searchlight removes this consistent BOLD difference, and so classification fails.

Schapiro 2013: response from Anna Schapiro

Anna Schapiro, first author of "Neural representations of events arise from temporal community structure", kindly provided additional details (even a diagram!), addressing the questions I asked about her research in my previous post. She gave me permission to share extracts from her reply here, which we hope will be useful to others as well.

She said that my descriptions and guesses were basically accurate, including that they used the searchmight toolbox to perform the searchlight analysis, and so used cubical searchlights by default.

My first question on the previous post was, "Were the item pairs matched for time as well as number of steps?" Anna replied that:

"We did at one point try to balance time as well as the number of steps between items. So we were averaging correlations between items that were both the same number of steps and the same amount of time apart. But I found that some of the time/step bins had very few items and that that made the estimates significantly more noisy. So I opted for the simpler step balance approach. Although it doesn't perfectly address the time issue, it also doesn't introduce any bias, so we thought that was a reasonable way to go."

about the node pairs

My second question was about how many correlation differences went into each average, or, more generally, which within-cluster and between-cluster pairs were correlated.

"Regarding choices of pairs, I think the confusion is that we chose one Hamiltonian path for each subject and used the forwards and backwards versions of that path throughout the scan (see the beginning of Exp 2 for an explanation of why we did this). Let's assume that that path is the outermost path through the graph, as you labeled in your post [top figure]. Then the attached figure [second figure] shows the within and between cluster comparisons that we used, with nodes indexed by rows and columns, and color representing the distance on the graph between the two nodes being compared. We performed the correlations for all node pairs that have the same color (i.e., were the same distance away on the graph) and averaged those values before averaging all the within or between cluster correlations."

Here's the node-numbered version of Figure 1 I posted previously, followed by the figure Anna created to show which node pairs would be included for this graph (the images should enlarge if you click on them).

Concretely, then, if the Hamiltonian path for a person was that shown in Figure 1 (around the outside of the graph, as the nodes are numbered), correlations would be calculated for the node pairs marked in the second figure on each of the approximately 20 path traversals. This is about 60 correlations for each path traversal, sorted into 30 within-cluster comparisons and 30 between-cluster comparisons ("about" 60 since some pairs on each traversal might be omitted depending on where the path started/ended).

The 30 pairs (correlations) for the within-cluster comparisons would be collapsed into a single number by averaging in two steps. First, the correlations of the same length would be averaged (e.g. three for length-4-within-cluster: 11-15, 6-10, 1-5), giving four averages (one for length-1 pairs, one for length-2 pairs, one for length-3 pairs, and one for length-4 pairs). Second, these four averages would be averaged, giving a single average for within- and between-cluster comparisons. This two-step averaging somewhat reduces the influence of path length imbalance (averaging all pairs together would include 12 length-1 pairs in the within-cluster comparison but only 3 length-1 pairs in the between-cluster comparison), though may not eliminate it completely. I wonder if picking three pairs of each length to average (i.e. all length-4 within-cluster but only a third of the length-1 within-cluster) would change the outcome?

group-level analysis

My third question was about how exactly the group analysis was performed. Anna replied that

"yes, we used a one-sample t-test on our within-between statistic values. In this scenario, randomise permutes the sign of the volumes for each subject. On each permutation, the entire volume for a particular subject is left alone or multiplied by -1. Then it looks for clusters of a certain size in this nonsense dataset. In the end it reports clusters in the true dataset are significantly more likely to be found than in these shuffled versions. I like using randomise because it preserves the spatial smoothness in every region of the brain in every subject. Searchlights may create smoothness that have a different character than other analyses, but we don't have to worry about it, since that smoothness is preserved in the null distribution in this permutation test."

I then asked, "Do you mean that it [randomise] permutes the sign of the differences for each person? So it is changing signs on the difference maps, not changing labels on the (processed) BOLD data then recalculating the correlations?", to which she replied that, "I feed randomise within-between maps, so yes, it's permuting the sign of the differences."


Thank you again, Anna Schapiro, for the helpful and detailed replies, and for allowing me to share our correspondence here!

searchlight analysis interpretation: worst case scenario

I like searchlight analysis. It avoids some curse-of-dimensionality headaches, covers the whole brain, makes pretty pictures (brain blob maps!), and can be easier for people accustomed to mass-univariate analyses.

But if I like searchlight analysis, why write a paper about it with "pitfalls" in the title? Well, because things can go badly wrong. I do not at all want to imply that searchlight analysis should be abandoned! Instead, I think that searchlight analyses need to be interpreted cautiously; some common interpretations do not always hold. Am I just picking nits, or does this actually matter in applications?

Here's an example of how searchlight analysis interpretation often is written up in application papers.

Suppose that the image at left is a slice of a group-level searchlight analysis results map for some task, showing the voxels surviving significance testing. These voxels form two large clusters, one posterior in region Y and the other a bit more lateral in region X (I'll just call them X and Y because I'm not great at anatomy and it really doesn't matter for the example - this isn't actually even searchlight data). We write up a paper, describing how we found significant clusters in the left X and Y which could correctly classify our task. Our interpretation is focused on possible contributions of areas X and Y to our task, drawing parallels to other studies talking about X and Y.

This is the sort of interpretation that I think is not supported by the searchlight analysis alone, and should not be made without additional lines of evidence.

Why? Because, while the X and Y regions could be informative for our task, the searchlight analysis itself does not demonstrate that these are more informative than other regions, or even that the voxel clusters themselves are informative (as implied in the interpretation). It is possible that the voxel clusters found significant in the searchlight analysis may not actually be informative outside the context of the particular searchlight analysis we ran (i.e. dependent upon our choice of classifier, distance metric, and group-level statistics).

The problem is the way my hypothetical interpretation shifts from the searchlight analysis itself to regions and clusters: The analysis found voxels which are at the center of significant searchlights, but the interpretation is about the voxels and regions, not the searchlights. Unfortunately, it is not guaranteed that if information is significantly present at one spatial scale (the searchlights) it will be present at smaller ones (the voxels) or larger (the regions).

Back to the hypothetical paper, it is possible that classifying the task using only the voxels making up cluster X could fail (i.e. instead of a searchlight analysis we make a ROI from a significant cluster and classify the task with just those voxels). This is one of my worst-case scenarios: the interpretation that regions X and/or Y have task information is wrong.

But that's not the only direction in which the interpretation could be wrong: X and Y could be informative, but a whole lot of other regions could also be informative, even more informative than X and Y. This is again a problem of shifting the scale of the interpretation away from the searchlights themselves: our searchlight analysis did not show that X and Y are the most significant clusters. One way of picturing this is if we did another searchlight analysis, this time using searchlights with the same number of voxels as Y: we could end up with a very different map (the center of Y will be informative, but many other voxels could also be informative, perhaps more than Y itself).

These are complex claims, and include none of the supporting details and demonstrations found in the paper. My goal here is rather to highlight the sort of searchlight analysis interpretations that the paper describe; the sorts of interpretations that are potentially problematic. But, note the "potentially problematic", not "impossible"! The paper (and future posts) describe follow-up tests that can support interpretations like in my scenario, ways to show we're not in one of the worst cases.

Schapiro 2013: "Neural representations of events arise from temporal community structure"

While not a methods-focused paper, this intriguing and well-written paper includes an interesting application of searchlight analysis which I'll explore a bit here. I'm only going to describe a bit of the searchlight - related analyses here, you really should take a look at the full paper.

First, though, they used cubical searchlights! I have an informal collection of searchlight shapes, and suspect that the authors used cubical searchlights from Francisco's legacy, though I couldn't find a mention of which software/scripts they used for the MVPA. (I don't mean to imply cubes are bad, just a less-common choice.)

a bit of background

Here's a little bit about the design relevant for the searchlight analysis; check the paper for the theoretical motivation and protocols. Briefly, the design is summarized in their Figure 1: Subjects watched long sequences of images (c). There were 15 images, not shown in random order, but rather in orders chosen by either random walks or Hamiltonian paths on the network in (a). I superimposed unique numbers on the nodes to make them easier to refer to later; my node "1" was not necessarily associated with image "1" (though it could have been).

Subjects didn't see the graph structure (a), just long (1,400 images) sequences of images (c). When each image appeared they indicated whether each image was rotated from its 'proper' orientation. The experiment wasn't about the orientation, however, but rather about the sequences: would subjects learn the underlying community structure?

The searchlight analysis was not a classification but instead rather similar to RSA (representational similarity analysis), though they didn't mention RSA. In their words,

"Thus, another way to test our prediction that items in the same community are represented more similarly is to examine whether the multivoxel response patterns evoked by each item come to be clustered by community. We examined these patterns over local searchlights throughout the entire brain, using Pearson correlation to determine whether activation patterns were more similar for pairs of items from the same community than for pairs from different communities."

Using Figure 1, the analysis is asking whether a green node (e.g. number 2) is more similar to other green nodes than to purple or orange nodes. It's not just a matter of taking all of the images and sorting them by node color, though - there are quite a few complications.

setting up the searchlight analysis

The fMRI session had 5 runs, each of which had 160 image presentations, during which the image orders alternated between random walks and Hamiltonian paths. They only wanted to include the Hamiltonian paths in the searchlight analysis (for theoretical reasons, see the paper), which I think would work out to around 5 eligible path-traversals per run (160/15 = 10.6/2 =~ 5); each node/image would have about 5 presentations per run. They didn't include images appearing at the beginning of a path-traversal, so I think there would be something less than 25 total possible image presentations to include in the analyses.

Hamiltonian paths in the graph mean that not all node orderings are possible: nodes of the same color will necessarily be visited sequentially (with the starting point's color potentially visited at the beginning and end of the path). For example, one path is to follow the nodes in the order of the numbers I gave them above: starting at 1 and ending at 15. Another path could be (1:5, 6,8,7,9,10, 11:15). But (1:5, 6,8,10,7,9, 11:15) is NOT possible - we'd have to got through 10 again to get out of the purple nodes, and Hamiltonian paths only visit each node once. Rephrased, once we reach one of the light-colored boundary nodes (1,5,6,10,11,15) we need to visit all the dark-colored nodes of that color before visiting the other boundary node of the same color.

This linking of order and group makes the searchlight analysis more difficult: they only want to capture same-cluster/different-cluster similarity differences due to cluster, not that the different-cluster images appeared separated by more time than the same-cluster images (since fMRI volumes collected closer together in time will generally be more similar to each other than fMRI volumes collected farther apart in time). They tried to compensate by calculating similarities for pairs of images within each path that were separated by the same number of steps (but see question 1 below). 

For example, there are three possible step-length-1 pairs for node 1: 15-1-2; 15-1-3; 15-1-4. The dark-colored nodes (2,3,4; 7,8,9; 12,13,14) can't be the "center" for any step-length-1 pairs since it takes at least 2 steps to reach the next cluster. Every node could be the "center" for a step-length-2 pair, but there are many more valid pairings for the dark-colored nodes than the light-colored ones.

The authors say that "Across these steps, each item participated in exactly four within-cluster pair correlations and exactly four across-cluster pair correlations.", but it's not clear to me whether this count means one correlation of each step-length or if only four pairings went into each average. It seems like there would be many more possible pairings at each step-length than four.

Once the pairings for each person have been defined calculating the statistic for each pairing on each searchlight would be relatively straightforward: get the three 27-voxel vectors corresponding to the item presentation, its same-cluster item presentation, and its different-cluster item presentation. Then, calculate the correlation between the item and its same-cluster and different-cluster items, Fisher-transform, and subtract. We'd then have a set of differences for each searchlight (one for each of the pairings), which are averaged and the average assigned to the center voxel.

I think this is an interesting analysis, and hopefully if you've read this far you found my description useful. I think my description is accurate, but had to guess in a few places, and still have a few questions:

1. Were the item pairs matched for time as well as number of steps? In other words, with the interstimulus interval jittering two steps could be as few as 3 seconds (1img + 1iti + 1img + 1iti + 1img) or as many as 11 seconds (1img + 5iti + 1img + 5iti + 1img).
2. How many correlation differences went into each average? Were these counts equal for step-lengths or in every subject?
3. How was the group analysis done? The authors describe using fsl's randomise on the difference maps; I guess a voxel-wise one-sample t-test, for difference != 0? What was permuted?

I'll gladly post any answers, comments, or clarifications I receive.

UPDATE, 24 May: Comments from Anna Schapiro are here.

Schapiro, A., Rogers, T., Cordova, N., Turk-Browne, N., & Botvinick, M. (2013). Neural representations of events arise from temporal community structure Nature Neuroscience, 16 (4), 486-492 DOI: 10.1038/nn.3331

supplemental information for "Searchlight analysis: Promise, pitfalls, and potential"

My commentary on searchlight analysis, "Searchlight analysis: Promise, pitfalls, and potential" is now up at NeuroImage . All of the supplemental information is there as well, but renamed ... which is rather annoying, since I refer to the R scripts by file name. I'd asked them to fix it (and make the .R files downloadable as text instead of zipped), but evidently they couldn't.

So, here is a mapping of the names, and a link to download all the files (named correctly) at once, saving some headaches. (FYI: .R files are plain text (.txt), suffixed .R to indicate that they contain R code.)

• "Supplementary material 1" is the supplementary text, downloaded as "mmc1.pdf".
• "Supplementary material 2", which downloads as "mmc2.txt", should be dataset.txt, the input data for the example dataset used in many of the examples in the supplementary material.
• "Supplementary material 3", downloads as mmc3.zip, containing the single file "NIM10275-mmc3.R", which is example1.R
• "Supplementary material 4", NIM10275-mmc4.zip, is dataset_exploration.R
• "Supplementary material 5", NIM10275-mmc5.zip is example3.R
• "Supplementary material 6", NIM10275-mmc6.zip is rareInformative.R
• "Supplementary material 7", NIM10275-mmc7.zip is searchlight_group.R
• "Supplementary material 8", NIM10275-mmc8.zip is searchlight.R
• "Supplementary material 9", NIM10275-mmc9.zip is sharedFunctions.R

I will post about some of the examples and suggestions in the paper; I added the new label SA:PPP to make it easier to find all the relevant posts. Feel free to submit questions or comments as well; anonymously or not as you prefer.

announcement: "Searchlight analysis: promise, pitfalls, and potential"

I'm pleased to announce that my recently-accepted Comment for NeuroImage, entitled "Searchlight analysis: promise, pitfalls, and potential" is now online (dx.doi.org/10.1016/j.neuroimage.2013.03.041). I plan to do a series of posts with expanded versions of some of the figures and supplemental material over the next weeks (it's that cool!).

It doesn't look like the supplemental material is officially online yet; let me know if you can't wait or have difficulty accessing the manuscript and I'll send it to you directly. I'm of course also interested in any feedback.

low-accuracy subjects' influence on group statistics

background

This post comes from something I saw in some of my actual searchlight results.The searchlight analysis was performed within a large anatomical area (~1500 voxels, svm within fairly small searchlights), separately in ~15 people. In most subjects, most all searchlights classified quite well, but in two subjects most searchlights classified quite poorly. As usual, there were a few low-accuracy searchlights in the "good" subjects (people in which most searchlights classified well), and a few high-accuracy searchlights in the "bad" subjects (people in which most searchlights did not classify accurately).

I made the group-level results by performing a t-test at each voxel: are the subjects' accuracies greater than chance? (I did not smooth the individual searchlight maps first.) I saw a worrisome pattern in the group maps: the peak (best t-value) areas tended to coincide where the "bad" subjects had their best (most accurate) searchlights. This outcome is not surprising (see below), but is not what we want: applying a strict threshold to the group map will identify voxels. But those particular voxels came out as peak because the "bad" subjects had accurate searchlights in those locations, not because the voxels were more accurate across subjects: the group map was overly influenced by the low-accuracy subjects.

simulation

Here's a little simulation to show the effect (R code is here). The code creates 15 searchlights, each of which contain 10 voxels. There are 12 people in the dataset, 10 of whom have signal in all searchlights, and 2 of whom have signal in all but the first two searchlights. It's a two-class classification, using a linear SVM, partitioning on the two "runs". Voxel values were sampled from a normal distribution, with different means for the voxels with signal but the same means for no signal.

The classification accuracies of the 12 people in each of the 15 searchlights are summarized in these boxplots (the figure should get larger if you click on it; run the code to get the underlying numbers). The accuracies vary a bit over the people and searchlights, as expected, given how the data was simulated. All of the "boxes" are well above chance, and all of the t-values (see figure labels) are above 5. So this seems reasonable.

The searchlights with the highest t-values are 1 and 2: the two searchlights which I assigned to have signal in the two "bad" subjects. You can see why in the boxplots: the bottom whisker in searchlights 1 and 2 only reach 0.6, while all the others have a whisker or outlier closer to chance. Some searchlights (like 8) have two outliers: the two "bad" subjects.

So the t-test didn't make an error: the first two searchlights should have the highest t-values, since they have the most individual accuracies the furthest above chance. But this could have led to an improper conclusion: If we had applied a high threshold (t > 8 in this simulation) we might think the first two searchlights are where the information falls, where in actuality it is evenly distributed.

what to do?

Follow-up testing or different group-level statistics can help to catch this type of situation. As often is the case, precise hypotheses are better (for example, if you want to find searchlights with significant classification in most subjects individually, test for that directly - don't assume that the best searchlights in a t-test will also have that property).

Here are a few suggestions for follow-up testing:

• Look at the individual subject's searchlight maps: Are a few subjects quite different than the rest (such as here, where two people had quite low accuracies compared to the others)? 
• Sensitivity testing can help: how much does the group-level map change when individual subjects are left out? 
• Do the group-level results align more closely with some subjects than others?

surface or volume searchlighting ... not a mixture

I saw a report with an analysis that used volumetric searchlights (spheres) within a grey matter mask (like the one at the right, though theirs was somewhat more dilated). The cortical thickness in the mask was usually around 10 voxels, and they ran a 3-voxel radius searchlight.

I do not think this is a good idea.

My short general suggestion is to do either volumetric searchlighting of the whole brain (or a large 3D subset of it - like the frontal lobes), or do surface-based searchlighting (e.g. Chen, Oosterhof), but not mix the two.

A fundamental problem is that a fairly small proportion of the searchlights will be fully in the brain using a volumetric grey matter mask.
In this little cartoon the yellow circles represent the spherical searchlights, and the two grey lines the extent of the grey matter mask. Many searchlights do not fully fall into the mask; the edges will be sampled differently than the center.

This cartoon tries to convey the same idea: only the strip in the middle of the mask is mapped by searchlights completely in the brain. If informative voxels are evenly distributed in the grey matter strip, searchlights in the middle (fully in the brain) could be more likely to be significant (depending on the exact classifier, information distribution, etc.), a distorted impression.

I've heard it suggested that it's better to use a grey matter mask because that's where the information should be. I don't think that's a good motivation. For one thing, running the searchlight analysis on the whole brain can serve as a control: if the "most informative areas" turn out to be the ventricles, something went wrong. For another, spatial normalization is not perfect. Depending on how things are run (searchlight on spatially-normalized brains or not, etc), there is likely to be some blurring of white and grey matter in the functional images. Letting the searchlights span both may capture more of the information actually present.

One final point. Volumetric searchlighting will put non-strip-wise-adjacent areas (areas touching across a sulcus) into the same searchlight. This might be ok, given the spatial uncertainties accepted in some fMRI. But if you want to avoid this, surface-based methods are the way to go, not a volumetric grey matter mask.

needles and haystacks: information mapping quirks

I really like this image and analogy for describing some of the distortions that can arise from searchlight analysis: a very small informative area ("the needle") can turn into a large informative area in the information map ("the haystack"), but the reverse is also possible: a large informative area can turn into a small area in the information map ("haystack in the needle").

I copied this image from the poster Matthew 
Cieslak, 
Shivakumar
 Viswanathan, 
and 
Scott 
T. 
Grafton
 presented last year at SfN (poster 626.16, Fitting and Overfitting in Searchlights, SfN2011). The current article covers some of the same issues as the poster, providing a mathematical foundation and detailed explanation.

They step through several proofs of information map properties, using reasonable assumptions. One result I'll highlight here is that the information map's representation of a fixed-size informative area will grow as searchlight radius increases (my phrasing, not theirs). Note that this (and the entire paper) is describing the  single-subject, not group level of analysis.

This fundamental 'growing' property is responsible for many of the strange things that can appear in searchlight maps, such as the edge effects I posted about here. As Viswanathan et al. point out in the paper, it also means that interpreting the number of voxels found significant in a searchlight analysis is fraught with danger: it is affected by many factors other than the amount and location of informative voxels. They also show that it is possible to have just 430 properly-spaced informative voxels create the entire brain to be marked as informative in the information map, using just 8 mm radius searchlights (that's not particularly large in the literature).

I recommend taking a look at this paper if you generate or interpret information maps via searchlight analysis, particularly if you have a mathematical bent. It nicely complements diagram- and description-based explanations of searchlight analysis (including, hopefully soon, my own). It certainly does not include all the aspects of information mapping, but provides a solid foundation for those it does include.


Shivakumar Viswanathan, Matthew Cieslak, & Scott T. Grafton (2012). On the geometric structure of fMRI searchlight-based information maps. arXiv: 1210.6317v1

permuting searchlight maps: Stelzer

Now to the proposals in Stelzer, not just their searchlight shape!

This is a dense methodological paper, laying out a way (and rationale) to carry out permutation tests for group-level classifier-based searchlight analysis (linear svm). This is certainly a needed topic; as pointed out in the article, the assumptions behind t-tests are certainly violated in searchlight analysis, and using the binomial is also problematic (they suggest that it is too lenient, which strikes me as plausible).

Here's my interpretation of what they propose:

1. Generate 100 permuted searchlight maps for each person. You could think of all the possible label (i.e. class, stimulus type, whatever you're classifying) rearrangements as forming a very large pool. Pick 100 different rearrangements for each person and do the searchlight analysis with this rearrangement. (The permuted searchlight analysis must be done exactly as the real one was - same cross-validation scheme, etc.)
2. Generate 100,000 averaged group searchlight maps. Each group map is made by picking one permuted map from each person (out of the 100 made for each person in step 1) and averaging the values voxel-wise. In other words, stratified sampling with replacement.
3. Do a permutation test at each voxel, calculating the accuracy corresponding to a p = 0.001 threshold. In other words, at each voxel you record the 100th biggest accuracy after sorting the 100,000 accuracies generated in step 2. (100/100000 = 0.001)
4. Threshold the 100,000 permuted group maps and the one real-labeled group map using the voxel-wise thresholds calculated in step 3. Now the group maps are binary (pass the threshold or not).
5. Apply a clustering algorithm to all the group maps. They clustered voxels only if they shared a face. I don't think they used a minimum cluster size, but rather called un-connected voxels clusters of size 1 voxel. (This isn't really clear to me.)
6. Count the number of clusters by size in each of the 100,000 permuted maps and 1 real map. (this gives counts like 10 clusters with 30 voxels in map #2004, etc.)
7. Generate the significance of the real map's clusters using the counts made in step 6. I think they calculated the significance for each cluster size separately then did FDR, but it's not obvious to me ("Cluster-size statistics" section towards end of "Materials and Methods").
8. Done! The voxels passing step 7 are significant at the cluster level, corrected for multiple comparisons (Figure 3F of paper). The step 4 threshold map can be used for uncorrected p-values (Figure 3E of paper).

Most of this strikes me as quite reasonable. I've actually previously implemented almost this exact procedure (minus the cluster thresholding) on a searchlight dataset (not linear svms).

The part that makes me twitch the most is step 2: turning the 100 maps for each person into 100,000 group-average maps. I've been wanting to post about this anyway in the context of my ROI-based permutation testing example. But in brief, what makes me uncomfortable is the way 100 maps turn into 100000. Why not just calculate 5 for each person?  5^12 >> 100,000 (they had 12 subjects in some of the examples).  Somehow 100 for each person feels more properly random than 5 for each person, but how many are really needed to properly estimate the variation? I will expand on this more (and give a few alternatives), hopefully somewhat soon.

 The other thing that makes me wonder is the leniency. They show (e.g. Figure 11) that many more voxels are called significant in their method than with a t-test, claiming that as closer to the truth. This relates to my concern about how to combine over subjects: using 100,000 group maps allows very small p-values. But if the 100,000 aren't as variable as they should be, the p-values will be inflated.


Stelzer, J., Chen, Y., & Turner, R. (2012). Statistical inference and multiple testing correction in classification-based multi-voxel pattern analysis (MVPA): Random permutations and cluster size control. NeuroImage DOI: 10.1016/j.neuroimage.2012.09.063



UPDATE (30 October): We discussed this paper in a journal club and a coworker explained that the authors do explain the choice of 100 permutations per person in Figure 8 and the section "Undersampling of the permutation space". They made a dataset with one searchlight and many examples (80, 120, 160), then varied the number of permutations they calculated for each individual (10, 100, 1000, 10,000). They then made 100,000 group "maps" as before (my step 2), drawing from each group of single-subject permutations. Figure 8 shows the resulting histograms: the curves for 100, 1000, and 10,000 individual permutations are quite similar, which they use as rationale for running 100 permutations for each person (my step 1).

I agree that this is a reasonable way to choose a number of each-person permutations, but I'm still not entirely comfortable with the way different permutation maps are combined. I'll explain and show this more in a separate post.

searchlight shapes: BrainVoyager

Rainer Goebel kindly provided a description and images of the searchlight creation used in BrainVoyager:

"BrainVoyager uses the "sphere" approach (as in our original PNAS paper Kriegeskorte, Goebel, Bandettini 2006), i.e. voxels are considered in a cube neighborhood defined by the radius and in this neighborhood only those voxels are included in the "sphere" that have an Euclidean distance from the center of less than (or equal to) the radius of the sphere. From your blog, I think the resulting shapes in BrainVoyager are the same as for pyMVPA.

Note, however, that in BrainVoyager the radius is a float value (not an integer) and this allows to create "spheres" where the center layer has a single element on each side at cardinal axes (e.g. with radius 1, 2, 3, 4... voxels, see snapshot below) but also "compact" spheres as you seem to have used by setting the radius, e.g. to 1.6, 1.8, 2.6, 2.8, 3.6...). "

At right is an image Rainer generated showing  radius 1.0 and 2.0 searchlights created in BrainVoyager.

I am intrigued by Rainer's comment that using non-integer radii will make more "compact" spheres; non-interger radii underscores the need to be explicit in describing searchlight shape in methods sections. It appears that pyMVPA requires integer radii, but the Princeton MVPA Toolbox does not.

random searchlight averaging

I just finished reading several papers by Malin Björnsdotter describing her randomly averaged searchlight method and am very excited to try it.

In brief, instead of doing an exhaustive searchlight analysis (a searchlight centered on every voxel), searchlights are positioned randomly, such that each voxel is in just one searchlight. The parcellation is repeated, so that each voxel is included in several different searchlights (it might be in the center of one, touching the center of another, in the edge of a third, etc.). You then create the final information map by assigning each voxel the average of all the searchlights it's a part of - not just the single searchlight it was the center of as in normal searchlight analysis.

Malin introduces the technique as a way to speed searchlight analysis: far fewer searchlights need to be calculated to achieve good results. That's certainly a benefit, as searchlight analyses can take a ridiculous amount of time, even distributed across a cluster.

But I am particularly keen on the idea of assigning each voxel the average accuracy of all the searchlights it's a part of, rather than the single centered searchlight accuracy. In this and other papers Malin shows that the Monte Carlo version of searchlight analysis produces more accurate maps than standard searchlighting, speculating that the averaging reduces some of the spurious results. I find this quite plausible: a single searchlight may be quite accurate at random; the chance high accuracy will be "diluted" by averaging with the overlapping searchlights. Averaging won't eliminate all distortions in the results of course, but might help mitigate some of the extremes.

Malin points out that averaged information maps can be constructed after doing an exhaustive searchlight mapping, not just a Monte Carlo one. I will certainly be calculating some of these with my datasets in the near future.


Björnsdotter M, Rylander K, & Wessberg J (2011). A Monte Carlo method for locally multivariate brain mapping. NeuroImage, 56 (2), 508-16 PMID: 20674749

searchlight shapes: searchmight

The searchmight toolbox was developed by Francisco Pereira and the Botvinick Lab and introduced in Information mapping with pattern classifiers. The toolbox is specialized for searchlight analysis, streamlining and speeding the process of trying different classification algorithms on the same datasets.

Francisco said that, by default, the toolbox uses cubical searchlights (other shapes can be manually defined). In the image of different one-voxel radius searchlights, this is the one that I labeled "edge, face, or corner": 26 surrounding voxels for a one-voxel radius searchlight. Larger radii are also cubes.


Pereira, Francisco, & Botvinick, Matthew (2011). Information mapping with pattern classifiers: A comparative study. NeuroImage. DOI: 10.1016/j.neuroimage.2010.05.026

pyMVPA searchlight shape: solved

Michael Hanke posted further description of the pyMVPA searchlight: it is the same as the Princeton MVPA toolbox searchlight I described earlier.

As Michael points out, it is possible to define other searchlight shapes in pyMVPA; this is the default for a sphere, not the sole option.

searchlight shapes: a collection

The topic of searchlight shapes is more complex than I'd originally expected. I'd like to compile a collection of what people are using. I added a "searchlight shapes" label to the relevant posts: clicking that label should bring up all of the entries.

I have contacted some people directly in the hopes of obtaining sufficient details to add their searchlight shapes to the collection. If you have details of a searchlight not included, please send them and I'll add it (or arrange for you to post them). My initial goal is to describe the main (i.e. default) spherical searchlights produced in the most widely used code/programs, but I also plan to include other types of searchlights, like the random ones created by Malin Björnsdotter.

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.