5 Statistical inference

After computing the parametric or robust statistical image for the test of the nonlinear interaction, we can perform different types of inference. Here, we perform inference two ways, using local maxima and cluster extents.

We suggest two methods for estimating the null distribution of the test statistic. The first is the permutation procedure (similar to randomise by FSL), which has been shown to work well in previous literature (and in our paper; REFS). The second is a Rademacher wild bootstrap procedure that we showed to work well in small samples. Note that, if using the permutation procedure, the test statistic image is no longer robust to heteroskedasticity, and p-values may not maintain the type 1 error rate. That said, in most cases, we have found in simulations that the permutation procedure works very well at controlling the type 1 error rate. For a paper, you only need to choose one of the two options. The results should be similar between the two methods. Naturally, we recommend the bootstrap procedure.

We first set some parameters for the analysis: * cft is the cluster forming threshold. This line converges the p-value thresholds to the chi-squared scale. * nboot is the number of bootstraps or permutations to run. These are used to compute adjusted and unadjusted p-values.

# cluster extent inference
## cluster forming threshold
cft = c(0.01, 0.001)
# number of bootstraps (or permutations) to run
nboot=5000

5.1 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 pbjInference 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.

if(!file.exists(pbjdatafile)){
  rcallResPBJ <- pbjInference(robustStatMap, rdata_rds=pbjdatafile, nboot = nboot, method='wild', mask = robustStatMap$mask, cft_p = cft)
  rcallResPBJ$ps_cpu_times()
}

if(!file.exists(permdatafile)){
  rcallResPerm <- pbjInferenceBG(robustStatMap, rdata_rds=permdatafile, nboot = nboot, method='permutation', mask = robustStatMap$mask, cft_p = 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_p = 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
}

5.2 pTFCE inference

This code performs probabilistic threshold free cluster enhancement (pTFCE) inference using the pTFCE package. It doesn’t use resampling and instead relies on GRF-based methods. We showed using simulations that it seems to work pretty well in large samples.

pbj::ptfce.statmap(robustStatMap)