Multiple Testing and Error Rate Control — A Guided Tour of Modern Statistics

Notes / Hypothesis Testing

3.2 Multiple Testing and Error Rate Control

3.2.1 Bonferroni’s Test and Fisher’s Test

We focus on global testing first with a number of null hypotheses \(H_{0,i}\). We may care about the global null hypothesis \[ H_0 = \bigcap _{i=1}^n H_{0,i}. \] We’ll assume that our p-values have the super-uniform property that (here \(p_i\) is a random variable coming out of the test) \[ {\mathbb{P}} (p_i\le t) = t \] under the null hypothesis.

Proof. The probability of rejecting the null under the null hypothesis is, by a union bound and the super-uniform property \begin{align*} \mathbb {P}(\text {reject}) &= \mathbb {P}\left (\bigcup _{i=1}^{n} \left \{p_i \leq \frac {\alpha }{n}\right \}\right ) \\ &\leq \sum _{i=1}^{n} \mathbb {P}\left (p_i \leq \frac {\alpha }{n}\right ) \\ &\leq n \cdot \frac {\alpha }{n} \\ &= \alpha . \end{align*} □

Now we need to discuss which test is better. Consider a particular “needle in a haystack" problem. Let \(Y_i \sim N(\mu _i, 1)\) independently for \(i = 1, \ldots , n\), and consider the global null \[ H_0 : \mu _1 = \cdots = \mu _n = 0 \quad \text {vs.} \quad H_1 : \exists \, i,\ \mu _i \neq 0. \] Under \(H_0\), Bonferroni’s statistic \(\max _i Y_i\) concentrates around \(\sqrt {2 \log n}\). So we have power only when some \(\mu _i > \sqrt {2 \log n}\) — the “needle in a haystack” threshold. Plotting power against \(\mu _i / \sqrt {2 \log n}\), there is a sharp phase transition at \(1\): power is \(\approx \alpha \) below and \(\approx 1\) above.

Assume the size of “needle" \(h = \sqrt {2r\log n}\). We will prove no \(\alpha \)-level test with some nontrivial power for \(r=1-\epsilon \). The point is to avoid having a composite alternative: instead, consider the Bayesian decision problem with null hypothesis \[ H_0 : \mu _i = 0 \text { for all } i \] and a “simple” alternative hypothesis \[ H_1 : \{\mu _i\} \sim \pi , \] where \(\pi \) selects a coordinate uniformly at random and sets its mean to \(\mu ^{(n)}\), keeping all other means to zero.

The most powerful test in such a setting where \(H_0, H_1\) are both simple hypotheses is the likelihood ratio test, and if we can show it has no power then everything else will have no power. Indeed, we have likelihoods \[ f_0(y) = \prod _{j=1}^{n} \frac {1}{\sqrt {2\pi }} \exp \left (-\frac {1}{2} y_j^2\right ), \quad f_1(y) = \frac {1}{n} \sum _{i=1}^{n} \frac {1}{\sqrt {2\pi }} \exp \left (-\frac {1}{2}(y_i - \mu )^2\right ) \prod _{j : j \neq i} \frac {1}{\sqrt {2\pi }} \exp \left (-\frac {1}{2} y_j^2\right ). \] and we want to reject when \(\frac {f_1}{f_0}\) is large. What’s nice is that all the terms cancel out except the shifted mean: \[ L = \frac {1}{n} \sum _{i=1}^{n} \exp \left (Y_i \mu - \frac {1}{2} \mu ^2\right ). \] Notice that this is different from Bonferroni — it’s a softmax instead of looking at the maximum, since if \(\mu = \infty \) this would be completely dominated by the maximum \(Y_i\). This is nice because it’s just a sum of iid terms of mean \(1\), and in fact if \(\mu \) is small enough this likelihood concentrates:

Under such proposition, if we set threshold \(T_n(\alpha )\) so that \({\mathbb{P}} _{{\mathbb{H}} _0}(L\ge T_n(\alpha )) = \alpha \), then \begin{align*} \mathbb {P}(\text {type II error}) = \mathbb {P}_{H_1}(L \leq T_n(\alpha )) &= \int \mathbf {1}\{L \leq T_n(\alpha )\} \, dP_{H_1} \\ &= \int L \cdot \mathbf {1}\{L \leq T_n(\alpha )\} \, dP_{H_0} \\ &= \int \mathbf {1}\{L \leq T_n(\alpha )\} \, dP_{H_0} + \int (L - 1) \mathbf {1}\{L \leq T_n(\alpha )\} \, dP_{H_0} \\ &\to (1-\alpha ) + 0, \end{align*}

