`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:

`::install_github("angeella/pARI") devtools`

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 \(0\) and variance \(0.05\) and the tests under the alternative come from a Normal distribution with mean \(10\) and variance \(0.05\). We simulate \(10\) tests under the null and \(10\) under the alternative considering \(10\) observations. Therefore, we have a matrix with dimensions \(10 \times 20\), where the rows represent the observations and the columns the variables, i.e., tests.

We expect that the lower bound of the number of true discoveries, considering the full set of hypotheses, equals to \(10\).

```
library(pARI)
<- 20 #number of tests
m <- 10 #number of observations
n <- matrix(rnorm(0.5*m*n, 0, 0.05),ncol=n,nrow=0.5*m) #tests under the null
X <- matrix(rnorm(0.5*m*n, 10, 0.05),ncol=n,nrow=0.5*m) #tests under the alternative
Y <- cbind(Y,X) #full set of datasets data
```

Then, we perform the sign-flipping test, using \(1000\) permutations, thanks to the function
`signTest`

(type `?pARI::signTest`

for more
details):

`<- signTest(t(data), 1000) pvalues `

and we plot it in \(-log_{10}\) scale:

`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.

`<- cbind(pvalues$pv,pvalues$pv_H0) pv `

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):

```
<- hommel(pv[,1], simes = TRUE)
hom 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)$discoveries`

We have at least \(10\) true discoveries considering the full set of hypotheses.

You can find several data sets in the fMRIdata package. So, first of all you need to install it by:

`::install_github("angeella/fMRIdata") devtools`

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-x`

where
`x`

is the number identifying the subjects, e.g.,
`sub-1`

, `sub-2`

and so on. - Write in
`max_sub`

the last `x`

number of your set of
subjects. - Write in `path`

the path where your set of nifti
files is

```
<- list()
Auditory_copes <- #write here your path where your set of nifti files is. Don't put the last /
path <- #write here you last x id subjects
max_sub <- sapply(c(1:max_sub),function(x) paste0(x))
sub_ids for (sid in 1:length(sub_ids)) {
<- RNifti::readNifti(paste0(path,"/sub-", sub_ids[sid] , ".nii.gz"))
Auditory_copes[[sid]]
}
```

**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\)
level value and the **threshold** in order to perform the
cluster map using a supra-threshold statistic rule:

```
= 0.1
alpha = 3.2 thr
```

then we can perform the Permutation-based ARI using the function
`pARIbrain`

(type `?pARI::pARIbrain`

for more
details):

`<- pARIbrain(Auditory_copes,thr=thr,mask=mask,alpha = alpha) out `

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_path
mask <- Statmap_path
Statmap
<- RNifti::readNifti(mask)
mask <- RNifti::readNifti(Statmap)
Statmap =mask!=0
mask!mask]=0
Statmap[::writeNifti(Statmap, file = "Statmap.nii.gz") RNifti
```

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
= get_array(Statmap)
Statmap = get_array(mask)
mask !mask]=0
Statmap[=cluster_threshold(Statmap>3.2)
clstr
=ARIbrain::ARI(Pmap = Pmap, clusters= clstr, mask=mask, Statmap = Statmap, alpha = alpha) res_ARI
```

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[]stat[]unipd[]it or insert a reproducible example using reprex on my issue github page.