{'fraction_non_alphanumeric': 0.04630825214227022, 'fraction_numerical': 0.021305317525531164, 'mean_word_length': 4.270674365848628, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 10, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Permutation-based true discovery proportions for fMRI cluster analysis

We propose a permutation-based method for testing a large collection of hypotheses simultaneously. Our method provides lower bounds for the number of true discoveries in any selected subset of hypotheses. These bounds are simultaneously valid with high confidence. It is particularly useful in functional Magnetic Resonance Imaging cluster analysis, where it provides a confidence statement on the percentage of truly activated voxels within clusters of voxels, avoiding the well-known spatial specificity paradox. We offer a user-friendly tool to estimate the percentage of true discoveries for each cluster while controlling the family-wise error rate for multiple testing and taking into account that the cluster was chosen in a data-driven way. The method adapts to the spatial correlation structure that characterizes functional Magnetic Resonance Imaging data, gaining power over parametric approaches.


Introduction
Functional Magnetic Resonance Imaging (fMRI) is the most frequently used technique to understand which regions of the human brain are activated as a consequence of a stimulus. Brain activation is measured as the correlation between a sequence of (cognitive) stimuli and the resulting blood oxygenation level dependent (BOLD) signal. The BOLD signal over the entire brain is measured in small cubes termed voxels, and for each of these, we test for significant BOLD activity. Typically, around 300,000 voxels are analyzed, so the resulting multiple testing problem has roughly 300,000 statistical tests.
Controlling Type I error at the voxel level usually negatively affects the power to detect activation (Nichols and Hayasaka, 2003). Therefore, cluster-extent based thresholding was developed to analyze the data at the level of clusters of contiguous voxels. This method is less conservative than voxel-wise inference since it exploits the spatial nature of the signal (Worsley et al., 1996). However, the assumptions behind this method require a very high Random Field Theory (RFT) based initial cluster-forming threshold, resulting in relatively small clusters left for significance testing (Eklund et al., 2016). Moreover, the method suffers from the spatial specificity paradox (Woo et al., 2014), meaning that the larger the cluster we find, the less we can say about the signal within that cluster. Since the method tests the hypothesis that none of the voxels in the cluster are active, rejecting this null hypothesis only allows the claim that there is at least one active voxel inside the cluster, hence, the larger the cluster, the less we can say about it. The number of active voxels and their spatial location remains unknown, and doing follow-up inference inside the cluster ("drilling down") leads to a double-dipping problem and inflated Type I error rate (Kriegeskorte et al., 2009).
These problems motivated Rosenblatt et al. (2018) to propose All-Resolution Inference (ARI), a method to compute the lower confidence bound for the true number of active voxels within a cluster (true discovery proportion (TDP)), simultaneously for all possible sets, e.g., all clusters of voxels. Simultaneous control permits users to drill down within clusters while maintaining error guarantees, thus resolving the spatial specificity paradox. ARI is based on the approach proposed by Goeman and Solari (2011) using closed testing (Marcus et al., 1976) with local Simes test (Simes, 1986) to control the family-wise error rate (FWER). The closed testing method has an exponential computational load in general; nevertheless, for this specific case Goeman et al. (2019) and Meijer et al. (2019) proposed a fast and exact linear time short-cut. ARI relies on the Simes inequality, assuming positive regression dependency on subsets (PRDS) (Sarkar, 2008). While the Simes inequality can be assumed to be valid for fMRI data (Nichols and Hayasaka, 2003), it can be conservative under strong positive dependence. This makes the method inefficient in the neuroimaging data framework since brain measurements have strong spatial dependence due to both physics and physiology.
Permutation tests assume only exchangeability under the null hypothesis (Pesarin and Salmaso, 2010) and can handle data having any correlation structure, adapting to that correlation structure both to keep type I error control and to gain power. Hemerik et al. (2019) proposed a permutation-based method, related to ARI, that adapts the procedure to the correlation structure of the p-values. However, this method finds TDP only for sets consisting of the smallest k p-values, simultaneously over k, and can therefore not handle spatially defined clusters. Moreover, the method allows much freedom in the choice of the shape of its rejection curve, and the optimal choice for fMRI data is not clear.
In this paper, we merge the strengths of ARI with the permutation-based method of Hemerik et al. (2019), adapting ARI to use permutations using the approach of Hemerik et al. (2019). The new method provides a lower bound for the TDP for all brain regions, allowing regions of interest to be chosen post-hoc, as in ARI, without compromising family-wise error control. By using permutation-based test statistics, the method gains in power compared to the parametric version of ARI because it adapts to the correlation structure. Moreover, permutation tests are robust, as widely demonstrated in the neuroimaging literature (Eklund et al., 2016;Winkler et al., 2016a,1) and can be used when the parametric assumptions of ARI are not satisfied. The permutation-based post-hoc method proposed here is similar to the one offered by Blanchard et al. (2020), which essentially generalized the approach of Hemerik et al. (2019) to arbitrary subsets of the hypotheses. However, we also propose here an iterative approach based on the idea presented by Hemerik et al. (2019) which uniformly improves Blanchard et al. (2020)'s method in most cases. This paper is organized as follows. Section 2 introduces the concept of closed testing based on a critical vector, revisiting the results from Goeman et al. (2021) and Rosenblatt et al. (2018). Then, in Section 3, we combine these results to obtain a permutation-based ARI, and its iterative version in Section 4. We discuss the families of critical curves to be used in Section 5 and which test and permutations we recommend to use in fMRI data in Section 6. Section 7 evaluates the performance of our method in comparison with the parametric version in fMRI data. We validate the method using the resting-state fMRI null data of Eklund et al. (2016) in Section 8. Finally, we perform some simulations in Section 9 in order to investigate the influence of the shape of the rejection curve in different scenarios. In Appendix H, a simple application to genomic data is presented.

Closed testing for true discovery proportions
In this section, we revisit some results from Rosenblatt et al. (2018) and Goeman et al. (2021) to introduce notation and to clarify the need for selective inference in fMRI data.
Suppose the brain B, with |B| = m, is composed of m voxels, and let 2 B be the collection of all subsets of the brain. Some of the voxels are truly active: let A ⊆ B be the unknown set of all truly active voxels. For a cluster of interest S ⊆ B we want to make inference on a(S) = |A ∩ S|, i.e., the number of truly active voxels in S, or equivalently the TDP, i.e., a(S)/|S|.
We assume that we have computed a test statistic for each voxel i, where i = 1, . . . , m, corresponding to the null hypothesis that the voxel is not active. Based on some knowledge or guess of the marginal null distribution of these test statistics, we may compute the corresponding (parametric) p-values where Ω is the sample space of the data X. In the parametric version, we will assume that these p-values will be valid, i.e., stochastically smaller than the uniform distribution of all inactive voxels. For the permutation-based method, we emphasize here that, to guarantee FWER control, we do not make any assumptions on the distribution of these p-values. The reason is that it will suffice that the p-values are computed in the same way for all permuted versions of the data Hemerik et al. (2019).
We will now revisit the simultaneous inference on TDP using closed testing and critical vectors. First, we define a critical vector.
where N = B \ A is the set of inactive voxels, and q (i) , 1 ≤ i ≤ |N |, are their sorted p-values.
The general parametric version of ARI assumes that, for a chosen error rate α ∈ [0, 1], there is a critical vector (l 1 , . . . , l m ), possibly random, expressed as Definition 1. If such a critical vector exists, then as a corollary to Lemma 8 from Goeman et al. (2021) we have the following proposition, which we prove in Appendix F. Theorem 1. Let l i satisfy (1). Then for every ∅ = S ⊆ B, is a lower (1 − α) confidence bound of a(S), simultaneously for all S ⊆ B, that is In ARI, the Simes-based critical vector is l i = iα/h, where h is a random variable that can be calculated using the short-cut defined by Goeman et al. (2019). It is the largest set size of a subset of the brain not rejected by the Simes test. Figure 1 illustrates computation ofā(S) as defined in Equation (2), where |S| = 1000. In the left part, the length of the dashed black segments are the z = 1 − u + |{i ∈ S : p i ≤ l u }| with u ∈ {1, . . . , |S|} described in Equation (2), while the solid red segment is the maximum value over u, i.e., the highest distance between the curve of observed p-values and critical vector l i , e.g., Simes-based. In the right part, we can see the trend ofā(S) over u. The maximum value of 1 − u + z = 232 is reached when u equals 97. This implies thatā(S) = max 1≤u≤|S| 1 − u + |{i ∈ S : p i ≤ l u }| = 232 is a lower confidence bound for the number of true discoveries in S.
The crucial assumption of Theorem 1 is that l i satisfies (1). In the case of the Simes test used by ARI, this follows from the PRDS assumption, commonly also adopted for the False Discovery Rate (FDR) Example of the computation ofā(S), i.e., the maximum 1 − u + z over u where u ∈ {1, . . . , 1000} and z = 1 − u + |{i ∈ S : p i ≤ l u }|.ā(S) is then the maximum value of 1 − u + z, represented by the red dotted line, attained when u equals 97.
controlling approach proposed by Benjamini and Hochberg (1995). Although this assumption is commonly accepted in neuroimaging Nichols and Hayasaka (2003), the critical values (l 1 , . . . , l m ) can be overly strict if p-values are positively correlated, leading to conservative results. Moreover, the Simes critical vector may also be too strict or too loose if the p-values are not well calibrated.

Permutation-based All-Resolutions Inference
To obtain a critical vector that leads to improved power, we propose a permutation procedure, based on results in Hemerik et al. (2019). The permutation method takes into account the dependence structure of the p-values and, therefore, often leads to a higher critical curve than parametric methods. Moreover, permutation methods not only adapt to the dependence structure but also to the marginal distributions of the p-values. This means that we do not require the null p-values to be uniformly distributed. Instead, we require that the null p-values are exchangeable with the corresponding post-permutation p-values (Assumption 1 in Hemerik et al. (2019)).
Following Hemerik et al. (2019), we consider a group of permutations or sign-flipping transformations or any other data transformation that preserves the distribution of the test statistics under the null hypothesis, such as rotations Hemerik et al. (2020). These are maps from the support of the data distribution to itself. Our method is based on w random permutations or sign-flipping transformations. Let p 1 1 , ..., p 1 m = p 1 , ..., p m be the p-values for the real data, and for every 2 ≤ j ≤ w, let p j 1 , ..., p j m be the p-values obtained for the j-th random permutation of the data.
Computing all possible permutations could be computationally infeasible, especially in the fMRI framework. However, Proposition 2 of Hemerik and Goeman (2018) states that if the permutation set has a group structure, and α ∈ [0/w, 1/w, . . . (w − 1)/w], the random permutations reach an exact α level. This means that the α level is exhausted if all hypotheses are true, and the error rate is at most α otherwise.
To obtain the permutation-based critical vector, the user must choose a family of candidate critical vectors. Examples of such candidate vectors are given in Section 5. We suppose that the candidate vectors are indexed by λ α ∈ Λ ⊆ R, so that l(λ α ) denotes the candidate vector corresponding to λ α . The family of candidate vectors is thus F = {l(λ α ) : λ α ∈ Λ}. We assume that the family of candidate vectors is monotone, in the sense that if λ 1 α , λ 2 α ∈ Λ and λ 1 α ≤ λ 2 α , then l i (λ 1 α ) ≤ l i (λ 2 α ) for every 1 ≤ i ≤ m. We define the permutation-based critical vector to be l(λ α ), where By Hemerik et al. (2019), the following holds, so that Theorem 1 applies.
The λ α -calibration permits to incorporate the unknown dependence structure of the data into the choice of the critical vector. As seen from (4), the vector λ α is computed in such a way that for at least (1 − α)100% of the permutations, all p-values lie above it. This is illustrated in Figure 2, considering a random sample of 100 permutation curves. The λ α parameter tends to be lower, if many null hypotheses are false, being "optimal" if Definition (4) considers only p i ∈ B \ A, i.e., the set of true null hypotheses. Figure 2: Example of l i (λ α ) computation using α = 0.10. The sorted p-values p (i) are plotted against their indexes i = 1, . . . , m with m = 150, 000. The dashed red line represents the highest critical curve, i.e., the optimal critical vector (l 1 (λ α ), . . . , l m (λ α )) than the gray ones, such that the α% p-values distribution (black curves plus red one) is below it.
We use the critical vector l i (λ α ) ∈ F in Theorem 1 instead of the Simes-based one employed in ARI to gain power in computingā(S). The computation time forā(S) is primarily related to the calculation of the permutation p-values distribution and on sorting them. Our algorithm, written in C++, takes around 0.3 seconds when considering 2000 variables and 20 observations for 1000 permutations. The computation time comes from a device with a processor having 1.8 GHz CPU and 16 GB of RAM. Finally, the method can be uniformly improved by its iterative version Hemerik et al. (2019), presented in the next section.

Iterative approach
We propose here an iterative method which uniformly improvesā(S) defined in (2) following the idea proposed by Hemerik et al. (2019). The confidence envelope is defined as the minimum confidence bound computed in the complementary set of the rejection set having a cardinality equal to the lower bound of the number of true discoveries found in the previous iteration. The improvement is then substantial only when the number of detectable false hypotheses is large. In short, the iterative method improvesā(S) sequentially, in each step using the bound obtained in the previous step. The method always converges after a finite (and usually small) number of steps.
We rephrase below Theorem 2 of Hemerik et al. (2019) to get an improvement of the calibration parameter λ α .
Theorem 3. Let λ 0 α = λ α as defined in (4), we define λ α (K) as, : p x ≤ c}, and we determine: Then λ 0 α ≤ λ 1 α ≤ · · · , and for a certain i ∈ N, Theorem (3) returns a confidence bound which uniformly improves the one defined in (1) with l(λ α ) defined in (2) as critical vector. Furthermore, if in every step we compute the improved bound for all c ∈ [0, 1] and take the best one, then we are always better than the method proposed by Blanchard et al. (2020) if the same family of curves F is used. There are multiple versions of the iterative method, and that one of them (the one where in each step we find the best c) is a uniform improvement of Blanchard et al. (2020)'s method.
The iterative method is uniformly more powerful than the single-step method defined in Section 3, and also is uniformly more powerful than the step-down approach presented by Blanchard et al. (2020) under certain conditions. However, the power gain has as a cost a high computational time. In fact, the calculation of λ it α can be computationally infeasible if a large number of hypotheses is considered as in the fMRI scenario. The iterative approach must compute the minimum of a set of size |S|!/(|S| −ā i (S))!ā i (S)!. Nevertheless, we suggest to use the approximated approach defined by Hemerik et al. (2019) which can be directly applied to our case. It simply calculates the minimum across sets randomly sampled from {K ∈ S : |K| =ā i (S)} for i ∈ N. The computation time, in this case, equals approximately 10 seconds analyzing 2000 hypotheses, 20 observations, and 1000 permutations. Finally, the approximated iterative approach provided valid inference in all simulations in Hemerik et al. (2019). Please see Appendix D for further details.