\[ \implies \mathbb {P}(\text {type I error})+\mathbb {P}(\text {type II error}) \to 1. \]

Our original model was that we had independent statistics \(X_i\) which are \(N(0, 1)\) under the null hypothesis and \(N(\mu _i, 1)\) under the alternative, but now we want to extend to a setting where we have a small fraction of non-null hypotheses. Thus we will now use a simple model where \[ H_0 : X_i \stackrel {\text {iid}}{\sim } N(0, 1), \quad H_1 : X_i \stackrel {\text {iid}}{\sim } (1 - \varepsilon ) N(0, 1) + \varepsilon N(\mu , 1). \] In the literature we historically parameterize \[ \varepsilon _n = n^{-\beta }, \quad \frac {1}{2} < \beta < 1, \] \[ \mu _n = \sqrt {2r \log n}, \quad 0 < r < 1. \] It turns out that we have a threshold curve \[ \rho ^*(\beta ) = \begin {cases} \beta - \frac {1}{2} & \frac {1}{2} < \beta \leq \frac {3}{4}, \\ (1 - \sqrt {1 - \beta })^2 & \frac {3}{4} \leq \beta \leq 1, \end {cases} \] such that Neyman-Pearson has full power for \(r > \rho ^*(\beta )\) (that is, we can adjust the test so that the sum of type I and type II error probabilities approaches \(0\)) and no power for \(r < \rho ^*(\beta )\) (that is, for any test, the limiting sum of type I and type II error is at least \(1\)). A crude calculation, we get power if \[ \max _{\text {non-null}} X_i \approx \sqrt {2r \log n} + \sqrt {2 \log n^{1-\beta }} > \sqrt {2 \log n}. \]

3.2.2 Higher Criticism

In general, the whole point is that in global testing we do not know \(\epsilon \) and \(\mu \) and thus cannot use the NP test in the first place. We introduce Tukey’s higher criticism. \[ HC_n^* = \max _{0 < \alpha \leq \alpha _0} \frac {F_n(\alpha ) - \alpha }{\sqrt {\alpha (1 - \alpha )/n}}. \]

Remark 3.26. The main problem is that near \(t = 0\), we are no longer approximately normal — \(B(n, p)\) is nicely approximated by a Gaussian if \(np\) is not too small, but if \(np\) is small we get something that’s Poisson, which has much heavier tails. The Burk-Jones statistic is an attempt to resolve this problem, but it is not fully effective. The idea is that for each \(t\), we can test whether \(n\hat {F}_n(t) < t\) using a likelihood ratio test \[ \log LR_n(t) = \begin {cases} nD(\hat {F}_n(t), t) & 0 \leq t \leq F_n(t), \\ 0 & \text {otherwise}, \end {cases} \] where \(D(p_0, p_1) = p_0 \log \frac {p_0}{p_1} + q_0 \log \frac {q_0}{q_1}\), \(q_i = 1 - p_i\) is the Kullback-Leibler divergence. We then define \[ BJ^+ = \max _{1 \leq i \leq n/2} nD\left (p_{(i)}, \frac {i}{n}\right ) \] that is, at what significance level we detect a divergence between what we see and what we expect.

3.2.3 False Discovery Rate

We now turn to identifying which hypotheses are non-null and to controlling the false discovery rate.

	accepted	rejected	total

true	\(U\)	\(V\)	\(n_0\)
false	\(T\)	\(S\)	\(n - n_0\)

total	\(n - R\)	\(R\)	\(n\)

Proof. Since \(V\) is nonnegative-integer-valued, we have \(\mathbb {P}(V \geq 1) \leq \mathbb {E}[V]\). And letting \(\mathcal {N}_0\) denote the set of null hypotheses, we have \[ \mathbb {E}[V] = \mathbb {E}\left [\sum _{i \in \mathcal {N}_0} \mathbf {1}\left \{p_i \leq \frac {\alpha }{n}\right \}\right ] = \sum _{i \in \mathcal {N}_0} \mathbb {P}\left (p_i \leq \frac {\alpha }{n}\right ) = n_0 \cdot \frac {\alpha }{n}, \] which completes the proof. □

