What's the formula for the Benjamini-Hochberg adjusted p-value?

Question

I understand the procedure and what it controls. So what's the formula for the adjusted p-value in the BH procedure for multiple comparisons?

---

Just now I realized the original BH didn't produce adjusted p-values, only adjusted the (non) rejection condition: https://www.jstor.org/stable/2346101. Gordon Smyth introduced adjusted BH p-values in 2002 anyways, so the question still applies. It's implemented in R as p.adjust with method BH.

Answer 1

The famous seminal Benjamini & Hochberg (1995) paper described the procedure for accepting/rejecting hypotheses based on adjusting the alpha levels. This procedure has a straightforward equivalent reformulation in terms of adjusted $p$-values, but it was not discussed in the original paper. According to Gordon Smyth, he introduced adjusted $p$-values in 2002 when implementing p.adjust in R. Unfortunately, there is no corresponding citation, so it has always been unclear to me what one should cite if one uses BH-adjusted $p$-values.

Turns out, the procedure is described in the Benjamini, Heller, Yekutieli (2009):

An alternative way of presenting the results of this procedure is by presenting the adjusted $p$-values. The BH-adjusted $p$-values are defined as $$p^\mathrm{BH}_{(i)} = \min\Big\{\min_{j\ge i}\big\{\frac{mp_{(j)}}{j}\big\},1\Big\}.$$

This formula looks more complicated than it really is. It says:

1. First, order all $p$-values from small to large. Then multiply each $p$-value by the total number of tests $m$ and divide by its rank order.
2. Second, make sure that the resulting sequence is non-decreasing: if it ever starts decreasing, make the preceding $p$-value equal to the subsequent (repeatedly, until the whole sequence becomes non-decreasing).
3. If any $p$-value ends up larger than 1, make it equal to 1.

This is a straightforward reformulation of the original BH procedure from 1995. There might exist an earlier paper that explicitly introduced the concept of BH-adjusted $p$-values, but I am not aware of any.

---

Update. @Zenit found that Yekutieli & Benjamini (1999) described the same thing already back in 1999:

Answer 2

First a to the point answer. Consider that $p_0$ is the (single test) $p$ value associated with value $z_0$ of the test statistic. The Benjamini-Hochberg FDR is computed in two steps ($N_0$ = # pvalues $\le$ $p_0$, $N$ = # pvalues):

• $\text{FDR }(p_0) = \frac{\quad p_0 \quad }{\frac{N_0}{N}}$

• $\text{FDR }(p_i) = \min (\text{FDR}(p_i), \text{FDR}(p_{i+1}))$

---

Now let's understand this. The (Bayesian) underlying idea is that observations come from a mixture of two distributions:

• $\pi_0 \: N$ observations from the null density $f_0(z)$
• $(1-\pi_0) \: N$ observations from alternative density $f_1(z)$.

What is observed is the mixture of those two:

• $f(z) = \pi_0 \cdot f_0(z) + (1-\pi_0) \cdot f_1(z)$

The (Bayesian) definitions are:

• $\text{Fdr} = \frac{\pi_0 \: (1-F_0(z_0))}{(1-F(z))}$ (a fraction of the tail areas)
• $\text{fdr} = \frac{\pi_0 \: f_0(z_0)}{f(z)}$ (a fraction of the tail densities)

As shown below, Fdr is equivalent to the Benjamini hocherg FDR when $\pi_0 \approx 1$ (which is the case in most bioinformatics studies)

(Based on Efron & Tibshirani's Computer Age Statistical Inference)