Bootstrap Analysis of CWAS

cwas_bootstrap.R

library(connectir)
library(boot)

# data = model
boot_mdmr <- function(formula, data, indices, sdist, factors2perm) {
    ###
    # Distances
    ###
    
    # We need to sample the distances based on the indices
    # that is we will create a new set of distances with subject indices based on indices
    # This will also create a local copy of the big matrix
    cat("Subset of subjects in distances\n")
    sdist <- filter_subdist(sdist, subs=indices)
    
    # Now we can gowerify
    cat("Gowerify\n")
    gmat <- gower.subdist2(sdist)
    
    # Info
    nvoxs <- ncol(gmat)
    nsubs <- sqrt(nrow(gmat))
    nperms <- 4999
    nfactors <- 1
    
    
    ###
    # Calculate memory demands
    ###
    opts <- list(verbose=TRUE, memlimit=20, blocksize=0, superblocksize=0)
    opts <- get_mdmr_memlimit(opts, nsubs, nvoxs, nperms, nfactors)
    # this will give opts$blocksize and opts$superblocksize that will
    # limit total RAM usage to 20GB
    
    ###
    # Get the model ready
    ###
    
    cat("Subset of subjects in model\n")
    model <- data.frame(data[indices,])
    
    
    ###
    # Call MDMR
    ###
    
    ret <- mdmr(gmat, formula, model, nperms, factors2perm, 
                 superblocksize=opts$superblocksize, blocksize=opts$blocksize)
    
    -log10(ret$pvals[,]) # or should I just return p-values?
}


###
# Bootstrap ROI-based IQ Results
###

# Set parallel processing
nthreads <- 8
set_parallel_procs(1, nthreads, TRUE)

# Read in the distances (using ROI-based distances and not voxelwise here)
#dpath <- "/home2/data/Projects/CWAS/nki/cwas/short/compcor_kvoxs_fwhm08_to_kvoxs_fwhm08/subdist.desc"
dpath <- "/home2/data/Projects/CWAS/nki/cwas/short/compcor_only_rois_random_k0800/subdist.desc"
sdist <- attach.big.matrix(dpath)

# Read in the model
mpath <- "/home2/data/Projects/CWAS/share/nki/subinfo/40_Set1_N104/subject_info_with_iq_and_gcors.csv"
model <- read.csv(mpath)
model <- subset(model, select=c("FSIQ", "Age", "Sex", "short_meanFD"))

# Set the formula
f     <- ~ FSIQ + Age + Sex + short_meanFD

# Call
results <- boot(data=model, statistic=boot_mdmr, R=500, formula=f, sdist=sdist, factors2perm="FSIQ")

This analysis is the first steps to dealing with a reviewer's suggestion for our CWAS paper.