Permutation testing using Finger, Foot, Lips data
This assumes that you have completed the fingerfootlips example
using the nlmemodel.R
model. You should have a directory called
nlmemodel/
with a debug.Rdata
file, pvalue maps, and tscore maps
for the finger, foot, and lips contrasts.
Here is an overview of the steps needed after fitting the model of interest :

Generate a permutation matrix using the R package
permute
, or your favorite tool (ensure that it respects your data structure and statistical question). 
Write a
processVoxel
function that returns a Zscore from the statistical test used for the model of interest. 
Compose the
readargs.R
file to use the permutation feature of Neuropointillist, and make sure the base directory has all the necessary files (at the moment, you have to copy the permutation matrix RDS file there). 
Run
npoint
to create the permutation job. 
Compose a batch script appropriate for your scheduler (we will use a SLURM batch in this example).

Run the N permutation jobs (this may take a little while).
To use AFNI’s Equitable Thresholding and Clustering (ETAC) method for cluster correction, follow these additional steps:
 Merge the permutation_Z_ maps into
.sdat
format.  Run
3dXClustSim
using the Zmaps as input.  Run
3dMultithresh
to apply the ETAC threshold to the output image for the statistical test from the observed data.
Caveat: If your test of interest is simple, AFNI already has great
tools that implement permutation testing—for example,
3dttest++
. FSL’s randomise can handle nested data for almost any
design you can generate with FEAT, which may include some
repeatedmeasures designs.
In longitudinal neuroimaging, the model that best captures the relationship of interest may not be one supported by existing software. In this case, you can use Neuropointillist to construct the model of interest and an appropriate permutation testing scheme.
Generating permutations
It’s important to understand the model you want to test, and specifically the null hypothesis of the specific term in the model you want to test. This is because, in complex designs, you must permute the data in a way that respects the structure of the data with regard to the null for the specific regression term you wish to test. For example, we might want to look at the association between a withinperson variable and an outcome that varies at the same period. The null hypothesis is that, after controlling for all of the other personlevel and timevarying covariates, the association between the target dependent variable and the outcome is no larger than what you might see if randomly ordered that dependent variable. In other words, the null is that you could order the observations withinperson willynilly, and you’d end up with associations that are about as big as you see when the observations are ordered just as they were collected. If the association is much bigger between the outcome and the predictor variable observations ordered justso, then you can decide to reject the null hypothesis.
There are many different ways of shuffling data within and between groupinglevels (exchangeability blocks; Winkler et al, 2014), and methods of developing appropriate permutations for this. In the fingerfootlips example, I will shuffle observations withinperson by shifting the observations by some random amount for each person, per permutation (this maintains the strongly autocorrelated structure of the data). I don’t shuffle observations at all between different participants, or between different runs, though you might have some reason to do so.
I say we shuffle observations, but really, what we do in this first step is create a matrix that tells us, for each permutation, how the observations should be ordered. It is thus a P*N matrix where P is the number of permutations (e.g., 1000) and N is the number of rows in the data. We have observations in these data, so our matrix is 1000 * 3660. Since you’ve already run the target model you can simply load the debug data and use that to create your permutation matrix. It’s important to be sure that the data set you use to generate permutations is exactly what will be used by the model. This means that any rows that would be automatically dropped by, e.g., nlme
, should be dropped when generating the permutation matrix. Here’s the code I’ve used to generate the permutation matrix:
#!/usr/bin/env Rscript
library(permute)
if(!file.exists('nlmemodel/debug.Rdata')){
stop('nlmemodel/debug.Rdata does not exist. Have you run the target model?')
} else {
load('nlmemodel/debug.Rdata')
}
attach(designmat)
nperm < 1000
set.seed(2322) #for reproducibility
ctrl.free < how(within = Within(type = 'series'),
plots = Plots(strata = idnum, type = 'none'),
blocks = time,
nperm = nperm)
perm_set < shuffleSet(n = idnum, control = ctrl.free)
permpath < file.path('permutations.RDS')
saveRDS(perm_set, permpath)
A much more well commented version of this code is in the example.fingerfootlips
directory. Run this as follows.
./04write_permutations.nlmemodel.R
This script outputs a single file called permutations.RDS
which needs to be copied into the permutation base directory once we’ve created it. But first we need to write the processVoxel funciton.
processVoxel for permutations
This processVoxel function will be run once for every voxel, for every permutation. The voxel is indexed with variable v
and the permutation with permutationNumber
. We want it to return a single value that is the test statistic estimated on data constructed so the null hypothesis is true. We will focus on the contrast for finger tapping.
There are many ways of generating permuted teststatistics in data more complex than simple experimental designs in which one can simply permute the labels of conditions. There are also different ways one can account for covariates (and not accounting for covariates leads to erroneous results). The Freedman and Lane (1983) procedure has been shown to be a robust method of accounting for covariates and applying a permutation matrix to data (Anderson & Legendre, 1999). This procedure, in brief, generates the permutation teststatistic for the ith permutation using the following steps:
1) Regress the dependent variable (Y) on any covariates (but not the variable of interest), saving the residuals and the predicted values of Y. 2) Permute the residuals according to a row of the permutation matrix; add the permuted residuals to the modelpredicted Y values to produce Y. 3) Regress Y on the variable of interest, X, and covariates; save the permutation test statistic for the association between X and Y*.
In the example below, we are regressing data from a voxel on a the convolved event impulses for the finger, foot, and lips task, as well as several covariates. In the part the implements the FreedmanLane procedure, I use notation from Winkler et al (2014). Note well: the variable permutationNumber
is supplied by Neuropointillist, so don’t worry that it is never defined; we will use it to pick out the correct row of the permutation matrix.
You can find the below code with additional comments in example.fingerfootlips/nlmemodel.permute.R
.
library(nlme)
#The code here makes it easy to test the function (uncomment and run
#these lines, and then you can run what is inside the function):
#
#load('nlmemodel/debug.Rdata')
#attach(designmat)
#v < 1e5
#permutationNumber < 1
processVoxel <function(v) {
#Get the brain data from voxel `v`
BRAIN < voxeldat[,v]
#Load the permutation matrix. `permutationNumber` is supplied by neuropointillist.
permutationMatrix < readRDS('permutations.RDS')
ithPermutation < permutationMatrix[permutationNumber, ]
#To implement FreedmanLane we first estimate the model without the target
#effect, extract the residuals, permute them, and add them to the
#modelpredicted y values; this new vector, y_star, then becomes the outcome
#variable for our effect of interest (Winkler, 2014).
#First, fit the residuals model and create Y*
p < try({
#exclude `finger` since that's what we want to generate null data for
residsModel < nlme::lme(BRAIN ~ 1 + Foot + Lips + WhiteMatter + X + Y + Z + RotX + RotY + RotZ,
random = ~1  idnum,
method = c("ML"),
na.action = na.omit,
corr = corAR1(form = ~ 1  idnum),
control = lmeControl(returnObject = TRUE, singular.ok = TRUE))
#get the residuals from the model without `finger`in other words, get
#all the variability that can't be explained by the covariates.
epsilon_z < resid(residsModel, level = 1)
#permute these residuals according the the scheme we established previously
#and loaded from the RDS file.
P_j.epsilon_z < epsilon_z[ithPermutation]
#get the modelexpected y values from the covariatesonly model
Zy < predict(residsModel, level = 1)
#compose our new y_star variable from the permuted residuals and the model
#predicted y values. Now we have a variable where the variance not explained
#by the covariates is, by construction, random with respect to the
#predictor, while we also retain all of the variance that is systematically
#related to the covariates.
y_star < P_j.epsilon_z + Zy
#return the model
residsModel
})
#Target model using the permuted data
e < try(mod < lme(y_star ~ 1 + Finger + Foot + Lips + WhiteMatter + X + Y + Z + RotX + RotY + RotZ,
random = ~ 1  idnum,
method = c("ML"),
na.action = na.omit,
corr = corAR1(form = ~ 1  idnum),
control = lmeControl(returnObject = TRUE, singular.ok = TRUE)))
if(inherits(e, "tryerror")) {
mod < NULL
}
if(!is.null(mod)) {
#We want to return the Z score that corresponds to the pvalue of the
#statistic we're interested in. To do so we compute the contrast for the
#Finger condition (there are also other ways to get these values):
#
# 1 + Finger + Foot + Lips + WhiteMatter + X + Y + Z + RotX + RotY + RotZ
contr < c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0)
out < anova(mod, L = contr)
finger.p < out[["pvalue"]]
finger.Z < qnorm(finger.p)
retvals < c(finger.Z)
} else {
# If we are returning 4 dimensional data, we need to be specify how long
# the array will be in case of errors
retvals < c(NULL)
}
names(retvals) < c("fingerZ")
return(retvals)
}
The code above is very similar to the processVoxel
function for the target model. In fact, except for the permutation steps, it should be essentially identical. The one major difference is that we only return a single statistic.
Using Neuropointillist
Almost all of the pieces are in place to run Neuropointillist. The final steps are to create a readargs.R
file, run npoint
, and copy permutations.RDS
to the permutations base directory.
To follow the example, first create the readargs.R
file. Let’s also copy a smaller mask into this directory and replace the full mask with this smaller mask to run this more quickly. Note that the original file that uses the full mask can be found in readargs.nlmemode.permute.R
.
cp ../example.rawfmri/smallermask.nii.gz .
cat > readargs.R << EOF
cmdargs < c("m","smallermask.nii.gz",
"set1", "setfilenames1.txt",
"set2", "setfilenames2.txt",
"setlabels1", "setlabels1.csv",
"setlabels2", "setlabels2.csv",
"model", "nlmemodel.permute.R",
"testvoxel", "500",
"output", "nlmemodel.perms/n.p.",
"debugfile", "debug.Rdata",
"slurmN", "1", #This is ignored
"permute", "1000")
EOF
Next run npoint
which creates the base directory, nlmemodel.perms
, and tests the model.
npoint
Now, copy our permutations into this new directory.
cp permutations.RDS nlmemodel.perms
We need to create a slurm job script, because permutation testing is life at the bleeding edge.
cat > permutations.slurm.sh << 'EOF'
#!/bin/bash
#Slurm submission options
#SBATCH o npointperm_%A_%a.out
#SBATCH mailtype=END
export OMP_NUM_THREADS=1
#You shouldn't need to edit below this
num=$(printf "%04d" $SLURM_ARRAY_TASK_ID)
dashm="n.p.0001.nii.gz"
model="$(pwd P)/nlmemodel.permute.R"
permfile="fingerZ.${num}permute.nii.gz"
design="n.p.designmat.rds"
echo running: srun c 1 /usr/bin/time verbose npointrun m ${dashm} model ${model} permutationfile ${permfile} d "${design}"
srun c 1 /usr/bin/time verbose npointrun m ${dashm} model ${model} permutationfile ${permfile} d ${design}
EOF
Now submit the slurm job.
sbatch array=11000 permutations.slurm.sh
When this completes, we’ll be left with 1,000 files called n.p.0001fingerZ.NNNN.nii.gz
, where NNNN is an index from 00011000.
Creating ETAC multithreshold
Convert to sdat format
After the permutations have run, you’ll be left with a number of nii.gz
files equal to the number of permutations you requested. These must first be transformed into the sdat
file format either individually, or combined into a single 4D file. To do so, simply run:
3dtoXdataset prefix n.p.0001fingerZ \
../smallermask.nii.gz \
n.p.0001fingerZ.*.nii.gz
Note that the sdat format requires you to specify a mask as the first argument to 3dtoXdataset
. The second argument is a list of files to be converted. This command results in a single 4D sdat
file.
Getting the thresholds
Now you can use the sdat
file to get the multithresh file. Bob Cox recommends using the global
option and not the local
option (personal communication). There are many other options one can adjust. See the help file. Below is a basic example command (be sure to set the environment variable to use the number of CPU cores available to you, the more the better):
export OMP_NUM_THREADS="1" #increase as much as you can
3dXClustSim inset smallermask.nii.gz \
n.p.0001fingerZ.sdat \
global \
prefix n.p.0001fingerZ.3dXClust
This will output three files:
 globalETAC.mthresh.n.p.0001fingerZ.3dXClust.A.5perc.niml
 n.p.0001fingerZ.3dXClust.mthresh.A.5perc+orig.BRIK
 n.p.0001fingerZ.3dXClustmthresh.A.5perc+orig.HEAD
