Mohsen Pirmoradiyan
← Statistical Inference

Benjamini–Hochberg: controlling the false discovery rate when testing thousands of features

Why testing every feature at 0.05 floods you with false positives, how the Benjamini–Hochberg step-up procedure controls the false discovery rate instead, and two easy ways to implement it wrong.

First explored December 2019·4 min read

University lab From my coursework in the MSc in Statistics and Machine Learning at Linköping University, rewritten as a note to review later.

multiple-testingfalse-discovery-ratehypothesis-testingfeature-selectionhigh-dimensional

In one of my labs I had a text dataset with 4,702 word features, and the task was to find which words really differ between two classes. My first move was to run a test per word and count the p-values below 0.05. Then the problem hit me. One test at the 5% level gives a false alarm 5% of the time. With 4,702 tests, you expect about 235 false alarms even if no word matters. So “how many p-values are below 0.05” is the wrong question. The Benjamini–Hochberg (BH) procedure answers a better one: which features can I call significant while keeping the share of false discoveries among them small?

What can go wrong when testing many hypotheses

Say we test mm null hypotheses and reject RR of them, and VV of those rejections are wrong (true nulls we rejected). Two ways to control errors:

  • Family-wise error rate: FWER=P(V≥1)\text{FWER} = P(V \ge 1), the chance of even one false rejection. Bonferroni controls it by testing each hypothesis at α/m\alpha/m. With mm in the thousands that’s extremely strict, and you find almost nothing.
  • False discovery rate: FDR=E[Vmax⁡(R,1)]\text{FDR} = E\left[\frac{V}{\max(R, 1)}\right], the expected fraction of rejections that are false. Controlling FDR at 5% means “among the features I flag, on average at most 5% are noise”. For screening many features this is usually what you want.

The BH step-up procedure

  1. Sort the p-values: p(1)≤p(2)≤⋯≤p(m)p_{(1)} \le p_{(2)} \le \dots \le p_{(m)}.
  2. Find the largest kk such that
p(k)≤km α.p_{(k)} \le \frac{k}{m}\,\alpha .
  1. Reject the hypotheses with the kk smallest p-values, p(1),…,p(k)p_{(1)}, \dots, p_{(k)}. All of them.

Benjamini and Hochberg proved this keeps FDR ≤m0mα≤α\le \frac{m_0}{m}\alpha \le \alpha (with m0m_0 the number of true nulls) when the tests are independent, and it’s known to hold under some kinds of positive dependence too.

The picture that helped me: plot the sorted p-values against their rank jj, and draw the line jα/mj\alpha/m. Find the last point that’s under the line. Everything to the left of it is rejected, even points that pop above the line before that.

Two ways I got it wrong

1. Counting raw p-values. My lab report said “281 null hypotheses are rejected if we use the rejection threshold 0.05”. That’s the count of unadjusted p-values below 0.05, with no correction at all. With 4,702 tests, a few hundred of those can be pure noise.

2. Comparing point by point. In a later version I sorted the p-values and labelled each one “reject” if p(j)<jα/mp_{(j)} < j\alpha/m. That’s not BH either. BH is a step-up rule: you look for the largest kk under the line and reject everything up to it. A p-value that sits slightly above its own line but before the cut-off point is still rejected. The point-by-point version can only reject fewer.

A related trap is computing BH-adjusted p-values by hand. They’re defined from the largest rank downward,

p~(j)=min⁡i≥j(mi p(i)) (capped at 1),\tilde p_{(j)} = \min_{i \ge j} \left( \frac{m}{i}\, p_{(i)} \right) \text{ (capped at 1)},

so each adjusted value is a running minimum of the adjusted values above it. When practising this for the exam I wrote the recursion with the previous raw p-value instead of the previous adjusted one, which gives wrong numbers. p.adjust(p, method = "BH") does it correctly; rejecting where p~≤α\tilde p \le \alpha is the same as the step-up rule.

bh_reject <- function(p, alpha = 0.05) {
  m <- length(p)
  o <- order(p)                                  # sort p-values, smallest first
  below <- which(p[o] <= (1:m) / m * alpha)      # sorted p_(j) under the line j * alpha / m
  k <- if (length(below)) max(below) else 0      # the LARGEST such j
  reject <- rep(FALSE, m)
  reject[o[seq_len(k)]] <- TRUE                  # reject p_(1), ..., p_(k), all of them
  reject
}
set.seed(6)
real <- rep(c(FALSE, TRUE), c(900, 100))         # 900 true nulls, 100 real effects
p <- ifelse(real, rbeta(1000, 0.3, 30), runif(1000))
rules <- list(raw = p < 0.05, bonferroni = p < 0.05 / 1000, bh = bh_reject(p))
sapply(rules, function(r) c(rejected = sum(r), false = sum(r & !real)))
#           raw  bonferroni  bh
# rejected  134          13  57
# false      40           0   3
all.equal(bh_reject(p), p.adjust(p, "BH") <= 0.05)   # TRUE

In this simulation the raw rule flags 134 features, but 40 of them are noise (30%). Bonferroni is clean but finds only 13. BH finds 57 with 3 false ones, about 5%, which is what it promises.

When to use it

  • Screening many features (genes, words, pixels, sensors) before a closer look.
  • When some false positives are acceptable but you want their share bounded.
  • Not when a single false positive is costly; then you want FWER control (Bonferroni or Holm).

Also keep in mind that each p-value is only as good as its test. A tt-test per word count assumes roughly normal, independent data, which sparse count features often aren’t.

Key takeaways

  • With mm tests at level α\alpha, expect about mαm\alpha false positives even when nothing is real.
  • FDR is the expected share of false discoveries among rejections; Bonferroni controls the much stricter FWER.
  • BH: sort p-values, find the largest kk with p(k)≤kα/mp_{(k)} \le k\alpha/m, reject the kk smallest.
  • It’s a step-up rule, not a point-by-point comparison; use p.adjust(..., "BH") to avoid hand-rolled mistakes.

Resources