Installing the pbj package

The latest stable version of the pbj package can be installed from github by specifying the master branch. Currently, the longitudinal branch includes some formatting changes and has edits to accommodate longitudinal data.

# install the latest versions of the packages to perform these analyses.
devtools::install_github('simonvandekar/pbj', ref='longitudinal')
#devtools::install_github('statimagcoll/NIsim')

Setup analysis files

In this section we walk through setting up the data for the group level image analysis. This first chunk of code loads the libraries needed to perform the group image analysis. The RNifti package is used to read and write Nifti images into R. The pbj package includes our functions needed to perform group level neuroimage analysis. The papayaWidget package is used to embed interactive visualizations into this html document. The ns package includes functions to fit natural cubic splines.

### LIBRARIES ###
library(pbj) # pbj package
library(RNifti) # Nifti I/O. Loaded with pbj
library(fslr) # for susan smoothing from FSL

## Loading required package: oro.nifti

## oro.nifti 0.11.0

## 
## Attaching package: 'oro.nifti'

## The following object is masked from 'package:RNifti':
## 
##     origin

## Loading required package: neurobase

library(papayaWidget) # image viewer
library(splines) # ns
library(magrittr) # %>%
library(callr) # used to run bootstraps in the background

After loading the libraries, we are ready to setup the files that we will need for the analysis. To perform the analysis, we need the following data

• The data file that includes

1. All covariates required for the analysis.
2. Directory paths to the .nii.gz VBM images, assumed to be spatially registered to template space.

• The mask to identify voxels to include in the analysis.
• A template image to use as the background for visualizing the results.

# sets path for data file
dbdatafile = '/media/disk2/pbj/data/rockland/demographic/RocklandBehavioral.csv'
# load in data file
dat = read.csv(dbdatafile, stringsAsFactors = TRUE)

# subset the data to relevant variables
dat = dat[,c('AnonymizedID', 'IDandSession', 'files', 'files2mm', 'files2mmsm8', 'age', 'sex')]

# View the top of the data file
# row.names, digits, and booktabs are optional arguments for displaying the
knitr::kable(head(dat[,c('AnonymizedID', 'IDandSession', 'age', 'sex')]), row.names = FALSE, digits=3, booktabs=T, caption='Covariates for first 10 participants.')

knitr::kable(head(dat[,'files2mm', drop=FALSE]), row.names = FALSE, digits=3, booktabs=T, caption='Nifti image paths for first 10 participants.')

# mask file (created May 12, 2021)
maskfile = "/media/disk2/pbj/data/rockland/neuroimaging/overlap_mask_2mm.nii.gz"

## CREATE MASK
# if(!file.exists(maskfile)){
#   imgs = readNifti(dat$files)
#   mask = imgs[[1]]
#   mask[,,] = 0
#   imgs = apply(simplify2array(imgs), 1:3, function(v) sum(v>0))
#   mask[sum(imgs)==nrow(dat)] = 1
#   writeNifti(mask, file=maskfile)
# }

The template image

The template image is included at this stage to make visualization seamless later on. In this case, all of the VBM data are in the MNI template space.

templatefile = '/usr/local/fsl/data/standard/MNI152_T1_2mm_brain.nii.gz'
# reads in image and creates lightbox view.
image(readNifti(templatefile))

Figure 1: Lightbox view of the template image. This is the default display for niftiImage objects.

# number of cores for parallel things
ncores = 16
# location to save results
pbjdatafile = '/media/disk2/pbj/pbj_ftest/nkirs_bootstrap_results.rdata'
permdatafile = '/media/disk2/pbj/pbj_ftest/nkirs_permutation_results.rdata'
# seed for randomization
set.seed(555)

Some image processing