Choice of family of curves
In the previous section, we consider a general family F of candidate vectors l(λ α ), λ α ∈ Λ. Here we will discuss several examples of such families, which we considered in the application later in the paper.
The first family F that we consider is inspired by Simes' probability inequality Simes (1986). The vectors are obtained by multiplying and shifting the Simes' critical vector. We denote the shift by δ ∈ {0, . . . , m − 1}. For every such δ, we have a different family, indexed by λ α ∈ R. The candidate critical vector l(λ α ) is defined by The shift parameter δ can be used to determine how sensitive the critical vector l(λ α ) will be to the smallest p-values. The parametric Simes-based approach corresponds to δ = 0 and λ α = α. We gain over that approach only if the λ α value is greater than α.
Regarding the choice of δ, note that l i (λ α ) ≤ 0 for i ≤ δ. As a consequence, we will findā(S) ≤ |S| − δ, andā(S) = 0 for all S with |S| ≤ δ. The value of δ therefore corresponds to the minimum size of a cluster that we are interested in detecting. To compensate, methods with large δ will often have a steeper slope λ α and consequently have more power for detecting large clusters. In addition, the lower confidence bound a(S) computed by the shifted version, i.e., δ > 0, can not reach the 100% true discovery proportion, since the maximum equals to (|S| − δ)/|S|. Further details about the shift parameter in computing bounds for the false discovery proportion can be found in Katsevich and Ramdas (2020).
The second example that we propose is a family F of candidate vectors that are derived from the asymptotically optimal rejection curves (AORC) considered in Finner et al. (2009) to control the FDR in an asymptotic Dirac uniform setting. Again we add a shift parameter δ ∈ {0, . . . , m} as above, and we have a different family of candidate vectors for each δ. The calibration parameters lie in Λ ⊆ R. The candidate critical vector l(λ α ) is defined by Our third example is related to the Higher Criticism method proposed by Donoho and Jin (2004). The candidate vectors are indexed by λ α ∈ Λ ⊆ R and are given by Finally, we, we consider the family of candidate vectors l(λ α ) defined as follows: Here λ α ∈ Λ = [0, 1] and F i (X) is the cumulative distribution function of the beta distribution Beta(i, m + 1 − i). This family was also considered in Hemerik et al. (2019). Further examples of candidate critical vectors can be found in Blanchard et al. (2020,0); Durand et al. (2020) and Hemerik et al. (2019), but we did not consider them here. The results obtained with our permutation method will depend on the critical vector l(λ α ) and hence on the choice of the family F = {l(λ α ) : λ α ∈ Λ}. However, one choice of a family will essentially never lead to a uniform improvement compared to another family, but only to improved TDP bounds for some sets of hypotheses and worse bounds for other sets of hypotheses. Thus, the most appropriate family will depend on which set of hypotheses we are most interested in, e.g., on whether we are interested in large or small clusters. Section 7 provides guidelines regarding a good choice of F for fMRI data.

