Wednesday, January 29, 2014

demo: R code to perform a ROI-based MVPA

In the previous post I shared R code to perform a searchlight analysis on a sample dataset; this post describes R code to perform a ROI-based analysis on the dataset. In addition to the code and dataset files linked in the previous sentence, the code requires a ROI mask, shown at left.

This code is essentially a simpler version of the searchlight code: instead of reading out which voxels to classify from the lookup and 3d files it gets the voxel coordinates from the ROI mask. The key part of this code is that the ROI mask and functional data must be fully aligned - as the images are read into R. This can easily go awry; I strongly suggest you spot-check several voxels in R and an outside program (such as MRIcroN) to make sure the laterality of all images agree.

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

Saturday, January 25, 2014

"tough love for fMRI"

I lifted the title from a great post on Chris Chambers' blog NeuroChambers: "Tough love for fMRI: questions and possible solutions". He has a very succinct summary of a lot of the problems the field has been struggling with: "MRI = v expensive method + chronically under-powered designs + intense publication pressure + lack of data sharing = huge fraud incentive." And not just an incentive to outright fraud, but also to practices that fall on the spectrum of questionable research practices.

I happened to read Chris Chambers' post the same day I read "Cluster-extent based thresholding in fMRI analyses: Pitfalls and recommendations", a new NeuroImage paper by Woo, Krishnan, and Wager, and I must say the combination is disheartening.

The paper describes mass-univariate analysis interpretation, but there are a lot of parallels between the problems in spatial specificity they point out and those that can occur in searchlight or ROI-based MVPA. In short: it is not appropriate to interpret sub-regions when an analysis produces a large cluster that spans several regions; some of the cluster voxels are significant, but you can't say which ones. You can't say that an area is "activated" simply because it is part of a larger activated area that passed a cluster threshold - it could be that only adjacent regions (other parts of the cluster) are actually activated.

The authors make several practical suggestions, one of which is to individually color-code each cluster that passed the cluster-extent thresholding (so you'd plot a group of blue voxels, another of red voxels, etc), rather than the common practice of plotting all surviving voxels with the same heatmap color scale. They further recommend "that figure legends and captions explicitly state that the true activation location and extent within each significant cluster cannot be determined." This pushes interpretation (properly) to the cluster level, paralleling  ROI-based MVPA, which presents descriptions of the properties of ROIs, not the voxels they contain.

This (and other advice in the paper) is very practical, but I suspect people will resist: it is a lot harder to come up with a compelling story (and get a blob in the "right spot") when clusters are properly interpreted. But this is the sort of "tough love" that is needed for fMRI.

ResearchBlogging.orgWoo CW, Krishnan A, & Wager TD (2014). Cluster-extent based thresholding in fMRI analyses: Pitfalls and recommendations. NeuroImage PMID: 24412399