19
$\begingroup$

The Benjamini-Hochberg procedure is a method that corrects for multiple comparisons and has a false discovery rate (FDR) equal to $\alpha$.

Or is it the family wise error rate, FWER? I am a bit confused about this. According to my below computations it seems to be the FWER that equals $\alpha$ and not the FDR.

Can we prove that this is true?

Let's assume that the multiple p-values for the different hypotheses are independent and the distribution of the p-value's (conditional on the null hypotheses being true) is uniform between $0,1$.


I can use a simulation to show that it get's close. With the below numbers $\alpha = 0.1$, and the number of times that I reject a hypothesis in this simulation is

$$\begin{array}{rcl} \alpha& =& 0.1\\ \text{observed FDR} &=& 0.100002 \pm 0.00030 \end{array}$$

with error based on $ \pm 2\sigma$ where $\sigma = \sqrt{\frac{0.1 \cdot 0.9}{ n}}$

set.seed(1)
m <- 10^6
n <- 10
a <- 0.1
k <- 1:n

sample <- function( plotting = F) {
  p <- runif(n)
  p <- p[order(p)]
  counts <- max(0,which(p<k/n*a))
  if (plotting) {
    plot(k,p, ylim = c(0,1) )
    lines(k,k/n*a)
  }
  counts
}

x <- replicate(m, sample())
s <- sum(x>0)/m
err_s <- sqrt(s*(1-s)/m)
c(s-2*err_s,s,s+2*err_s) 
$\endgroup$
1
  • $\begingroup$ I understand now that the FWER is not just the case of type I error when all the hypotheses are correct, but also when one or more of them are correct and we do not want a type I error for the non-correct hypotheses. $\endgroup$ Commented Nov 5, 2020 at 17:30

3 Answers 3

19
+100
$\begingroup$

Intro

Let us start with some notation: We have $m$ simple hypotheses we test, with each null numbered $H_{0,i}$. The global null hypothesis can be written as an intersection of all the local nulls: $H_0=\bigcap_{i=1}^{m}{H_{0,i}}$. Next, we assume that each hypothesis $H_{0,i}$ has a test statistic $t_i$ for which we can compute the p-value $p_i$. More specifically, we assume that for each $i$ the distribution of $p_i$ under $H_{0,i}$ is $U[0,1]$.

For example (Chapter 3 of Efron), consider a comparison of 6033 genes in two groups: For each gene in $1\le i\le 6033$, we have $H_{0,i}:$ "no difference in the $i^{th}$ gene" and $H_{1,i}:$ "there is a difference in the $i^{th}$ gene". In this problem, $m=6033$. Our $m$ hypotheses can be divided into 4 groups:

$H_{0,i}$ Accepted Rejected Sum
Correct $U$ $V$ $m_0$
Incorrect $T$ $S$ $m-m_0$
Sum $m-R$ $R$ $m$

We do not observe $S,T,U,V$ but we do observe $R$.

FWER

A classic criterion is the familywise error, denoted $FWE=I\{V\ge1\}$. The familywise error rate is then $FWER=E[I\{V\ge1\}]=P(V\ge1)$. A comparison method controls the level-$\alpha$ $FWER$ in the strong sense if $FWER\le\alpha$ for any combination $\tilde{H}\in\left\{H_{0,i},H_{1,i}\right\}_{i=1,...,m}$; It controls the level-$\alpha$ $FWER$ in the weak sense if $FWER\le\alpha$ under the global null.

Example - Holm's method:

  1. Order the p-values $p_{(1)},...,p_{(m)}$ and then respectively the hypotheses $H_{0,(1)},...,H_{0,(m)}$.
  2. One after the other, we check the hypotheses:
    1. If $p_{(1)}\le\frac{\alpha}{m}$ we reject $H_{0,(1)}$ and continue, otherwise we stop.
    2. If $p_{(2)}\le\frac{\alpha}{m-1}$ we reject $H_{0,(2)}$ and continue, otherwise we stop.
    3. We keep rejecting $H_{0,(i)}$ if $p_{(i)}\le\frac{\alpha}{m-i+1}$
  3. We stop the first time we get $p_{(i)}>\frac{\alpha}{m-i+1}$ and then reject $H_{0,(1)},...,H_{0,(i-1)}$.

Holm's method example: we've simulated $m=20$ p-values, then ordered them. Red circles indicate rejected hypotheses, green circles were not rejected, the blue line is the criterion $\frac{\alpha}{m-i+1}$: A visual example of Holm's method

The hypotheses rejected were $H_{0,(1)},H_{0,(2)},H_{0,(3)}$. You can see that $p_{(4)}$ is the smallest p-value larger than the criterion, so we reject all hypotheses with smaller p-values.

