# Control of the False Discovery Proportion

Authors: Yongchao Ge; Xiaochun Li
Journal: Journal of Probability and Statistics
 (2012)

---

## Abstract

Consider the multiple testing problem of testingm null hypotheses H1,…,Hm, among which m0 hypotheses are truly null. Given the P-values for each hypothesis, the question of interest is
how to combine the P-values to find out which hypotheses are false nulls and possibly to make
a statistical inference on m0. Benjamini and Hochberg proposed a classical procedure that can
control the false discovery rate (FDR). The FDR control is a little bit unsatisfactory in that it
only concerns the expectation of the false discovery proportion (FDP). The control of the actual
random variable FDP has recently drawn much attention. For any level 1−α, this paper proposes
a procedure to construct an upper prediction bound (UPB) for the FDP for a fixed rejection
region. When 1−α=50%, our procedure is very close to the classical Benjamini and Hochberg
procedure. Simultaneous UPBs for all rejection regions' FDPs and the upper confidence bound
for the unknown m0 are presented consequently. This new proposed procedure works for finite
samples and hence avoids the slow convergence problem of the asymptotic theory.


---

## Body

## 1. Introduction

In this paper, we consider the problem of testingm null hypotheses H1,…,Hm, among which m0 hypotheses are truly null. We shall assume that P-values are available for individual hypotheses. In a seminal paper, Benjamini and Hochberg [1] proposed the false discovery rate (FDR) as an alternative to the classically defined family-wise error rate (FWER). The proposed FDR achieves a good balance between the P-value itself and the FWER correction [2]; the former may give too many false positives, and the latter may give too many false negatives. However, the control of the FDR is a little bit unsatisfactory in that it only concerns the expectation of the false discovery proportion (FDP). In practice, researchers may be interested in more detailed statistical inference on the actual random variable FDP, not just its expectation. The goal of this paper is to provide a simple procedure to control the FDP.Let us first introduce some notation. Givenm hypotheses H1,…,Hm, let the complete index set be M={1,…,m}, M0 the unknown subset of M for which the null hypotheses are true, and M1=M∖M0 the subset for which null hypotheses are false. Denote that  m0=|M0|,m1=|M1|, where |·| denotes the cardinality of a set. The P-values for testing the m hypotheses are P1,…,Pm. A fixed rejection region for the P-values can conveniently be taken as [0,t] (0<t<1). The value of t could be 0.05, for example. Define the number Rt of all rejected hypotheses and the number Vt of falsely rejected hypotheses, respectively,(1.1)Rt=∑i=1mI(Pi≤t),Vt=∑i∈M0I(Pi≤t).
Following the notation of Korn et al. [3], and Genovese and Wasserman [4], Lehmann and Romano [5], the false discovery proportion is defined to be the proportion of falsely rejected null hypotheses among the rejected ones,(1.2)Qt=VtRt,
where the ratio is defined to be zero when the denominator is zero. For a given fixed rejection region [0,t], Rt, Vt, and Qt are random variables. R, V, and Q will be shorthand for Rt, Vt, and Qt respectively, when the rejection region [0,t] is clear from the context. The false discovery rate of Benjamini and Hochberg [1] is(1.3)FDR=E(Q).A good understanding ofQ will lead investigators to pick an appropriate rejection region [0,t] of P-values. As Q is an unobservable random variable depending on the observed P-values and the rejection region [0,t], the quantity FDR just describes the expectation of this random variable Q. One way to have a more detailed statistical inference on the random variable Q is to derive its distribution, which is very difficult unless a strong assumption can be imposed on the P-values from the false null hypotheses. A conservative approach is to compute an upper prediction bound for Q so that we can safeguard against excessive type I errors. In Section 2, for a fixed rejection region [0,t], we can compute an upper prediction bound (UPB) Q¯1-α(t) for Qt such that(1.4)pr(Qt≤Q¯1-α(t))≥1-α.IfQ¯1-α(t) had been a nonrandom variable, then it should be always no less than the 1-α quantile of the random variable Qt. When 1-α=50%, our procedure is very close to the classical BH procedure of Benjamini and Hochberg [1]. In other words, the BH procedure gives us an approximate 50% upper prediction bound (UPB) for Q. With different degrees of being conservative, one should take 1-α at 0.9, 0.95, and 0.99 to ensure high coverage of the false discovery proportion. We also describe how to compute an upper confidence bound (UCB) for m0, the number of true null hypotheses. The UCB for m0 can be used to improve the estimate Q¯1-α(t). In practice, the rejection region [0,t] needs to be adapted to the actual dataset. In Section 3, we give a procedure to construct an upper prediction band for Qt for all t∈(0,1), and this upper prediction band can be used to pick a data-defined rejection region [0,τ] of P-values such that the false discovery proportion Q can be controlled at target level γ with prediction accuracy 1-α, that is,(1.5)pr(Qτ≤γ)≥1-α.Thus with probability at least1-α, the value of Q is γ or less. For the independent true null P-values, Genovese and Wasserman [4], Meinshausen and Rice [6] also worked on the control of the FDP in the sense of the above equation. However, their results are based on asymptotic theory, while our focus is on the finite-sample results and avoids the slow convergence problem of the asymptotic theory. Other works such as Lehmann and Romano [5], Romano and Shaikh [7, 8],  and van der Laan et al. [9] proposed procedures that allow dependence in the P-values but have potentially lost statistical power as the dependence information is not exploited. Section 4 presents a focused statistical inference by restricting the rejection regions onto {[0,t]:t∈[t0,t0′]}, which unifies the results of Sections 2 and 3. Section 5 generalizes the results from independent data to less-independent situations, and Section 6 gives our discussion.


## 2. Finding a1-α UPB for the False Discovery Proportion for a Fixed Rejection Region

