Equivalence testing of MRI-derived morphometric measures

email: heath.pardoe@nyulangone.org

www: https://sites.google.com/site/hpardoe

A manuscript describing these analyses has been published in the Journal of Neuroimaging [link]. If you use the data from this study, please acknowledge our work using the following citation:

Pardoe, H.R., Cutter, G., Alter, R.A., Kucharsky Hiess, R., Semmelroch, M., Parker, D., Farquharson, S., Jackson, G. and Kuzniecky, R. (2016) Pooling morphometric estimates: a statistical equivalence approach, Journal of Neuroimaging, 26(1): 109-115

Analyses are carried out using the R software package.

Equivalence testing for univariate analyses.

An equivalence analysis of Freesurfer-derived hippocampal volumes estimated from different head coils.

1. Install and load the equivalence R package. The "install.packages" command installs the library so only needs to be run once.

> install.packages("equivalence")

> library("equivalence")

2. Load in your data; the example csv file I've provided are hippocampal volume estimates obtained by processing MPRAGE acquisitions obtained from 16 healthy controls, using two different head coils (20 channel and 32 channel) and the Freesurfer longitudinal processing stream

> hv <- read.csv("two_coils_hv.csv")

3. Define an equivalence margin. For this example I've set it to 5% of the mean hippocampal volume.

> equivalence.margin <- 0.05*mean(c(hv$ch.20,hv$ch.32))

4. Do the two-one sided test (TOST).

> tost(x = hv$ch.20, y = hv$ch.32, alpha = 0.05, epsilon = equivalence.margin, paired = T)

Output should look similar to that shown below; p = 4.28e-15, which is less than 0.05, therefore these estimates can be considered equivalent.

Equivalence testing for vertex-based analyses

In this example I will show how to carry out a vertex-wise equivalence analysis of Freesurfer-derived cortical thickness estimates from different head coils. I've provided a zip archive that contains some example thickness surfaces from people that were scanned twice using the same sequence but with different head coils [zip].

1. Install R package "equivalence", and download load.mgh.R, save.mgh.R, read.label.R, and vertex.wise.equivalence.R from https://github.com/hpardoe/struct.mri.

2. Load functions

source("vertex.wise.equivalence.R")

3. The vertex.wise.equivalence function expects a dataframe containing the subject ID for each paired measurement. An easy way to make this dataframe is to import it from a text file (text file subject.list is provided with the zip archive).

head.coil.subject.df <- read.csv("subject.list",h=F)

4. Do vertex wise equivalence analysis (takes a couple of minutes on my machine)

out <- vertex.wise.equivalence(head.coil.subject.df, equivalence.margin.percent = 5, surface.file = "lh.thickness.fwhm10.fsaverage.mgh", subjects.dir="./equivalence.coil.thickness.surfaces")

The p-values are now stored in a list; out$tost.p stored two one-sided test p-values and p values from an old-fashioned paired T-test are in out$ttest.p. You can do things like:

Calculate median TOST p-value over hemisphere

median(out$tost.p,na.rm=T)

Plot distribution of p-values

plot(density(out$tost.p, na.rm=T))

Correct for multiple comparisons

tost.p.adjusted <- p.adjust(p = out$tost.p, method = "fdr")

Save p-value as a surface; for this it is sometimes more convenient to save the log-transform to make it easier to display; in this case you can get thresholds by running -log10(alpha) eg. -log10(0.05) = 1.3, which means to threshold the surface image at alpha = 0.05 you should use an actual threshold of 1.3

tost.surface <- out$template.surface

tost.surface$x <- -log10(tost.p.adjusted)

save.mgh(vol = tost.surface, fname = "my.surface.mgh")