Multiple testing with anytime-valid Monte-Carlo p-values

Besag and Clifford [7] introduced a sequential permutation p-value:

where $h$ is some predefined parameter, $L_{B}:=\sum_{t=1}^{B}\mathbbm{1}\{Y_{t}\geq Y_{0}\}$ is the random number of “losses” after $B$ permutations and $\gamma(h,B)=\min(\inf\{t\in\mathbb{N}:L_{t}=h\},B)$ is a stopping time. In case of $B=\infty$, we just write $\gamma(h)$. The Besag-Clifford (BC) p-value only allows to stop at the particular time $\gamma(h,B)$ and not at any other stopping time. Hence, it allows to stop early when $h$ losses occurred before all $B$ permutations were sampled and therefore only saves computations under the null hypothesis, but not when the evidence against the null is strong. For this reason, it is not well-suited to adapt the number of permutations to the number of rejections when multiple hypotheses are considered. However, the BC p-value permits a simple anytime-valid generalization which solves this issue.

For simplicity, assume $B=\infty$, meaning the Besag-Clifford method would only stop after observing $h$ losses (the case $B<\infty$ is also possible [12], however, we do not gain much and the notation is more complicated). Suppose we plan to use the Besag-Clifford p-value $P^{\mathrm{BC}}_{\gamma(h)}$, are currently at some time $t<\gamma(h)$, and since the evidence against $H_{0}$ already seems very strong, we would like to stop the process to save computational time. What (non-trivial) p-value could we report at this time while ensuring (4) is fulfilled? The answer is very simple, just report the p-value under the most conservative future, meaning if all upcoming test statistics are losses. Since we need to sample $h-L_{t}$ test statistics until we hit the stopping time $\gamma(h)$, this leads to the following generalization of the Besag-Clifford p-value

Since $P_{\gamma(h)}^{\mathrm{avBC}}=P^{\mathrm{BC}}_{\gamma(h)}$, its anytime-validity (4) follows immediately by the validity of the Besag-Clifford p-value at time $\gamma(h)$

where $\tau$ is an arbitrary stopping time. Apart from Fischer and Ramdas [12] mentioning the anytime-valid BC method as special instance of their general approach (see next subsection), we have not seen it being exploited in the literature before. The simple idea of stopping early by assuming a conservative future was previously explored for the classical permutation p-value, but not found to be very useful [21, 18]. However, it is much more efficient for the BC method and particularly useful in the multiple testing setting, because it allows to adapt the stopping time to the individual significance level. To see this, note that for a fixed level $\alpha$, the anytime-valid BC method rejects the same hypotheses as the classical permutation p-value for $B=\lceil h/\alpha\rceil-1$, although it might stop earlier [42]. Hence, if the significance level an individual hypothesis is tested at is not predefined, as it is the case in multiple testing, the anytime-valid BC method naturally adapts the number of permutations to this data-dependent level.

In case of $h=1$, Fischer and Ramdas [12] referred to the anytime-valid BC p-value as the aggressive strategy, as it already stops after the first loss. It also requires the least possible number of permutations, which can be useful in scenarios with enormous computational demand, and therefore deserves a special mention.

Fischer and Ramdas [12] introduced a general approach to anytime-valid permutation testing based on a testing by betting technique [19, 41]. In particular, the anytime-valid BC method is a concrete instance of their algorithm. Furthermore, they were able to reduce the number of permutations compared to the classical permutation p-value and the BC method substantially, while achieving a similar power.

Their approach works as follows. At each step $t=1,2,\ldots$, the statistician bets on the outcome of the loss indicator $I_{t}=\mathbbm{1}\{Y_{t}\geq Y_{0}\}$ by specifying a betting function $B_{t}:\{0,1\}\rightarrow\mathbb{R}_{\geq 0}$. If $I_{t}=r$, the wealth of the statistician gets multiplied by $B_{t}(r)$, $r\in\{0,1\}$. Starting with an initial wealth of $W_{0}=1$, after $t$ rounds of gambling the wealth of the statistician is given by

