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.
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
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.
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
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.
| 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.
m includes the tests you ran and did not report. A
correction applied to the six you liked is not a correction.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