The input VBM data are at 1mm isotropic resolution and unsmoothed. It is common to downsample and apply smoothing before group level analysis to speed computation, accommodate small regional differences in registration, and reduce noise. This code downsamples the images to 2mm isotropic resolution and applies smoothing at 8mm FWHM. It first checks whether the downsampled, smoothed output files exist; if not, it downsamples and applies smoothing using functions from the fslr package.

The second to last line creates an interactive viewer embedded into the html document using papayaWidget ot visualize the effect of the smoothing in the first participant. The last line sets the files variable to point to the downsampled smoothed data.

if(!all(file.exists(dat$files2mmsm8))){
  ## CREATE DOWNSAMPLED DATA
  invisible(mcmapply(flirt, infile=dat$files, outfile=dat$files2mm, opts = '-applyxfm', MoreArgs=list(reffile=templatefile, retimg=FALSE), mc.cores = ncores ))
  ## SMOOTH DOWNSAMPLED DATA
  # convert FWHM to sigma
  sigma = 8/2.355
  invisible(mcmapply(susan, file=dat$files2mm, outfile=dat$files2mmsm4, MoreArgs=list(sigma=sigma, dimg='3', n_usans='0'), mc.cores=ncores)) 
  # sigma = 8: using 8mm smoothing
}

# view files
# papaya(c(template, dat$files2mm[1], dat$files2mmsm4[1]))
# set downsampled smoothed files as the main file
papaya(c(templatefile, dat$files2mm[1], dat$files2mmsm8[1]))

Model specification

In order to study sex differences in the mean adult life-span curve of gray matter volume we fit the model \[ Y_i(v) = \alpha_0(v) + \alpha_1(v)\text{sex}_i + f_1(\text{age}_i, v) + \text{sex}_i \times f_2(\text{age}_i, v) + E_i(v), \] where \(Y_i(v)\) denotes the VBM image for participant \(i\) and voxel location \(v\). \(\text{sex}_i\) is an indicator variable that is equal to 1 if participant \(i\) is a male. We choose to fit \(f_1\) and \(f_2\) with natural cubic splines on 4 degrees of freedom. The unknown parameters \(\alpha_0(v)\), \(\alpha_1(v)\), \(f_1\), and \(f_2\) quantify the association between each variable and the VBM image at location \(v\).

To study sex differences in the life-span curve, we test the nonlinear age by sex interaction which corresponds to parameters that control the shape of the function \(f_2\). The null hypothesis is that \(f_2(\text{age}_i) = 0\) for all values of \(\text{age}_i\) and all voxel locations \(v\), versus the alternative hypothesis that this is not true, \[\begin{align*} H_0: & f_2(\text{age}, v) = 0 \\ H_1: & f_2(\text{age}) \ne 0. \end{align*}\] Because \(f_2\), is modeled with 4 degrees of freedom the test for \(H_0\) is an F-test at each voxel location, \(v\), on 4 degrees of freedom. This model specification is used to test the nonlinear age by sex interaction. The null hypothesis at each location is that males and females follow the same developmental trajectory.

pbj uses the full/reduced model specification to test variables in the data set. The full model, lmfull, specifies a formula that includes all model terms. The star indicates to add all main effects and their interactions between the two variables. The ns function comes from the splines package and is used to specify natural cubic splines to model the nonlinear effect of age. There is no left side to the equation (the dependent/outcome variable) because it is determined by the imaging variables. The reduced model, lmred specifies to include main effects of sex and the nonlinear age effect, but not their interaction. Comparing lmfull to lmred will test the terms of the full model that are related to the interaction. The third line subsets the data to rows that have no missing variables (complete case analysis).

dat$files = dat$files2mmsm8
# Testing the interaction between nonlinear age x sex
lmfull =  ~ ns(age, df=4) * sex
lmred = ~ sex + ns(age, df = 4)
dat = dat[apply(!is.na(dat[,all.vars(as.formula(lmfull))]), 1, all), ]

# add a demographic summary table using tangram

Robust test statistics

