function for calculating continuous FWER correction for LSM (see http://dx.doi.org/10.1016/j.neuropsychologia.2017.08.025)

cFWER.R

#Thresholding for VLSM using permutation based statistics 

# ANTsR needs to be installed and loaded before the function can be run

#cFWER takes the following arguments 
#   lesmat - a matrix containing the lesion data from the patient scans 
#   behaveData - a vector containing only the behavioral data
#   numPerm - a single integer argument that tells the function how many permutations to run
#   voxThreshs - a vector containing the voxel numbers for thresholding 
#   flip - a bolean flag telling the function whether or not to flip the t values output from the 
#         vlsm regression 
#   sortOrder - boolean flag that tells function to sort output tVals from high to low or low to high

cFWER <- function(lesMat, behaveData, numPerm=1000, voxThreshs=c(1,10,100,1000), flip=FALSE, sortHtoL=TRUE){
  
  #run a basic vlsm regression 
  base_vlsm <- bigLMStats(lm(lesMat ~ behaveData))
  
  #if flip is set to TRUE it flips the values from the vlsm
  if(flip == TRUE){
    base_vlsm$beta.t = -(base_vlsm$beta.t)
  }
  
  #preallocate space for the t val matrix
  tValMatrix <-  matrix(nrow=length(voxThreshs), ncol=numPerm)
  
  # run permutations, get voxThreshs t-values from each one
  for (i in 1:numPerm){
    #permute the behavioral scores
    scores <- sample(behaveData)
    
    #run the vlsm
    perm_vlsm <- bigLMStats(lm(lesMat ~ scores))
    
    if(flip == TRUE){
      #flip the t values 
      perm_vlsm$beta.t = -(perm_vlsm$beta.t)
    }
    
    #sort the t values in the order specified by sortHtoL 
    sortedTvals <- sort(perm_vlsm$beta.t, decreasing = sortHtoL)
    
    #get the t-values values for the given voxel numbers and store them
    tValMatrix[,i] <- sortedTvals[voxThreshs] 
    
    print(paste("Finished permutation:",i))
  }
  
  #Compute the 95% t-threshold value for each value of voxThresh
  # i.e., each row in the tValMatrix 
  #threshVals <- aaply(tValMatrix, 1, .fun=function(x){quantile(x, 0.95)})
  threshVals <- apply(tValMatrix, 1, function(x) quantile(x, 0.95))
  
  #count how many voxels in original VLSM passed each threshold 
  baseVLSM_numVoxPass <- matrix(nrow=length(voxThreshs))
  for(j in 1:length(voxThreshs)){
    if(sortHtoL == TRUE){
      baseVLSM_numVoxPass[j] <- sum(base_vlsm$beta.t > threshVals[j])
    }else{
      baseVLSM_numVoxPass[j] <- sum(base_vlsm$beta.t < threshVals[j])
    }
    
  }
  
  #creates base threshData table  
  threshData <- data.frame(
    numVoxels = voxThreshs,
    permThresh.95 = threshVals,
    numVoxPassThresh = baseVLSM_numVoxPass)
  
  return(threshData)
}