where the concrete structure of the wealth depends on the betting strategy used. The wealth process $(W_{t})_{t\in\mathbb{N}}$ is a test martingale with respect to the filtration generated by the indicators $(\mathcal{I}_{t})_{t\in\mathbb{N}}$, if

By the optional stopping theorem, this implies that $(W_{t})_{t\in\mathbb{N}}$ is an anytime-valid e-value, meaning $\mathbb{E}_{H_{0}}[W_{\tau}]\leq 1$ for any stopping time $\tau$. In addition, by Ville’s inequality,

This shows that a hypothesis could be rejected at level $\alpha$ as soon as the wealth is larger or equal than $1/\alpha$. Ville’s inequality is often exploited by defining an anytime-valid p-value as

However, we introduce a more powerful approach in Section 3, which is particularly useful for multiple testing.

As a special instance of their general approach, Fischer and Ramdas [12] proposed the binomial mixture strategy, which particularly leads to desirable asymptotic behavior. For a predefined parameter $c\in[0,1]$, the wealth of the binomial mixture strategy is given by

where $\mathrm{Bin}(L_{t};t+1,c)$ is the CDF of a binomial distribution with size parameter $t+1$ and probability $c$. Fischer and Ramdas [12] showed that

for $t\to\infty$, where $P_{\mathrm{lim}}$ is the limiting permutation p-value $P_{\mathrm{lim}}=\lim\limits_{B\to\infty}P^{\mathrm{perm}}_{B}$. Therefore, if we choose $c<\alpha$ for some predefined significance level $\alpha$, the binomial mixture strategy rejects almost surely after a finite number of permutations if $P_{\mathrm{lim}}<c$. Hence, we can make the loss compared to the limiting permutation p-value arbitrarily small by choosing $c<\alpha$ arbitrarily close to $\alpha$. In addition to this desirable asymptotic behavior, an advantage of the binomial mixture strategy compared to the Besag-Clifford method is that it automatically adapts to the difficulty of the decision. If the evidence against the null is strong, the wealth will go fast to $1/c$, and if the evidence for the null hypothesis is strong, it will go fast to $0$. Consequently, we obtain a fast decision if the data are unambiguous, but if the decision is close, meaning the limiting permutation p-value is close to $\alpha$, the binomial mixture strategy will take its time and we can always opt to continue sampling.

In multiple testing scenarios the level of an individual hypothesis might not be predefined and depends on the rejection or non-rejection of the other hypotheses. In the following section, we will show how we can still use the binomial mixture strategy and all other $\alpha$-dependent betting strategies in multiple testing procedures.

For each null hypothesis $H_{0}^{i}$ we choose a betting strategy as described in Section 2.2 we would use if we test at level $\alpha^{\prime}\in(0,1)$. Let $W_{t}^{i,\alpha^{\prime}}$ be the wealth of the strategy for hypothesis $H_{i}$ and level $\alpha^{\prime}$ after $t$ permutations. Going forward, we assume that for all $\alpha_{1}^{\prime}<\alpha_{2}^{\prime}$,

This means that if our strategy for level $\alpha_{1}^{\prime}$ rejects at level $\alpha_{1}^{\prime}$, our strategy for some larger level $\alpha_{2}^{\prime}>\alpha_{1}^{\prime}$, needs to reject at level $\alpha_{2}^{\prime}$ as well. In particular, this is satisfied if we use for each $\alpha^{\prime}\in(0,1)$ the binomial mixture strategy with parameter $c=b\alpha^{\prime}$ for some constant $b\in(0,1)$. To see this, note that

and $(1-\mathrm{Bin}(\ell;t+1,b\alpha^{\prime}))/b$ is monotone increasing in $\alpha^{\prime}$. In the following, we therefore also use the parameter $b$ instead of $c$ to parameterize the binomial mixture strategy in an $\alpha$-independent way.