To test the parameters corresponding to the \(f_2\), we use a robust test statistic. The model is fit using the lmPBJ function and computes the coefficients, the robust test statistic, and objects required to perform inference in the next step. The robust test statistic ensures that the standard errors are consistently estimated, even if homoskedasticity is violated, i.e., it does not assume the spatial covariance of the image for each subject is the same. Unlike the parametric Gaussian random field methods implemented in SPM, the methods used in pbj yield cluster extent p-values that are valid regardless of data smoothness, or cluster-forming threshold (CFT) and perform well in small samples. pbj also supports the permutation procedure similar to what is implemented in FSL’s randomise REF.

The first argument to lmPBJ is the character vector of outcome images. The two subsequent arguments are the full and reduced formulas. The mask argument takes a niftiImage object or character vector to identify which voxels to include in the analysis. template is an optional argument that can be specified to make visualization easier after the model is fit. data is a data frame including the variables referenced by the two model formulas. The last three arguments are default values and adjust how the statistics are computed from the analysis.

• robust is a logical specifying whether robust standard errors should be used, which makes the test statistics valid even under some types of model misspecification.
• HC3 is a small sample adjustment to the robust standard error estimates that improves coverage performance (it makes the test statistics appear closer to a normal distribution).
• transform uses an inverse CDF transformation to make the test statistics appear closer to a normal distribution in small sample sizes. This is a sensible transformation if the null hypothesis is true, but might not be a good idea if the null is false.

# lmPBJ -- fits the model and computes the statistical image for the covariate of interest
## using robust test statistics (`robust = TRUE`)
robustStatMap <- lmPBJ(dat$files, form=lmfull, formred=lmred, mask=maskfile, template = templatefile, data=dat, robust = TRUE, HC3 = TRUE, transform = 'none')

The result from this command is an object of class statMap that contains the following objects:

• stat is a vector the test statistics values for all voxels in the mask.
• coef is a matrix of the parameters tested for all voxels in the mask.
• sqrtSigma is a list of objects that contain everything required to perform statistical inference in the next step. It includes settings that the user chose when running the lmPBJ command and features of the test statistic, such as the degrees of freedom and residual degrees of freedom.
• formulas is a list of the formulas used to run lmPBJ.
• data is the data set covariates used to fit the model.

Visualizing the results

When multiple parameters are being tested the test statistic is computed as a chi-squared statistical image. We can visualize the image to get an idea of what the results look like before computing p-values for topological features, such as cluster extent p-values. The image command can be used to create static images at this stage to create figures for papers or quickly visualize the results. Because the test statistic is approximately chi-square distributed, we can compute uncorrected image thresholds for visualizing the results. If we want to view where the uncorrected p-value \(p<0.01\), then we can use the command qchisq(0.01, df=robustStatMap$sqrtSigma$df, lower.tail=FALSE)= to obtain an uncorrected threshold for the image. For convenience, this can be done with the image function by setting the argument pCFT=0.01. Running image without any arguments defaults to a lightbox view. Alternatively, the user can specify another plane to visualize, or a set of specific slices.

image(robustStatMap, pCFT=0.01)

Figure 2: Lightbox view of results for the test of the age by sex interaction, with a CFT of \(p=0.01\).

image(robustStatMap, pCFT=0.01, plane = 'sagittal')

Figure 3: Lightbox view of results for the test of the age by sex interaction, with a CFT of \(p=0.01\).

The base graphics function layout can be used to visualize multiple images in one figure in a custom layout.

# custom layout visualizing inferior frontal cluster
layout(matrix(1:6, byrow=TRUE, nrow=1)); image(robustStatMap, pCFT=0.01, plane = 'axial', slice=25:30)

Figure 4: An example of a custom layout for the test of the age by sex interaction, with a CFT of \(p=0.01\).

# Visualize the chi square image
statmapFolder = tempfile() # generate a random file folder name
dir.create(statmapFolder) # create the directory
files = write.statMap(robustStatMap, outdir=statmapFolder) # write the statistical image