Choice of permutations and t-statistic
The correlation between the sequence of cognitive stimuli and BOLD expresses brain activation. The changing hemodynamic impacts the local intensity of the magnetic resonance signal, i.e., the voxel intensity. Therefore, the intensity of each voxel becomes the unit of interest. The inference is then characterized by high spatial correlation and heteroscedasticity across voxels and across subjects if the between-subjects analysis is considered, due to different backgrounds, e.g., quality of the recording. The permutation approach is useful in this situation, where parametric tests fail due to violations of assumptions.
In this section, we review which permutation test is valid and powerful to perform fMRI group analysis consisting of multi-subject studies to explore the differences in BOLD response recorded under two experimental conditions Helwig (2019); Holmes et al. (1996); Winkler et al. (2014). Group fMRI data are widely analyzed using a two-stage summary statistics approach within a mixed model. This approach uses ordinary least squares (OLS) methods Mumford and Nichols (2009), in particular one-or two-sample t-tests, using within-subject parameter estimates as observations.
Let the first level within-subject model for each voxel i ∈ {1, . . . , m} and each subject j ∈ {1, . . . , J}: where Y ij ∈ R n is the brain signal of subject j in voxel i, n is the total number of time points, J is the total number of subjects, m is the total number of voxels, X j ∈ R n×p is the design matrix, where p regressors of interest, β ij ∈ R p is the vector of parameters, and ij ∈ R n is the vector of autocorrelated and non-independent error terms. Let β 1ij be the parameter relative to the first experimental condition, while β 2ij to the second experimental condition for the subject j, we then assume for simplicity p = 2. We make inference on the contrasts of parameter estimates involving brain activation differences, i.e., D ij =β 1ij −β 2ij , so: where µ i is the unknown parameter of interest representing the between-subject mean activation in voxel i, and ij are the error terms ∼ N (0, Σ). To make inference on µ i , the one-sample t-test is performed for each voxel i: . So, we have m statistical tests to analyze, one for each voxel i, i.e., H 0i : µ i = 0, that create a statistical brain mapping.
Nevertheless, we need valid permutations to have a valid permutation testing procedure. It needs a null-invariant transformation of the data, i.e., the joint distribution of the p-values under H 0i does not change Hemerik and Goeman (2018). In this case, H 0i : The compound symmetry is weaker than normality and also allows for heteroskedasticity. It can be justified by subtraction of two sample means with the same (arbitrary) distribution. Therefore, under H 0i , we can flip the sign at random of each D ij Winkler et al. (2014), always taking the identity permutation as the first transformation to have an exact α method Hemerik and Goeman (2018); Pesarin and Salmaso (2010).
The same approach can be used in the case of two-sample t-test. Let G j = {1, 2} expresses the group label for the j-th subject, the null hypothesis is then defined as H 0i : µ 1i = µ 2i . The exchangeability assumption implies (D ij |G j = 1) d = (D ij |G j = 2) for each voxel i, we can just shuffle the subject-group labels at random to compute the p-values null distribution. Permutation-based tests can be applied in various hypothesis testing's situations, e.g., tests for linear models even in the presence of nuisance effects Helwig (2019); Solari et al. (2014); Winkler et al. (2014), and tests for generalized linear models Hemerik et al. (2020).

