Causal Inference I(g) — A/B Testing: Multiple Comparisons and FDR (R companion)

p.adjust — Bonferroni, Holm, and Benjamini-Hochberg

This companion reproduces the Python notebook with base R's p.adjust: the multiplicity problem (testing many metrics inflates the family-wise error rate to ~1), FWER control via bonferroni / holm (safe but low power), and FDR control via BH (Benjamini-Hochberg — far more power at a controlled false-discovery rate). Same simulation of 200 metrics with 40 true movers.

1. The multiplicity problem

200 metrics, 40 truly moved, 160 null. Under naive 0.05 testing the chance of at least one false positive among 160 nulls is essentially 1, and about 5% of the nulls (~8) cross 0.05 by luck every time.

set.seed(0); K<-200; m1<-40; m0<-K-m1
trial<-function(seed, effect=3){ set.seed(seed)
  z<-c(rnorm(m1,effect,1), rnorm(m0,0,1)); p<-2*(1-pnorm(abs(z))); truth<-c(rep(TRUE,m1),rep(FALSE,m0)); list(p=p,truth=truth) }
cat("FWER (prob of >=1 false positive) under naive 0.05 testing:\n")
for(k in c(1,5,10,50,160)) cat(sprintf("  %3d null metrics: FWER = %.3f\n", k, 1-(1-0.05)^k))
tr<-trial(0); cat(sprintf("\nExample: %d false positives among %d nulls (expected %.0f)\n", sum(tr$p[!tr$truth]<0.05), m0, 0.05*m0))
options(repr.plot.width=13, repr.plot.height=4.2); par(mfrow=c(1,2))
ks<-1:200; plot(ks, 1-(1-0.05)^ks, type="l", col="#c53030", lwd=2, xlab="number of null tests", ylab="P(>=1 false positive)", main="FWER explodes with the number of metrics"); abline(h=0.05, col="#2f855a", lty=2)
hist(tr$p[!tr$truth], breaks=20, col="grey", xlim=c(0,1), main="Null p-values uniform -> ~8 below 0.05", xlab="p-value"); abline(v=0.05, col="#c53030", lty=2)
par(mfrow=c(1,1))

FWER (prob of >=1 false positive) under naive 0.05 testing:

    1 null metrics: FWER = 0.050
    5 null metrics: FWER = 0.226
   10 null metrics: FWER = 0.401
   50 null metrics: FWER = 0.923
  160 null metrics: FWER = 1.000

Example: 5 false positives among 160 nulls (expected 8)

2. FWER and FDR control via p.adjust

p.adjust implements every standard correction. method="bonferroni" and "holm" control the family-wise error rate (Holm uniformly more powerful); method="BH" controls the false-discovery rate. We evaluate each over many simulated experiments, tracking realized FWER, FDR, and power (fraction of the 40 true movers detected).

evaluate<-function(method, nsim=500){ fdp<-pw<-fw<-numeric(nsim)
  for(s in 1:nsim){ tr<-trial(s)
    rej<- if(method=="naive") tr$p<0.05 else p.adjust(tr$p, method=method)<0.05
    R<-sum(rej); V<-sum(rej & !tr$truth); S<-sum(rej & tr$truth)
    fdp[s]<-V/max(R,1); pw[s]<-S/m1; fw[s]<-(V>=1) }
  c(FWER=mean(fw), FDR=mean(fdp), power=mean(pw)) }
cat(sprintf("%-12s %6s %6s %6s\n","method","FWER","FDR","power"))
for(m in c("naive","bonferroni","holm","BH")){ r<-evaluate(m)
  cat(sprintf("%-12s %6.3f %6.3f %6.3f\n", m, r["FWER"], r["FDR"], r["power"])) }
cat("\nBonferroni/Holm hold FWER ~0.05 (power ~0.26); BH holds FDR ~0.05 with ~2x the power.\n")

method         FWER    FDR  power

naive         1.000  0.186  0.852
bonferroni    0.036  0.003  0.255
holm          0.036  0.003  0.260
BH            0.612  0.037  0.610

Bonferroni/Holm hold FWER ~0.05 (power ~0.26); BH holds FDR ~0.05 with ~2x the power.

3. Visualizing Benjamini-Hochberg

BH sorts the p-values and rejects all up to the largest $i$ with $p_{(i)}\le\frac{i}{K}q$ — the points below the BH line. The bar chart contrasts power and realized error across methods: FWER methods are safe but weak; BH recovers far more true movers at a controlled FDR.

tr<-trial(3); n<-length(tr$p); ord<-order(tr$p); ps<-tr$p[ord]; ts<-tr$truth[ord]; line<-0.05*(1:n)/n
below<-which(ps<=line); cut<- if(length(below)) max(below) else 0
options(repr.plot.width=13, repr.plot.height=4.4); par(mfrow=c(1,2))
plot(1:80, line[1:80], type="l", col="#c53030", lwd=2, ylim=c(0,0.06), xlab="rank of p-value", ylab="p-value", main="Benjamini-Hochberg: reject below i/K*q")
points(1:80, ps[1:80], col=ifelse(ts[1:80],"#2f855a","grey"), pch=19, cex=.6); if(cut>0) abline(v=cut, col="#2b6cb0", lty=3)
res<-sapply(c("naive","bonferroni","holm","BH"), evaluate)
bp<-barplot(rbind(power=res["power",], FDR=res["FDR",]), beside=TRUE, col=c("#2f855a","#c53030"), names.arg=colnames(res), main="BH: best power at a controlled error rate", ylab="rate", legend.text=c("power","FDR"), args.legend=list(x="topright",bty="n")); abline(h=0.05, lty=2)
par(mfrow=c(1,1))

4. Summary

p.adjust reproduced the multiplicity story: naive testing of 200 metrics gave a ~100% family-wise error rate and ~8 false positives per experiment; Bonferroni/Holm held the FWER at 5% but detected only ~26% of true movers; Benjamini-Hochberg held the FDR at ~5% while recovering ~61% — roughly double the power.

Guidance: pre-register a few primary metrics (FWER control or none) and apply BH to the wider secondary/guardrail dashboard; a raw scan of many metrics is hypothesis-generating, not confirmatory. Cross-links: the across-metrics twin of the across-time peeking problem (1d); the threshold/selection logic connects to conformal prediction and model selection in the ML arc; pre-registering primaries echoes the pre-analysis discipline of honest RCTs (1a).