We define the $p$-value for the $i$-th hypothesis at step $t\in\mathbb{N}$ as

For example, the anytime-valid version of the BC method (6) is the result of such an $\alpha$-dependent calibration, however, Fischer and Ramdas [12] just used this implicitly and did not write this approach down in its general form. Note that if we use the same strategy for each $\alpha$, then (13) is trivially satisfied and the p-value in (14) reduces to the one in (8).

The proof follows by combining Ville’s inequality (7) and Condition (13) to note that for every stopping time $\tau_{i}$ we have

It might seem computationally challenging to calculate the p-value (14) at each step $t$. However, we do not need to evaluate $W_{t}^{i,\alpha}$ for each $\alpha\in[0,1]$, but in order to compare $P_{t}^{i}$ with some individual significance level $\alpha^{i}$ we can just check whether $W_{t}^{i,\alpha^{i}}\geq 1/\alpha^{i}$, saving a lot of computational effort. Sometimes the p-value (14) also has a closed form, as it is the case for the anytime-valid BC method (6) [12].

Our general algorithm requires two ingredients, an anytime-valid permutation p-value $(P_{t}^{i})_{t\in\mathbb{N}}$ for each hypothesis $H_{0}^{i}$, $i\in\{1,\ldots,M\}$, and a monotone multiple testing procedure $d$. A multiple testing procedure $d:[0,1]^{M}\to\{0,1\}^{M}$ mapping the p-values to the decisions ($1=\text{reject; }0=\text{accept}$) is called monotone, if $d$ is coordinatewise nonincreasing. Hence, if some of the p-values become smaller, all previous rejections remain and additional rejections might be obtained. Recall the two properties of anytime-valid p-values from Section 2: (A), $(P_{t}^{i})_{t\in\mathbb{N}}$ is nonincreasing in $t$, and (B), $(P_{t}^{i})_{t\in\mathbb{N}}$ is valid at all data-dependent random times. Consequently, due to (B), we can stop the sampling process for $H_{0}^{i}$ based on all available data, and particularly as soon as $H_{0}^{i}$ can be rejected by $d$, without violating the validity of $(P_{t}^{i})_{t\in\mathbb{N}}$. In addition, due to (A) and the monotonicity of $d$, we can already report the rejection of $H_{0}^{i}$ at that time, as the decision won’t change if we continue testing. Our general method is summarized in Algorithm 1. It should be noted that in its most general case, error control is only provided if the multiple testing procedure works under arbitrary dependence. This is captured in the following theorem, whose proof follows immediately from the explanations above.

For FWER control, the class of monotone multiple testing procedures under arbitrary dependence includes all Bonferroni-based procedures such as the Bonferroni-Holm [25], the sequentially rejective graphical multiple testing procedure [8] and the fixed sequence method [34]. For FDR control, the Benjamini-Yekuteli (BY) procedure [6] is most common. In order to control the desired error rate under arbitrary dependence, the procedures usually need to protect against the worst case distribution. Therefore, improvements can be derived under additional assumptions. A typical assumption is positive regression dependence on a subset (PRDS), under which Hochbergs’s [24] and Hommel’s [26] procedure provide uniform improvements of Holm’s procedure and the famous Benjamini-Hochberg (BH) procedure uniformly improves BY. However, if PRDS is required, caution with respect to the choice of betting strategy and stopping time needs to be taken. We address this in the next section by the example of the BH procedure.

When PRDS is required, one needs to be careful with choosing the stopping time, as it could possibly induce some kind of negative dependence. For example, suppose we stop sampling for hypothesis $H_{0}^{i}$ as soon as $P_{t}^{i}\leq m_{t}^{*}\alpha/M$. The smaller the other p-values, the sooner we could stop, which potentially induces some negative dependence between the p-values, even if the data used for calculating the different p-values is independent. However, the following proposition shows that stopping early for rejection with the BH procedure cannot inflate the FDR, which is very important as it allows to adapt the number of permutations to the number of rejections.