fMRI data Application
In this section, the permutation-based ARI method is evaluated using fMRI data. Two datasets from https://openfmri.org/ are analyzed. Both datasets have the same experimental design, i.e., a block design with two stimuli. Pre-processing and first-level data analysis were performed using FSL Jenkinson et al. (2012). Registration to MNI space was done using FLIRT Jenkinson et al. (2012); Jenkinson and Smith (2001), motion correction using MCFLIRT Jenkinson et al. (2002), and brain-extraction using BET Smith (2002). We applied spatial smoothing using a Gaussian kernel of 6mm FWHM. Finally, we applied a high-pass filter to the time-series data (Gaussian-weighted least-squares straight-line fitting, with sigma=64.0s). The parameter estimates (copes), i.e., D j ∈ R m , were used as input in the pARI Andreella (2022) package developed in R R Core Team (2017). These parameter estimates are instead downloadable by installing the fMRIdata R package Andreella (2020).
For all the analyses, the α level is taken as 0.05 for a two-sided alternative hypothesis. We use 1000 permutations: 999 random permutations plus the identity. The approximated iterative approach (100 random combinations) presented in Section 4 is then applied. The results using the single-step method, i.e., λ α computed on the full set of hypotheses, are reported in Appendix A and B.
We chose δ as 0, 1, 9, and 27, to account for signal spreading out in clusters with size at least equals 0, 1, 9 and 27 voxels. The third powers were considered to exploit the three-dimensional structure of the voxels.

