---
title: "Which multiple-comparison correction do you actually need?"
subtitle: "Bonferroni, Holm, and Benjamini-Hochberg on one real dataset"
output:
  html_document:
    toc: true
    toc_float: true
    theme: readable
    df_print: kable
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
options(scipen = 999)
```

## 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.

```{r load}
p <- read.csv("wine_correlation_pvalues.csv", stringsAsFactors = FALSE)
p <- p[order(p$p_value), ]
m <- nrow(p)
alpha <- 0.05
m
```

## What "significant" means before you correct

```{r raw}
raw_hits <- sum(p$p_value < alpha)
raw_hits
```

`r raw_hits` of `r m` 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:

```{r expected-false}
alpha * m
```

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

```{r corrections}
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)
)
```

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

```{r verify-against-base-r}
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))
```

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

**Holm and Bonferroni return the same answer here: `r holm`.**

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:

```{r near-miss}
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)
)
```

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

**The choice that actually moved the answer was FWER versus FDR**: `r bh - holm` extra
discoveries, `r bh` against `r holm`.

## 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 | `r holm` survivors | `r bh` survivors |

The FDR promise is weaker on purpose. It says: of the `r bh` 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 `r bh - holm` 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 `r bh`.
- **Very small p-values are still exactly representable here, contrary to what you may
  have been told about this dataset.** The smallest is `r sprintf("%.3g", min(p$p_value))`,
  and there are `r sum(p$p_value == 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.