In the following subsections, we discuss how to choose the stopping times $\tau_{1}^{\prime},\ldots,\tau_{M}^{\prime}$ in Theorem 4.1 such that the p-values $P_{\tau_{1}^{\prime}}^{1},\ldots,P_{\tau_{M}^{\prime}}^{M}$ are PRDS.

We differentiate between two different stops; stop for rejection and stop for futility. By stop for futility, we mean that we stop early because it is unlikely that we would reject the hypothesis if we continue sampling. Thus, provided that the wealth of the used strategy is nonincreasing in the number of losses, one could write a stop for futility quite generally as

where the parameters $z_{t}^{i}$ can either be constant or dependent on the number of losses of the other hypotheses based on the multiple testing procedure used. For example, the stopping time of the Besag-Clifford method $\gamma_{i}(h,B)$ is obtained by $z_{1}^{i}=\ldots=z_{B-1}^{i}=h-1$ and $z_{B}=-1$. A stop for futility in combination with a stop for rejection includes most reasonable stopping times. Theorem 4.1 shows that stopping for rejection is never a problem. Therefore, it remains to show that stopping for futility does not violate the PRDS condition as well. One obvious sufficient condition is if the data for different hypotheses is independent and the stopping time $\tau_{i}^{\prime}$ only depends on the data for the corresponding hypothesis $H_{0}^{i}$, which is, for example, the case if the $z_{t}^{i}$ in (17) are constants. This is captured in the following proposition.

In Proposition 4.2 we use that $\boldsymbol{P}_{\mathrm{lim}}$ is independent on $I_{0}$ instead of the assumption of independent data for the different hypotheses. This makes sense since if the data is independent, then $\boldsymbol{P}_{\mathrm{lim}}$ is independent as well, thus it is a less restrictive and mathematically better graspable assumption. In addition, $P_{\mathrm{lim}}^{i}$ can also be interpreted as the p-value we would choose if we knew the distribution of $Y_{0}^{i}$ under $H_{0}^{i}$. Hence, it is a reasonable goal to ensure that our approach provides valid FDR control in those cases where the BH procedure with the limiting permutation p-values provides FDR control.

If there is some positive dependence between the data for the different hypotheses, proving PRDS for the anytime-valid p-values becomes more difficult. The reason for this is that $P_{t}^{i}$ is not only a function of the losses at step $t$ but might depend on the losses of all previous steps $L_{1}^{i},\ldots,L_{t}^{i}$. This can lead to situations in which a larger p-value $P_{t}^{i}$ yields more evidence against a large $P_{\mathrm{lim}}^{i}$, which makes it hard to use that $\boldsymbol{P}_{\mathrm{lim}}$ is PRDS for proving that $\boldsymbol{P}_{t}$ is PRDS.

We expect the effect from the above phenomenon, that larger bounds for the number of losses can lead to a lower probability for large values of the limiting permutation p-value, to be very minor and not inflating the FDR. If we choose a betting strategy with nonincreasing wealth in the number of losses and the $z_{t}^{i}$ in (17) are constant, then the stopped anytime-valid p-values $P_{\tau_{i}^{\prime}}^{i}$, $i\in\{1,\ldots,M\}$, are coordinatewise nondecreasing and nonrandom functions of the number of losses $(L_{t}^{i})_{t\in\mathbb{N}}$. Due to De Finetti’s theorem, $\boldsymbol{P}_{\mathrm{lim}}$ being strongly PRDS implies that $(L_{t_{1}}^{1},\ldots,L_{t_{M}}^{M})$ is PRDS for any $t_{1},\ldots,t_{M}\in\mathbb{N}$ (see also the proof of Proposition 4.3). Hence, in this case $P_{\tau_{1}^{\prime}}^{1},\ldots,P_{\tau_{M}^{\prime}}^{M}$ potentially have some kind of positive association as well and therefore we argue that it is reasonable to replace the classical permutation p-values with our anytime-valid ones in the BH procedure. Even if the $z_{t}^{i}$ depend in a nonincreasing way on the losses of the other hypotheses, we do not see any reason why this should hurt the PRDS of the stopped p-values. Note that in practice the classical permutation p-values can usually not be checked for PRDS either and therefore PRDS is simply assumed if there is no explicit evidence against it. Since we can additionally stop for rejection (Theorem 4.1), this provides a large class of possible betting strategies and stopping times.

