Thursday, September 18, 2014

demo: R code to perform a voxelwise t-test

It's easy to perform a voxelwise t-test (t-test at each voxel individually). Programs for mass-univariate analysis (like SPM and fsl) of course do this (and much more), but sometimes you just want to do a simple t-test across subjects at each voxel.

The demo code linked in this post does a voxelwise t-test in R. It takes as input a set of 3d NIfTI files, where each NIfTI is assumed to come from a different person, each voxel of which contains a statistic describing effect strength (for example, accuracy resulting from a searchlight analysis). The code reads the 3d NIfTI images into a 4d array (people in the fourth dimension), then performs a t-test at each voxel, saving the t-values for each voxel in a new NIfTI image. This figure shows the t-value image produced by the demo code.

The demo R code together with the input images to run the demo are available here, and the R code alone here.

Here is the key part of the code. The 4d array big contains the 3d statistic images for each person (subjects are the 4th dimension of the array). Then the plyr function aaply calls my little getT function (below and in the code file), which calculates the t-value at each voxel (ie over the subjects), creating a 3d array of t-values.

 # function to calculate the t-test at each voxel and return the t value  
 getT <- function(x) {    
  # can't do a t-test if variance is zero, so check before trying.  
  if (var(x) == 0) {   
   stat <- NA;   
  } else {   
   stat <- t.test(x, alternative="greater", mu=0.5)$statistic;   
  }  
  return(stat)  
 }  
 # plyr function aaply calls getT at each voxel in the 4d array big,  
 # creating a 3d output array of t-values (what getT returns)  
 t.img <- aaply(big, c(1,2,3), getT);   

No comments:

Post a Comment