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.


2 Answers 2


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:

enter image description here

  • $\begingroup$ That's the answer I was expecting, +1. I remember reading about Gordon Smyth implementation of the adjusted p value as well and not knowing who to cite, cool to see there's a "canon" citation to this. $\endgroup$
    – Firebug
    Apr 10, 2019 at 17:52
  • 1
    $\begingroup$ I believe an even earlier reference exists: Yekutieli and Benjamini (1999) (pdf version available here). Definition 2.4 describes how the original 1995 FDR procedure can be rephrased in terms of adjusted p-values. Credit to this blog post where I found about this. $\endgroup$
    – Zenit
    Oct 15, 2019 at 20:46
  • $\begingroup$ @Zenit Oh wow! Great find! I should update my answer. $\endgroup$
    – amoeba
    Oct 15, 2019 at 21:33
  • 4
    $\begingroup$ @Firebug That's because I didn't want to divert citations away from the original BH (1995) paper. I wasn't aware of the p-value formula in Yekutieli & Benjamini (1999). When I talked with Yoav Benjamini about it in 2004, he said that they had decided not to express the procedure in the form of adjusted p-values. I would still tend to cite the BH 1995 paper instead of the 1999 paper. $\endgroup$ Apr 18, 2022 at 6:43
  • 1
    $\begingroup$ Interesting @GordonSmyth ! Thank you for the input, that's a piece of statistical history right there :) $\endgroup$
    – Firebug
    Apr 18, 2022 at 9:19

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)$

enter image description here

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)

enter image description here

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

  • $\begingroup$ But in practice... if I fit an lme4 model or a cox model choosing the best variables with the stepwise method... How do I applythe FDR correction to the p-values of each variable? $\endgroup$
    – skan
    Jun 15, 2023 at 11:05

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.