Fischer and Ramdas [12] proposed to stop for futility if the wealth drops below $\alpha$. This can easily be transferred to a condition of the form (17) with constant parameters $z_{t}^{i}$. In the multiple testing case the level an individual hypothesis is tested at is not fixed in advance. For this reason, we propose to stop sampling for $H_{0}^{i}$ at time $t$, if $W_{t}^{i,\alpha_{t}^{\max}}<\alpha_{t}^{\max}$, where $\alpha_{t}^{\max}:=\alpha(|A_{t}|+m_{t}^{*})/M$ and $A_{t}$ is the index set of hypotheses for which the testing process did not stop before step $t$. Hence, $\alpha_{t}^{\max}$ can be interpreted as the maximum level a hypothesis will be tested at by the BH procedure according to the information up to step $t$. When using the binomial mixture strategy, the combined stopping time of this stop for futility with a stop for rejection is almost surely finite (see (10)). The detailed algorithm of the entire BH procedure using the binomial mixture strategy is illustrated in Appendix D. Of course, this is only one possible stop for futility and there are many other reasonable choices. For example, in situations where the sampling process takes several weeks or months it would also be possible to choose the stopping time interactively, meaning to revisit the data at some point and decide for which hypotheses to continue sampling based on the study interests and the evidence gathered so far. Our simulations in Section 5.3 confirm that the binomial mixture strategy with this stop for futility does not inflate the FDR if the limiting permutation p-values are PRDS.

First, we evaluate the power and average number of permutations per hypothesis using the aggressive strategy, the anytime-valid generalization of the BC method with $h=10$ (6), the binomial mixture strategy (9) with $b=0.9$, the AMT algorithm by Zhang et al. [48] with $\delta=0.01$ and the classic permutation p-value (1). The anytime-valid BC method was applied with its incorporated stop for futility and the binomial mixture strategy with the stop described in Section 4.3, but all methods were stopped after a maximum number of $B=\numprint{10000}$ permutations. Both anytime-valid methods were stopped for rejection as shown in (16). In addition to applying the classical permutation p-value for $B=\numprint{10000}$, we also evaluated it for $B=200$ as a reference. It should be noted that the AMT algorithm was applied at an overall level of $\alpha-\delta$ such that FDR control at level $\alpha$ is provided [48].

The results in Figure 1 were obtained by averaging over $10$ independently simulated trials. Similarly as in [48], we set the standard values of the simulation parameters to $\pi_{A}=0.4$, $\mu_{A}=2.5$, $\alpha=0.1$ and $M=1000$, while one of them was varied in each of the plots.

It can be seen that the anytime-valid BC method and binomial mixture strategy lead to a similar power as the classical permutation p-value for $B=\numprint{10000}$ permutations, while being able to reduce the number of permutations by orders of magnitude. The anytime-valid BC method performs slightly better than the binomial mixture one, however, the binomial mixture strategy provides more flexibility with respect to the stopping time and one could, for example, continue sampling for the hypotheses where the data looks promising but did not lead to a rejection after the first $\numprint{10000}$ permutations. In addition, it has advantages with respect to early reporting of rejections as described in the next section.

When the number of permutations of the classical permutation p-value is reduced to $200$, the power reduces substantially, particularly if the proportion of false hypotheses, strength of the alternative or significance level is small. Since the anytime-valid BC and binomial mixture method also need approximately $200$ permutations on average, this shows that the performance of these sequential methods cannot be accomplished by the permutation p-value with a fixed number of permutations. The AMT algorithm was outperformed by the anytime-valid BC method and the binomial mixture strategy in terms of power and number of permutations in all considered scenarios. The use of the aggressive strategy can be reasonable when the main goal is to reduce the number of permutations, while a power loss is acceptable.

