Purpose

You ran a lot of tests at once. Some came back significant. How many of those do you actually believe?

This document works one real dataset end to end and answers the question the textbooks usually skip: not “should I correct” (yes) but “correct for what”. Every number below is computed in this file from the committed CSV. Nothing is asserted.

The data

UCI Wine Quality (red): 1,599 wines, 12 chemistry variables, every pairwise Pearson correlation. That is choose(12, 2) = 66 hypothesis tests run on one dataset, which is exactly the situation multiplicity corrections exist for.

This is also the corpus the live standard_multiple_comparisons tool is verified against, so this document, the tool and the video all read one set of numbers.

p <- read.csv("wine_correlation_pvalues.csv", stringsAsFactors = FALSE)
p <- p[order(p$p_value), ]
m <- nrow(p)
alpha <- 0.05
m
## [1] 66

What “significant” means before you correct

raw_hits <- sum(p$p_value < alpha)
raw_hits
## [1] 55

55 of 66 tests come back below .05.

That number is not trustworthy and the reason is arithmetic, not opinion. If every one of the 66 nulls were true, the expected number of false positives at alpha = .05 is:

alpha * m
## [1] 3.3

So you would expect roughly three “findings” from pure noise. The corrections below differ in what they promise to do about that.

The three corrections, computed

bonferroni <- sum(p$p_value < alpha / m)

# Holm: step DOWN the sorted p-values, comparing each to alpha/(m-i+1),
# and stop at the first failure. Everything from there on is retained.
holm <- 0
for (i in seq_len(m)) {
  if (p$p_value[i] < alpha / (m - i + 1)) holm <- holm + 1 else break
}

# Benjamini-Hochberg: the LARGEST i whose p is under i*alpha/m.
bh <- 0
for (i in seq_len(m)) {
  if (p$p_value[i] <= i * alpha / m) bh <- i
}

data.frame(
  method    = c("none (raw p < .05)", "Bonferroni", "Holm", "Benjamini-Hochberg"),
  controls  = c("nothing", "FWER", "FWER", "FDR"),
  survivors = c(raw_hits, bonferroni, holm, bh)
)
method controls survivors
none (raw p < .05) nothing 55
Bonferroni FWER 43
Holm FWER 43
Benjamini-Hochberg FDR 54

Cross-check against R’s own implementation, because a hand-rolled loop is exactly the kind of thing that is quietly wrong:

c(bonferroni = sum(p.adjust(p$p_value, "bonferroni") < alpha),
  holm       = sum(p.adjust(p$p_value, "holm")       < alpha),
  BH         = sum(p.adjust(p$p_value, "BH")         < alpha))
## bonferroni       holm         BH 
##         43         43         54

The finding, which is not the one most explanations lead with

Holm and Bonferroni return the same answer here: 43.

That is not a bug and it is not unusual. Holm is uniformly at least as powerful as Bonferroni, so it can never do worse, but it only does better when a p-value falls in the narrow band between alpha/m and alpha/(m-i+1). Real p-value distributions tend to be bimodal: genuine effects are tiny, everything else is large, and that band is empty.

Here it was very nearly not empty:

i <- holm + 1                       # the first test Holm rejected
data.frame(
  rank          = i,
  comparison    = p$comparison[i],
  p_value       = p$p_value[i],
  holm_bar      = alpha / (m - i + 1),
  missed_by     = p$p_value[i] - alpha / (m - i + 1)
)
rank comparison p_value holm_bar missed_by
44 volatile acidity ~ total sulfur dioxide 0.0022139 0.0021739 0.0000399

It missed gaining one discovery by about four hundred-thousandths.

The choice that actually moved the answer was FWER versus FDR: 11 extra discoveries, 54 against 43.

What the two guarantees actually promise

Bonferroni / Holm Benjamini-Hochberg
controls family-wise error rate false discovery rate
the promise P(one or more false positives among ALL tests) <= alpha the expected proportion of your rejections that are false <= alpha
you should want it when one false positive is expensive: a drug label, a published claim, a shipped decision you are screening, and a known fraction of dead ends is an acceptable cost of finding the live ones
on this data 43 survivors 54 survivors

The FDR promise is weaker on purpose. It says: of the 54 correlations you are about to chase, about 5% may be noise. If chasing a dead end is cheap, that is a good trade for the 11 extra leads. If each one costs a clinical trial, it is not.

The chooser

  1. Are you screening or confirming? Screening (which of these 66 is worth a second look?) points at BH. Confirming (which of these can I publish?) points at FWER.
  2. If FWER: use Holm, never Bonferroni. Same guarantee, never fewer discoveries, occasionally more. There is no scenario where Bonferroni is the better of the two; it survives on habit and on being one line of mental arithmetic.
  3. How many tests are you actually running? m includes the tests you ran and did not report. A correction applied to the six you liked is not a correction.

Honest limitations

  • These 66 tests are not independent. Wine chemistry variables are correlated with each other, so the tests share information. Bonferroni and Holm are valid regardless (they make no independence assumption). BH as used here assumes independence or positive regression dependence; under arbitrary dependence the conservative variant is Benjamini-Yekutieli, which would return fewer than 54.
  • Very small p-values are still exactly representable here, contrary to what you may have been told about this dataset. The smallest is 4.06e-220, and there are 0 exact zeros against a double-precision floor near 2.2e-308. Nothing underflowed. (The tool’s registry entry carries an “underflow floored at 1e-300” note; that describes an earlier computation of this corpus, not this one, and I checked rather than repeating it.) What these tiny values do mean is that the correlation is unmistakable, not that it is infinitely certain: they are conditional on the linearity and normality assumptions of a Pearson test, which no p-value re-examines.
  • This is one dataset. The finding “Holm equals Bonferroni” is a property of this p-value distribution, not a general law. The general law is only that Holm is never worse.

References

  • Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics 6(2), 65-70.
  • Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate. JRSS B 57(1), 289-300.
  • Benjamini, Y. and Yekutieli, D. (2001). The control of the false discovery rate under dependency. Annals of Statistics 29(4), 1165-1188.
  • Cortez, P. et al. (2009). Wine Quality Data Set. UCI Machine Learning Repository.

Your turn

Bring your own data and the question you actually need answered.

CympleData Scientist Send me your data and question, I’ll send you the analytics. ds@mcpanalytics.ai