Runs cluster correction on the output from MDMR analysis using connectir (http://github.com/czarrar/connectir). This can also convert ROI data to voxelwise data.

inline_connectir.R

# This file will be called the `le_correcter.R` file shown below
# and both files must be in the same directory.
# This file contains the low-level functions to find the clusters
# as well as convert ROIs to voxelwise data.

suppressPackageStartupMessages(library(connectir))
suppressPackageStartupMessages(library(inline))

plugin_bigmemory <- function() {     
   l <- getPlugin("RcppArmadillo")     
      
   l$includes <- paste(l$includes, '   
 #include "bigmemory/BigMatrix.h"     
 #include "bigmemory/MatrixAccessor.hpp"     
 #include "bigmemory/bigmemoryDefines.h"     
 #include "bigmemory/isna.hpp"     
    ')     
    
    l$LinkingTo <- c("bigmemory", "BH", l$LinkingTo)     
      
    l$Depends <- c("bigmemory", l$Depends)     
      
    return(l)     
} 
registerPlugin("bigmemory", plugin_bigmemory)     

plugin_connectir <- function() {     
   l <- getPlugin("bigmemory")     
      
   l$includes <- paste(l$includes, '     
 #include "connectir/connectir.h"     
    ')     
    
    l$LinkingTo <- c("connectir", l$LinkingTo)     
      
    l$Depends <- c("connectir", l$Depends)     
      
    return(l)     
} 
registerPlugin("connectir", plugin_connectir)     

# Convert roi-level data to be voxelwise
rois_to_voxelwise <- cxxfunction( signature(sA = "object", sRow = "numeric", sRoiInds='List', sNvoxs = "numeric"), 
' 
  try{
    BM_TO_ARMA_ONCE(sA, A)
    double row = DOUBLE_DATA(sRow)[0] - 1;
    Rcpp::List list_roi_inds(sRoiInds);
    double nvoxs = DOUBLE_DATA(sNvoxs)[0];
    
    arma::vec voxs(nvoxs);
    for (size_t i = 0; i < A.n_cols; ++i) {
        SEXP s_roi_inds = list_roi_inds[i];
        Rcpp::NumericVector roi_inds(s_roi_inds);
        int ninds = roi_inds.size();
        for (int j = 0; j < ninds; ++j)
            voxs[roi_inds[j]-1] = A(row,i);
    }
    
    return Rcpp::wrap( voxs );
  } catch( std::exception &ex ) {
    forward_exception_to_r( ex );
  } catch(...) { 
    ::Rf_error( "c++ exception (unknown reason)" ); 
  }
  return R_NilValue; // -Wall
', plugin = "connectir")

get_row <- cxxfunction( signature(sA = "object", sRow = "numeric"), 
' 
  try{
    BM_TO_ARMA_ONCE(sA, A)
    double row = DOUBLE_DATA(sRow)[0] - 1;
    double nvoxs = A.n_cols;
    
    arma::vec voxs(nvoxs);
    for (size_t i = 0; i < A.n_cols; ++i) {
        voxs[i] = A(row,i);
    }
    
    return Rcpp::wrap( voxs );
  } catch( std::exception &ex ) {
    forward_exception_to_r( ex );
  } catch(...) { 
    ::Rf_error( "c++ exception (unknown reason)" ); 
  }
  return R_NilValue; // -Wall
', plugin = "connectir")


# Cluster Table
cluster_table_src = '
using namespace arma;

vec all_inds = Rcpp::as<vec>(R_all_inds);
vec voxs_to_use = Rcpp::as<vec>(R_voxs_to_use);
vec clust_vals = Rcpp::as<vec>(R_clust_vals);
vec raw_vals = Rcpp::as<vec>(R_raw_vals);
vec offset_inds = Rcpp::as<vec>(R_offset_inds);
Rcpp::NumericVector C_last_ind(R_last_ind);
double last_ind = C_last_ind[0];

double clust_i, clust_num, clust_size, clust_mass = 0;
int c_ind, ind;

size_t n_voxs = all_inds.n_elem;
size_t n_inds = offset_inds.n_elem;
vec nei_inds(n_inds);
mat center_inds(0,1);

mat sizes(0,1); mat masses(0,1);

for(size_t i = 0; i < n_voxs; ++i)
{
    c_ind = all_inds(i);
    
    if (voxs_to_use(c_ind) == 0) continue;
    
    //printf("\\nnum: %.2f\\n", clust_num);
    
    // One less voxel
    voxs_to_use(c_ind) = 0;
    
    // Initialize number, size, and mass
    clust_num = clust_num + 1;
    clust_size = 1;
    clust_mass = raw_vals(c_ind);
    
    //printf("size: %.2f\\n", clust_size);
    
    // Save center index
    center_inds.insert_rows(0, 1, false);
    center_inds(0,0) = c_ind;
    
    // Update cluster voxel val
    clust_vals(c_ind) = clust_num;
    
    while (center_inds.n_elem > 0) 
    {
        // Get neighbors
        c_ind = center_inds(0);
        nei_inds = c_ind + offset_inds;
    
        // Loop through each neighbor
        for(size_t j = 0; j < n_inds; ++j)
        {
            ind = nei_inds(j);
    
            // TODO: test speed improvement with [i]
            if (ind > 0 && ind < last_ind) {
                if (voxs_to_use(ind) == 1)
                {
                    // Save index to later look at its neighors
                    center_inds.insert_rows(1, 1, false);
                    center_inds(1,0) = ind;
            
                    // One less voxel
                    voxs_to_use(ind) = 0;
                
                    // Update size, mass, and cluster voxel val
                    clust_size = clust_size + 1;
                    clust_mass = clust_mass + raw_vals(ind);
                    clust_vals(ind) = clust_num;
                }
            }
        }

        // Get rid of the center ind that started this
        center_inds.shed_row(0);
    }
    
    clust_i = clust_num - 1;
    
    sizes.insert_rows(clust_i, 1, false);
    masses.insert_rows(clust_i, 1, false);
    
    sizes(clust_i,0) = clust_size;
    masses(clust_i,0) = clust_mass;
    
    //printf("size: %.1f\\n", clust_size);
}

return Rcpp::List::create(Rcpp::Named("n") = Rcpp::wrap( clust_num ),
                          Rcpp::Named("sizes") = Rcpp::wrap( sizes ),
                          Rcpp::Named("masses") = Rcpp::wrap( masses ), 
                          Rcpp::Named("clust") = Rcpp::wrap( clust_vals ));
'

inline_cluster_table <- cxxfunction( 
    signature(R_all_inds = "numeric", R_voxs_to_use = "numeric", 
              R_clust_vals = "numeric", R_raw_vals = "numeric", 
              R_offset_inds = "numeric", R_last_ind = "numeric"), 
    body = cluster_table_src, 
    plugin = "RcppArmadillo"
)


# New cluster table making use of c++ inline function
cpp.cluster.table <- function(x, vox.thr=0, dims=NULL, mask=NULL, 
                              nei=1, nei.dist=3, pad=1) 
{
    # TEMPORARY
    #x <- vox.fstats[,1]; vox.thr <- 1.5; dims <- hdr$dim; nei <- 1; nei.dist <- 3; pad <- 1
    #orig_mask <- mask
    # TEMPORARY
    
    if (is.null(mask))
        mask <- rep(T, length(x))
    if (is.null(dims)) {
        if (length(dim(x)) != 3)
            stop("If dims isn't provided, then x must be a 3D array")
        dims <- dim(x)
    }
    
    nx <- dims[1]; ny <- dims[2]; nz <- dims[3]
    
    # This pads the mask with zero-voxels around the image
    # note that x is only the non-zero voxels in the mask
    if (pad > 0) {
        mask.pad <- array(F, c(nx+pad*2, ny+pad*2, nz+pad*2))
        mask.pad[(pad+1):(pad+nx), (pad+1):(pad+ny), (pad+1):(pad+nz)] <- T
        
        # New dimensions
        nx <- dim(mask.pad)[1]; ny <- dim(mask.pad)[2]; nz <- dim(mask.pad)[3]
        
        # Get new mask
        mask.pad <- as.vector(mask.pad)
        tmp <- rep(F, length(mask.pad))
        tmp[mask.pad] <- mask
        mask <- tmp
        rm(tmp)
    }
    
    # Vector version of the data, mask, and clusters
    xfull <- vector("numeric", nx*ny*nz)
    mfull <- vector("logical", nx*ny*nz)
    cfull <- mfull*0
    
    # Get neighbours to check for clusters
    # default is 27 (face, edge, corner touching)
    nmat <- expand.grid(list(i=-nei:nei, j=-nei:nei, k=-nei:nei))
    voxdists <- rowSums(abs(nmat))
    nmat <- nmat[voxdists<=nei.dist,]
    offset.inds <- nmat$k*nx*ny + nmat$j*nx + nmat$i    # neighbors relative to center node
    rm(nmat, voxdists)  # only keep offset.inds
    
    # Threshold data and save as new mask
    tmp.mask <- x > vox.thr
    nvoxs <- sum(tmp.mask)
    mfull[mask] <- tmp.mask
    xfull[mask][tmp.mask] <- x[tmp.mask]
    rm(tmp.mask)
    
    last.ind <- length(mfull) + 1   # index w/ end of the data
    
    all.inds <- which(mfull)
    
    ct <- inline_cluster_table(all.inds, mfull, cfull, xfull, 
                               offset.inds, as.double(last.ind))
    
    cfull <- ct$clust
    if (pad > 0)
        cfull <- cfull[mask.pad]
    dim(cfull) <- dims
    
    if (length(ct$sizes) == 0) {
        ct$n <- 0
        ct$sizes <- matrix(0, 1, 1)
        ct$masses <- matrix(0, 1, 1)
    }
    
    list(
        nclust=ct$n, 
        max.size=max(ct$sizes), 
        max.mass=max(ct$masses), 
        size=ct$sizes[,1],
        mass=ct$masses[,1],
        clust=cfull
    )
}

le_correcter.R

#!/usr/bin/env Rscript

# This script will
# - take the subject distance and mdmr directory as input
# - if needed, convert roi data to voxelwise
#
# - create the unthresh logp maps (15k perms)
# - create the unthres fdr-corrected logp maps
# - run cluster correction
#   - determine fstat threshold (5k perms)
#   - determine cluster thresholds (5k perms)
#   - create the clusterized results
#
# It must be in the same directory as the above script `inline_connectir.R`
#
# Usage:
# ./le_correcter.R ${distance_directory} ${mdmr_directory}
#
# Where distance_directory is the full path to the subject distances directory 
# and mdmr_directory is the full path to the MDMR directory.
#
# This will create a new directory called cluster_correct_v05_c05 
# (voxel level threshold of p < 0.05 and cluster level threshold of p < 0.05) 
# within the mdmr directory you input.
# This will have a clust_logp_X.nii.gz file that is the thresholded -logp file with clusters corrected based on a permutation test.
# There is also an easythresh directory that runs everything with easythresh as well. It will run easythresh first since it is very fast. 
# In both permutation and easythresh outputs, the voxel values represent significance estimates from permutation testing in the case of the permutation output and using GRF for easythresh. 
#
# For now to change the threshold, you need to manually modify the file and specifically the following lines:
# vox.thresh <- 0.05
# clust.thresh <- 0.05

cat("Loading necessary libraries and functions\n")
source("inline_connectir.R")    # loads connectir

# Take the subject distance and mdmr directories as input
# the roi file is optional

args <- commandArgs(trailingOnly = TRUE)

if (length(args) < 2) {
    stop("usage: correct.R distance-dir mdmr-dir [roi-file]")
}

dist.dir <- args[1]
mdmr.dir <- args[2]
roi.file <- args[3] # optional

roify <- !is.na(roi.file)

## fixed but extra options
### overwrite output if it exists?
overwrite <- FALSE
### voxel and cluster thresholds
vox.thresh <- 0.05
clust.thresh <- 0.05
### # of perms for cluster parts
### should include non-permuted data
nperms4clust <- 5000

# Output Directory
vstr <- sub("0.", "", as.character(vox.thresh))
cstr <- sub("0.", "", as.character(clust.thresh))
outdir <- file.path(mdmr.dir, sprintf("cluster_correct_v%s_c%s", vstr, cstr))

if (file.exists(outdir)) {
    vstop("output directory '%s' already exists", outdir)
} else {
    dir.create(outdir)
}


# Checks
vcat(T, "\nChecks")

## check if ROI file is needed
vcat(T, "...are rois needed")
if (is.na(roi.file) && file.exists(file.path(dist.dir, "mask.txt"))) {
    stop("ROI file needed but not given")
}


# Setup
vcat(T, "\nSetup")

## vector of permuted factor names
vcat(T, "...factor names")
load(file.path(mdmr.dir, "modelinfo.rda"))
factors <- names(attr(modelinfo$qrhs, "factors2perm"))
nfactors <- length(factors)

## p-value matrix from all permutations
vcat(T, "...p-values from all permutations")
f <- file.path(mdmr.dir, "pvals.desc")
pvals.mat <- attach.big.matrix(f)

## list of permuted matrices
vcat(T, "...load permutations")
files <- file.path(mdmr.dir, sprintf("fperms_%s.desc", factors))
list.fperms <- lapply(files, function(f) attach.big.matrix(f))

## determine number of permutations
vcat(T, "...determine number of permutations for clustering")
nperms <- nrow(list.fperms[[1]])
if (nperms < nperms4clust) {
    vcat(T, "...changing nperms4cluster to %i", nperms)
    nperms4clust <- nperms
} else {
    vcat(T, "...will use %ik permutations for clustering", nperms4clust)
}

## subset of permuted matrices for clustering
vcat(T, "...take subset of permutations for clustering")
tmp <- lapply(list.fperms, function(fperms) {
    deepcopy(fperms, rows=1:nperms4clust, shared=FALSE)
})
rm(list.fperms); invisible(gc(F,T))
list.fperms <- tmp

## new p-values for subset
vcat(T, "...calculate p-values for new subset of permutations")
pvals4clust.mat <- mdmr.fstats_to_pvals(list.fperms)


# ROI stuff
if (roify) {
    vcat(T, "ROI Setup")
    
    ## read in rois and create mask
    vcat(T, "...read in ROIs and create mask")
    rois <- read.nifti.image(roi.file)
    hdr  <- read.nifti.header(roi.file)
    mask <- as.vector(rois!=0)
    rois <- rois[mask]
    nvoxs <- sum(mask)

    ## get unique rois and related indices
    vcat(T, "...determine unique rois and indices")
    urois <- sort(unique(rois))
    urois <- urois[urois!=0]
    nrois <- length(urois)
    rois.inds <- lapply(urois, function(ur) which(rois==ur))

    ## simple method for rois => voxelwise 
    simple.rois2voxelwise <- function(roi.data, vox.rois) {
        vox.data <- vector("numeric", length(vox.rois))
    
        urois <- sort(unique(vox.rois))
        urois <- urois[urois!=0]
        nrois <- length(urois)
    
        for (ri in 1:nrois)
            vox.data[vox.rois==urois[ri]] <- roi.data[ri]
    
        return(vox.data)
    }
} else {
    vcat(T, "...read mask")
    
    maskfile <- file.path(dist.dir, "mask.nii.gz")
    hdr <- read.nifti.header(maskfile)
    mask <- read.mask(maskfile)
    
    rois.inds <- NULL
}

fstat_from_pval <- function(fstats, pvals, pthr) {
    df <- data.frame(fstats=fvals, logps=-log10(pvals))
    model <- lm(fstats ~ logps, data=df)
    fthr <- predict(model, data.frame(logps=-log10(pthr)))
    return(fthr)
}

permute_max_sizes <- function(fperms, fthr, hdr, mask, rois.inds=NULL) {
    roify <- !is.null(rois.inds)
    nperms <- nrow(fperms)
    nvoxs <- sum(mask)
    
    max.sizes <- laply(1:nperms, function(i) {
        if (roify) {
            img <- rois_to_voxelwise(fperms, as.double(i), rois.inds, as.double(nvoxs))
        } else {
            img <- get_row(fperms, as.double(i))
        }
        ct <- cpp.cluster.table(img, fthr, hdr$dim, mask)
        ct$max.size
    }, .progress="text")
    
    return(max.sizes)
}


# The Meat
for (fi in 1:nfactors) {
    factor <- factors[fi]
    vcat(T, "\nFactor: %s", factor)
    
    
    ## Create unthresholded logp maps
    vcat(T, "...unthresholded logp and fdr logp maps")
    
    pvals       <- -log10(pvals.mat[,fi])
    fdr.pvals   <- -log10(p.adjust(pvals.mat[,fi], "fdr"))
    
    if (roify) {
        pvals     <- simple.rois2voxelwise(pvals, rois)
        fdr.pvals <- simple.rois2voxelwise(fdr.pvals, rois)
    }
        
    outfile <- file.path(outdir, sprintf("logp_%s.nii.gz", factor))
    write.nifti(pvals, hdr, mask, outfile=outfile, overwrite=overwrite)
    
    outfile <- file.path(outdir, sprintf("fdr_logp_%s.nii.gz", factor))
    write.nifti(fdr.pvals, hdr, mask, outfile=outfile, overwrite=overwrite)
    
    
    ## Comparison Easy Thresh Clustering
    vcat(T, "...comparison with easythresh")
    
    easy.dir <- file.path(outdir, "easythresh")
    dir.create(easy.dir, FALSE)
    
    # Convert p-value to zstats
    vcat(T, "...pvals to zstats")
    
    zstats <- qt(pvals.mat[,fi], Inf, lower.tail=F)
    if (roify) zstats <- simple.rois2voxelwise(zstats, rois)
    
    outfile0 <- file.path(easy.dir, sprintf("zstat_%s_tmp.nii.gz", factor))
    write.nifti(zstats, hdr, mask, outfile=outfile0, overwrite=overwrite)
    
    # Fix to remove any infinite values
    curdir <- getwd()
    setwd(easy.dir)
    outfile <- file.path(easy.dir, sprintf("zstat_%s.nii.gz", factor))
    cmd <- sprintf("3dcalc -a %s -expr a -prefix %s", outfile0, outfile)
    vcat(T, cmd)
    system(cmd)
    file.remove(outfile0)
    setwd(curdir)
    
    # Brain mask
    if (roify) {
        # Save mask file
        vcat(T, "...saving mask file")
        maskfile <- file.path(outdir, "mask.nii.gz")
        write.nifti(mask, hdr, outfile=maskfile, overwrite=overwrite)
    } else {
        maskfile <- file.path(dist.dir, "mask.nii.gz")
    }
    
    # Background (underlay) image
    bgfile <- file.path(dist.dir, "bg_image.nii.gz")
    
    # Run it
    curdir <- getwd()
    setwd(easy.dir)
    zstatfile <- outfile
    zthr <- qt(vox.thresh, Inf, lower.tail=F)
    cmd <- sprintf("easythresh %s %s %.4f %.4f %s zstat_%s --mm", 
                    zstatfile, maskfile, zthr, clust.thresh, bgfile, factor)
    vcat(T, cmd)
    system(cmd)
    setwd(curdir)
    
        
    # Cluster correct
    vcat(T, "...Cluster Correcting")
    
    fperms <- list.fperms[[fi]]
    
    # Determine pseudo-F given a particular p-value
    # this would be the voxel threshold before clustering
    vcat(T, "...determine pseudo-F threshold with a p=%.4f", vox.thresh)
    
    fvals <- fperms[1,]
    pvals <- pvals4clust.mat[,fi]
    fthr <- fstat_from_pval(fvals, pvals, vox.thresh)
    
    
    # Determine clusters from non-permuted data
    vcat(T, "...clusters in non-permuted data")
    
    if (roify) {
        img <- simple.rois2voxelwise(fvals, rois)
    } else {
        img <- fvals
    }
    orig.ct <- cluster.table(img, fthr, hdr$dim, mask)
    
    # Maximum cluster sizes across permutations
    vcat(T, "...max clust sizes across permutations")
    
    max.sizes <- permute_max_sizes(fperms, fthr, hdr, mask, rois.inds)
    
    # Cluster significance
    vcat(T, "...determine significant clusters")
    
    clust.sig <- sapply(orig.ct$size, function(ocs) {
        sapply(ocs, function(s) sum(s<max.sizes)/nperms4clust)
    })
    
    clust <- orig.ct$clust[mask]
    w.clusts <- which(rev(clust.sig<clust.thresh))
    clust.keep <- clust*0    # empty vector
    for (i in 1:length(w.clusts)) clust.keep[clust==w.clusts[i]] <- 1
    
    logps <- -log10(pvals.mat[,fi])
    if (roify) logps <- simple.rois2voxelwise(logps, rois)
    clust.logps <- logps * clust.keep
    
    # Save
    vcat(T, "...saving")
    
    outfile <- file.path(outdir, sprintf("clust_%s.nii.gz", factor))
    write.nifti(clust.keep, hdr, mask, outfile=outfile, overwrite=overwrite)
    
    outfile <- file.path(outdir, sprintf("clust_logp_%s.nii.gz", factor))
    write.nifti(clust.logps, hdr, mask, outfile=outfile, overwrite=overwrite)    
    
    outfile <- file.path(outdir, sprintf("clustinfo_%s.rda", factor))
    save(orig.ct, max.sizes, clust.sig, fthr, fvals, pvals, img, logps, clust.keep, 
         file=outfile)
}