It should be noted that the behavior of the methods does not change much with the number of hypotheses, since the other parameters remain fixed, which implies that $m_{t}^{*}/M$, and thus the significance level of BH procedure, remain approximately constant. Lastly, we would like to highlight that the results for all these different constellations of the simulation parameters were obtained with the same hyperparameters for the sequential methods, which thus seem to be universally applicable choices.

In the previous subsection, we have shown that a lot of permutations can be saved by our sequential strategies, while the power remains similar to the classical permutation p-value. However, reducing the total number of permutations is not the only way of increasing efficiency with sequential permutation tests. We can also report already made decisions, and particularly rejections, before the entire process has stopped. This might not be important if the entire procedure only takes some hours or one day to run. However, in large-scale trials generating all required permutations and performing the tests can take up to several months or even longer. This can slow down the research process, as the results are crucial for writing a scientific paper or identifying follow up work. Therefore, it would be very helpful to be able to already report the unambiguous decisions at an earlier stage of the process. Note that this is not possible with the AMT algorithm by Zhang et al. [48], since no decisions can be obtained until the entire process has finished.

In Figure 2 we show the distribution of the rejection times in an arbitrary simulation run using the anytime-valid BC and the binomial mixture strategy in the standard setting described in Section 5.1. In this case, the binomial mixture strategy rejected $305$ and the anytime-valid BC method $304$ hypotheses, while the former needed a mean number of $459$ and the latter of $328$ permutations to obtain a rejection. However, the distribution of the stopping times looks very different. More than $50\%$ of the rejections made by the binomial mixture strategy were obtained after $147$ permutations and more than $75\%$ after less than $180$ permutations, while the rest is distributed quite broadly up to $\numprint{10000}$ permutations. In contrast, all rejections by the anytime-valid BC procedure were made at time $328$ such that the entire process was stopped at that time.

Basically, this is because for a fixed level $\alpha$ the rejections made by the anytime-valid BC procedure are not obtained in a sequential manner, but we can just reject all hypotheses with $L_{t_{\alpha}}\leq h-1$, where $t_{\alpha}=\lceil h/\alpha\rceil-1$ (see Section 2). Hence, with the BH procedure we sample until there is a $m^{*}\in\{1,\ldots,M\}$ such that $L_{t_{\alpha m^{*}/M}}\leq h-1$ for at least $m^{*}$ hypotheses. Since $L_{t_{\alpha m^{*}/M}}\leq h-1$ for all hypotheses that have not been stopped for futility yet, we can stop the entire process and reject all the remaining hypotheses. However, if the sampling process is stopped before reaching the time $t_{\alpha m^{*}/M}$ such that $L_{t_{\alpha m^{*}/M}}\leq h-1$ for $m^{*}$ hypotheses (in our example time $328$), the anytime-valid BC method would not reject any hypothesis, while the binomial mixture strategy could have already achieved a lot of rejections.

Consequently, while the anytime-valid BC method needed less permutations in total, the binomial mixture strategy would be more useful if early reporting of the rejections is desired, since the majority of decisions is obtained much faster. In practice, the binomial mixture strategy could be used to start interpreting the results while keeping to generate further permutations in order to increase the total number of discoveries. Indeed, if the data for some of the undecided hypotheses looks promising, sampling could be continued after the first $\numprint{10000}$ permutations to possibly make even more rejections. In this case, one could also consider increasing the precision of the binomial mixture strategy by choosing the parameter $b$ closer to $1$, since a larger average number of permutations might be acceptable when the majority of decisions is obtained fast.