It is quite easy to prove that this method controls level-$\alpha$ $FWER$ in the strong sense. A counterexample would be the Simes method, which only controls the level-$\alpha$ $FWER$ in the weak sense.

FDR

The false discovery proportion ($FDP$) is a softer criterion than the $FWE$, and is defined as $FDP=\frac{V}{\max\{1,R\}}=\left\{\begin{matrix} \frac{V}{R} & R\ge1\\ 0 & o.w\end{matrix}\right.$. The false discovery rate is $FDR=E[FDP]$. Controlling level-$\alpha$ $FDR$ means that The false if we repeat the experiment and rejection criteria many times, $FDR=E[FDP]\le\alpha$.

It is very easy to prove that $FWER\ge FDR$: We start by claiming $I\{V\ge1\}\ge\left\{\begin{matrix} \frac{V}{R} & R\ge1\\ 0 & o.w\end{matrix}\right.$.

If $V=0$ then $I\{V\ge1\}=0=\left\{\begin{matrix} \frac{V}{R} & R\ge1\\ 0 & o.w\end{matrix}\right.$. In the above table, $R\ge V$ so for $V>0$ we get $\frac{V}{\max\{1,R\}}\le1$ and $I\{V\ge1\}\ge\frac{V}{\max\{1,R\}}$. This also means $E\left[I\{V\ge1\}\right]\ge E\left[\frac{V}{\max\{1,R\}}\right]$, which is exactly $FWER\ge FDR$.

Another very easy claim if that controlling level-$\alpha$ $FDR$ implies controlling level-$\alpha$ $FWER$ in the weak sense, meaning that under the global $H_0$ we get $FWER=FDR$: Under the global $H_0$ every rejection is a false rejection, meaning $V=R$, so $$FDP=\left\{\begin{matrix} \frac{V}{R} & R\ge1\\ 0 & o.w\end{matrix}\right.=\left\{\begin{matrix} \frac{V}{V} & V\ge1\\ 0 & o.w\end{matrix}\right.=\left\{\begin{matrix} 1 & V\ge1\\ 0 & o.w\end{matrix}\right.=I\{V\ge1\}$$ and then $$FDR=E[FDP]=E[I\{V\ge1\}]=P(V\ge1)=FWER.$$

B-H

The Benjamini-Hochberg method is as follows:

  1. Order the p-values $p_{(1)},...,p_{(m)}$ and then respectively the hypotheses $H_{0,(1)},...,H_{0,(m)}$
  2. Mark as $i_0$ the largest $i$ for which $p_{(i)}\le \frac{i}{m}\alpha$
  3. Reject $H_{0,(1)},...,H_{0,(i_0)}$

Contrary to the previous claims, it is not trivial to show why the BH method keeps $FDR\le\alpha$ (to be more precise, it keeps $FDR=\frac{m_0}{m}\alpha$). It isn't a short proof either, ,his is some advanced statistical courses material (I've seen it in one of my Master of Statistics courses). I do think that for the extent of this question, we can simply assume it controls the FDR.

BH example: again we've simulated $m=20$ p-values, then ordered them. Red circles indicate rejected hypotheses, green circles were not rejected, the blue line is the criterion $\frac{i\cdot\alpha}{m}$:

A visual example of the BH method

The hypotheses rejected were $H_{0,(1)}$ to $H_{0,(10)}$. You can see that $p_{(11)}$ is the largest p-value larger than the criterion, so we reject all hypotheses with smaller p-values - even though some of them ($p_{(6)},p_{(7)}$) are larger than the criterion. Compare this (largest p-value larger than the criterion) and Holm's method (smallest p-value larger than the criterion).

Proving $FDR\le\alpha$ for BH

For $m_0=0$ (which means each $p_i$ is distributed under $H_{1,i}$) we get $FDR\le\alpha$ because $V=0$, so assume $m_0\ge1$. Denote $V_i=I\{H_{0,i}\text{ was rejected}\}$ and $\mathcal{H}_0\subseteq\{1,...,m\}$ the set of hypotheses for which $H_{0,i}$ is correct, so $FDP=\sum_{i\in\mathcal{H}_0}{\frac{V_i}{\max\{1,R\}}}=\sum_{i\in\mathcal{H}_0}{X_i}$. We start by proving that for $i\in\mathcal{H}_0$ we get $E[X_i]=\frac{\alpha}{m}$:

$$X_i=\sum_{k=0}^{m}{\frac{V_i}{\max\{1,R\}}I\{R=k\}}=\sum_{k=1}^{m}{\frac{I\{H_{0,i}\text{ was rejected}\}}{k}I\{R=k\}}=\sum_{k=1}^{m}{\frac{I\left\{ p_i\le\frac{k}{m}\alpha \right\}}{k}I\{R=k\}}$$

Let $R(p_i)$ the number of rejections we get if $p_i=0$ and the rest of the p-values unchanged. Assume $R=k^*$, if $p_i\le\frac{k^*}{m}\alpha$ then $R=R(p_i)=k^*$ so in this case, $I\left\{ p_i\le\frac{k^*}{m}\alpha \right\}\cdot I\{R=k^*\}=I\left\{ p_i\le\frac{k^*}{m}\alpha \right\}\cdot I\{R(p_i)=k^*\}$. If $p_i>\frac{k^*}{m}\alpha$ we get $I\left\{ p_i\le\frac{k^*}{m}\alpha \right\}=0$ and again $I\left\{ p_i\le\frac{k^*}{m}\alpha \right\}\cdot I\{R=k^*\}=I\left\{ p_i\le\frac{k^*}{m}\alpha \right\}\cdot I\{R(p_i)=k^*\}$, so overall we can deduce

$$X_i=\sum_{k=1}^{m}{\frac{I\left\{ p_i\le\frac{k}{m}\alpha \right\}}{k}I\{R=k\}}=\sum_{k=1}^{m}{\frac{I\left\{ p_i\le\frac{k}{m}\alpha \right\}}{k}I\{R(p_i)=k\}},$$

and now we can calculate $E[X_i]$ conditional on all p-values except $p_i$. Under this condition, $I\{R(p_i)=k\}$ is deterministic and we overall get:

$$E\left[I\left\{ p_i\le\frac{k}{m}\alpha \right\}\cdot I\{R(p_i)=k\}\middle|p\backslash p_i\right]=E\left[I\left\{ p_i\le\frac{k}{m}\alpha \right\}\middle|p\backslash p_i\right]\cdot I\{R(p_i)=k\}$$ Because $p_i$ is independent of $p\backslash p_i$ we get

$$E\left[I\left\{ p_i\le\frac{k}{m}\alpha \right\}\middle|p\backslash p_i\right]\cdot I\{R(p_i)=k\}=E\left[I\left\{ p_i\le\frac{k}{m}\alpha \right\}\right]\cdot I\{R(p_i)=k\}\\=P\left(p_i\le\frac{k}{m}\alpha\right)\cdot I\{R(p_i)=k\}$$

We assumed before that under $H_{0,i}$, $p_i\sim U[0,1]$ so the last expression can be written as $\frac{k}{m}\alpha\cdot I\{R(p_i)=k\}$. Next,

$$E[X_i|p\backslash p_i]=\sum_{k=1}^m\frac{E\left[I\left\{ p_i\le\frac{k}{m}\alpha \right\}\cdot I\{R(p_i)=k\}\middle|p\backslash p_i\right]}{k}=\sum_{k=1}^m\frac{\frac{k}{m}\alpha\cdot I\{R(p_i)=k\}}{k}\\=\frac{\alpha}{m}\cdot\sum_{k=1}^m{I\{R(p_i)=k\}}=\frac{\alpha}{m}\cdot 1=\frac{\alpha}{m}.$$

Using the law of total expectation we get $E[X_i]=E[E[X_i|p\backslash p_i]]=E\left[\frac{\alpha}{m}\right]=\frac{\alpha}{m}$. We have previously obtained $FDP=\sum_{i\in\mathcal{H}_0}{X_i}$ so

$$FDR=E[FDP]=E\left[\sum_{i\in\mathcal{H}_0}{X_i}\right]=\sum_{i\in\mathcal{H}_0}{E[X_i]}=\sum_{i\in\mathcal{H}_0}{\frac{\alpha}{m}}=\frac{m_0}{m}\alpha\le\alpha\qquad\blacksquare$$

Summary

We've seen the differences between strong and weak sense of $FWER$ as well as the $DFR$. I think that you can now spot yourself the differences and understand why $FDR\le\alpha$ does not imply that $FWER\le\alpha$ in the strong sense.

$\endgroup$
7
  • $\begingroup$ Nice explanation of FWER (weak and strong) and FDR. I guess that my view in the footnote is about FWER in the weak sense, assuming that all hypotheses are true. $\endgroup$ Commented Nov 24, 2021 at 15:08
  • $\begingroup$ @SextusEmpiricus True, that's the global null. Also, see the visual examples I've added. $\endgroup$
    – Spätzle
    Commented Nov 24, 2021 at 15:37
  • $\begingroup$ They are nice explanations, but the question is about the derivation/proof of the FDR for the BH procedure. $\endgroup$ Commented Nov 24, 2021 at 15:46
  • 1
    $\begingroup$ Thank you @Spätzle that looks like nice work (that I still have to digest) more than worthy of the bounty. I was myselve hoping that I could find an easy solution by the ordered p-values $p_1, p_2, \cdots p_n$ and the conditions $p_k \geq \frac{k}{n} \alpha$ by computing the integral $$\int_{p_1 = \frac{1}{n} \alpha}^1 \int_{p_2 = max(p_1,\frac{2}{n} \alpha)}^1 \dots \int_{p_n = max(p_{n-1},\frac{n}{n} \alpha)}^1 n! dp_n \dots dp_2 dp_1$$ but I still have to write it out as a piecewise polynomial and express it generally. I will see how it compares to your solution. $\endgroup$ Commented Nov 25, 2021 at 8:42
  • 1
    $\begingroup$ We've studied this proof during one of my MA courses and I took it directly from my notes $\endgroup$
    – Spätzle
    Commented Aug 31, 2023 at 5:33
8
$\begingroup$

Caveat: this is just a set of notes for the proof in Spätzle's answer, and nothing more. Please, read his answer for the main notations. Additional notations: $N_0\subset\{1,\dots,m\}$ are the indexes of the true null hypotheses, and $q_i=(i/m)q$.

Suppose that $i\in N_0$. Key facts:

  1. The corresponding (true) null hypothesis is (wrongly) rejected by BH if and only if $p_i\leq q_R$.

  2. The random variables $p_i$ and $R$ are not independent. For instance, suppose that $q<1$ and observe that $$ \mathbb{P}(p_i\leq q\mid R=m)=1\ne P(p_i\leq q)=q. $$

  3. Define $R_i$ as the number of BH rejections among the $p$-values $(\{p_j\}_{j=1}^m\setminus\{p_i\})\cup\{0\}$.

  4. Crucial step: $\{p_i\leq q_r,R=r\}=\{p_i\leq q_r,R_i=r\}$. In words: if a null hypothesis has been rejected by BH, replacing the corresponding $p_i$ by zero doesn't change the total number of BH rejections.

  5. Since BH suppose that all the $p$-values are independent, we have that $p_i$ is independent of $R_i$, by the very definition of $R_i$.

The remaining steps are routine:

\begin{align*} \text{FDR} &= \mathbb{E}\!\left[\sum_{i\in N_0} \frac{\mathbb{I}(p_i\leq q_R)}{\max\{1,R\}}\right] = \sum_{i\in N_0}\sum_{r=1}^m \frac{1}{r} \mathbb{P}(p_i\leq q_r,R=r) \\ &= \sum_{i\in N_0}\sum_{r=1}^m \frac{1}{r} \mathbb{P}(p_i\leq q_r)\mathbb{P}(R_i=r)=\frac{m_0}{m}q\leq q. \end{align*}

$\endgroup$
6
$\begingroup$

Geometric interpretation

The values of the different p-values $p_1,p_2,\dots, p_n$ are distributed in a hypercube and rejection occurs when the point falls inside a region.

The case of 2 variables

For this case we can see easily that the rejection rate is $\alpha$ by adding the area's in the below figure together

plot

Algebraic computation for more variables

We can represent the above area's by the following product where each $x_k$ represents whether for a p-value we have $\alpha \frac{k-1}{n} < p<\alpha \frac{k}{n} $ and the last $x_{n+1}$ represents that the $p>\alpha$

$$(x_1+x_2+ \dots +x_{n+1})^n$$

... to be continued

$\endgroup$
3
  • $\begingroup$ Personally, I like the approach described by Schweder & Spjøtvoll (1982) where you plot (1 - p-value) on the x-axis vs. index number=k on the y-axis. A regression line can then be fit from the origin to the y-axis intercept representing estimated # of true null hypotheses, with estimated rejected H0s straying above the line. Question: is there a more rigorous way to fit the line aside from eyeballing? Or is this the procedure Benjamini & Hochberg suggested? $\endgroup$
    – RobertF
    Commented Apr 27, 2023 at 14:50
  • $\begingroup$ @RobertF it is a long time ago that this question was on my mind, and the other answer helped me a lot to let it rest (I still want to figure out how to finish the problem from my original thought, although that has become a different issue). Your proposed plot of $k$ as function of versus $1-p$ is not directly clear to me how this relates to my answer here. Especially, I am not sure what you mean by regression line. In the plot from spätzli's answer we use a plot of $p$ as function of $k$ and the line $k\alpha/m$ is fixed (there's no eyeballing). $\endgroup$ Commented Apr 27, 2023 at 15:30
  • $\begingroup$ I could be wrong, it may not be a true regression line but an approximate fit. I need to find a copy of Schweder & Spjøtvoll (1982) and reread, has been a while since I looked at this topic as well. When I read explanations of the B-H Method, p-value plots or histograms similar to your own are often presented, but the intuitive connection to the actual procedure isn't explained. I suspect behind the math, the B-H procedure is simply a methodical way of doing the S&S p-value plot and rejecting the null hypothesis for tests that deviate from the line. $\endgroup$
    – RobertF
    Commented Apr 27, 2023 at 16:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.