## Writing stat and coef images.

## Writing sqrtSigma object.

papaya(c(templatefile, files$stat )) # view the statistical image

Running statistical inference for the robust test statistic images

There are many different options for performing inference for neuroimaging data and the pbjInference function can accommodate a wide range of possibilities. Below are several examples of ways to run pbjInference. It is best to include any inference methods you want to consider, so that you don’t have to rerun the bootstrap or permutation procedure each time. Here are some key arguments to the pbjInference: * statMap is the object returned by lmPBJ. * statistic is a function used to compute common topological features of the image. * nboot is the number of bootstraps to perform. * method controls what resampling procedure to use. * ... arguments passed to the statistic function. * runMode controls what output is returned. Default is “cdf”, which returns the cumulative distribution functions used to compute adjusted and unadjusted p-values for each topological feature.

The mmeStat function is the default statistic function for pbjInference and includes options to run inference on the most common types of topological features, including local maxima, cluster mass inference, and cluster extent inference. To perform cluster extent inference (the default), the mmeStat function requires that the user specify a vector cft, which are the thresholds to use to form contiguous clusters in the statistic image. This is passed to pbjInference by the user. The final command saves the results in an Rdata file for future use.

In the code below, we use the pbjInferenceBG command to run the boostrapping in the background. This uses the r_bg function from the callr package to run the command in the background, so that we can continue to use this session.

# NEED TO EDIT THIS STILL.
if(!file.exists(pbjdatafile)){
  rcallResPBJ <- pbjInferenceBG(robustStatMap, rdata_rds=pbjdatafile, nboot = nboot, method='wild', mask = robustStatMap$mask, cft = cft)
}

if(!file.exists(permdatafile)){
  rcallResPerm <- pbjInferenceBG(robustStatMap, rdata_rds=permdatafile, nboot = nboot, method='permutation', mask = robustStatMap$mask, cft = cft)
}

Below is the standard called to pbjInference that runs in this session (not in the background). This run will perform inference on local maxima and cluster mass statistics.

# computes cluster mass inference, and inference for maxima. Does not compute CEI
maximaCEIstatMap <- pbjInference(robustStatMap, nboot = nboot, method='wild', mask = robustStatMap$mask, cft = cft, max=TRUE, CMI=TRUE, CEI=FALSE)

The following code loads in the results of the call to pbjInferenceBG. The resulting object is the statMap object from the call to lmPBJ, with the inference results added.

if(file.exists(permdatafile)){
  load(permdatafile)
  permStatMap <- statMap
}
if(file.exists(pbjdatafile)){
  load(pbjdatafile)
  robustStatMap <- statMap
}

Cluster tables

permTable = table.statMap(permStatMap)
robustTable = table.statMap(robustStatMap)
# This creates a table using the larger cluster forming threshold specified above.
stringentTable = table.statMap(robustStatMap, cft=cft[2])

knitr::kable(permTable[1:5,], row.names = FALSE,
             digits=3, booktabs=T, caption = "Cluster summary table for permutation test results of the age by sex interaction.")

Table 2: Cluster summary table for permutation test results of the age by sex interaction.

knitr::kable(robustTable[1:5,], row.names = FALSE,
             digits=3, booktabs=T, caption = "Cluster summary table for wild bootstrap test results of the age by sex interaction.")

Table 2: Cluster summary table for wild bootstrap test results of the age by sex interaction.

Writing output and interactive visualization

pbjOutdir = tempdir()
result = write.statMap(robustStatMap, pbjOutdir)
# visualizes template with statistical map and adjusted p-values overlayed
papaya(c(templatefile, result$stat, result$CEI1[1]))

#papaya(c(templatefile, result$stat, result$CMI1[1]))

Creating figures

#image(robustStatMap, plane = 'axial', cft=cft[1], alpha=0.06, roi=1)

