Benjamini-Hochberg Procedure: A Comprehensive Technical Analysis of FDR Control

1. Introduction

1.1 The Multiple Testing Problem

A test at the 0.05 level has a 5% chance of a false positive when there is no effect. Run 20 independent tests of true nulls and the chance of at least one false positive is about 64%; run 1,000 and you should expect around 50 false positives. Testing many metrics, segments, products or variants therefore turns routine significance into a source of false leads.

1.2 What This Paper Covers

• Two error rates: the family-wise error rate (FWER) and the false discovery rate (FDR)
• The BH algorithm, with a small worked example you can run
• The conditions under which its guarantee holds
• A reproducible simulation comparing no correction, Bonferroni and BH
• Guidance on when to use which

2. Background

2.1 Family-Wise Error Rate and Bonferroni

The family-wise error rate is the probability of making at least one false rejection among all tests. The Bonferroni correction tests each hypothesis at level alpha / m, where m is the number of tests, and guarantees FWER of at most alpha whatever the dependence between tests. Holm's step-down procedure keeps the same guarantee and is never less powerful. The price is that the per-test threshold shrinks as m grows, so moderate real effects are often missed.

2.2 False Discovery Rate

Benjamini and Hochberg (1995) proposed controlling a different quantity. If R is the number of rejections and V the number of false rejections, the false discovery proportion is V / R (taken as 0 when R = 0), and the FDR is its expected value. Controlling FDR at 0.10 means that, on average over repeated studies, no more than about 10% of the results you flag are false. That is a weaker guarantee than FWER control, and it allows more discoveries.

FDR is an average over repetitions. It does not say that exactly 10 of any particular 100 flagged results are false, and it does not replace follow-up of the results that matter.

3. Approach

3.1 The Benjamini-Hochberg Algorithm

The threshold rises with rank, from q / m for the smallest p-value to q for the largest. A small p-value can therefore be rejected even when it is above q / m, provided enough other p-values are small too. A small example with five p-values:

import numpy as np
from statsmodels.stats.multitest import multipletests

p = np.array([0.001, 0.008, 0.039, 0.041, 0.600])
m, q = len(p), 0.05
thresholds = np.arange(1, m + 1) / m * q          # (i / m) * q for the sorted p-values
below = p <= thresholds
k = np.max(np.nonzero(below)[0]) + 1 if below.any() else 0
print("thresholds:", np.round(thresholds, 3))
print("rejected the smallest", k, "p-values")
print("statsmodels fdr_bh:", multipletests(p, alpha=q, method="fdr_bh")[0])

thresholds: [0.01 0.02 0.03 0.04 0.05]
rejected the smallest 2 p-values
statsmodels fdr_bh: [ True  True False False False]

The third p-value, 0.039, is above its threshold of 0.03 and the fourth, 0.041, is above 0.04, so only the two smallest are rejected.

3.2 When the Guarantee Holds

Benjamini and Hochberg proved that for independent test statistics the procedure gives FDR at most (m₀ / m) * q, where m₀ is the number of true null hypotheses. Benjamini and Yekutieli (2001) showed that the same bound holds under a form of positive dependence (positive regression dependence on a subset), and gave a more conservative variant, which divides q by the sum 1 + 1/2 + ... + 1/m, that is valid under any dependence. The guarantee also assumes the p-values themselves are valid: tests that violate their assumptions give p-values that BH cannot repair.

3.3 Choosing an Error Rate

The number of tests does not decide the error rate; the cost of mistakes does. If each flagged result triggers inexpensive follow-up, a 10% false share may be fine. If one false positive would be acted on directly, such as a safety or regulatory decision, a family-wise method is the safer choice. Writing down the cost of a false positive and of a miss before testing makes the choice explicit.

4. A Reproducible Simulation

The simulation runs 1,000 two-sided z-tests, 100 of which have a real shift of 3 standard errors; the other 900 are true nulls. It compares no correction at 0.05, Bonferroni at 0.05 and BH at q = 0.10, averaged over 500 repetitions. In the second setting all test statistics share a correlation of 0.5, a positive dependence.

import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

def simulate(rho, reps=500, m=1000, n_true=100, effect=3.0, seed=11):
    """m two-sided z-tests; n_true of them have a real shift of `effect` standard errors.
    rho is the correlation shared by all test statistics (0 = independent)."""
    rng = np.random.default_rng(seed)
    truth = np.zeros(m, dtype=bool)
    truth[:n_true] = True
    out = {"uncorrected": [], "Bonferroni": [], "BH": []}
    for _ in range(reps):
        shared = rng.normal()
        z = np.sqrt(rho) * shared + np.sqrt(1 - rho) * rng.normal(size=m) + effect * truth
        p = 2 * stats.norm.sf(np.abs(z))
        calls = {"uncorrected": p < 0.05,
                 "Bonferroni": multipletests(p, alpha=0.05, method="bonferroni")[0],
                 "BH": multipletests(p, alpha=0.10, method="fdr_bh")[0]}
        for name, hit in calls.items():
            false_disc = (hit & ~truth).sum()
            out[name].append((false_disc / max(hit.sum(), 1), (hit & truth).sum() / n_true, false_disc))
    return {k: np.array(v).mean(axis=0) for k, v in out.items()}

for rho in (0.0, 0.5):
    print(f"correlation between test statistics = {rho}  (1000 tests, 100 real effects, 500 repetitions)")
    print(f"  {'method':12s} {'false share of discoveries':>27s} {'real effects found':>19s} {'false discoveries':>18s}")
    for name, (fdp, power, fd) in simulate(rho).items():
        print(f"  {name:12s} {fdp:27.3f} {power:19.3f} {fd:18.1f}")

Output:

correlation between test statistics = 0.0  (1000 tests, 100 real effects, 500 repetitions)
  method        false share of discoveries  real effects found  false discoveries
  uncorrected                        0.343               0.852               44.7
  Bonferroni                         0.003               0.145                0.0
  BH                                 0.090               0.614                6.2
correlation between test statistics = 0.5  (1000 tests, 100 real effects, 500 repetitions)
  method        false share of discoveries  real effects found  false discoveries
  uncorrected                        0.241               0.852               39.5
  Bonferroni                         0.004               0.134                0.0
  BH                                 0.050               0.563                6.9

With uncorrected tests about a third of the discoveries (independent case) were false. Bonferroni made almost no false discoveries and found 13% to 15% of the real effects. BH found about 56% to 61% of the real effects; its false share was 0.090 with independent tests, close to q times the fraction of true nulls (0.10 times 0.9), and 0.050 with the shared correlation (with other random seeds this second figure came out near 0.05 to 0.06). These figures describe this design, with a fairly strong effect and 10% real effects; with weaker or rarer effects every method finds fewer.

5. Practical Guidance

6. Conclusion

BH trades a controlled share of false positives for much better power than Bonferroni, which suits exploratory screening where results are followed up. In the simulation here it found about 60% of the real effects at a false share near its target, where Bonferroni found under 15%. The procedure is only as good as the p-values it receives and the conditions behind its guarantee, so state the error rate, the level and the family of tests.