Auditory Data
We analyzed data from 140 subjects passively listening to vocal (i.e., speech), and non-vocal sounds, collected by Pernet et al. (2019), available at https://openneuro.org/datasets/ds000158/versions/1.0.0. We estimated the statistics map regarding the contrast that describes the difference of neural activation during vocal and non-vocal stimuli for each participant, i.e., D j . The hypothesis testing is then constructed considering H 0i : µ i = 0 with two-sided alternative, where µ i is the mean J j=1βvocal j −β non-vocal j /J computed for each voxel i = 1, . . . , m, as described in Equation (9).
In concordance with results from earlier studies Calvert et al. (1997) . While our method allows any method for forming clusters, we started from a map computed using RFT with a cluster-forming-threshold equalling |T i | > 3.2. This threshold is quite liberal, and therefore we will make additional inferences inside these clusters with a threshold of |T i | > 4. Table 1 includes the lower bounds of the proportion of active voxels (π(S) =ā(S)/|S|), the size of the cluster (|S|), the FWER-corrected p-values (p F W ER ) from classical cluster analysis and the mm coordinates of the maximum. The FWER p-values based on the clusterwise RFT are reported only for the first cluster-forming-threshold equals to |T i | > 3.2, since this method does not allow double-dipping. The results are computed using the permutation-based ARI with the Simes and AORC family using δ = 1. We compared these methods with the original parametric ARI calculated using the R package ARIbrain  As can be seen, the permutation-based ARI, i.e., columns 10 for the Simes family and 11 for the AORC family in Table 1, has a better performance overall than the parametric approach, i.e., column 12 in Table  1. However, the two families of candidate curves return very similar results; this likely reflects the similar structure of these two families of critical vectors. We also applied the shifted versions with δ > 1 (see Appendix A for the results), but found that the loss of power in small clusters is not sufficiently offset by the gain in power in the larger clusters. We believe that this is due to the conservativeness of the null p-values as shown in Figure 3. The family of critical vectors based on the Higher Criticism provide lower TDP than the ones given by the Simes and AORC families, and we put the results in Appendix A. The family based on the Beta quantile instead does not work on fMRI data due to the large number of variables that make the beta parameters unmanageable in terms of numerical precision. In addition, we believe that the weakness of both these two families is due also to mismatch between the design of the curves based on independent p-values that contrasts with a high correlation in the actual data. Figure 4 shows the TDP bounds as a cluster brain map using the results using the Simes family confidence bound. In these maps, the user can directly interpret activation as the proportion of truly active voxels inside a cluster.

Rhyme Data
Subsequently, we analyzed data from 13 subjects making rhyming judgments for pairs of either words or pseudo-words, collected by Xue and Poldrack (2020) and available at https://openneuro.org/datasets/ ds000003/versions/1.0.0. The analysis follows directly the one performed in Section 7.1, but the neural activation during the word stimulus was analyzed.
We  Xue and Poldrack (2007). The cluster map is thresholded the same as the previous dataset: using cluster-wise RFT with a threshold of |T i | > 3.2. We then drilled down π(S) using a threshold of |T i | > 4.    Figure 5 shows the null distribution of the p-values from the one-sample t-test (two-sided alternative) for the contrast regarding the word stimulus. As in Section 7.1, the Table 2 represents the results using the Simes family in column 13, and the AORC family in column 14 with δ equalling 27. As can be seen, the power improvement over the parametric method, i.e., column 15 of Table 2, is striking in this data set. We decided to take δ equals 27 since, in this case, all clusters based on RFT have large sizes. The results using δ ∈ {0, 1, 9} are shown in Appendix B. Once again, the critical vectors based on Higher Criticism and the beta distribution do not work well due to the high correlation in the data.
Finally, Figure 6 shows the TDP represented in Table 2 as cluster brain map.