rang = range(robustStatMap$stat)
layoutMat = cbind(matrix(1:6, nrow=2, byrow=TRUE) %x% matrix(1, nrow=3, ncol=3) , 7)
layout(layoutMat)
image(robustStatMap, plane = 'axial', cft=cft[1], alpha=0.06, slice=26:31, clusterID = FALSE)
#mtext('Probability', side=3, outer = TRUE, cex=1*cex, font=2)
# these margins adjust the height and width of the color bar bottom,left,top,right
fgcol='white'
par(mar=c(8,4,8,0.5), mgp=c(3,0.6,0), fg=fgcol, col.axis=fgcol, col.lab=fgcol, col.main = fgcol, col.sub=fgcol)
colorBar(pbj:::redyellow(64), min=cft[1], max=rang[2], nticks=4, ylab = 'Chi-sq')

plotResults = function(statMap,  emForm=NULL, method='CEI', roiInds=NULL, ind=1, data=NULL){
  pbjObj = statMap$pbj
  if(is.null(pbjObj)) stop('Must run pbjInference to plot results.')
  st = table.statMap(statMap, method, cft)
  ind = inferenceIndex(statMap$pbj$obsStat, method=method, cft=cft)
  rois = pbjObj$ROIs[[ind]]
  obsstat = pbjObj$obsStat[[ind]]
  data = statMap$data
  
  
  # load data
  imgs = RNifti::readNifti(statMap$images)
  # fit model on average data
  full = as.formula(statMap$formulas$full)
  red = as.formula(statMap$formulas$reduced)
  fullT = terms(full)
  redT = terms(red)
  # formula elements to condition on
  term = attr(fullT, 'term.labels')[!attr(fullT, 'term.labels') %in% attr(redT, 'term.labels') ]
  condVars = sapply(all.vars(full), function(x) grepl(x, term) )
  condVars = names(condVars)[condVars]
  robust = statMap$sqrtSigma$robust
  data = get_all_vars(full, data=data)
  for(roiInd in roiInds){
    data$y = sapply(imgs, function(img) mean(img[ rois==roiInd]) )
    
    fullModel = lm(update.formula(full, y ~ .), data=data)
    #redModel = lm(update.formula(red, y ~ .), data=data)
    
    # emmeans argument using condVars
    # known warning about as.numeric on characters
    suppressWarnings(ats <- apply(data[,condVars], 2, function(x) sort(unique(ifelse(is.na(as.numeric(x)), x, as.numeric(x) ) ) ) ))
    emg = emmeans::ref_grid(fullModel, at=ats, data=data)
    plotdf = as.data.frame(emmeans::emmeans(emg, ~ age | sex))
    
    
    cex=1.5
    # graphical parameters
    fgcol = 'white'
    bgcol = 'black'
    par(mgp=c(1.7,.7,0), lwd=1.5, lend=2,
        cex.lab=cex, cex.axis=0.8*cex, cex.main=1*cex,
        mar=c(3,3,2.2,0), bty='l', oma=c(0,0,0,0), bg=bgcol, fg=fgcol, col.axis=fgcol, col.lab=fgcol, col.main = fgcol, col.sub=fgcol)
    
    # plot predictions with raw data
    cols=c('#fb9a99', '#a6cee3', '#e31a1c', '#1f78b4')
    plot(data[,'age'], data[,'y'], col=cols[2+as.numeric(factor(data[,'sex']))], pch=16, ylab='Gray matter volume (AU)', xlab='Age (Years)')
    by(plotdf, plotdf$sex, function(df){
      polygon(c(df$age, rev(df$age)), c(df$lower.CL, rev(df$upper.CL)), col=scales::alpha(cols[as.numeric(df$sex)], alpha=0.8), border=NA)
      points(df$age, df$emmean, type='l', col=cols[as.numeric(df$sex)+2])
    } )
  }
}
plotResults(robustStatMap, roiInds = 1:4)