To evaluate the FDR control of BH with our anytime-valid p-values, we generate the vector $(Y_{0}^{1},\ldots,Y_{0}^{M})$ as multivariate normal with mean $0$ and $\mathrm{cov}(Y_{0}^{i},Y_{0}^{j})=\rho$ for $i\neq j$, $\rho\in\{0,0.1,0.3,0.5,0.7,0.9\}$. In this case, the limiting permutation p-values are strongly PRDS [6]. The other simulation parameters and parameters for the permutation testing strategies were chosen as in our standard setting introduced in Section 5.1. The FDR is obtained by averaging over $1000$ independent trials. All procedures controlled the FDR at the prespecified level $\alpha=0.01$ for all considered parameters $\rho$ (see Figure 3). This is not surprising, since the BH procedure with the aggressive strategy, the anytime-valid BC method and the classical permutation p-value control it provably and we gave reasonable arguments why it should with the binomial mixture strategy as well. Also note that in this case the BH procedure controls the FDR at level $0.06=0.6\alpha$, since the proportion of false hypotheses was set to $\pi_{A}=0.4$ [5].

Henderson et al. [23] evaluated the neural responses in the higher visual cortex to natural scene images using fMRI data from the Natural Scenes Dataset [1]. For each considered voxel in each of the eight participants, resulting in a total number of more than $\numprint{150000}$ voxels, they modeled the voxel response by texture statistics [36] in a ridge regression, where the regularization parameter was obtained by cross-validation. To evaluate model accuracy, they calculated the coefficient of determination $R^{2}$ on a held-out validation set and accessed significance using a permutation test. For this, they permuted the image labels in the training and validation data and performed the entire fitting and evaluation process for each of these permuted datasets. The final decisions were obtained by applying the BH procedure to the classical permutation p-values at level $\alpha=0.01$.

Due to the large number of hypotheses and refitting step for each permutation, this problem requires an enormous computational effort. Henderson et al. [23] drew $1000$ permutations for each hypothesis and made the resulting $R^{2}$ publicly available at https://osf.io/8fsnx/, which we also use in our analysis. Note that $1000$ permutations is arguably somewhat low for testing such a large number of hypotheses at level $\alpha=0.01$, but it turned out to not be an issue because there were many very small p-values, and thus the rejection threshold of the BH procedure was not too small. The proportion of rejections and average number of permutations per hypothesis obtained by the different strategies are illustrated in Table 1. We applied the anytime-valid BC method and the binomial mixture strategy with parameters $h=3$ and $b=0.6$, respectively, to keep the number of permutations low and stopped at the latest when all $1000$ permutations were drawn. For the AMT algorithm [48] we changed $\delta$ to $0.001$ due to the lower $\alpha$ in this case.

The results show that, although Henderson et al. [23] have already chosen a rather low number of permutations, our sequential methods were able to further reduce it significantly while rejecting a similar number of hypotheses. Also note that the binomial mixture strategy could hypothetically have continued sampling permutations for only those hypotheses where the data looks promising after generating the first $1000$ permutations, in order to possibly obtain further rejections. If the proportion of low p-values was smaller, e.g., if less than $10\%$ of the p-values were smaller than $\alpha$, then no hypothesis could have been rejected by Henderson et al. [23] no matter how strong the evidence was. In contrast, our sequential strategies with parameters $h=3$ and $b=0.6$ would still have reasonable power, but might have sampled more than $1000$ permutations for some of the hypotheses. In this case the AMT algorithm (for a maximum number of $B=1000$ permutations) would be powerless as well. The R code for reproducing this analysis is available at github.com/fischer23/MC-testing-by-betting/.

Until this point we have focused only on the total number of permutations, ignoring running time. For example, we could not study running time in our brain imaging example (above) because the test statistics were pre-computed. In this section we apply our sequential permutation testing method to analyze an RNA sequencing (RNA-seq) dataset, carefully profiling the running time of our method and comparing it to competitors.

Typically, RNA-seq data are analyzed using a parametric regression method, such as DESeq2 [31] or Limma [38]. However, in a recent, careful analysis of 13 RNA-seq datasets, Li et al. [29] found that popular tools for RNA-seq analysis — including DESeq2 — did not control type-I error on negative control data (i.e., data devoid of signal). This likely was because the parametric assumptions of these methods broke down. Li et al. instead recommended use of the Mann-Whitney (MW) test [32] — a classical, nonparametric, two-sample test — for RNA-seq data analysis.¹¹1According to Li et al. [29], the RNA-seq data should contain at least 16 samples for the MW test to have sufficient power.

Li et al. applied the MW test to analyze several datasets, including a dataset of human adipose (i.e., fat) tissue generated by the Genotype-Tissue Expression (GTEx) project [30]. Subcutaneous (i.e., directly under-the-skin) adipose tissue and visceral (i.e., deep within-the-body) adipose tissue samples were collected from 581 and 469 subjects, respectively. RNA sequencing was performed on each of these $n=1,050$ tissue samples to measure the expression level of each of $m=54,591$ genes. This experimental procedure yielded a tissue-by-gene expression matrix $Y\in\mathbb{N}^{n\times m}$ and a binary vector $X\in\{0,1\}^{n}$ indicating whether a given tissue was subcutaneous ($X_{i}=0$) or visceral ($X_{i}=1$). Li et al. sought to determine which genes were differentially expressed across the subcutaneous and visceral tissue samples. To this end, for each gene $j\in\{1,\dots,m\}$, Li et al. performed a MW test to test for association between the gene expressions $Y_{1,j},\dots,Y_{n,j}$ and the indicators $X_{1},\dots,X_{n}$, yielding p-values $p_{1},\dots,p_{m}$. Finally, Li et al. subjected the p-values to a BH correction to produce a discovery set. Li et al. used the asymptotic version of the MW test, as it is much more computationally efficient than the finite sample version of the MW test in high-multiplicity settings. The finite sample MW test generally is calibrated via permutations.

We sought to explore whether our sequential permutation testing procedure would enable computationally efficient application of the finite-sample MW test to the RNA-seq data. We thus conducted an experiment in which we applied three methods to analyze the RNA-seq data with the BH procedure: (i) the asymptotic MW p-value (implemented via the wilcox.test() function in R); (ii) a permutation p-value based on the MW statistic using a fixed number $B=5m/\alpha$ of permutations across hypotheses; and (iii) the anytime-valid BC p-value based on the MW test statistic with Algorithm 1. We implemented our anytime-valid BC method in C++ for maximum speed (see Appendix C for more information). We set the tuning parameter $h$ to $15$ in the latter method, and we set the nominal FDR across methods to $\alpha=0.1$. We compared the methods with respect to their running time and number of discoveries.

The classical permutation p-value was slowest, taking over three days and 14 hours to complete. The asymptotic p-value and the anytime-valid permutation p-value, on the other hand, were much faster, running in about two minutes and three and a half minutes, respectively. All three methods rejected a similar proportion of the null hypotheses ($56.4-56.9$% across methods). We concluded on the basis of this experiment that BH with the anytime-valid permutation p-value inherited the strenghts of both the asymptotic p-value and the classical permutation p-value: like the asymptotic p-value, the anytime-valid permutation p-value was fast, and like the classical permutation p-value, the anytime-valid permutation test avoided asymptotic assumptions. Additional information about the RNA-seq data analysis is relayed in Appendix C.

It is important to note that a significant limitation of the MW test is that it does not allow for the straightforward inclusion of covariates, such as the sex of the subject from which a given tissue sample was derived. Our goal in this section was to reproduce the analysis of Li et al., swapping their asymptotic MW test for a computationally efficient, finite-sample alternative. How to test for differential expression under minimal assumptions — while adjusting for covariates — is an ongoing topic of research in statistical genomics [35].

The authors thank Leila Wehbe for stimulating practical discussions. LF acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project number 281474342/GRK2224/2. AR was funded by NSF grant DMS-2310718.