Validating Permutation-based ARI
In the fMRI framework, the null data consists of resting-state fMRI data, where the null hypothesis of mean zero activation between groups is true, due to not containing consistent BOLD activity. Eklund et al. (2016) found that many software programs such as FSL and SPM Penny et al. (2011), do not control  properly the probability or the average proportion of the false positives in the cluster-wise inference of brain activity, caused by the assumption of spatial autocorrelation having squared exponential shape required for the Gaussian RFT. However, ARI does not make this assumption; for that, we want to analyze the false positive rate of the permutation-based ARI using resting-state fMRI data with no signal from Oulu data provided by the 1000 Functional Connectomes Projects Biswal et al. (2010). The pre-processing pipeline follows the one used in Eklund et al. (2016). In particular, we analyzed the Oulu data from https://tinyurl.com/clusterfailure considering fMRI images pre-processed by FSL Jenkinson et al. (2012) with a level of smoothness equals 6mm FWHM and 6 different first level designs (four event activity paradigms, and two block activity paradigms). The Oulu data consists of 103 subjects; however, to estimate the false positive rate, this set of subjects is not sufficient. In addition, Eklund et al. (2016) found asymmetric errors in the case of permutation test for the one-sample t-test using the Oulu data; therefore, we validate the permutation-based and parametric ARI, performing the two-sample t-tests. We select two groups of 20 subjects by randomly permuting 100 times the subject numbers and selecting the first 40 of this permuted data set. Eklund et al. (2016) underlines that the estimate of the false positive rate is unbiased, even if these random datasets are not independent. Finally, the set of voxels used as a cluster map is the whole-brain mask. Please see Appendix G for the results the one-sample t-tests. Figure 7 shows the FWER estimate considering six different first level designs (i.e., two-block activity paradigms: boxcar10 (10-s on-off), boxcar30 (30-s on-off) and four event activity paradigms: E1 (single event of 2-s activation, 6-s rest), E2 (single event 1-to 4-s activation, 3-to 6-s rest, randomized), E3 (13 events of 3-6 s for each task), and E4 (13 events of 3-6 s for each task, randomized). See Eklund et al. (2016) and Eklund et al. (2019) for more details about tested parameter combinations.
To sum up, in most cases, the parametric-based ARI returns a false positive rate equal to 0, while the permutation-based ARI with the Simes family returns false positive rates greater than 0 and inside the FWER confidence interval (dotted red lines). Considering the boxcar10, E1, and E2 designs, the families with lower shifts, i.e., 0 and 1, are more powerful than imposing the shift equals 3, 9, and 27. Therefore, both methods (i.e., parametric-based and permutatin-based ARI) control the FWER as expected.

Simulation study
We simulate data considering the simple following model: where Y i ∈ R m , with i = 1, . . . , n, n is the number of independent observation and m is the total number of variables. The noise i ∈ R m follows the standard normal distribution with mean 0 and equi-correlation variance structure, i.e., i ∼ N (0, Σ ρ 2 ), where ρ is the level of equi-correlation between pairs of variables. The signal µ is computed considering the difference in means having power of the one-sample t-test equals 0.8. The signal µ is equal to 0 under the null hypothesis.
First of all, we want to understand how the improvement of the nonparametric TDP lower bound changes concerning ρ and the proportion of null hypotheses π 0 . Let n = 50, m = 1000, ρ 2 ∈ {0, 0.01, . . . , 0.99, 1} and π 0 ∈ {0.6, 0.7, 0.8, 0.9}, we simulate data 1000 times and the mean ofπ(S m ) over simulation is represented. However, high values of ρ 2 are unrealistic in real applications. The Simes family of confidence bound without shift is taken into account to compare with the parametric approach directly. Having no prior knowledge about the structure of the set of hypotheses to analyze, we consider the full set of hypotheses, i.e., S m . Figure 8 shows the difference ofπ(S m ) computed using the permutation and parametric methods over the ρ 2 and π 0 values. As expected, the permutation approach gets some power with respect to the   parametric one in the case of correlation between pairs of variables. It can handle any type of dependence structure of the p-values. Secondly, we want to examine why certain families of critical curves do not provide good results in Section 7. The Higher Criticism critical vector (7), the Beta critical vector (8) and the Simes critical vector (5) are then used to computedā(S m ) using simulated data with π 0 = 0.9, m = 1000 and n = 50. As previously, we repeat the simulations 1000 times for each framework, and the mean value ofā(S m ) is computed. Figure  9 shows the behaviour of these three families of critical vectors respect to ρ 2 ∈ {0, 0.01, . . . , 0.99, 1}. In Section 5, we said that the Higher Criticism and Beta families could be problematic in the case of a strong correlation between tests. As expected, the Beta critical vector does not work in the case of a strong correlation between variables. This is due also to computational numerical difficulties, i.e., when the dashed line in Figure 9 disappears. The Higher Criticism family seems to work, but it loses power with an increase in correlation.
Thirdly, we want to analyze how the Simes family of critical curves (5) works if anti-conservative p-values distribution is considered. Let n = 50, m = 1000, ρ = 0 and π 0 = 0.9, we computeā(S m ) for each 1000 simulations, and once again the mean over simulations is reported. Figure 10 showsā(S m ) considering the Simes family using δ ∈ {0, . . . , 30}. We can note that the shifted version works well in the case of anti-conservative p-values if the corrected value for the tuning parameter δ is chosen, described by the red dotted line, i.e., δ = 8. Therefore, we explore how the Simes family of critical curves (5) works with different values of ρ and S size. The left part of Figure 11 shows the mean of the lower bounds for the TDP considering the full set of hypotheses, i.e., S m , and δ ∈ {0, 5, 10, 15, 20} over 1000 simulations. We can note that in almost all scenarios, the shifted version outperforms the unshifted ones. The difference gets smaller if ρ 2 increases. However, the situation changes if we compute the TDP for a smaller set of hypotheses than S m as shown in the right part of Figure 11. In this case, we randomly sample 40 hypotheses from the false null ones, i.e., S 40 .
Finally, the performance of the iterative approach proposed in Section 4 is compared with respect to the single-step one presented in Section 3. Let consider directly S π 1 the set of true discoveries, therefore   π(S π 1 ) = 1. Figure 12 shows the true discovery proportionπ(S π 1 ) computed on 1000 simulated datas with π 0 equals 0.9, ρ 2 ∈ {0, 0.2, 0.4} and different levels of power used to simulate the data. In this case, we consider m = 50 so that we can use the exact iterative method. First of all, we can see how the approximated iterative version equals the exact one, and, more important, how both of them uniformly improve the single-step approach.
To sum up, we suggest using the Higher Criticism and Beta families if the correlation across the variables is supposed to be low. Besides, we recommend considering the shifted version of the Simes or AORC family if the interest is in large sets of hypotheses rather than in small ones. and if the distribution of the p-value is expected to be anti-conservative, with a reasonable prior value of δ with respect to the data analyzed.
We include some simulation analyses to examine the power of the iterative approach presented in Section 4 and the influence of the number of combinations chosen for its approximated version in Appendix D.