For the sake of simplicity, we will first assume that theP-values from the true null hypotheses are following mutually independently uniform distribution U[0,1]. We have no further assumptions on the P-values from false null hypotheses. This assumption is the same as in Benjamini  and Hochberg [1]. In Section 5 we will generalize the result to less independent situations. For a fixed rejection region [0,t] of the P-values, we would like to find the 1-α upper prediction bound (UPB) for the false discovery proportion Qt. As we mentioned in Section 1, the distribution of Qt is unknown. However, for any given experimental data, the total number of rejections, Rt, can be easily obtained by (1.1). Under the assumption that true null P-values are independently distributed as U[0,1], Vt has a binomial distribution Bin(m0,t). Let Ui,i=1,…,m be random variables mutually independently distributed as U[0,1], and Nm0,t=∑i=1m0I(Ui≤t) distributed as Bin(m0,t), hence,(2.1)Vt=dNm0,t.The1-α quantile for Vt is the 1-α quantile C1-α(m0,t) of the distribution Bin(m0,t). Here C1-α(m0,t) is defined as(2.2)C1-α(m0,t)=min⁡{k:pr(Nm0,t≤k)≥1-α}.
As Rt can be computed from the observed data, a 1-α UPB for Qt can be estimated by(2.3)Q¯1-α(m0,t)=C1-α(m0,t)Rt.Lemma 2.1.
For any given0≤m1≤m2≤m, (a)
C1-α(m1,t)≤C1-α(m2,t),(b)
m1-C1-α(m1,t)≤m2-C1-α(m2,t), and(c)
letg(k)=C1-α(k,t) and h(k)=k-C1-α(k,t). The values that g(k+1)-g(k) and h(k+1)-h(k) take can only be zero or one.Proof.
By noting thatNm1,t≤Nm2,t when m1≤m2, we have pr(Nm1,t≤k)≥pr(Nm2,t≤k). Applying the definition of C1-α, we obtain the result for part (a).
Note that(2.4)C1-α(m1,t)=min⁡{k:pr(Nm1,t≤k)≥1-α}=min⁡{k:pr(m1-Nm1,t≥m1-k)≥1-α}=m1-max⁡{k′:pr(m1-Nm1,t≥k′)≥1-α}(usek′toreplacem1-k).
Therefore,
(2.5)m1-C1-α(m1,t)=max⁡{k′:pr(m1-Nm1,t≥k′)≥1-α}.
Similarly, we can obtain that
(2.6)m2-C1-α(m2,t)=max⁡{k′:pr(m2-Nm2,t≥k′)≥1-α}.
By noting that
(2.7)m1-Nm1,t=∑i=1m1I(Ui>t)≤∑i=1m2I(Ui>t)=m2-Nm2,t,
we have
(2.8)pr(m1-Nm1,t≥k′)≤pr(m2-Nm2,t≥k′). Combining this with (2.5) and (2.6) leads to the result for part (b).
Parts (a) and (b), respectively, say thatg(k) and h(k) are both increasing functions of k. Simple algebra can establish that {g(k+1)-g(k)}+{(h(k+1)-h(k)}=1. Both {g(k+1)-g(k)} and {h(k+1)-h(k)} are nonnegative due to the increasing property of functions g(k) and h(k), and hence 0≤g(k+1)-g(k)≤1 and 0≤h(k+1)-h(k)≤1. The only values that {g(k+1)-g(k)} and {h(k+1)-h(k)} take can only be zero and one as functions g(k) and h(k) only take integer values. Thus, we complete the proof for part (c).Lemma 2.2.
For any given t,0<t<1, Q¯1-α(m0,t) of (2.3) is a 1-α UPB for the false discovery proportion Qt, that is,
(2.9)pr(Qt≤Q¯1-α(m0,t))≥1-α.The proof is straightforward by using the fact thatVt=dNm0,t. We have(2.10)pr(Qt≤Q¯1-α(m0,t))=pr(VtRt≤C1-α(m0,t)Rt,Rt>0)+pr(Rt=0)=pr(Vt≤C1-α(m0,t),Rt>0)+pr(Rt=0)=pr(Vt≤C1-α(m0,t))≥1-α.
In the third line, we have used the fact that {Vt≤C1-α(m0,t)andRt=0} is the same as the set {Rt=0}, which is obtained by noting that Vt must be zero when Rt=0. Following this proof, we can easily see that(2.11){Qt≤Q¯1-α(m0,t)}={Vt≤C1-α(m0,t)}.
The basic construction of Q¯1-α(m0,t) in (2.3) is the idea central to formulating prediction inference for Qt. In practice, m0 is an unknown parameter. The most conservative approach is to replace m0 with m, in which case we obtain a conservative1-α UPB for Qt. The independence assumption among true null P-values can be used to give a confidence inference for m0; thus, we can find a better estimate of the UPB for Qt. For any given 0<λ<1, a 1-α UCB for m0 is given by(2.12)m¯0,1-α(λ)={max⁡k=0,…,m-1{k:h(k)=m-Rλ}ifh(m)<m-Rλ,motherwise,
where h(k)=k-C1-α(k,λ) as defined in Lemma 2.1(c). Since h(0)=0 and h(k+1)-h(k) takes value of only zero and one, there exists at least one k,k∈{0,…,m-1} such that h(k)=m-Rλ when h(m)<m-Rλ. Therefore, m¯0,1-α(λ) in (2.12) is well defined. The parameter λ in (2.12) is used to construct a UCB for m0; more discussion about it can be seen in  Remark 2.6 of the following theorem.Theorem 2.3.
(a)m¯0,1-α(λ) is a conservative 1-α UCB for m0, that is,
(2.13)pr(m0≤m̅0,1-α(λ))≥1-α.
(b) Especially, ifλ takes the same value as t in the P-value rejection region, then
(2.14)pr(Qt≤Q¯1-α(m¯0,1-α(t),t),m0≤m¯0,1-α(t))≥1-α.Proof.
Usem¯0 as a shorthand of m¯0,1-α(λ) in this proof. We want to establish that
(2.15){m0≤m¯0}={h(m0)≤h(m¯0)}.
The fact that function h(k) is increasing in k leads to {m0≤m¯0}⊆{h(m0)≤h(m¯0)}. On the other hand, if m0>m¯0, then m¯0 is strictly less than m, and we must have h(m¯0)=m-Rλ according to (2.12). m¯0 is the maximum of k such that h(k)=m-Rλ, and hence h(m0)≠m-Rλ=h(m¯0) as m0>m¯0. The increasing property of h(k) leads to h(m0)≥h(m¯0). Combining this with h(m0)≠h(m¯0), we obtain that h(m0)>h(m¯0); therefore, we conclude
(2.16){m0>m¯0}⊆{h(m0)>h(m¯0)}
and complete the proof of (2.15).
Note that(2.17){h(m0)≤h(m¯0)}={h(m0)≤m-Rλ,m¯0<m}∪{h(m0)≤h(m),m¯0=m}={h(m0)≤m-Rλ,m¯0<m}∪{m¯0=m}(h(m0)≤h(m)always holds)⊇{h(m0)≤m-Rλ},
we have
(2.18)pr(m0≤m¯0)=pr(h(m0)≤h(m¯0)≥pr(m0-C1-α(m0,λ)≤m-Rλ))≥pr(m0-C1-α(m0,λ)≤m0-Vλ)(usingm-Rλ≥m0-Vλ)=pr(Vλ≤C1-α(m0,λ))≥1-α(NotethatVλ=dNm0,λ). Hence, we have the proof of part (a). When λ is set to be t, we have that the set {m0≤m¯0,1-α(t)} contains {Vt≤C1-α(m0,t)} from the above derivation, and that {Qt≤Q¯1-α(m0,t)}={Vt≤C1-α(m0,t)} from the derivation of Lemma 2.2. Therefore,
(2.19){Qt≤Q¯1-α(m¯0,1-α(t),t),m0≤m¯0,1-α(t)}⊇{Qt≤Q¯1-α(m0,t),m0≤m¯0,1-α(t)}⊇{Vt≤C1-α(m0,t)}.
Thus, we have the proof of part (b) of this theorem.Remark 2.4.
Theorem2.3 gives researchers a good sense of the total number m0 of true null hypotheses. Other papers, for example, Storey et al. [10], Benjamini and Hochberg [11], and Langaas et al. [12], gave only point estimates of m0 or π0=m0/m. Part (a) gives a confidence inference for m0, and part (b) gives a simultaneous statement for the Qt and m0, which is more interesting. Meinshausen [13] gives a confidence for m0 by using resampling methods, while ours exploited the independence information so that it works for finite samples.Remark 2.5.
Theorem 2.3 of Göb [14] implies that for a binomial distribution, the difference between the median and mean is less than 1, that is, |C0.5(m0,t)-m0t|<1. From (2.3), we know the 50% UPB for the Qt can be estimated by C0.5(m0,t)/Rt. This UPB for Qt is very close to m0t/Rt with a difference smaller than 1/Rt. Replacing m0 by m in m0t/Rt is equivalent to the classical BH procedure. For a very large Rt, the term 1/Rt can be ignored, and the BH procedure offers an approximate estimate of the 50% UPB for Qt.Remark 2.6.
Whenk is large, the distribution Bin(k,λ) can be closely approximated by N(kλ,kλ(1-λ)). Let z1-α be the 1-α quantile of a standard normal distribution. After some algebraic manipulations, we obtain a 1-α UCB for m0(2.20)m¯0,1-α(λ)≈{12(1-λ)(z1-αλ(1-λ)+z1-α2λ(1-λ)+4(m-Rλ)(1-λ))}2.
Taking 1-α=0.5, we have (m-Rλ)/(1-λ), which is equivalent to (2.3) of Storey et al. [10]. For most practical applications, one can set the value of λ=0.5. Fine tuning of the parameter λ will be discussed in Section 3.3.
When the rejection region[0,t] is small, the UCB for m0 obtained from part (b) of Theorem 2.3 may be too conservative. It may be advantageous to have separate values for λ and t. Part (a) of Theorem 2.3 implies the following.Corollary 2.7.
Replacingm0 by its the upper confidence bound m¯0,1-α2(λ) and α by α1 in (2.3), we define
(2.21)Q¯1-α1,1-α2(t,λ)=C1-α1(m¯0,1-α2(λ),t)Rt.
Then, Q¯1-α1,1-α2(t,λ) is a conservative 1-α(α=α1+α2) UPB for the false discovery proportion Qt.


## 3. Upper Prediction Bounds and Simultaneous Inferences

### 3.1. The Setup

In Section2, the UPBs for Q are only valid for a fixed rejection region [0,t] of P-values. In practice, researchers will not fix the rejected region [0,t] but adapt it to the actual data. The logic is the same as with single hypothesis testing. In single hypothesis testing with nested rejection regions {Γα:0<α<1}, for an observed statistic T, one will find the rejection region that contains the observed statistic with the smallest type I error α, that is,(3.1)P-value(T)=min⁡{α:T∈Γα}.The same logic can be applied to our false discovery proportion. In this case, we will try to find the largest rejection region[0,t] such that the false discovery proportion Q is not more than γ, say 10%, with probability 1-α. Define(3.2)τ=max⁡{t:Q¯1-α1,1-α2(t,λ)≤γ}.
We then reject any hypothesis whose P-value is no greater than τ. If τ is independent of Q and Q¯, then we can expect that(3.3)pr(Qτ≤γ)≥pr(Qτ≤Q¯1-α1,1-α2(τ,λ))≥1-α.Asymptotically,τ and (Q,Q¯) may be independent: this question is open for future research. To overcome the independence assumption of τ and (Q,Q¯), we seek an alternative approach: to find simultaneous UPBs for all rejection regions [0,t], t∈(0,1), that is, to find an upper prediction bandQ¯t such that(3.4)pr(Qt≤Q¯tfort∈(0,1))≥1-α.Hence we have the simultaneous inferences onQt for each rejection region [0,t],t∈(0,1). Following the definition of C1-α(n,t) in (2.2) to construct the UPB for Qt, we want to define the simultaneous critical values of Nn,t. Using the distribution of maxt∈(0,1)Nn,t is unwise as maxt∈(0,1)Nn,t=Nn,1, which takes value n with probability one. A better approach is to center Nn,t, that is,(3.5)sup⁡t∈(0,1)(Nn,t-nt).This leads to a test statistic related to the Kolmogorov-Smirnov test statistic, which gives an upper confidence band for a cumulative distribution functionF(x). It turns out that  this method leads to very high UPBs when t is close to zero or one. Therefore, we normalize Nn,t, that is,(3.6)Z̃n=sup⁡t∈(0,1)(Nn,t-nt)nt(1-t).
Note that Z̃n is continuously distributed even though each Nn,t is discretely distributed. Let z̃1-α(n) be the 1-α quantile of Z̃n, that is,(3.7)pr(Z̃n≤z̃1-α(n))=1-α.
We can then redefine Q¯ as(3.8)Q̃1-α(m0,t)=m0t+z̃1-α(m0)m0t(1-t)Rt.
Corresponding to Lemma 2.2 and Corollary 2.7, we have similar results below.Corollary 3.1.
For any given0<t<1, Q̃1-α(m0,t) of (3.8) is an exactly 1-α upper prediction band for the false discovery proportion Qt, that is,
(3.9)pr(Qt≤Q̃1-α(m0,t)∀t∈(0,1))=1-α.Corollary 3.2.
Denote thatm¯0=m¯0,1-α2(λ). Define
(3.10)Q̃1-α1,1-α2(t,λ)=m¯0t+z̃1-α1(m¯0)m¯0t(1-t)Rt.
Let α=α1+α2. Then Q̃1-α1,1-α2(t,λ) is a conservative 1-α upper prediction band for the false discovery proportion Qt, that is,
(3.11)pr(Qt≤Q̃1-α1,1-α2(t,λ)∀t∈(0,1))≥1-α.Remark 3.3.
Using the same idea as in the proof of Lemma2.2, the proof of the above corollaries is straightforward after converting the comparison between Q and Q̃ to the comparison between Vt and m0t+z̃1-α(m0)m0t(1-t). This conversion provides a powerful tool for understanding the false discovery proportion.Remark 3.4.
The formulation ofQ̃1-α(m0,t) in (3.8) is motivated by the normal approximation of Nn,t. But our definition of Q̃1-α(m0,t) gives exact UPBs simultaneously for all t∈(0,1) due to the exactness of the quantile z̃1-α.Remark 3.5.
Meinshausen and Rice [6] and Donoho and Jin [15] also utilize the empirical process Z̃n. However, they focus on the asymptotic theory for Z̃n, which may face the slow convergence problem described in the next section. Our focus is on the finite sample control.Remark 3.6.
Letf(t)=(k-nt)/nt(1-t). After some simplifications, we have f′(t)·(nt(1-t))3=-n2t(1-t)-1/2·(k-nt)[n(1-2t)]=-n/2·[k(1-t)+(n-k)t], and then f(t) is a decreasing function in t,  0<t<1 for 0≤k≤n. Equation (3.6) can be simplified to be
(3.12)max⁡k=1,…,nsup⁡t∈[U(k),U(k+1))k-ntnt(1-t)=max⁡k=1,…,nk-nU(k)nU(k)(1-U(k)),where U(k) is the kth smallest ordered one among the n samples of U[0,1] distribution. This formulation facilitates the computation of the distribution of Z̃n by Monte Carlo methods. The standard error associated with the Monte Carlo simulations in computing the probability in (3.7) is no greater than α(1-α)/B, where B is the number of simulations.


### 3.2. Computing the Distribution ofZ̃n

In order to make simultaneous inferences, we need to know the distribution ofZ̃n defined in (3.6). Example 1 of Jaeschke [16] showed that asymptotically, for any x,(3.13)lim⁡n→∞pr(Z̃n≤x+2lnlnn+(1/2)lnlnlnn-(1/2)ln⁡π2lnlnn)=exp⁡[-exp⁡(-x)].
This implies that Z̃n/2lnlnn converges to 1 in probability as n goes to ∞. Jaeschke [16] claimed that this probability convergence is of almost no practical use. This is where we need to be cautious using asymptotic results. Figure 1 shows the poor approximation of the asymptotic result, even for a very large n=105. Noe and Vandewiele [17] gave an iterative algorithm to compute the exact probability pr(Z̃n≤z̃). Their algorithm is only good for very small n due to the computational time and propagation of precision errors in representing real numbers in computer. Equation (24) of their paper gives an approximate formula for n=1,…,100,(3.14)pr(Z̃n≤z̃)≈1-(z̃)-2-(2-3n-1)(z̃)-4-(10-57n-1+48n-2)(z̃)-6-(74-1021n-1+2743n-2-1797n-3)(z̃)-8-(706-19123n-1+111905n-2-213619n-3+120132n-4)(z̃)-10.
This approximation is very good for z̃≥4 but is away from the true probability when z̃<4. For our applications, the 50% quantile (median) of Z̃ is very useful, but the approximation of (3.14) is poor there.Figure 1
Plot of the probabilitypr(Z̃n≤z̃) for 1≤z̃≤20 with n=105. The x-axis is plotted on a log scale. The black solid curve is computed from 106 Monte Carlo simulations. The red-dashed curve is based on (3.13). The green dot-dashed curve is computed from (3.14).In order to overcome the above poor approximation, we propose to use the Monte Carlo method to obtain the probabilitypr(Z̃≤z̃). Figures 2 and 3 give the probability for z̃∈[1,20] with 106 simulations for n=1,10,102,103,104, and 105. The Monto Carlo method generated quantiles z̃ for n=1,…,100 are almost the same as those quantiles that were able to be computed by the exact algorithm in Noe and Vandewiele [17]. Our two figures show that the distribution of Z̃n does not change dramatically from n=1,…,105. This property is beneficial for a multiple testing problem with large number of hypotheses, as it will not be overpenalized. Table 1 gives the quantiles of Z̃n with n=105.Table 1
The quantiles ofZ̃n of Figure 3, where n=105. This table is estimated by 106 Monto Carlo simulations. The column 1-α gives the probabilities,  and the column z̃ gives the quantiles of Z̃n.
1-αz̃1-αz̃1-αz̃1-αz̃0.52.370.692.760.883.520.96755.740.512.390.72.780.893.610.975.950.522.410.712.810.93.70.97256.190.532.420.722.840.9053.750.9756.470.542.440.732.860.913.810.97756.780.552.460.742.890.9153.880.987.230.562.480.752.920.923.960.98257.730.572.50.762.950.9254.050.9858.310.582.520.772.990.934.150.98759.040.592.540.783.020.9354.260.9910.040.62.560.793.060.944.40.992511.560.612.580.83.090.9454.550.99514.240.622.60.813.130.954.730.995514.870.632.620.823.180.95254.840.99615.750.642.640.833.220.9554.950.996516.60.652.660.843.270.95755.080.99717.970.662.690.853.320.965.220.997519.690.672.710.863.390.96255.360.682.730.873.450.9655.55Figure 2
Plot of the probabilitypr(Z̃n≤z̃) for 1≤z̃≤20. The probability is computed with 106 Monte Carlo simulations. The curves from the top to the bottom correspond to n=1,10,102,103,104, and 105.Figure 3
A blowup of Figure2. This shows the part of the probability pr(Z̃n≤z̃) that is close to 1 as 1−P is drawn on a log scale for the y-axis.


### 3.3. More about the Upper Confidence Bound form0

In computing the UCB form0 and consequently the UPB for Qt, we rely on the unspecified parameter λ. A conventional choice of λ is 0.5. It is tempting to use minλ∈(0,1)m¯0,1-α(λ) as the best UCB for m0. This approach should be avoided as it may lead to an overoptimistic UCB. We can use the same idea in computing the simultaneous upper prediction bounds for Qt to find an UCB for m0. Equation (2.20) motivates to the following theorem.Theorem 3.7.
Definem̃0,1-α(λ) as
(3.15){12(1-λ)(z̃¯1-α(m)λ(1-λ)+(z̃¯1-α(m))2λ(1-λ)+4(m-Rλ)(1-λ))}2,
where
(3.16)z̃¯1-α(m)=max⁡⁡n=1,…,mz̃1-α(n).
Let
(3.17)m̃0,1-α=min⁡λ∈(0,1)m̃0,1-α(λ).
Using m̃0,1-α to replace m0 in (3.8) results in Q̃1-α(m̃0,1-α,t). We have
(3.18)pr(m0≤m̃0,1-α,Qt≤Q̃1-α(m̃0,1-α,t)∀t∈(0,1))≥1-α.
Thus simultaneously m̃0,1-α is a 1-α UCB for m0 and Q̃ is a 1-α upper prediction band.Proof.
Note that whenx>0,
(3.19)x(1-λ)-xz̃¯1-α(m)λ(1-λ)-(m-Rλ)≤0
if and only if
(3.20)x≤{12(1-λ)(z̃¯1-α(m)λ(1-λ)+(z̃¯1-α(m))2λ(1-λ)+4(m-Rλ)(1-λ))}2.
Therefore,
(3.21){m0≤m̃0,1-α}={m0≤m̃0,1-α(λ)∀λ∈(0,1)}={m0(1-λ)-m0z̃¯1-α(m)λ(1-λ)-(m-Rλ)≤0∀λ∈(0,1)}={max⁡λ∈(0,1)m0(1-λ)-(m-Rλ)m0λ(1-λ)≤z̃¯1-α(m)}⊇{max⁡λ∈(0,1)Vλ-m0λm0λ(1-λ)≤z̃1-α(m0)}.
The last step follows: (i) z̃¯1-α(m) is no less than z̃1-α(n) for any n≤m, and (ii) m-Rλ≥m0-Vλ. The fact (ii) gives
(3.22)m0(1-λ)-(m-Rλ)≤m0(1-λ)-(m0-Vλ)=Vλ-m0λ. Following the same idea as in the proof of Theorem 2.3 part (b), we can show that the set {m0≤m̃0,1-αandQt≤Q̃1-α(m̃0,1-α,t)forallt∈(0,1)} is a superset of {maxt∈(0,1)(Vt-m0t)/m0t(1-t)≤z̃1-α(m0)}. Therefore,
(3.23)pr(m0≤m̃0,1-α,Qt≤Q̃1-α(m̃0,1-α,t)∀t∈(0,1))≥pr(max⁡t∈(0,1)Vt-m0tm0t(1-t)≤z̃1-α(m0))=1-α(NotethatVt=dNm0,t).Readers should note that the maximum quantile z̃¯1-α(m) defined in (3.16) is only used to construct m̃0,1-α, The construction of Q̃1-α(m̃0,1-α,t) itself does not use the maximum but the quantile z̃1-α(m), while Q̃1-α(m̃0,1-α,t) still depends on the maximum quantiles indirectly through m̃0,1-α.


### 3.4. The Algorithm

Putting all these pieces together, we describe the procedure to compute the upper prediction band forQt and the UCB for m0 in Algorithm 1. Note that we have to compute the quantile z̃¯1-α(m) and z̃1-α(m̃0,1-α). This is very time consuming for large m, which is typically from thousands to tens of thousands. The computationally time can be reduced by the following strategies.(1)
After careful study of the two equations (3.8) and (3.15), we find that if we replace all z̃1-α(n) and z̃¯1-α(n) by z̃¯1-α(N), where N≥n, the conclusions of Corollaries 3.1 and 3.2 and Theorem 3.7 still hold.(2)
The quantilez̃1-α(n) is an increasing function of n, as shown by the Monte Carlo simulations in Figures 2 and 3. The rigorous mathematical proof of this finding is open to future research. For practical applications, we can first use Monte Carlo simulations to verify this property for the range of n that is related to the project and then replace all z̃¯1-α(n) by z̃1-α(n).(3)
Figure2 shows that z̃1-α(n) is very close to z̃1-α(N) if n is close to a large N, say more than 100. Therefore, in practical computations, we can first compute and store a representative sequences of the quantiles z̃1-α(n) for n=n1,…,nI, and consequently we can get an upper bound for z̃1-α(n) for n=1,…,m. In computing the quantiles z̃1-α(n), we recommend to have at least 104 Monte Carlo simulations in order to get an accurate quantile computation for tail part. Even with 106 simulations, we still see a small amount of random noise in the tail part in Figure 3.Algorithm 1:Compute the simultaneous UPBs for the false discovery proportion and the UCB for m0.
For any givenα and γ,(1) Compute m̃0,1-α(λ) of (3.15) for some pre-specified λi’s,  say λi=i/1000, for i=1,…,999.(2) Compute m̃0,1-α=minim̃0,1-α(λi). This m̃0,1-α is the 1-α UCB for m0.Ifm̃0,1-α exceeds m, replace it by m.(3) Sort the observed P-values such that P(1)≤⋯≤P(m), and use (3.8) to compute the 1-αsimultaneous UPBs for the false discovery proportionQ, that is, for i=1,…,m,Q̃1-α(P(i))=(1/i)(m̃0,1-αP(i)+z̃1-α(m̃0,1-α)m̃0,1-αP(i)(1-P(i))).IfQ̃1-α(P(i)) exceeds 1, replace it by 1.(4) Compute τ=max{P(i):Q̃1-α(P(i))≤γ},reject the hypotheses whoseP-values are no greater than τ,which ensures that the false discovery proportionQ is not exceeding γ with probability 1-α.


## 3.1. The Setup

In Section2, the UPBs for Q are only valid for a fixed rejection region [0,t] of P-values. In practice, researchers will not fix the rejected region [0,t] but adapt it to the actual data. The logic is the same as with single hypothesis testing. In single hypothesis testing with nested rejection regions {Γα:0<α<1}, for an observed statistic T, one will find the rejection region that contains the observed statistic with the smallest type I error α, that is,(3.1)P-value(T)=min⁡{α:T∈Γα}.The same logic can be applied to our false discovery proportion. In this case, we will try to find the largest rejection region[0,t] such that the false discovery proportion Q is not more than γ, say 10%, with probability 1-α. Define(3.2)τ=max⁡{t:Q¯1-α1,1-α2(t,λ)≤γ}.
We then reject any hypothesis whose P-value is no greater than τ. If τ is independent of Q and Q¯, then we can expect that(3.3)pr(Qτ≤γ)≥pr(Qτ≤Q¯1-α1,1-α2(τ,λ))≥1-α.Asymptotically,τ and (Q,Q¯) may be independent: this question is open for future research. To overcome the independence assumption of τ and (Q,Q¯), we seek an alternative approach: to find simultaneous UPBs for all rejection regions [0,t], t∈(0,1), that is, to find an upper prediction bandQ¯t such that(3.4)pr(Qt≤Q¯tfort∈(0,1))≥1-α.Hence we have the simultaneous inferences onQt for each rejection region [0,t],t∈(0,1). Following the definition of C1-α(n,t) in (2.2) to construct the UPB for Qt, we want to define the simultaneous critical values of Nn,t. Using the distribution of maxt∈(0,1)Nn,t is unwise as maxt∈(0,1)Nn,t=Nn,1, which takes value n with probability one. A better approach is to center Nn,t, that is,(3.5)sup⁡t∈(0,1)(Nn,t-nt).This leads to a test statistic related to the Kolmogorov-Smirnov test statistic, which gives an upper confidence band for a cumulative distribution functionF(x). It turns out that  this method leads to very high UPBs when t is close to zero or one. Therefore, we normalize Nn,t, that is,(3.6)Z̃n=sup⁡t∈(0,1)(Nn,t-nt)nt(1-t).
Note that Z̃n is continuously distributed even though each Nn,t is discretely distributed. Let z̃1-α(n) be the 1-α quantile of Z̃n, that is,(3.7)pr(Z̃n≤z̃1-α(n))=1-α.
We can then redefine Q¯ as(3.8)Q̃1-α(m0,t)=m0t+z̃1-α(m0)m0t(1-t)Rt.
Corresponding to Lemma 2.2 and Corollary 2.7, we have similar results below.Corollary 3.1.
For any given0<t<1, Q̃1-α(m0,t) of (3.8) is an exactly 1-α upper prediction band for the false discovery proportion Qt, that is,
(3.9)pr(Qt≤Q̃1-α(m0,t)∀t∈(0,1))=1-α.Corollary 3.2.
Denote thatm¯0=m¯0,1-α2(λ). Define
(3.10)Q̃1-α1,1-α2(t,λ)=m¯0t+z̃1-α1(m¯0)m¯0t(1-t)Rt.
Let α=α1+α2. Then Q̃1-α1,1-α2(t,λ) is a conservative 1-α upper prediction band for the false discovery proportion Qt, that is,
(3.11)pr(Qt≤Q̃1-α1,1-α2(t,λ)∀t∈(0,1))≥1-α.Remark 3.3.
Using the same idea as in the proof of Lemma2.2, the proof of the above corollaries is straightforward after converting the comparison between Q and Q̃ to the comparison between Vt and m0t+z̃1-α(m0)m0t(1-t). This conversion provides a powerful tool for understanding the false discovery proportion.Remark 3.4.
The formulation ofQ̃1-α(m0,t) in (3.8) is motivated by the normal approximation of Nn,t. But our definition of Q̃1-α(m0,t) gives exact UPBs simultaneously for all t∈(0,1) due to the exactness of the quantile z̃1-α.Remark 3.5.
Meinshausen and Rice [6] and Donoho and Jin [15] also utilize the empirical process Z̃n. However, they focus on the asymptotic theory for Z̃n, which may face the slow convergence problem described in the next section. Our focus is on the finite sample control.Remark 3.6.
Letf(t)=(k-nt)/nt(1-t). After some simplifications, we have f′(t)·(nt(1-t))3=-n2t(1-t)-1/2·(k-nt)[n(1-2t)]=-n/2·[k(1-t)+(n-k)t], and then f(t) is a decreasing function in t,  0<t<1 for 0≤k≤n. Equation (3.6) can be simplified to be
(3.12)max⁡k=1,…,nsup⁡t∈[U(k),U(k+1))k-ntnt(1-t)=max⁡k=1,…,nk-nU(k)nU(k)(1-U(k)),where U(k) is the kth smallest ordered one among the n samples of U[0,1] distribution. This formulation facilitates the computation of the distribution of Z̃n by Monte Carlo methods. The standard error associated with the Monte Carlo simulations in computing the probability in (3.7) is no greater than α(1-α)/B, where B is the number of simulations.


## 3.2. Computing the Distribution ofZ̃n

In order to make simultaneous inferences, we need to know the distribution ofZ̃n defined in (3.6). Example 1 of Jaeschke [16] showed that asymptotically, for any x,(3.13)lim⁡n→∞pr(Z̃n≤x+2lnlnn+(1/2)lnlnlnn-(1/2)ln⁡π2lnlnn)=exp⁡[-exp⁡(-x)].
This implies that Z̃n/2lnlnn converges to 1 in probability as n goes to ∞. Jaeschke [16] claimed that this probability convergence is of almost no practical use. This is where we need to be cautious using asymptotic results. Figure 1 shows the poor approximation of the asymptotic result, even for a very large n=105. Noe and Vandewiele [17] gave an iterative algorithm to compute the exact probability pr(Z̃n≤z̃). Their algorithm is only good for very small n due to the computational time and propagation of precision errors in representing real numbers in computer. Equation (24) of their paper gives an approximate formula for n=1,…,100,(3.14)pr(Z̃n≤z̃)≈1-(z̃)-2-(2-3n-1)(z̃)-4-(10-57n-1+48n-2)(z̃)-6-(74-1021n-1+2743n-2-1797n-3)(z̃)-8-(706-19123n-1+111905n-2-213619n-3+120132n-4)(z̃)-10.
This approximation is very good for z̃≥4 but is away from the true probability when z̃<4. For our applications, the 50% quantile (median) of Z̃ is very useful, but the approximation of (3.14) is poor there.Figure 1
Plot of the probabilitypr(Z̃n≤z̃) for 1≤z̃≤20 with n=105. The x-axis is plotted on a log scale. The black solid curve is computed from 106 Monte Carlo simulations. The red-dashed curve is based on (3.13). The green dot-dashed curve is computed from (3.14).In order to overcome the above poor approximation, we propose to use the Monte Carlo method to obtain the probabilitypr(Z̃≤z̃). Figures 2 and 3 give the probability for z̃∈[1,20] with 106 simulations for n=1,10,102,103,104, and 105. The Monto Carlo method generated quantiles z̃ for n=1,…,100 are almost the same as those quantiles that were able to be computed by the exact algorithm in Noe and Vandewiele [17]. Our two figures show that the distribution of Z̃n does not change dramatically from n=1,…,105. This property is beneficial for a multiple testing problem with large number of hypotheses, as it will not be overpenalized. Table 1 gives the quantiles of Z̃n with n=105.Table 1
The quantiles ofZ̃n of Figure 3, where n=105. This table is estimated by 106 Monto Carlo simulations. The column 1-α gives the probabilities,  and the column z̃ gives the quantiles of Z̃n.
1-αz̃1-αz̃1-αz̃1-αz̃0.52.370.692.760.883.520.96755.740.512.390.72.780.893.610.975.950.522.410.712.810.93.70.97256.190.532.420.722.840.9053.750.9756.470.542.440.732.860.913.810.97756.780.552.460.742.890.9153.880.987.230.562.480.752.920.923.960.98257.730.572.50.762.950.9254.050.9858.310.582.520.772.990.934.150.98759.040.592.540.783.020.9354.260.9910.040.62.560.793.060.944.40.992511.560.612.580.83.090.9454.550.99514.240.622.60.813.130.954.730.995514.870.632.620.823.180.95254.840.99615.750.642.640.833.220.9554.950.996516.60.652.660.843.270.95755.080.99717.970.662.690.853.320.965.220.997519.690.672.710.863.390.96255.360.682.730.873.450.9655.55Figure 2
Plot of the probabilitypr(Z̃n≤z̃) for 1≤z̃≤20. The probability is computed with 106 Monte Carlo simulations. The curves from the top to the bottom correspond to n=1,10,102,103,104, and 105.Figure 3
A blowup of Figure2. This shows the part of the probability pr(Z̃n≤z̃) that is close to 1 as 1−P is drawn on a log scale for the y-axis.


## 3.3. More about the Upper Confidence Bound form0

In computing the UCB form0 and consequently the UPB for Qt, we rely on the unspecified parameter λ. A conventional choice of λ is 0.5. It is tempting to use minλ∈(0,1)m¯0,1-α(λ) as the best UCB for m0. This approach should be avoided as it may lead to an overoptimistic UCB. We can use the same idea in computing the simultaneous upper prediction bounds for Qt to find an UCB for m0. Equation (2.20) motivates to the following theorem.Theorem 3.7.
Definem̃0,1-α(λ) as
(3.15){12(1-λ)(z̃¯1-α(m)λ(1-λ)+(z̃¯1-α(m))2λ(1-λ)+4(m-Rλ)(1-λ))}2,
where
(3.16)z̃¯1-α(m)=max⁡⁡n=1,…,mz̃1-α(n).
Let
(3.17)m̃0,1-α=min⁡λ∈(0,1)m̃0,1-α(λ).
Using m̃0,1-α to replace m0 in (3.8) results in Q̃1-α(m̃0,1-α,t). We have
(3.18)pr(m0≤m̃0,1-α,Qt≤Q̃1-α(m̃0,1-α,t)∀t∈(0,1))≥1-α.
Thus simultaneously m̃0,1-α is a 1-α UCB for m0 and Q̃ is a 1-α upper prediction band.Proof.
Note that whenx>0,
(3.19)x(1-λ)-xz̃¯1-α(m)λ(1-λ)-(m-Rλ)≤0
if and only if
(3.20)x≤{12(1-λ)(z̃¯1-α(m)λ(1-λ)+(z̃¯1-α(m))2λ(1-λ)+4(m-Rλ)(1-λ))}2.
Therefore,
(3.21){m0≤m̃0,1-α}={m0≤m̃0,1-α(λ)∀λ∈(0,1)}={m0(1-λ)-m0z̃¯1-α(m)λ(1-λ)-(m-Rλ)≤0∀λ∈(0,1)}={max⁡λ∈(0,1)m0(1-λ)-(m-Rλ)m0λ(1-λ)≤z̃¯1-α(m)}⊇{max⁡λ∈(0,1)Vλ-m0λm0λ(1-λ)≤z̃1-α(m0)}.
The last step follows: (i) z̃¯1-α(m) is no less than z̃1-α(n) for any n≤m, and (ii) m-Rλ≥m0-Vλ. The fact (ii) gives
(3.22)m0(1-λ)-(m-Rλ)≤m0(1-λ)-(m0-Vλ)=Vλ-m0λ. Following the same idea as in the proof of Theorem 2.3 part (b), we can show that the set {m0≤m̃0,1-αandQt≤Q̃1-α(m̃0,1-α,t)forallt∈(0,1)} is a superset of {maxt∈(0,1)(Vt-m0t)/m0t(1-t)≤z̃1-α(m0)}. Therefore,
(3.23)pr(m0≤m̃0,1-α,Qt≤Q̃1-α(m̃0,1-α,t)∀t∈(0,1))≥pr(max⁡t∈(0,1)Vt-m0tm0t(1-t)≤z̃1-α(m0))=1-α(NotethatVt=dNm0,t).Readers should note that the maximum quantile z̃¯1-α(m) defined in (3.16) is only used to construct m̃0,1-α, The construction of Q̃1-α(m̃0,1-α,t) itself does not use the maximum but the quantile z̃1-α(m), while Q̃1-α(m̃0,1-α,t) still depends on the maximum quantiles indirectly through m̃0,1-α.


## 3.4. The Algorithm

Putting all these pieces together, we describe the procedure to compute the upper prediction band forQt and the UCB for m0 in Algorithm 1. Note that we have to compute the quantile z̃¯1-α(m) and z̃1-α(m̃0,1-α). This is very time consuming for large m, which is typically from thousands to tens of thousands. The computationally time can be reduced by the following strategies.(1)
After careful study of the two equations (3.8) and (3.15), we find that if we replace all z̃1-α(n) and z̃¯1-α(n) by z̃¯1-α(N), where N≥n, the conclusions of Corollaries 3.1 and 3.2 and Theorem 3.7 still hold.(2)
The quantilez̃1-α(n) is an increasing function of n, as shown by the Monte Carlo simulations in Figures 2 and 3. The rigorous mathematical proof of this finding is open to future research. For practical applications, we can first use Monte Carlo simulations to verify this property for the range of n that is related to the project and then replace all z̃¯1-α(n) by z̃1-α(n).(3)
Figure2 shows that z̃1-α(n) is very close to z̃1-α(N) if n is close to a large N, say more than 100. Therefore, in practical computations, we can first compute and store a representative sequences of the quantiles z̃1-α(n) for n=n1,…,nI, and consequently we can get an upper bound for z̃1-α(n) for n=1,…,m. In computing the quantiles z̃1-α(n), we recommend to have at least 104 Monte Carlo simulations in order to get an accurate quantile computation for tail part. Even with 106 simulations, we still see a small amount of random noise in the tail part in Figure 3.Algorithm 1:Compute the simultaneous UPBs for the false discovery proportion and the UCB for m0.
For any givenα and γ,(1) Compute m̃0,1-α(λ) of (3.15) for some pre-specified λi’s,  say λi=i/1000, for i=1,…,999.(2) Compute m̃0,1-α=minim̃0,1-α(λi). This m̃0,1-α is the 1-α UCB for m0.Ifm̃0,1-α exceeds m, replace it by m.(3) Sort the observed P-values such that P(1)≤⋯≤P(m), and use (3.8) to compute the 1-αsimultaneous UPBs for the false discovery proportionQ, that is, for i=1,…,m,Q̃1-α(P(i))=(1/i)(m̃0,1-αP(i)+z̃1-α(m̃0,1-α)m̃0,1-αP(i)(1-P(i))).IfQ̃1-α(P(i)) exceeds 1, replace it by 1.(4) Compute τ=max{P(i):Q̃1-α(P(i))≤γ},reject the hypotheses whoseP-values are no greater than τ,which ensures that the false discovery proportionQ is not exceeding γ with probability 1-α.


## 4. A Focused Inference onQ and m0: A Unified Approach

In many applications, it may be unnecessary to compute the simultaneous UPBs forQt for all t in (0,1) and using m̃1-α(λ) for all λ in (0,1) to derive a 1-α UCB for m0. In most applications, it may be reasonable to restrict the rejections onto {[0,t]:t∈[t0,t0′]}. The t0 can take value of 0.01/m based on Bonferroni FWER control at level 0.01. It is rare to consider a smaller rejection region than this. The t0′ can take value of 0.05 as it is rare to consider a larger rejection region than [0,0.05] even in a single hypothesis testing problem. For the same reason, we can also restrict λ onto [λ0,λ0′] in (3.17). The interval [λ0, λ0′] can be taken as a region close to one as the minimum of m̃0,1-α(λ) is reached when λ is close to 1 [18], but if λ is too close to 1, m̃0,1-α(λ) is not stable. One good choice of λ0 can be 0.8, and λ0′ can be 0.95.The above scenario is a focused inference onQ and m0. The z̃ in Section 3.2 will be a little bit more conservative for us. We can redefine Z̃n* as in the following:(4.1)Z̃n*=max⁡(sup⁡t∈[t0,t0′](Nn,t-nt)nt(1-t),sup⁡λ∈[λ0,λ0′](Nn,λ-nλ)nλ(1-λ)).From thisZ̃n* we can define the 1-α quantile z̃1-α*(n), and derive results similar to Theorem 3.7. Figure 4 shows quantiles z̃1-α*(n) for n=1,…,105 for [t0,t0′]=[0.01/n,0.05] and [λ0,λ0′]=[0.8,0.95]. Table 2 gives the numerical values of z̃1-α*(n) for n=105. It clearly shows that z̃1-α*(n) is around 10% smaller than the unrestricted quantiles z̃1-α*(n). For small values of α, say that  α≤0.01, the former is at least 25% smaller than the latter.Table 2
The quantiles ofZ̃n* of Figure 4, where n=105. The column 1-α gives the probabilities,  and the column z̃* gives the quantiles of Z̃n*.
1-αz̃*1-αz̃*1-αz̃*1-αz̃*0.52.10.692.550.883.380.96755.150.512.120.72.580.893.470.975.280.522.140.712.610.93.570.97255.430.532.160.722.640.9053.620.9755.580.542.180.732.670.913.680.97755.780.552.210.742.70.9153.740.985.990.562.230.752.730.923.820.98256.20.572.250.762.770.9253.890.9856.480.582.270.772.80.933.970.98756.820.592.30.782.840.9354.060.997.230.62.320.792.880.944.180.99257.730.612.350.82.920.9454.30.9958.370.622.370.812.960.954.430.99558.520.632.390.823.010.95254.520.9968.630.642.420.833.060.9554.610.99658.770.652.440.843.110.95754.70.9978.940.662.470.853.170.964.790.99759.120.672.50.863.230.96254.890.682.520.873.30.9655.02Figure 4
The probability distribution ofZ̃n* is computed with 106 Monte Carlo simulations. The interval [t0,t0′] takes value of [0.01/n,0.05], and interval [λ0,λ0′] takes value of [0.80,0.95]. The curves from the top to the bottom correspond to n=1,10,102,103,104, and 105.Corollary 4.1.
Definem̃0,1-α*(λ) as
(4.2){12(1-λ)(z̃¯1-α*(m)λ(1-λ)+(z̃¯1-α*(m))2λ(1-λ)+4(m-Rλ)(1-λ))}2,
where z̃¯1-α*(m)=maxn=1mz̃1-α*(n). Let
(4.3)m̃0,1-α*=min⁡(min⁡λ∈[λ0,λ0′]m̃0,1-α*(λ),min⁡λ∈[t0,t0′]m̃0,1-α*(λ)). Define Q̃*(4.4)Q̃1-α*(m0,t)=m0t+z̃1-α*(m0)m0t(1-t)Rt.
Replacing m0 by m̃0,1-α* results in Q̃1-α*(m̃0,1-α*,t). We have that m̃0,1-α* is a 1-α UCB for m0, and Q̃* is a 1-α upper prediction band for Q for t∈[t0·t0′], that is,
(4.5)pr(m0≤m̃0,1-α*,Qt≤Q̃1-α*(m̃0,1-α*,t)∀t∈[t0,t0′])≥1-α.Note that the1-α UCB for m0 takes not only the minimum of m̃0,1-α*(λ) for λ∈[λ0,λ0′], but also the minimum of m̃0,1-α*(t) for t∈[t0,t0′]. This advantage is due to the construction of the Z̃n*, which takes maximum over these two intervals. The details of the calculation are summarized in Algorithm 2. The proof of this corollary is the same as that in Theorem 3.7.Algorithm 2:Focused simultaneous inferences on the UPBs for the false discovery proportion and the UCB for m0.
For any givenα and γ, choose t0=0.01/m, t0′=0.05, λ0=0.8,λ0′=0.95.(1) Compute m̃0,1-α*(λ) of (4.2) for some pre-specified λi’s in the region (λ0,λ0′)and pre-specifiedti’s in the region (t0,t0′),sayλi=λ0+(λ0′-λ0)i/1000, ti=t0+(t0′-t0)i/1000 for i=0,…,1000.(2)Compute m̃0,1-α*=min(minim̃0,1-α*(λi),minim̃0,1-α*(ti)).Thism̃0,1-α* is the 1-α UCB for m0. If m̃0,1-α* exceeds m, replace it by m.(3)Sort the observed P-values such that P(1)≤⋯≤P(m),and use (4.4) to compute the 1-α UPB for the false discovery proportion Q,that is, forP(i)∈[t0,t0′]Q̃1-α*(P(i))=(1/i)(m̃0,1-α*P(i)+z̃1-α*(m̃0,1-α*)m̃0,1-α*P(i)(1-P(i))).IfQ̃1-α*(P(i)) exceeds 1, replace it by 1.(4)Compute τ=max{P(i)∈[t0,t0′]:Q̃1-α*(P(i))≤γ},reject the hypotheses whoseP-values are no greater than τ,which ensures that the false discovery proportionQ is not exceeding γ withprobability1-α.By settingλ0=λ0′=t and t0=t0′=t, this corollary is equivalent to Theorem 2.3 through some algebra manipulations, while Theorem 2.3 uses the exact confidence bound from the binomial distribution without relying on the quantiles of Z̃n*. Furthermore, Theorem 3.7 is exact a special case of Corollary 4.1 by setting λ0=0,λ0′=1, and t0=0,t0′=1 and by considering open intervals rather close intervals. The focused inference thus unifies both the fixed rejection approach and simultaneous approach. We should be cautious of selecting [t0,t0′] and [λ0,λ0′] based on the observed data, which may result in overoptimistic false discovery proportions. These settings have to be decided before the data are generated. A careful study of choosing appropriate values for [t0,t0′] and [λ0,λ0′] is open for future research.


## 5. Generalizing the Results to Less-Independent Situations

The results of Sections2, 3, and 4 are based on the assumption that the true null P-values are independently distributed as U[0,1]. Given this, we need no further assumptions concerning the false null P-values. This independence assumption can be weakened as in the following:Binomial Dominant Condition: One has Vt≤dNm0,tfor0<t<1.The notationX≤dY means that random variable X is stochastically no greater than random variable Y, that is, pr(X≤x)≥pr(Y≤x) for any x. Replacing the independence assumption by the binomial dominant condition, the results corresponding to Lemma 2.2, Theorem 2.3, and Corollary 2.7 in Section 2 still hold for a fixed rejection region. For the simultaneous UPBs in Sections 3 and 4, we need a stronger assumption than the binomial dominant condition as the joint distribution of {Vt,t∈(0,1)} needs to be specified. We can replace the binomial dominant condition by the following.Joint Binomial Dominant Condition: (Vt1,…,Vtk)≤d(Nm0,t1,…,Nm0,tk) for any k=1,2,…, and t1,…,tk∈(0,1). Here Nm0,t=∑i=1m0I(Ui≤t), and Ui,i=1,…,m0 are mutually independently distributed as distribution U[0,1]. The notation (X1,…,Xk)≤d(Y1,…,Yk) means for any x1,…,xk, pr(X1≤x1,…,Xk≤xk)≥pr(Y1≤x1,…,Yk≤xk).Replacing the independence assumption by this joint binomial dominant condition, the results in Sections3 and 4 are still valid for the upper prediction band for Q and the UCB for m0. A special case for this joint binomial dominant condition is that when the true null P-values are independent with distribution stochastically no smaller than U[0,1]. This happens when the null hypothesis is composite or the statistic to test the null hypothesis is not a continuous random variable.More generally, we would like the construction of upper prediction band forQ not to rely on the independence assumption or any kind of weak dependence assumption (the binomial dominant condition or the joint binomial dominant condition). The method of Romano and Shaikh [7, 8] can be applied without any assumptions on the dependence, but may potentially have lost power due to that the correlation structure of the data has not been exploited. A resampling procedure [13, 19] has been proposed to address this limitation.


## 6. Discussion

The method of this paper applies to data where true nullP-values are independent, or to slightly dependent data where the joint binomial dominant condition is satisfied. This assumption does not rely on any specification for the false null P-values. In this paper we used the idea of considering a fixed rejection region to construct a UPB for Qt and a UCB for m0. By utilizing the normalized empirical process Z̃n=supt∈(0,1)(Nn,t-nt)/nt(1-t), we find simultaneous UPBs for Qt for all t∈(0,1) and can further modify the construction of the UCB for m0. The result of Theorem 3.7 gives the joint statement about the UCB for m0 and the simultaneous UPBs for the false discovery proportions Q. A focused approach in Corollary 4.1 unifies the result of the fixed rejection region method and the simultaneous approach.The method in this paper is based on finite samples and avoids the slow convergence problem of the asymptotic theory for the empirical processZ̃n. The Monte Carlo simulations give very accurate estimates of the quantiles for Z̃n. The standard error associated with the Monte Carlo simulations in computing the probability in (3.7) is no greater than α(1-α)/B, where B is the number of simulations.In the dataset where the test statistics are not independent or do not satisfy joint binomial dominant condition, the method in this paper may not be guaranteed to work. One can alternatively use the methods proposed in Romano and Shaikh [7, 8], Meinshausen [13], Ge et al. [19]. The method proposed in this paper can be potentially extended to dependent data by using resamplings, and this work is open for future research.



---

*Source: 320425-2012-06-07.xml*