pARI is the package developed to compute the permutation-based All-Resolution Inference (ARI) method. Therefore, this method does not assume any distribution of the null distribution of the p-values. It needs to satisfy the exchangeability assumption as all permutation-based methods. For further details, please refers to Pesarin et al. (2010).
Like parametric ARI, this method aims to compute simultaneous lower confidence bounds for the number of true discoveries, i.e., active voxels, in the fMRI framework. The function takes as input the list of copes, i.e., contrast maps, one for each subject, given by neuroimaging tools as FSL, SPM, etc.
With these data, you can insert a cluster map that can be the high-level output from some neuroimaging software, a region of interest (ROI), etc. If you want to construct these cluster maps using a supra-threshold statistic rule, you can specify the threshold into the function's argument thr.
Therefore, the function pARIbrain returns the lower bounds of true discoveries, i.e., active voxels, for each cluster coming from the cluster map inserted.
You can insert these cluster maps as many times as you want, because the Permutation-based ARI, as the parametric version, allows for circular analysis, still controlling for the multiplicity of inferences.
You can install the released version of pARI with:
devtools::install_github("angeella/pARI")for the developed version or
install.packages("angeella/pARI")for the CRAN version.
Here, you can perform a toy example, using simulated data where the tests under the null hypotheses come from a Normal distribution with mean
We expect that the lower bound of the number of true discoveries, considering the full set of hypotheses, equals to
library(pARI)
m <- 20 #number of tests
n <- 10 #number of observations
X <- matrix(rnorm(0.5*m*n, 0, 0.05),ncol=n,nrow=0.5*m) #tests under the null
Y <- matrix(rnorm(0.5*m*n, 10, 0.05),ncol=n,nrow=0.5*m) #tests under the alternative
data <- cbind(Y,X) #full set of datasets
Then, we perform the sign-flipping test, using signTest(type ?pARI::signTest for more details):
pvalues <- signTest(t(data), 1000)and we plot it in
plot(-log(pvalues$pv,base = 10), pch = 20)We create the p-values matrix where the rows represent the permutations and the columns the variables, i.e., the first row represents the observed p-values; the remain rows represent the p-values under the null distribution.
pv <- cbind(pvalues$pv,pvalues$pv_H0)Then, we use the parametric approach considering the full set of hypotheses, i.e., ix equals c(1:5), using the function hommel and discoveries from the hommel package (type ?hommel::hommel and hommel::discoveries for more details):
hom <- hommel(pv[,1], simes = TRUE)
discoveries(hom,ix = c(1:20),alpha = 0.05)
and the permutation-based one using the function pARI (type ?pARI::pARI for more details)
pARI(t(data),ix = c(1:20),alpha = 0.1,family = "simes", B= 1000)$discoveriesWe have at least
You can find several data sets in the fMRIdata package. So, first of all you need to install it by:
devtools::install_github("angeella/fMRIdata")for the developed version or
install.packages("angeella/fMRIdata")for the CRAN version.
This is a basic example using fMRI data from the Auditory dataset. We need the following data:
1. The list of copes, one for each subject, in nifti file format:
data(Auditory_copes)
str(Auditory_copes)Alternatively, you can construct the list of copes using your data in this way:
- Rename the nifti files as
sub-xwherexis the number identifying the subjects, e.g.,sub-1,sub-2and so on. - Write in
max_subthe lastxnumber of your set of subjects. - Write in
paththe path where your set of nifti files is
Auditory_copes <- list()
path <- #write here your path where your set of nifti files is. Don't put the last /
max_sub <- #write here you last x id subjects
sub_ids <- sapply(c(1:max_sub),function(x) paste0(x))
for (sid in 1:length(sub_ids)) {
Auditory_copes[[sid]] <- RNifti::readNifti(paste0(path,"/sub-", sub_ids[sid] , ".nii.gz"))
}
2. the mask, which is a 3D array of logicals (i.e. TRUE/FALSE means in/out of the brain). Alternatively, it may be a (character) nifti file name. If omitted, all voxels are considered.
data(Auditory_mask)
str(Auditory_mask)3. the
alpha = 0.1
thr = 3.2then we can perform the Permutation-based ARI using the function pARIbrain(type ?pARI::pARIbrain for more details):
out <- pARIbrain(Auditory_copes,thr=thr,mask=mask,alpha = alpha)if you prefer to insert some cluster map, you can just add the argument cluster instead of thr. The argument cluster accepts the map as nifti file or as character name (path where the cluster map is). You can create the Random Field Theory based cluster map using FSL. Let the Statmap, that you can compute using
Statmap(copes,alternative = "two.sided",path = "your path",mask = mask)P.S: If you are using your copes images, before using the function Statmap you need to perform an affine registration using the command flirt of FSL to have the copes images with dimensions 91 x 109 x 91 instead of 64 x 64 x 30:
flirt -in /your_feat_directory/stats/cope3.nii -ref /your_feat_directory/reg/standard.nii -applyxfm -init /your_feat_directory/reg/example_func2standard.mat -out cope_flirt
the new cope image in this case is called cope_flirt.
After using the Statmap function, you need to type in R:
mask <- mask_path
Statmap <- Statmap_path
mask <- RNifti::readNifti(mask)
Statmap <- RNifti::readNifti(Statmap)
mask=mask!=0
Statmap[!mask]=0
RNifti::writeNifti(Statmap, file = "Statmap.nii.gz")So, you have now the Statmap.nii.gz that you will use in FSL to perform the cluster analysis, typing in the shell the following commands:
fslmaths Statmap.nii.gz -mas mask.nii.gz Statmap_mask.nii.gz
smoothest -d X -r res4d -m mask
where X is the number of subjects. Then, you have the smoothness estimate value (DLH) and the number of voxels in the mask (VOLUME) that you insert in the following command:
cluster -i Statmap_mask -t 3.2 -p 1 --dlh=DLH --volume=VOLUME --mm -o cluster.nii > cluster.txt
Finally, the cluster.nii is the cluster map that you can use in ARI.
For help about this FSL code, see the FSL book code.
Finally, you can produce also the True Discovey Proportion brain map (type ?pARI::map_TDP for more details):
map_TDP(out,path= getwd(), name = "tdp", mask)Then, you can compare it with the parametric method ARI using the ARI package:
data(Auditory_Pmap)
data(Auditory_mask)
data(Auditory_Statmap)
#Create Clusters using a threshold equal to 3.2
Statmap = get_array(Statmap)
mask = get_array(mask)
Statmap[!mask]=0
clstr=cluster_threshold(Statmap>3.2)
res_ARI=ARIbrain::ARI(Pmap = Pmap, clusters= clstr, mask=mask, Statmap = Statmap, alpha = alpha)
Andreella, A., Hemerik, J., Finos, L., Weeda, W., & Goeman, J. (2023). Permutation-based true discovery proportions for functional magnetic resonance imaging cluster analysis. Statistics in Medicine.
Rosenblatt, J. D., Finos, L., W., W. D., Solari, A., and Goeman, J. J. (2018). All-resolutions inference for brain imaging. NeuroImage, 181:786-796.
Hemerik, J., Solari, A., and Goeman, J. J. (2019). Permutation-based simultaneous confidence bounds for the false discovery proportion. Biometrika, 106(3):635-649.
Eklund, A., Nichols, E. T., and Knutsson, H. (2016). Cluster failure: Why fMRI inferences for spatial extent have inflated false-positive rates. Pnas, 113(28):7900-7905.
Please write to angela.andreella[\at]unive[\dot]it or insert a reproducible example using reprex on my issue github page.