Conclusions
Our proposed method finds simultaneous lower bounds for the TDP over all possible hypothesis subsets using the permutation theory in a computationally efficient way. As a simultaneous method, it allows the decision of which hypothesis sets to analyze to be entirely flexible and post-hoc, that is, the user can choose it after seeing the data and revise the choice as often as he/she wants. It is particularly useful in fMRI and MRI single and multi-subjects analysis to infer inside the clusters, resolving the so-called spatial specificity paradox, without falling into the double-dipping problem.
Permutation-based methods are recommended whenever they can be used, gaining power over the parametric approaches, especially when the p-values are strongly dependent, as for fMRI data. In this work, we used permutation theory to calculate the critical vector needed for ARI. Our method adapts to the correlation structure of the data in an exact way, by means of the calibration of the parameter λ α . In this way, our method remedies the existing issues of anti-conservativeness and conservativeness. Indeed, we found that the permutation-based method has more power than the parametric approach both in simulated and real data and confirmed FWER control using resting-state null data. Permutation methods are not assumption-free but require exchangeability of the test statistics under the null hypothesis. We suggest to use the OLS one-sample t-test for fMRI group analysis, having a fast and straightforward computation of the permutation null distribution, randomly flipping the sign of each subject's contrast. The exchangeability assumption needed to perform permutation-based methods is satisfied, i.e., the error terms of the model need to be symmetric around 0. Even if permutations are employed to perform the method proposed, the computation time remains low, e.g., around 250 seconds using the single-step method, while around around 1 hour using the iterative approach with 10 combinations, having 150, 000 hypotheses and 1000 permutations. The computation time is related to a device with a processor having 1.8 GHz CPU and 16 GB of RAM, finally, the R package pARI available on CRAN Andreella (2022) based on the C++ language was used.
The proposed method is general, allowing different families of confidence bounds. The choice of the family is critical since it directly influences the bounds for the true discovery proportions, and thus the power properties of the method. We analyzed several such families and found that those based on Simes and AORC are preferred to analyze the group difference between two stimuli in the fMRI framework. Shifted versions seem attractive but lose power for small clusters. We suggest using δ > 1 if the interest is in analyzing large clusters, while δ = 1 if the interest is in examining small clusters. Finally, drilling down may increase or decrease the lower bound for the TDP of some subclusters if a large cluster is analyzed, e.g., the first cluster found in Section 7.2. This suggestion is also confirmed by the simulation analysis. We found that the shifted versions gain power if the raw p-values are anti-conservative. Other types of families that we analyzed, based on Higher Criticism Donoho and Jin (2004) and Beta quantiles, do not seem to perform well in fMRI data analysis due to the strong correlation among the voxels, as also illustrated in the simulation study of Section 9. Generally, we suggest a family of critical vectors more concentrated on small p-values if the number of rejected hypotheses may be low, and a family of critical vectors more diffuse if the number of rejected hypotheses may be high.
Finally, we implement the iterative approach proposed by Hemerik et al. (2019) in the simultaneous post-hoc inference scenario, which uniformly improves the Blanchard et al. (2020)'s bounds in most cases. There is a power gain here, but for fMRI data with sparse signal, the gain is small and comes at a large computational cost.
Our presently used method provides a useful and practical selective inference for fMRI data that exploits the advantages of permutation theory and the closed-testing procedure, resolving the spatial specificity paradox with quite fast computation time. The proposed method would be applicable not only for the fMRI data but more for any other data types that may yield multiple testing problems and cluster-wise inference (e.g., electroencephalography data and genomic data).

Author contributions
All authors have directly participated in the planning and execution of the presented work.

Conflict of interest
The authors declare no potential conflict of interests.

DATA AVAILABILITY STATEMENT
The fMRI data used are available from http://github.com/angeella/fMRIdata, and the approach is developed as R package available on CRAN (https://CRAN.R-project.org/package=pARI). Table 3 contains the results using the single-step approach presented in Section 3 with δ = 1 as Table 1 of Subsection 7.1. Instead Table 4 shows the results using the Simes family of confidence bounds, considering δ equal respectively to 0, 9 and 27. In the same way, Table 5 proposes the results using the AORC family. Finally, Table 6 shows the results using the family of critical vectors based on the higher criticism. In all analysis, we consider α equals 0.05. Table 7 includes the results applying the single-step method presented in Section 3 imposing δ = 27, following the structure of Table 2 of Subsection 7.2. The performance having δ ∈ {0, 1, 9} is proposed in Table 8 using the Simes family of confidence bounds, and in Table 9 using the AORC family. Finally, Table  10 shows the results using the family of critical vectors based on the Higher Criticism. In all the analysis, the α level equals 0.05.

C. Word object analysis
We analyze the dataset provided by Duncan et al. (2009). It consists of 48 subjects looking 4 different visual stimuli: written words, pictures of objects, scrambled pictures of the same objects, and consonant letter strings. It is a block design where each block contains 16 stimuli from a single category using a one-back task, i.e., two runs are performed. Therefore, in this case, a third-level analysis was carried out. Table 11 reports the results regarding the contrast of the activation difference between word and consonant string Table 3: Auditory data: Clusters S identified with threshold |T| > 3.2 and Active Proportion Percentagē π(S) using Simes and AORC families (δ = 1) and parametric ARI, "drill down" clusters at |T| > 4. The single-step method is applied. The size of the clusters |S| and the voxel coordinates (x, y, z) Figure 13 shows the TDP as a cluster brain map regarding the results using the Simes family confidence bound.
If you are interested in analyzing other possible contrasts, e.g., words versus scrambled pictures, you can find the full dataset in Andreella (2020).