Suppose we have \(n\) hypotheses \(H_{(1)}, \ldots , H_{(n)}\) corresponding to the ordered \(p\)-values \(p_{(1)} \leq \cdots \leq p_{(n)}\) (so we choose \(H_{(1)}\) to be the hypotheses with the most surprising \(p\)-value, and so on). Now we will compare \(p\)-values with an adaptive threshold based on what we’ve seen so far.

What we do here is called a Holm’s procedure procedure:

First, compare with Bonferroni’s threshold: if \(p_{(1)} \leq \frac {\alpha }{n}\), then reject \(H_{(1)}\) and move to the next step. Otherwise, reject nothing (accept all null hypotheses).
Now in general for step \(i\), if \(p_{(i)} \leq \frac {\alpha }{n - i + 1}\), then reject \(H_{(i)}\) and go to the next step (\(i + 1\)). Otherwise, accept \(H_{(i)}, H_{(i+1)}, \ldots , H_{(n)}\) and stop.
Finally, if \(p_{(n)} \leq \alpha \), then we reject \(H_{(n)}\); otherwise we accept it.

Notice that \[ \{V \geq 1\} = \{\text {procedure reached } i_0\} \cap \left \{p_{(i_0)} \leq \frac {\alpha }{n - i_0 + 1}\right \}, \] and \(\mathbb {P}(V \geq 1) \leq \alpha \).

However, FWER is so stringent that we often return nothing if we require FWER control.

We’ll now consider a procedure generally more powerful than Holm’s procedure – instead of Hochberg’s procedure, we consider a step-up procedure with critical values \(\alpha _i = \frac {\alpha i}{n}\), which is far less conservative than \(\frac {\alpha }{n - i + 1}\).

Proof. Let \(V_i = \{H_i \text { rejected}\}\) be the indicator function for hypothesis \(i\) being rejected. By definition, we have \[ \mathrm {FDP} = \sum _{i \in \mathcal {H}_0} \frac {V_i}{R \vee 1}. \] Now, it suffices to show that for any null \(i\), we have \(\mathbb {E}\left [\frac {V_i}{R \vee 1}\right ] = \frac {\alpha }{n}\). (This is somehow “the only answer we can get” because the nulls are uniform and thus the random variables are exchangeable.) To prove this claim, notice that we can do casework over the value of \(R\) and write \[ \frac {V_i}{R \vee 1} = \sum _{k=1}^{n} \frac {V_i \mathbf {1}\{R = k\}}{k} = \sum _{k=1}^{n} \frac {\mathbf {1}\left \{p_i \leq \frac {\alpha k}{n}\right \} \mathbf {1}\{R = k\}}{k}, \] since assuming \(R = k\), we know the threshold for rejection is \(\frac {\alpha k}{n}\). Notice that on the event \(p_i \leq \frac {\alpha k}{n}\), changing \(p_i\) to zero doesn’t change the threshold, meaning whenever we reject \(H_i\), the number (and identity) of rejections is the same. So we can write the above expression as \[ \sum _{k=1}^{n} \frac {\mathbf {1}\left \{p_i \leq \frac {\alpha k}{n}\right \} \mathbf {1}\{R(p_i \to 0) = k\}}{k}, \] where this notation means that we set this null \(p\)-value to zero. Now we can take the expectation of this quantity conditioned on all other \(p\)-values – the only randomness is in \(p_i\) here, so \begin{align*} \mathbb {E}\left [\frac {V_i}{R \vee 1} \,\bigg |\, p_1, \ldots , p_{i-1}, p_{i+1}, \ldots , p_n\right ] &= \sum _{k=1}^{n} \frac {\frac {\alpha k}{n} \mathbf {1}\{R(p_i \to 0) = k\}}{k} \\ &= \sum _{k=1}^{n} \frac {\alpha }{n} \mathbf {1}\{R(p_i \to 0) = k\} \\ &= \frac {\alpha }{n}. \end{align*}

Then finally taking the expectation over the last \(p\)-value yields the result. □

Consider dependency,

The following proof is also by Professor Candés and a former student:

Proof. Let \(\alpha _i = \frac {i\alpha }{n}\); much like in the proof before, it suffices to show that \(\mathbb {E}\left [\frac {V_i}{R \vee 1}\right ] = \frac {\alpha }{n} S(n)\). We again have \[ \frac {V_i}{R \vee 1} = \sum _{k=1}^{n} \frac {\mathbf {1}\{p_i \leq \alpha _k\} \mathbf {1}\{R = k\}}{k}, \] and we look at where \(p_i\) can fall. Summing over the possible ranks it can take on, we have \[ \sum _{k=1}^{n} \sum _{\ell =1}^{k} \frac {\mathbf {1}\{\alpha _{\ell -1} \leq p_i \leq \alpha _\ell \} \mathbf {1}\{R = k\}}{k} = \sum _{\ell =1}^{k} \sum _{k \geq \ell } \frac {\mathbf {1}\{\alpha _{\ell -1} \leq p_i \leq \alpha _\ell \} \mathbf {1}\{R = k\}}{k} \] just by swapping the order of summation. But now if we do the \(k\)-sum first, we’re just looking at the probability of getting a particularly high number of rejections, so this simplifies to \[ \sum _{\ell =1}^{n} \frac {\mathbf {1}\{R \geq \ell \}}{R} \mathbf {1}\{p_i \in [\alpha _{\ell -1}, \alpha _\ell ]\}. \] Everything so far has been an equality, so “nothing interesting” has happened yet. But now we can simplify the first fraction to be bounded by \(\frac {1}{\ell }\), \[ \sum _{\ell =1}^{n} \frac {1}{\ell } \mathbf {1}\{p_i \in [\alpha _{\ell -1}, \alpha _\ell ]\} = \sum _{\ell =1}^{n} \frac {1}{\ell } \frac {\alpha }{n} = S(n) \frac {\alpha }{n}, \] and then the rest of the proof proceeds as before. What’s surprising is that the result of Guo and Rao shows that there are distributions of \(p\)-values for which this inequality is indeed tight! □

3.2.4 E-Values

Another bad news is in many cases we can not get p-values easily. The e-value is a generic replacement of the p-value which will handle this problem of optional continuation.

Definition 3.34. Recall that a null hypothesis \(\mathcal {H}_0\) is basically a collection of probability measures. An \(e\)-variable \(E\) for testing \(\mathcal {H}_0\) is a nonnegative random variable such that \[\sup _{P_0 \in \mathcal {H}_0} \mathbb {E}_{P_0}[E] \leq 1.\] A realization of an \(e\)-variable is called an \(e\)-value. Meanwhile, a \(p\)-variable for testing \(\mathcal {H}_0\) is a nonnegative random variable that satisfies \[\sup _{P_0 \in \mathcal {H}_0} \mathbb {P}_{P_0}(P \leq \alpha ) \leq \alpha \] for all \(\alpha \in (0,1)\), and a realized value of a \(p\)-variable is a \(p\)-value.

Proof. Indeed, we have \[\mathbb {P}\left (\frac {1}{E} \leq \alpha \right ) = \mathbb {P}\left (E \geq \frac {1}{\alpha }\right ) \leq \frac {\mathbb {E}[E]}{1/\alpha } \leq \alpha .\] □

We’ll establish an \(e\)-value \(E_1\) on the first batch. From there, we evaluate an \(e\)-value \(E_2\) on the next batch, but only if the outcome is in a certain range (promising but not conclusive) and the external factors take on certain values (things that we cannot plan) – otherwise we stop early. Then depending on the outcomes and external factors up until the second batch, we decide whether or not to compute \(E_3\), and so on. But the point is that after \(\tau \) total data batches, the final result we report is the product \[ V_\tau = \prod _{i=1}^{\tau } E_i. \] In particular, we’re allowed to choose whether to continue on depending on whether each individual \(E_i\) is above some threshold of our choice, and the point is that we’ll still be able to control the type I error:

Theorem 3.37. Regardless of the stop-continue rule, as long as \(\tau \) is a stopping time, \(V_\tau \) is itself an \(e\)-value. More formally, let \(\mathcal {F}_t\) be a filtration. Suppose that for all \(t\), the conditional \(e\)-variable \(E_t\) is a nonnegative random variable which is \(\mathcal {F}_t\)-measurable, and such that for all \(P_0 \in \mathcal {H}_0\) we have \(\mathbb {E}_{P_0}[E_t | \mathcal {F}_{t-1}] = 1\). (This is easy to check in practice.) Then \(V_t = \prod _{i \leq t} E_i\) is a nonnegative supermartingale under the null. Thus by the optional stopping theorem, for any stopping time \(\tau \) with respect to the filtration, \(V_\tau = \prod _{t=1}^{\tau } E_t\) is an \(e\)-value. (In particular, \(V_t\) is an \(e\)-value for each fixed \(t\).)

Now we focus on multiple testing.