Instantly share code, notes, and snippets.

Embed
What would you like to do?
A function for "continuous" family-wise error rate correction
#Thresholding for VLSM using permutation based statistics
#Packages that need to be installed and loaded before the function can be run
# ANTsR
# gtools
# plyr
#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 <- permute(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)})
#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)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment