Multiple testing with anytime-valid Monte-Carlo p-values
Lasse Fischer Aaditya Ramdas


Security Analysis of WiFi-based Sensing Systems: Threats from Perturbation Attacks

    Hangcheng Cao, 
        Wenbin Huang, 
        Guowen Xu, 
        Xianhao Chen, 
        Ziyang He, 
        Jingyang Hu, 
        Hongbo Jiang, Senior Member, IEEE,
        and Yuguang Fang, Fellow, IEEE

    April 29, 2024
==================================================================================================================================================================================================================


In contemporary problems involving  genetic or neuroimaging data,  thousands of hypotheses need to be tested. Due to their high power, and finite sample guarantees on type-1 error under weak assumptions, Monte-Carlo permutation tests are often considered as gold standard for these settings. However, the enormous computational effort required for (thousands of) permutation tests is a major burden. 
Recently, <cit.>  constructed a permutation test for a single hypothesis in which the permutations are drawn sequentially one-by-one and the testing process can be stopped at any point without inflating the type I error. They showed that the number of permutations can be substantially reduced (under null and alternative) while the power remains similar. We show how their approach  can be modified to make it suitable for a broad class of multiple testing procedures. In particular, we discuss its use with the Benjamini-Hochberg procedure and illustrate the application on a large dataset.






§ INTRODUCTION


In a classical Monte-Carlo permutation test, we observe some test statistic Y_0, generate B additional test statistics Y_1,…, Y_B and calculate the permutation p-value by 

    _B=1+∑_t=1^B 1{Y_t≥ Y_0}/1+B.

  In order to be able to reject the null hypothesis at level α, we would need to generate at least B≥ 1/α -1 permuted datasets with corresponding test statistics. This can already lead to considerably computational effort when the data generation process is demanding. However, it is even much higher when M hypotheses are tested instead of a single one. First, we need to perform a permutation test for each hypothesis, leading to a minimum number of M(1/α -1) permutations in total. Second, when multiple testing corrections are performed, each hypothesis is tested at a lower individual level than the overall level α. For example, if the Bonferroni correction is used, each hypothesis is tested at level α/M leading to the requirement of at least M(M/α -1)≈  M^2/α permutations. Even if we use more powerful multiple testing procedures such as the Benjamini-Hochberg (BH) procedure <cit.>, this lower bound for the total number of permutations usually remains. The problem is that B need to be chosen in a manner that protects against the worst case in which all p-values are large except for one that is then essentially tested at level α/M. In this paper, we consider applying multiple testing procedures to anytime-valid permutation p-values which allows to stop resampling early and reject (or accept) a hypothesis as soon as the corresponding p-value is rejected (or accepted) by the multiple testing procedure. In this way, we can adapt the number of permutations to the number of rejections and potentially save a lot of permutations.

Despite our terminology and notation focusing on permutation tests, everything in this paper identically applies to other Monte Carlo methods such as conditional randomization tests as well, which are also common in the causal inference literature, and have been applied to genetic settings under the so-called model-X assumption <cit.>.

There are many existing approaches considering sequential Monte-Carlo tests for multiple testing — in particular for the BH procedure <cit.>. All these works can be categorized into the following three approaches, from which only the first and third one are useful for saving computational resources. 

    
  * The first idea is to use early stopping rules that ensure obtaining the same decisions as applying the multiple testing procedure to the permutation p-value _B for a fixed B∈ℕ with high probability <cit.>.
    
  * Another line of work provides bounds for the probability of obtaining a different decision than the limiting permutation p-value <cit.>. Here, the main goal is not to reduce the number of permutations but to minimize the resampling risk <cit.>.
    
  * As a third approach, <cit.>, <cit.> and <cit.> consider applying multiple testing procedures to the sequential permutation p-value by <cit.>.
   
The problem with the third approach is that the Besag-Clifford p-value only allows to stop early for accepting the null hypothesis because the evidence is weak, and it cannot be stopped early for rejection. For this reason, the previous works either did not stop for rejection <cit.> or use heuristics that do not provide provable error control <cit.>. To the best of our knowledge, the only existing method allowing to adapt the number of permutations to the number of rejections while guaranteeing FDR control, is the AMT algorithm by <cit.>. They fix the number of permutations B∈ℕ in advance and consider the  generation process of the test statistics as a multi-armed bandit (MAB) problem, and their guarantee falls into the first category above. 


We will follow a strategy similar to the third path, however, instead of the Besag-Clifford p-value, we consider the anytime-valid permutation p-values recently constructed by <cit.>. These reduce the number of permutations under both the null and the alternative, without compromising power or type I error. With this, our approach allows to adapt the number of permutations to the number of rejections while providing provable error control. If the used multiple testing procedure provides the desired error control under arbitrary dependence of the p-values, our approach works with all anytime-valid permutation p-values and allows to stop fully data-adaptively. If a specific dependence structure is required, some care regarding the stopping time and used anytime-valid permutation p-value needs to be taken. For brevity, we restrict to the BH procedure in this case and prove that FDR control is guaranteed, if the limiting permutation p-values are independent, or if the limiting permutation p-values are PRDS and we use the anytime-valid BC method. In addition, we argue heuristically why we also do not expect the FDR to be inflated, if the limiting permutation p-values are PRDS and we use other strategies than the anytime-valid BC method. Consequently, in contrast to the AMT algorithm <cit.>, our method can be applied with a large variety of p-value based multiple testing procedures, allows to stop for rejection data-adaptively at any time and does not require a prespecified maximum number of permutations. Furthermore, our method outperformed the AMT algorithm significantly in all considered experiments (see Sections <ref> and <ref>).   




















In Section <ref>, we begin with a recap of the betting approach for sequential permutation tests by <cit.>. Afterward, we extend this methodology to multiple testing by using α-dependent e-to-p calibration (Section <ref>). The proposed technique is straightforward for all p-value based multiple testing procedures that do not require assumptions on the dependence structure. In Section <ref>, we discuss what needs to be considered when applying multiple testing procedures for which the p-values must be independent or PRDS with the example of the Benjamini-Hochberg procedure. Finally, we demonstrate the application of the BH procedure with sequential permutation p-values to simulated data and real neuroimaging data in Sections <ref> and <ref>, respectively. 




§ RECAP: SEQUENTIAL MONTE-CARLO TESTING BY BETTING


In this section, we recap the method by <cit.> on Monte Carlo testing of a single hypothesis via betting, since this is a key module in the multiple testing methods that we return to in the following sections.

Suppose we observed some data X_0 and have the ability to generate additional data X_1,X_2,…, e.g., by permuting the treatment labels of X_0 in a treatment vs. control trial. Let Y_i=T(X_i) for some test statistic T. If the same generating mechanism is used independently for each X_i, i≥ 1, then the test statistics Y_1, Y_2, … are always exchangeable conditional on Y_0, which means that the joint distribution of the sequence does not change for any finite permutation of the test statistics. We consider testing the null hypothesis 

    H_0: Y_0,Y_1, … are exchangeable

against the alternative 

    H_1: Y_1, Y_2, … are exchangeable conditional on  Y_0,  and  Y_0  is stochastically larger than  Y_1, Y_2,… .

Note that these hypotheses implicitly assume that we sample the permutations with replacement, since there would be an upper bound for the number of permutations otherwise. However, this is the usual state of affairs in practice, since remembering which permutations have already been drawn is computationally intensive and if the data size is large, not all permutations can be drawn anyway <cit.>. Thus, the standard approach to this testing problem uses Monte Carlo sampling, that is the permutation p-value (<ref>). In general, one would want to choose B as large as possible, since the power increases with B. However, the computational effort can be too high for large B. In order to reduce the computational cost in permutation testing, <cit.> introduced a sequential permutation p-value:

    _γ(h,B)=
                h/γ(h,B),     K_γ(h,B)=h 
    
                (K_γ(h,B)+1)/(B+1),    otherwise,

where h is some predefined parameter, K_B:=∑_t=1^B 1{Y_t≥ Y_0} is the random number of “losses” after B permutations and γ(h,B)=min( inf{t∈ℕ:K_t=h}, B) is a stopping time. In case of B=∞, we just write γ(h). The Besag-Clifford (BC) p-value only allows to stop at the particular time  γ(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. 


To solve this issue, <cit.> recently introduced an anytime-valid permutation test based on a testing by betting technique <cit.>. This makes it possible to look at the data at any time point and then decide whether to stop and make a decision or to continue sampling based on the evidence gathered so far. They showed that this approach is fundamental for the construction of (sequential) permutation tests, meaning that any rank-based permutation test can be obtained by betting. 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,…, the statistician bets on the outcome of the relative rank R_t=∑_i=0^t 1{Y_i≥ Y_t} by specifying a non-negative t+1-dimensional betting vector B_t. If R_t=r, the wealth of the statistician gets multiplied by B_t(r). Starting with an initial wealth of W_0=1, after t rounds of gambling the wealth of the statistician is given by

    W_t=∏_s=1^t B_s(R_s).

The wealth process (W_t)_t∈ℕ is a test martingale with respect to the filtration generated by the relative ranks (ℛ_t)_t∈ℕ, if 

    ∑_r=1^t+1 B_t(r)=t+1      (t∈ℕ).

By the optional stopping theorem, this implies that (W_t)_t∈ℕ is an anytime-valid e-value, meaning 𝔼_H_0[W_τ]≤ 1 for any stopping time τ. In addition, by Ville's inequality (𝔭_t)_t∈ℕ, where

    𝔭_t:=1/max_s=1,…,tW_s,

is an anytime-valid p-value, meaning ℙ_H_0(𝔭_τ≤α)≤α for all α∈ [0,1]. In particular, this shows that a hypothesis could be rejected as soon as the wealth is larger or equal than 1/α. The concrete structure of the wealth depends on the betting strategy used. 




  
An anytime-valid generalization of the Besag-Clifford method.
The wealth of the aggressive strategy after t permutations is given by 

    _t= (t+1) 1{Y_1< Y_0,…,Y_t< Y_0}.

This is the most aggressive strategy, as the wealth hits zero as soon as one of the generated test statistic is larger or equal than the observed one. However, it also requires the least possible number of permutations. At the stopping time γ(1,B), the p-value of the aggressive strategy equals the BC p-value with parameter h=1. Since the aggressive strategy is also valid at any stopping time τ < γ(1,B), it can be seen as a generalization of the BC method for h=1. <cit.>  showed that the BC method for arbitrary h can also be generalized by their strategy, which results in the p-value

    𝔭_t^avBC=h/t+h-K_t

at time t ≤γ(h), if no maximum number of permutations is prespecified. For t > γ(h), we just have 𝔭_t^avBC=𝔭_γ(h)^avBC. This anytime-valid BC method might seem very simple, since 𝔭_t^avBC is just the classical BC p-value (with B=∞) in the case that only losses occur after step t. However, we have not seen it being explored in the literature before and it seems particularly useful in the multiple testing setting. Note that for a fixed level α, the anytime-valid BC method is equivalent to the classical permutation p-value for B=⌈ h/α⌉ -1, although it might stop earlier <cit.>. 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. 

This simple idea of stopping early because the decision will remain the same if we continue sampling was previously explored for the classical permutation p-value, but not found to be very useful <cit.>. However, it is much more efficient for the BC method. First, because it automatically adapts to the individual significance level. Second, since usually h≪ B, it is much less conservative to assume that all permutations are losses until hitting h losses than B permutations. For example, suppose h=10 and B=999. If α=0.01, the anytime-valid generalizations of the BC method and the permutation p-value are equivalent. They can both stop for accepting the null hypothesis as soon as the number of losses equals 10 or stop for rejecting the null if the number of losses after 990+k permutations is k, k∈{0,1,…,9}. Now suppose α increases to 0.05. Suddenly, the anytime-valid BC method allows to stop for rejection after 190+k permutations, if k losses are observed, while at least 950 permutations would need to be sampled such that a rejection can be obtained by the permutation p-value.




  
The binomial mixture strategy.
A more sophisticated betting strategy that particularly leads to desirable asymptotic behavior is the binomial mixture strategy <cit.>. For a predefined parameter c∈ [0,1], the wealth of the binomial mixture strategy is given by 

    _t^c=(1-Bin(K_t;t+1,c))/c,

where Bin(K_t;t+1,c) is the CDF of a binomial distribution with size parameter t+1 and probability c. <cit.> showed that 

    _t^c|{=}a.s.→
            1/c,     if ∈ [0,c) 
    
            0,     if ∈ (c,1]
   for t→∞,  where  is the limiting permutation p-value =lim_B→∞_B. Therefore, if we choose c<α for some predefined significance level α, the binomial mixture strategy rejects almost surely after a finite number of permutations if <c. Hence, we can make the loss compared to the limiting permutation p-value arbitrarily small by choosing c<α arbitrarily close to α. 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 α, 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 α-dependent betting strategies in multiple testing procedures.  
    




§ MULTIPLE TESTING WITH ANYTIME-VALID PERMUTATION TESTS BY Α-DEPENDENT P-VALUE CALIBRATION


Suppose we have M null hypotheses H_0^1,…, H_0^M of the form (<ref>) with corresponding alternatives H_1^1,…, H_1^M as in (<ref>). We denote by Y_0^i, Y_1^i,… the sequence of test statistics for H_0^i and by K_B^i the random number of losses for hypothesis H_0^i after B permutations. Furthermore, let I_0 ⊆{1,…,M} be the index set of true hypotheses, R⊆{1,…,M} the set of rejected hypotheses and V=I_0∩ R the set of falsely rejected hypotheses. Two of the most common error rates considered in multiple testing are the familywise error rate (FWER) and the false discovery rate (FDR). The FWER is defined as probability of rejecting any true null hypothesis

    FWER:=ℙ(|V|>0).

The FDR is the expected proportion of falsely rejected hypotheses among all rejections 

    FDR:=𝔼[|V|/|R| 1].

The aim is to control either the FWER or FDR at some prespecified level α∈ (0,1). Although we focus on FWER and FDR in this paper for concreteness, our approach also works with all other common error metrics (FDP tail probabilities, k-FWER, etc.).
Control of the FWER implies control of the FDR <cit.> and is therefore more conservative. Hence, FDR controlling procedures often lead to more rejections and are more appropriate in large-scale exploratory hypothesis testing. In non-exploratory validation studies, FWER control is usually the norm. 

Many multiple testing methods take p-values 𝔭^1,…,𝔭^M as inputs and reject H_0^i, if 𝔭^i≤α_i for some potentially data-dependent individual significance level  α_i∈ [0,1). Therefore, we do not know in advance at which individual significance level the hypothesis H_0^i will be tested at. However, most betting strategies derived by <cit.>, including the binomial mixture strategy defined in Section <ref>, require the significance level as input and thus it need to be fixed in advance. In the following, we show how this can be circumvented by calculating an anytime-valid e-value for each level α∈ [0,1] and then calibrating these into an anytime-valid p-value by taking the smallest level α at which the corresponding e-value can reject the null hypothesis. 

For each null hypothesis H_0^i we define a strategy we would use if we test at level α' ∈ (0,1). Let W_t^i,α' be the wealth of the strategy for hypothesis H_i and level α' after t permutations. We suppose that for all α_1' < α_2', 

    {W_t^i,α_1'≥ 1/α_1'}⊆{W_t^i,α_2'≥ 1/α_2'}.

This means that if our strategy for level α_1' rejects at level α_1', our strategy for some larger level α_2'>α_1', needs to reject at level α_2' as well. In particular, this is satisfied if we use for each α'∈ (0,1) the binomial mixture strategy with parameter c=bα' for some constant b∈ (0,1). To see this, note that 

    _t^c ≥ 1/α' ⟺1-Bin(K_t;t+1,b α')/b≥ 1