The .niml
file is something that contains the thresholds, and it does not appear this is required for applying the threshold. It’s plain text, so feel free to inspect it.
Applying the thresholds
The 3dMultiThresh
command creates a version of your grouplevel model that is appropriately thresholded according to the above ETAC mthresh files. In addition to this thresholded map, it may also produce, when passed the allmask
option, a multivolume dataset where each volume is a binary mask of voxels that pass one of the tests.
The following code presumes that the grouplevel model output for the variable of interest is a file called ../nlmemodel/n.finger.tstat.nii.gz
, and that it has been set up so that AFNI knows that it is a map of tstatistics with specific degrees of freedom. It may be simpler to ensure that the grouplevel map is output in terms of Zscores. It still may be necessary to run, e.g., 3drefit fizt group_stats_map.nii.gz
.
Note: I set prefix
and allmask
below to create the new files in the directory of the target model, ../nlmemodel/
.
3dMultiThresh mthresh n.p.0001fingerZ.3dXClust.mthresh.A.5perc+orig \
input ../nlmemodel/n.finger.tstat.nii.gz \
prefix ../nlmemodel/n.finger.tstat.multithreshed.nii.gz \
allmask ../nlmemodel/n.finger.tstat.multithreshed.allmask.nii.gz \
nozero
The nozero
option avoids creating new files if no clusters survive correction.
You are now able to visualize your thresholded statistical map (n.finger.tstat.multithreshed.nii.gz
in this example) in AFNI or whatever your favorite program is. If you do want to use AFNI, it helps to have a standard anatomical map for use as an underlay. You can copy one from FSL by executing cp $AFNI_HOME/MNI_avg152T1+tlrc* ./
. It can also be helpful to set the space for the multithreshed images so that AFNI knows it can overlay it on a standard anatomical map. Do so by running 3drefit view tlrc space MNI n.finger.tstat.multithreshed.*nii.gz
in whatever directory you saved those files.
References
Anderson, M. J., & Legendre, P. (1999). An empirical comparison of permutation methods for tests of partial regression coefficients in a linear model. Journal of Statistical Computation and Simulation, 62(3), 271–303. https://doi.org/10.1080/00949659908811936
Freedman, D., & Lane, D. (1983). A Nonstochastic Interpretation of Reported Significance Levels. Journal of Business & Economic Statistics, 1(4), 292–298. JSTOR. https://doi.org/10.2307/1391660
Winkler, A. M., Ridgway, G. R., Webster, M. A., Smith, S. M., & Nichols, T. E. (2014). Permutation inference for the general linear model. NeuroImage, 92, 381–397. https://doi.org/10.1016/j.neuroimage.2014.01.060