D. Iterative approach
In this section, the iterative version, defined in Theorem 3, is examined following the simulation analysis proposed in Section 9. The method is applied directly on S π 1 the set of true discoveries, therefore π(S π 1 ) = 1. Figure 14 illustrates the behavior of the approximated iterative method using a different number of combinations. The approximation version becomes exact when the number of combinations goes to infinity. However, as we can see in Figure 14, the results using only 10 combinations are nearly equal to the results using 1000 combinations. In addition, looking at Figure 15, we can see that the method is robust if a different number of variables are considered, i.e., the lines in Figure 15 are below 1 (true discovery proportion).

E. Proof of Theorem 2
Proof. In Lemma 8 by Goeman et al. (2021), define l i:n = l i for every i ≥ 1 and n ≥ 0. This lemma then implies that: are valid simultaneous bounds, as was to be shown.

F. Proof of Theorem 3
We do not report a formal proof of Theorem 3 since you can directly refer to (Hemerik et al., 2019, p. 643). Indeed, the relationship between the definition of the λ α calibration parameter and the power of the  Figure 14: Simulated true discovery lower bounds for S π 1 over different values of ρ 2 and power using the approximated iterative version with 10, 50, 100, and, 1000 random combinations. The solid black line represents the true discovery proportion π(S π 1 ) = 1.  Figure 15: Simulated true discovery lower bounds for S π 1 over different values of m variables using the approximated iterative version with 10, . . . , 100 random combinations. The solid red line represents the true discovery proportion π(S π 1 ) = 1. Permutation−based Simes δ=1 Permutation−based Simes δ=3 Permutation−based Simes δ=9 Permutation−based Simes δ=27 Figure 16: Estimated FWER considering six different first level designs, i.e., two-block activity paradigms: boxcar10 (10-s on-off), boxcar30 (30-s on-off), and four event activity paradigms, i.e., E1 (single event of 2-s activation, 6-s rest), E2 (single event 1-to 4-s activation, 3-to 6-s rest, randomized), E3 (13 events of 3-6 s for each task), and E4 (13 events of 3-6 s for each task, randomized), and six different methods to compute the TDP's lower bound (parametric Simes and permutation-based Simes considering five different values of the shift parameter, i.e., δ ∈ {0, 1, 3, 9, 27}). The solid line represents the nominal FWER equals 0.05, while the dotted ones stand for the lower bound and the upper bound since FWER is distributed as a Binomial with mean equals α = 0.05 and variance α(1 − α)/D, where D is the number of simulations, i.e., D = 100. method is evident, i.e., a large value of λ α leads to more power. We mainly have rephrased Hemerik et al.
(2019)'s theorem based on the concept of upper bound for the false discovery proportion in terms of the λ α calibration parameter. To sum up, Theorem 3 gets an improvement of λ α which gives an improved critical vector, and then using Theorem 2 we get an improved lower (1 − α) confidence bound of a(S) simultaneously for all S ⊆ B.

G. Validating Permutation-based ARI
We propose here the results of performing the one-sample t-tests instead of the two-sample t-tests in the Oulu data Eklund et al. (2016) following the same procedure of Section 8. Figure 16 is structured as Figure  7 presented in Section 8. Again, we can note how the parametric-based ARI returns in most cases an estimated FWER equals 0, while the permutation-based ARI gain power controlling the FWER at the same time.

H. Gene cluster analysis
We show here a simple example using Montgomery et al. (2010) and Pickrell et al. (2010) data, i.e., an expression set that combines two studies transcriptome genetics using second generation sequencing HapMap. Figure 17: Gene cluster analysis: we compute the lower bound for the true discovery percentage for each cluster coming from hierarchical cluster analysis. The first cut returns two clusters represented by two rectangles at the bottom (the blue and the green ones). So, using a second cut, we split the first cluster (blue rectangle) into two sub-clusters colored by blue scales and the second cluster (green rectangular) into three sub-clusters colored by green scales.
The data comprises 52580 genes and 129 samples (i.e., 60 unrelated Caucasian individuals of European descent and 69 unrelated Nigerian individuals). After pre-processing steps, we perform a two-sample t-test for each gene, and we define the sets of interest the ones computed by the hierarchical cluster analysis. Figure 17 represents the cluster dendrogram with the relative clusters cutting the three two times. The lower bound for the true discovery percentage is computed for each cluster of genes using the Simes family with δ = 0.
Alternatively, we can retrieve the gene sets for an organism under investigation from databases such as Gene Ontology (GO, http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/catalog/org_list.html). In this case we consider the KEGG pathways for human genes. The lower confidence bounds for the proportion of differentially expressed genes (TDP) for each pathway are represented in Table 12 using the Simes family with δ = 0 (single-step). Finally, Figure 18 shows the lower confidence bounds forπ(S) for each pathway with corresponding size |S|. Table 12: Gene cluster analysis: Clusters S identified by KEGG pathways for human genes and active proportion percentageπ(S) using Simes family (δ = 0). The single-step method is applied. The size of the cluster |S| is reported for each pathway.

PATHWAY NAME
|S|π ( Figure 18: Gene cluster analysis: Lower confidence bounds forπ(S) for each pathway with corresponding size |S|. The permutation-based ARI is applied with Simes family and δ = 0 (single-step version).