and (1-Bin(k;t+1,b α'))/b is monotone increasing in α'. In the following, we therefore also use the parameter b instead of c to parameterize the binomial mixture strategy in an α-independent way.

We define the p-value for the i-th hypothesis at step t∈ℕ as 
    𝔭_t^i=inf{α∈ (0,1)|∃ s∈{1,…,t}: W_s^i,α≥ 1/α}.

For example, the anytime-valid version of the BC method (<ref>) is the result of such an α-dependent calibration, however, <cit.> 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 α, as it is for example the case with the aggressive strategy, then (<ref>) is trivially satisfied and the p-value in (<ref>) reduces to the one in (<ref>). 


    If (<ref>) holds, then (𝔭_t^i)_t∈ℕ in (<ref>) is an anytime-valid p-value.


Due to Ville's inequality and Condition (<ref>), for every stopping time τ_i we have 

    ℙ_H_0^i(𝔭_τ_i^i≤α) =ℙ_H_0^i(∃ t∈{1,…,τ_i}: W_t^i,α≥ 1/α}) ≤α.


It is important to note that (𝔭_t^i)_t∈ℕ is not only valid at any stopping time τ_i, but also at any random time T_i  which is not required to be adapted to the filtration generated by the sequential ranks for H_0^i. This means, we can use any (data-dependent) information to stop the testing process without violating the validity of the stopped p-value. In particular, we can stop and reject H_0^i as soon as 𝔭_t^i is less or equal than α_i for some t≥ 1, or equivalently, if W_t^i,α_i≥ 1/α_i, where α_i is the (possibly) data-dependent threshold obtained by a multiple testing procedure. Nevertheless, since random and stopping times are equivalent for anytime-valid p-values <cit.>, we will stick to the term stopping time in the following.

Consequently, we can choose any betting strategy from the very general class provided by (<ref>), stop at any time we want and then apply the multiple testing procedure to the stopped p-values. This is the first multiple testing method for sequential permutation p-values providing such flexibility.  We are mainly interested in monotone multiple testing procedures <cit.>. That means, if d:[0,1]^M→{0,1}^M is the multiple testing procedure mapping the p-values to the decisions (1=reject;  0=accept), then d should be coordinatewise nonincreasing. Hence, if some of the p-values become smaller, all previous rejections remain and additional rejections might be obtained. A useful property of our anytime-valid p-value (𝔭_t^i)_t∈ℕ is that it is nonincreasing in time (see (<ref>)). Together with the monotonicity of the multiple testing procedure, this implies that as soon as  𝔭_t^i can be rejected by the multiple testing procedure, we can stop the sampling process for H_0^i and report a rejection, since the decision would not change if we continue sampling. Our general method is summarized in Algorithm <ref>. This also shows that we do not need to evaluate W_t^i,α for each α∈ [0,1] in order to calculate 𝔭_t^i, but to compare 𝔭_t^i with α_i we can just check whether W_t^i,α_i≥ 1/α_i, saving a lot of computational effort.
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.


If d is a monotone multiple testing procedure with FWER or FDR control under arbitrary dependence of the p-values, then Algorithm <ref> controls the FWER or FDR respectively at level α.





For FWER control, the class of monotone multiple testing procedures under arbitrary dependence includes all Bonferroni-based procedures such as the Bonferroni-Holm <cit.>, the sequentially rejective graphical multiple testing procedure <cit.> and the fixed sequence method <cit.>. For FDR control, the Benjamini-Yekuteli (BY) procedure <cit.> 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 <cit.> and Hommel's <cit.> 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.














§ THE BENJAMINI-HOCHBERG PROCEDURE WITH ANYTIME-VALID PERMUTATION P-VALUES



The Benjamini-Hochberg (BH) procedure <cit.> rejects all hypotheses H_0^i with p-value 𝔭^i≤ m^* α /M, where

    m^*=max{ m∈{1,…,M}: ∑_i=1^M 1{𝔭^i≤ mα /M}≥ m }

with the convention max(∅)=0. In the following we write m^*_t for the BH threshold at time t∈ℕ when the p-values in (<ref>) are replaced by anytime-valid p-values at time t. The BH procedure controls the FDR when the p-values are positive regression dependent on a subset (PRDS). In order to define PRDS we need the notion of an increasing set. A set D⊆ℝ^M is increasing, if z∈ D implies y∈ D for all y≥ z.  

[PRDS <cit.>]
    A random vector of p-values 𝔭 is weekly PRDS on I_0 if for any null index i ∈ I_0 and increasing set D⊆ℝ^M, the function x↦ℙ(𝔭∈ D 
|𝔭^i ≤ x) is nondecreasing in x for any x∈ℝ. We call 𝔭 strongly PRDS on I_0, if “𝔭^i ≤ x" is replaced with “𝔭^i = x".

Strong PRDS implies weak PRDS, but for controlling the FDR with the BH procedure, weak PRDS on I_0 is sufficient <cit.>. When we just write PRDS in the remainder of this paper, we always mean weak PRDS. 


For example, the PRDS condition holds when the null p-values are independent from each other and the non-nulls. However, it also holds when there is some kind of positive dependence, which is an appropriate assumption in many trials. To guarantee FDR control by applying the BH procedure to our anytime-valid p-values, we need to ensure that our stopped p-values are PRDS. In the following we discuss when this is a reasonable assumption.



 §.§ Stopping early for rejection

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 𝔭_t^i≤ m_t^*α/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.



  Let τ_i' be a stopping time for each i∈{1,…,M} such that the stopped anytime-valid p-values  𝔭_τ_1'^1,…,𝔭_τ_M'^M are PRDS and define 
    τ_i=inf{t≥ 1:τ_i'=t 𝔭^i_t≤ m_t^* α /M}.
 Then the BH procedure applied on 𝔭_τ_1^1,…,𝔭_τ_M^M rejects the same hypotheses as if applied on 𝔭_τ'_1^1,…,𝔭_τ'_M^M and therefore controls the FDR.  


    Note that 𝔭_t^i, i∈{1,…,M}, is nonincreasing and m_t^* is nondecreasing in t. Hence, if 𝔭_t^i≤ m_t^*α /M for some t∈ℕ, then 𝔭_s^i≤ m_s^*α /M for all s ≥ t.


In the following subsections, we discuss how to choose the stopping times τ_1',…,τ_M' in Theorem <ref> such that the p-values 𝔭_τ_1'^1,…,𝔭_τ_M'^M are PRDS.











 §.§ Stopping early for futility with independent p-values

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 
    τ_i'=inf{t≥ 1: K_t^i>z_t^i},
 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 γ_i(h,B) is obtained by z_1^i=…=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 <ref> 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 τ_i' only depends on the data for the corresponding hypothesis H_0^i, which is, for example, the case if the z_t^i in (<ref>) are constants. This is captured in the following proposition.


    Suppose we generate the permuted samples independently for all hypotheses and the vector of limiting p-values 𝔭_lim is independent on I_0. Furthermore, for each i∈{1,…,M}, let the stopping time  τ_i' be given by τ_i'=inf{t∈ℕ: K_t^i > z_t^i}, where z_t^i are constant parameters. Then the stopped p-values 𝔭_τ'=(𝔭_τ_1'^1,…,𝔭_τ_M'^M) are independent on I_0. In particular, applying BH to 𝔭_τ=(𝔭_τ_1^1,…,𝔭_τ_M^M), where τ_i is defined as in (<ref>), controls the FDR.


    Due to De Finetti's theorem K_t^i, i∈{1,…,M} and t≥ 1, follows a mixture binomial distribution with size parameter t and random probability ^i. Since 𝔭_lim is independent on I_0,  𝔭_τ'
 is independent on I_0 as well.



 §.§ Stopping early for futility with PRDS p-values


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 𝔭_t^i is not only a function of the losses at step t, but might depend on the losses of all previous steps K_1^i,…, K_t^i. This can lead to situations in which a larger p-value 𝔭_t^i yields more evidence against a large ^i, which makes it hard to use that 𝔭_lim is PRDS for proving that 𝔭_t is PRDS. 
To illustrate this, consider a simple example in which for some x<x^*: 𝔭_3^i≤ x is equivalent to K_1^i≤ 0 ∪ K_2^i≤ 0 ∪ K_3^i≤ 0, which reduces to K_1^i≤ 0, and 𝔭_3^i≤ x^* is equivalent to K_1^i≤ 0 ∪ K_2^i≤ 0 ∪ K_3^i≤ 1, which reduces to K_1^i≤ 0 ∪ K_3^i≤ 1.  Recall that K_t^i follows a binomial distribution with size parameter t and random probability ^i. With this, under the assumption that ^i follows a uniform distribution, one can check that ℙ(^i> 0.9| K_1^i≤ 0)=0.01>0.0091=ℙ(^i> 0.9| K_1^i≤ 0 ∪ K_3^i≤ 1). To understand why this is the case, note that K_1^i≤ 0 yields the following set of possible loss sequences {(0,1,1),(0,0,1),(0,1,0),(0,0,0)}, while K_1^i≤ 0 ∪ K_3^i≤ 1 only adds (0,0,1) to this set. For this reason, K_1^i≤ 0 ∪ K_3^i≤ 1 provides more evidence against small values for ^i, but also more evidence against very large values of ^i. Therefore, ^i is not PRDS on 𝔭_t^i in general. 

Nevertheless, we expect the effect from the above phenomenon, that larger bounds for the number of losses can lead to a lower probability for very 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 (<ref>) are constant, then the stopped anytime-valid p-values 𝔭_τ_i'^i, i∈{1,…,M}, are coordinatewise nondecreasing and nonrandom functions of the number of losses (K_t^i)_t∈ℕ. Due to De Finetti's theorem, 𝔭_lim being PRDS implies that (K_t^1,…, K_t^M) is PRDS for any t∈ℕ. Hence, in this case 𝔭_τ_1^1,…, 𝔭_τ_M^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. Since we can additionally stop for rejection (Theorem <ref>), this provides a large class of possible betting strategies and stopping times. <cit.> proposed to stop for futility if the wealth drops below α. This can easily be transferred to a condition of the form (<ref>) 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,α_t^max<α_t^max, where α_t^max:=α (|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, α_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 if ^i is continuously distributed (see (<ref>)). The detailed algorithm of the entire BH procedure using the binomial mixture strategy is illustrated in Appendix <ref>. 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. 

We believe that for most applications this heuristic argumentation should be reasonable enough to replace the classical permutation p-values with our sequential ones in the BH procedure. Note that the PRDS of the classical permutation p-values cannot usually be checked in practice either and is just assumed if there is no explicit evidence against this. However, if this heuristic argumentation for PRDS of the stopped p-values is not allowed, one could still use the anytime-valid version of the BC method (<ref>) with its incorporated stop for futility γ(h) as shown by the next proposition.
 









    Suppose we generate the permuted samples independently for all hypotheses and the vector of limiting p-values 𝔭_lim is strongly PRDS on I_0. Then the anytime-valid BC p-values  𝔭_γ(h)^avBC=(𝔭_γ_1(h)^avBC,1,…,𝔭_γ_M(h)^avBC,M), where γ_i(h)=inf{t≥ 1: K_t^i=h},  are weakly PRDS on I_0. In particular, applying BH to 𝔭_τ^avBC=(𝔭_τ_1^avBC,1,…,𝔭_τ_M^avBC,M), where τ_i is defined as in (<ref>) for τ_i'=γ_i(h), controls the FDR.

The proof can be found in the appendix.
















     

















 











































§ SIMULATED EXPERIMENTS

In this section, we aim to characterize the behavior of applying the BH procedure with the proposed sequential permutation p-values using simulated data. We consider an independent Gaussian mean multiple testing problem. That means, each hypothesis is given by H_0^i:𝔼[Y_0^i]=0, i∈{1,…,M}, where Y_0^i follows a standard normal distribution under the null hypothesis and a shifted standard normal distribution with mean μ_A under the alternative. The probability for a hypothesis being false was set to π_A∈ [0,1] and the generated test statistics Y_1^i,…, Y_B^i were sampled from N(0,1). Note that this matches the setting described in Section <ref>. 




 §.§ Power and average number of permutations

First, we evaluate the power and average number of permutations per hypothesis using the aggressive strategy (<ref>), the anytime-valid generalization of the BC method with h=10 (<ref>), the binomial mixture strategy (<ref>) with b=0.9, the AMT algorithm by <cit.> with δ=0.01 and the classic permutation p-value (<ref>). The early stops for our sequential methods were chosen as described in Section <ref>, but all methods were stopped after a maximum number of B=10000 permutations. In addition to applying the classical permutation p-value for B=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 α-δ such that FDR control at level α is provided <cit.>.




The results in Figure <ref> were obtained by averaging over 10 independently simulated trials. Similarly as in <cit.>, we set the standard values of the simulation parameters to π_A=0.4, μ_A=2.5, α=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=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 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.




 §.§ Early reporting of decisions

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 <cit.>, since no decisions can be obtained until the entire process has finished. 

In Figure <ref> 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 <ref>. 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 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 α the rejections made by the anytime-valid BC procedure are not obtained in a sequential manner, but we can just reject all hypotheses with K_t_α≤ h-1, where t_α=⌈ h/α⌉ -1 (see Section <ref>). Hence, with the BH procedure we sample until there is a m^*∈{1,…,M} such that K_t_α m^*/M≤ h-1 for at least m^* hypotheses. Since K_t_α m^*/M≤ 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_α m^*/M such that K_t_α m^*/M≤ 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 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.










§ REAL DATA ANALYSIS



<cit.> evaluated the neural responses in the higher visual cortex to natural scene images using fMRI data from the Natural Scenes Dataset <cit.>. For each considered voxel in each of the eight participants, resulting in a total number of more than 150000 voxels, they modeled the voxel response by texture statistics <cit.> 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 α=0.01. 

Due to the large number of hypotheses and refitting step for each permutation, this problem requires an enormous computational effort. <cit.> 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 very low for testing such a large number of hypotheses at level α=0.01, but it worked out since there were many very small p-values. The proportion of rejections and average number of permutations per hypothesis obtained by the different strategies are illustrated in Table <ref>. 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 <cit.> we changed δ to 0.001 due to the lower α in this case.



The results show that, although <cit.> have already chosen a rather low number of permutations, our sequential methods were able to further reduce it significantly while rejecting a similar amount of hypotheses. Also note that the binomial mixture strategy could have continued sampling for hypotheses where the data looks promising after generating the first 1000 permutations to possibly obtain further rejections. Moreover, <cit.> were lucky that their risky approach of sampling just 1000 permutations worked out and the brain regions showed the expected response. If the proportion of low p-values was smaller, e.g., if less than 10% of the p-values were smaller than α, then no hypothesis could have been rejected 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. Note that in this case the AMT algorithm (for a maximum number of B=1000 permutations) would be powerless as well. 






§ CONCLUSION

There are several advantages of using the proposed sequential permutation p-values in large-scale multiple testing problems:

    
  * They automatically adapt the number of permutations to the number of rejections and thus the proportion of false hypotheses and strengths of the alternatives. For this reason, there is no need to adapt the parameter choice to the unknown data generating process. We found h=10 for the anytime-valid BC method and b=0.9 for the binomial mixture strategy to be universally applicable. Depending on whether the main focus lies in increasing precision or reducing the number of permutations, h and b could also be chosen larger or smaller, respectively.
    
  * Due to their ability to adapt to the actual significance level obtained by the multiple testing procedure and stopping early as soon as a decision can be made, they reduce the required number of permutations by orders of magnitude. In particular, simulations showed that there is no fixed number of permutations such that the classical permutation p-value performs similar as the sequential permutation p-values. The anytime-valid version of the Besag-Clifford method performed best in terms of power and required number of permutations.
    
  * They allow to report decisions early such that one can start with writing a paper or identifying follow up work before the entire testing process is finished. In particular, the binomial mixture strategy was found to be useful for such a proceeding.
    
  * During computing times of several months or more, the initial goals and ideas might change based on external information or the data observed so far. The sequential permutation tests allow to adapt the stopping time and even the betting strategy interactively based on the data. It should be noted that the betting strategy can only be changed for the upcoming permutations and not for the permutations that have already been observed.



 In this paper, we showed how anytime-valid permutation tests can be used with p-value based multiple testing procedures. Another type of multiple testing procedures does not rely on individual p-values and explicitly uses the permutation tests to adapt to the unknown dependence structure of the data. This includes the famous MaxT approach for FWER control by Westfall & Young <cit.>, which is particularly powerful when there is a large positive correlation between the test statistics. In Appendix <ref>, we derive sequential versions of such permutation based multiple testing procedures for FWER and simultaneous FDP control. These permutation based multiple testing procedures are based on the closure principle <cit.> and therefore only consist of tests at level α, which reduces the need for sequential permutation tests. Still, we believe that these can be useful in certain applications.
 
 
 






 §.§ Acknowledgments

The authors thank Leila Wehbe and Timothy Barry 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.



plainnat






§ SEQUENTIAL PERMUTATION BASED MULTIPLE TESTING PROCEDURES


 We consider the same setting as in Section <ref>. However, instead of defining the p-values as in (<ref>), we construct a level-α permutation test ϕ_I for each intersection hypothesis H_0^I=⋂_i∈ i H_0^i, I⊆{1,…,M}, and then use the closure principle to obtain decisions for individual hypotheses.  There are two versions of the closure principle. The initial version was proposed by <cit.> who showed that the rejection set defined by 
 
    R={i∈{1,…,M}: ϕ_I=1  for all  I  with  i∈ I}

controls the FWER. However, <cit.> showed that the closure principle can also be used for simultaneous true discovery control. In particular, they proved that

    ℙ(d(S)≤ |S∩ I_1|  for all  S⊆{1,…,M}) ≥ 1-α,

where d(S)min{|S∖ I| : I⊆ℕ, ϕ_I=0 } and I_1={1,…,M}∖ I_0.
In an earlier paper, <cit.> proposed an equivalent approach without using the closure principle. 

Most permutation based multiple testing procedures are based on choosing some combination function C_I for each I⊆ M, defining the intersection tests ϕ_I by the classical permutation p-value applied on C_I((Y_0^i)_i∈ I), C_I((Y_1^i)_i∈ I), …, C_I((Y_B^i)_i∈ I) and then applying the closure principle with the intersection tests ϕ_I. To make sure that the ϕ_I are indeed level-α tests such that this yields a valid procedure, the following assumption is usually made.









 
     The vectors of test statistics corresponding to true hypotheses (Y_0^i)_i∈ I_0, (Y_1^i)_i∈ I_0, …, are jointly exchangeable.


This general approach encompasses many existing permutation based multiple testing methods.
For example, the MaxT approach for FWER control <cit.> is obtained by C_I((Y_j^i)_i∈ I)=max_i∈ I Y_j^i. <cit.> focus on procedures providing true discovery control with C_I((Y_j^i)_i∈ I)=∑_i∈ I Y_j^i. Furthermore, the method by <cit.>, which uniformly improves the popular significance analysis of microarrays (SAM) procedure <cit.>, is obtained by C_I((Y_j^i)_i∈ I)=#{i∈ I:Y_j^i ∈ D^i} for some prespecified sets D^i.

It is straightforward to derive sequential versions of these multiple testing methods by replacing the classical permutation tests ϕ_I by the anytime-valid ones introduced in Section <ref>. We summarized this general approach in Algorithm <ref>.



For a large number of hypotheses M this general approach is computationally infeasible. For this reason, short cuts have been proposed for specific choices of the combination function C_I <cit.>. In case of the MaxT approach, the entire closed test can be performed with a maximum number of M intersection tests <cit.>. It might be the case that using these short cuts only works for specific betting strategies. For example, the MaxT short cut can be used with all betting strategies with nonincreasing wealth for an increasing number of losses. The detailed procedure for applying the MaxT approach with the binomial mixture strategy is provided in Algorithm <ref>.






























§ OMITTED PROOFS





Let Y, Z and X be real valued random variables on (Ω, 𝒜, 𝒫).  If Y is weakly PRDS on X, Z strongly PRDS on Y, and Z and X independent conditional on Y, then Z is weakly PRDS on X.


Let D be an increasing set and x^*≥ x, then
    
    ℙ(Z ∈ D| X≤ x)    = 𝔼_Y| X≤ x[ℙ(Z∈ D| Y, X≤ x)] 
       = 𝔼_Y| X≤ x[ℙ(Z∈ D| Y)] 
       = ∫_Ωℙ(Z∈ D| Y=y) dF_Y| X≤ x(y) 
       ≤∫_Ωℙ(Z∈ D| Y=y) dF_Y| X≤ x^*(y) 
       = ℙ(Z ∈ D| X≤ x^*),

    where the inequality follows from the fact that P(Z∈ D| Y=y) is increasing in y and F_Y| X≤ x^*(y)≤ F_Y| X≤ x(y) for all y.




    Let x,x^*∈ [0,1] with x^*≥ x, D⊆ [0,1]^M be an increasing set and i∈ I_0 be arbitrary but fixed. In the following, we just write 𝔭 and  𝔭^i instead of 𝔭^avBC and 𝔭^avBC,i, respectively. We want to show that 
    ℙ(
    𝔭_γ(h)^-i∈ D|𝔭_γ_i(h)^i≤ x)≤ℙ(𝔭_γ(h)^-i∈ D|𝔭_γ_i(h)^i≤ x^*),

where 𝔭_γ(h)^-i=
(𝔭_γ_1(h)^1,…,𝔭_γ_i-1(h)^i-1, 𝔭_γ_i+1(h)^i+1,…, 𝔭_γ_M(h)^M).
    The proof mainly consists of showing the following two claims.
     
         
  * ^i is weakly PRDS on 𝔭_γ_i(h)^i.
         
  * 𝔭_γ(h)^-i is strongly PRDS on 𝔭_lim.
     
 Since we assumed that 𝔭_lim is strongly PRDS on I_0, the first claim and Lemma <ref> imply that 𝔭_lim is weakly PRDS on 𝔭_γ_i(h)^i. Together with the second claim and the fact that 𝔭_γ(h)^-i is independent of 𝔭_γ_i(h)^i conditional on 𝔭_lim (since we sample independently for all hypotheses), the final proposition follows by Lemma <ref>.
     
 We start with proving the first claim. Note that 𝔭_γ_i(h)^i≤ x iff K_t_x≤ h-1,
where t_x=⌈ h/x ⌉ -1.  With this and the fact that K_t_x^i|^i=p follows a binomial distribution with size parameter t_x and probability p, we obtain


    ℙ(^i≤ p^* |𝔭_γ_i(h)^i≤ x)    = ∫_0^p^*ℙ( 𝔭_γ_i(h)^i≤ x |^i= p  ) dp /∫_0^1ℙ( 𝔭_γ_i(h)^i≤ x |^i= p  ) dp 
       = ∑_k=0^h-1t_x  k∫_0^p^* p^k (1-p)^t_x-k dp /∑_k=0^h-1t_x  k∫_0^1 p^k (1-p)^t_x-k dp  
       = ∑_k=0^h-1  (1-Bin(k;t_x+1,p^*))
       /h
       ,

where Bin(k;t_x+1,p^*) is the CDF of a binomial distribution with size parameter t_x+1 and probability p^*. Since t_x is decreasing in x, ℙ(^i≤ p^* |𝔭_γ_i(h)^i≤ x) is decreasing in x as well.
 

For the second claim, note that 𝔭_γ_j(h)^j is independent of 𝔭_lim^-j and 𝔭_γ(h)^-j conditional on ^j for all j∈{1,…, M}. Therefore, it only remains to show that 𝔭_γ_j(h)^j is strongly PRDS on ^j. Since K_t_u^j|^j=p follows a binomial distribution with probability p, it is immediately implied by the fact that the CDF of a binomial distribution is decreasing in its probability parameter.









§ DETAILED ALGORITHMS




























































