ARIbrain is the package for All-Resolution Inference in neuroscience.
Here we show how to compute lower bound for proportion of active voxels (or any other spatially located units) within given clusters.
The clusters can be defined a priori, on the basis of previous knowledges or on the basis of anatomical regions. Clusters of such a kind are usually called ROIs. There are no limitations to the number of ROIs that can be evaluated in the same analysis; the lower bounds for each ROI is valid simultaneously for all estimates (i.e. corrected for multiplicity).
Even more interestigly, the clusters can be defined on the basis of the same data. This is true, since the ARI allows for circular analysis, still controlling for multiplicity of inferences.
In the following we show an analysis where clusters are defined by a supra-threshold-statistic rule. This is the typical case of cluster-wise analysis followed by a multiplicity correction based on Random Field Theory. Here we follow an alternative way: we provide lower bound for proportion for the estimate of active voxels.
The sintax of the function is (type ?ARIbrain::ARI for more details)
ARI(Pmap, clusters, mask = NULL, alpha = 0.05, Statmap = function(ix) -qnorm(Pmap[ix]), summary_stat = c("max", "center-of-mass"), silent = FALSE)
The main input parameters of ARI() are:
Pmap: the map of p-values andclusters: the map of cluster index.The function accepts both character file names and 3D arrays. Therefore the minimal sintax is
ARI(Pmap, clusters)
Others maps (parameters) are:
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.Statmap which is a 3D array of statistics (usually t-values) on which the summaries are based. File name is also accepted.In order to perfom the analysis you need:
zstat.nii.gz containing the test statistic used in the analysismask.nii.gz (not mandatory, but usefull)cluster.nii.gz image of cluster index.You simply need to run on the shell:
cluster -z zstat1.nii.gz -t 3.2 -o cluster.nii.gz
This will create the cluster.nii.gz that you need.
hint: In case it retun an error message like
cluster: error while loading shared libraries: libutils.so: cannot open shared object file: No such file or directory,
type into the shell (replacing the path with your own path of the file fsl.sh):
source /etc/fsl/5.0/fsl.sh
and try again.
Get a complete help for FSL at
https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Cluster
library(ARIbrain)
pvalue_name <- system.file("extdata", "pvalue.nii.gz", package="ARIbrain")
cluster_name <- system.file("extdata", "cluster_th_3.2.nii.gz", package="ARIbrain")
zstat_name <- system.file("extdata", "zstat.nii.gz", package="ARIbrain")
mask_name <- system.file("extdata", "mask.nii.gz", package="ARIbrain")
res_ARI=ARI(Pmap = pvalue_name, clusters= cluster_name,
mask=mask_name, Statmap = zstat_name)## A hommel object for 145872 hypotheses.
## Simes inequality is assumed.
## Use p.adjust(), discoveries() or localtest() to access this object.
##
## With 0.95 confidence: at least 10857 discoveries.
## 3387 hypotheses with adjusted p-values below 0.05.
##
## Size FalseNull TrueNull ActiveProp dim1 dim2 dim3 Stat
## cl18 6907 5179 1728 0.74981902 17 57 38 7.826075
## cl17 4607 3409 1198 0.73996093 76 53 39 7.512461
## cl16 385 0 385 0.00000000 75 71 52 4.536705
## cl15 249 15 234 0.06024096 20 65 63 4.882405
## cl14 168 0 168 0.00000000 55 60 32 4.590014
## cl13 108 0 108 0.00000000 67 78 36 4.653047
## cl12 32 0 32 0.00000000 46 92 31 3.532471
## cl11 30 0 30 0.00000000 71 59 61 3.958764
## cl10 28 0 28 0.00000000 51 57 41 3.572852
## cl9 23 0 23 0.00000000 39 51 34 4.069620
## cl8 18 0 18 0.00000000 52 50 35 3.909603
## cl7 13 0 13 0.00000000 32 54 36 3.654027
## cl6 11 0 11 0.00000000 43 49 37 3.496598
## cl5 9 0 9 0.00000000 15 36 46 3.762685
## cl4 7 0 7 0.00000000 46 53 40 3.504706
## cl3 2 0 2 0.00000000 46 44 35 3.330204
## cl2 1 0 1 0.00000000 19 87 41 3.237725
## cl1 1 0 1 0.00000000 52 58 46 3.376935
## cl0 133273 0 133273 0.00000000 42 56 43 3.199741str(res_ARI)## num [1:19, 1:8] 6907 4607 385 249 168 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:19] "cl18" "cl17" "cl16" "cl15" ...
## ..$ : chr [1:8] "Size" "FalseNull" "TrueNull" "ActiveProp" ...library(RNifti)
Tmap = readNifti(system.file("extdata", "zstat.nii.gz", package="ARIbrain"))
# compute p-values from Test statistic (refering to Normal distribution, right-side alternative)
Pmap=pnorm(-Tmap)
#Read the mask file.
mask = RNifti::readNifti(system.file("extdata", "mask.nii.gz", package="ARIbrain"))
# Make sure that it is a logical map by: ()!=0
mask=mask!=0
#Create Clusters using a threshold equal to 3.2
Tmap[!mask]=0
clstr=cluster_threshold(Tmap>3.2)
table(clstr)## clstr
## 0 1 2 3 4 5 6 7 8 9
## 890030 1 1 2 7 9 11 13 18 23
## 10 11 12 13 14 15 16 17 18
## 28 30 32 108 168 249 385 4607 6907res_ARI=ARI(Pmap,clusters = clstr,mask = mask,Statmap = Tmap)## A hommel object for 145872 hypotheses.
## Simes inequality is assumed.
## Use p.adjust(), discoveries() or localtest() to access this object.
##
## With 0.95 confidence: at least 10857 discoveries.
## 3387 hypotheses with adjusted p-values below 0.05.
##
## Size FalseNull TrueNull ActiveProp dim1 dim2 dim3 Stat
## cl18 6907 5179 1728 0.74981902 17 57 38 7.826075
## cl17 4607 3409 1198 0.73996093 76 53 39 7.512461
## cl16 385 0 385 0.00000000 75 71 52 4.536705
## cl15 249 15 234 0.06024096 20 65 63 4.882405
## cl14 168 0 168 0.00000000 55 60 32 4.590014
## cl13 108 0 108 0.00000000 67 78 36 4.653047
## cl12 32 0 32 0.00000000 46 92 31 3.532471
## cl11 30 0 30 0.00000000 71 59 61 3.958764
## cl10 28 0 28 0.00000000 51 57 41 3.572852
## cl9 23 0 23 0.00000000 39 51 34 4.069620
## cl8 18 0 18 0.00000000 52 50 35 3.909603
## cl7 13 0 13 0.00000000 32 54 36 3.654027
## cl6 11 0 11 0.00000000 43 49 37 3.496598
## cl5 9 0 9 0.00000000 15 36 46 3.762685
## cl4 7 0 7 0.00000000 46 53 40 3.504706
## cl3 2 0 2 0.00000000 46 44 35 3.330204
## cl2 1 0 1 0.00000000 52 58 46 3.376935
## cl1 1 0 1 0.00000000 19 87 41 3.237725
## cl0 133273 0 133273 0.00000000 42 56 43 3.199741hom=hommel::hommel(Pmap[mask])
(thr_p=hommel::concentration(hom))## [1] 0.0008627815(thr_z=-qnorm(thr_p))## [1] 3.133804Tmap[!mask]=0
clstr=cluster_threshold(Tmap>thr_z)
table(clstr)## clstr
## 0 1 2 3 4 5 6 7 8 9
## 889443 1 2 3 9 16 21 35 38 39
## 10 11 12 13 14 15 16
## 66 119 194 257 410 4748 7228res_ARI_conc=ARI(Pmap,clusters = clstr,mask = mask,Statmap = Tmap)## A hommel object for 145872 hypotheses.
## Simes inequality is assumed.
## Use p.adjust(), discoveries() or localtest() to access this object.
##
## With 0.95 confidence: at least 10857 discoveries.
## 3387 hypotheses with adjusted p-values below 0.05.
##
## Size FalseNull TrueNull ActiveProp dim1 dim2 dim3 Stat
## cl16 7228 5179 2049 0.71651909 17 57 38 7.826075
## cl15 4748 3409 1339 0.71798652 76 53 39 7.512461
## cl14 410 0 410 0.00000000 75 71 52 4.536705
## cl13 257 15 242 0.05836576 20 65 63 4.882405
## cl12 194 0 194 0.00000000 55 60 32 4.590014
## cl11 119 0 119 0.00000000 67 78 36 4.653047
## cl10 66 0 66 0.00000000 39 51 34 4.069620
## cl9 39 0 39 0.00000000 51 57 41 3.572852
## cl8 38 0 38 0.00000000 46 92 31 3.532471
## cl7 35 0 35 0.00000000 71 59 61 3.958764
## cl6 21 0 21 0.00000000 52 50 35 3.909603
## cl5 16 0 16 0.00000000 15 36 46 3.762685
## cl4 9 0 9 0.00000000 46 53 40 3.504706
## cl3 3 0 3 0.00000000 19 87 41 3.237725
## cl2 2 0 2 0.00000000 46 44 35 3.330204
## cl1 1 0 1 0.00000000 59 70 26 3.154989
## cl0 132686 0 132686 0.00000000 31 71 52 3.133518