Causal Inference I(g) — A/B Testing: Multiple Comparisons and False Discovery Rate

When you test many metrics (or many variants), 5% per test is not 5% overall

The peeking notebook showed that looking many times in sequence inflates the error rate. The same disaster strikes when you look many times in parallel — testing a treatment against a dashboard of dozens of metrics, or comparing many variants (A/B/C/D/…) at once. Each test at $\alpha=0.05$ has a 5% false-positive chance, so a dashboard of 200 metrics in which 160 are truly null yields about 8 "significant" movers by chance alone, and the probability of at least one false positive races toward certainty.

This notebook builds the two families of corrections and the trade-off between them:

1. The multiplicity problem — how the family-wise error rate (FWER) explodes with the number of tests.
2. FWER control — Bonferroni and Holm: bound the probability of any false positive. Safe, but so conservative that real effects are missed.
3. FDR control — Benjamini-Hochberg: bound the expected proportion of false discoveries. The right tool for large metric screens — far more power at a controlled error rate.

Simulation-first with a known set of true/null metrics. Python leads (from-scratch Bonferroni/Holm/BH); the R companion uses p.adjust. This is the across-metrics companion to the across-time Sequential Testing notebook (1d) — same inflation, different axis.

1. The multiplicity problem

We simulate an experiment scored on 200 metrics, of which 40 truly moved and 160 are null. Each metric yields one z-test p-value. Testing each at $\alpha=0.05$, the family-wise error rate (FWER) — the chance of at least one false positive among the 160 nulls — is $1-(1-0.05)^{160}\approx 1$: essentially every experiment produces spurious "winners." With large metric dashboards, uncorrected testing guarantees false discoveries.

import numpy as np, matplotlib.pyplot as plt, warnings
from scipy.stats import norm
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
K=200; m1=40; m0=K-m1
def trial(seed, effect=3.0):
    r=np.random.default_rng(seed)
    z=np.r_[r.normal(effect,1,m1), r.normal(0,1,m0)]     # first m1 real, rest null
    p=2*(1-norm.cdf(np.abs(z))); truth=np.r_[np.ones(m1,bool),np.zeros(m0,bool)]
    return p,truth
ks=np.arange(1,201); fwer_theory=1-(1-0.05)**ks
print(f"Family-wise error rate (prob of >=1 false positive) under naive 0.05 testing:")
for k in [1,5,10,50,160]: print(f"  {k:3d} null metrics: FWER = {1-(1-0.05)**k:.3f}")
p,truth=trial(0)
print(f"\nOne example experiment ({m0} null metrics): {(p[~truth]<0.05).sum()} false positives by chance (expected {0.05*m0:.0f})")
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
ax[0].plot(ks,fwer_theory,color=RED,lw=2); ax[0].axhline(0.05,color=GREEN,ls="--",label="nominal 0.05")
ax[0].set_xlabel("number of (null) tests"); ax[0].set_ylabel("P(at least one false positive)"); ax[0].set_title("FWER explodes with the number of metrics"); ax[0].legend()
ax[1].hist(p[truth],bins=20,color=GREEN,alpha=.7,label="truly moved (40)"); ax[1].hist(p[~truth],bins=20,color=GREY,alpha=.6,label="null (160)")
ax[1].axvline(0.05,color=RED,ls="--"); ax[1].set_xlabel("p-value"); ax[1].set_ylabel("metrics"); ax[1].set_title("Null p-values are uniform -> ~8 fall below 0.05 by luck"); ax[1].legend()
plt.tight_layout(); plt.show()
print("Null p-values are uniform on [0,1], so ~5% of the 160 nulls land under 0.05 every time. Reporting those as wins is the")
print("multiple-comparisons trap -- the parallel-testing twin of peeking.")

Family-wise error rate (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

One example experiment (160 null metrics): 6 false positives by chance (expected 8)

Null p-values are uniform on [0,1], so ~5% of the 160 nulls land under 0.05 every time. Reporting those as wins is the
multiple-comparisons trap -- the parallel-testing twin of peeking.

2. Family-wise error control — Bonferroni and Holm

The strictest goal is to bound the probability of any false positive (the FWER) at $\alpha$. Bonferroni simply tests each metric at $\alpha/K$ — guaranteed FWER $\le\alpha$, but brutally conservative: with $K=200$ the per-test threshold is 0.00025, so only the very strongest effects survive and power collapses. Holm's step-down is uniformly more powerful (and still controls FWER): sort the p-values and compare the $i$-th smallest to $\alpha/(K-i+1)$, stopping at the first failure. Both keep the FWER at the nominal level; both pay for it in missed true effects.

def holm(p,a=0.05):
    n=len(p); order=np.argsort(p); rej=np.zeros(n,bool)
    for i,idx in enumerate(order):
        if p[idx]<=a/(n-i): rej[idx]=True
        else: break
    return rej
def evaluate(method,nsim=500):
    fdp=[];pw=[];fwer=0
    for s in range(nsim):
        p,truth=trial(s)
        if method=="naive": rej=p<0.05
        elif method=="bonferroni": rej=p<0.05/K
        elif method=="holm": rej=holm(p)
        elif method=="bh": rej=bh(p)
        R=rej.sum();V=(rej&~truth).sum();S=(rej&truth).sum()
        fdp.append(V/max(R,1));pw.append(S/m1);fwer+=(V>=1)
    return np.mean(fdp),np.mean(pw),fwer/nsim
def bh(p,q=0.05):
    n=len(p);order=np.argsort(p);thr=q*np.arange(1,n+1)/n;passed=p[order]<=thr
    k=np.where(passed)[0];rej=np.zeros(n,bool)
    if len(k): rej[order[:k.max()+1]]=True
    return rej
for meth in ["naive","bonferroni","holm"]:
    fdr,power,fw=evaluate(meth)
    print(f"  {meth:12s}: FWER {fw:.3f}, power {power:.3f}, false-discovery rate {fdr:.3f}")
print("  Bonferroni/Holm hold FWER at ~0.05 but detect only ~1 in 4 of the true movers.")

  naive       : FWER 1.000, power 0.857, false-discovery rate 0.189

  bonferroni  : FWER 0.062, power 0.255, false-discovery rate 0.006
  holm        : FWER 0.062, power 0.260, false-discovery rate 0.006
  Bonferroni/Holm hold FWER at ~0.05 but detect only ~1 in 4 of the true movers.

3. False discovery rate control — Benjamini-Hochberg

Controlling the FWER is often the wrong goal for a metric screen: with hundreds of metrics you don't need zero false positives, you need the proportion of your reported winners that are spurious to be small. That is the false discovery rate (FDR) — the expected fraction of rejections that are false — and Benjamini-Hochberg (BH) controls it. Sort the p-values $p_{(1)}\le\dots\le p_{(K)}$, find the largest $i$ with $p_{(i)}\le \frac{i}{K}q$, and reject all metrics up to it. BH controls the FDR at $q$ under independence (and positive dependence), and delivers far more power than FWER methods — the standard choice for experimentation platforms scoring many metrics.

for meth in ["bonferroni","holm","bh"]:
    fdr,power,fw=evaluate(meth)
    print(f"  {meth:12s}: FDR {fdr:.3f}, power {power:.3f}")
print(f"  BH controls FDR at 0.05 while recovering ~2x the true movers of Bonferroni.")
# visualize the BH procedure on one experiment
p,truth=trial(3); n=len(p); order=np.argsort(p); ps=p[order]; ts=truth[order]
line=0.05*np.arange(1,n+1)/n; kmax=np.where(ps<=line)[0]
cut=kmax.max() if len(kmax) else -1
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].plot(np.arange(1,n+1),line,color=RED,lw=2,label="BH line  i/K * q")
ax[0].scatter(np.arange(1,n+1),ps,c=[GREEN if t else GREY for t in ts],s=12)
if cut>=0: ax[0].axvline(cut+1,color=BLUE,ls=":",label=f"reject the smallest {cut+1}")
ax[0].set_xlim(0,80); ax[0].set_ylim(0,0.06); ax[0].set_xlabel("rank of p-value"); ax[0].set_ylabel("p-value"); ax[0].set_title("Benjamini-Hochberg: reject below the i/K*q line"); ax[0].legend(fontsize=8)
res={m:evaluate(m) for m in ["naive","bonferroni","holm","bh"]}
x=np.arange(4); w=0.38; names=list(res)
ax[1].bar(x-w/2,[res[m][1] for m in names],w,color=GREEN,label="power (true movers found)")
ax[1].bar(x+w/2,[res[m][0] for m in names],w,color=RED,label="false-discovery rate")
ax[1].axhline(0.05,color="k",ls="--",lw=1); ax[1].set_xticks(x); ax[1].set_xticklabels(names,fontsize=8); ax[1].set_title("BH: best power at a controlled error rate"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print("BH accepts a small, controlled fraction of false discoveries (~5%) in exchange for detecting far more real effects than")
print("the all-or-nothing FWER methods. For a dashboard of many metrics or many variants, BH is the default.")

  bonferroni  : FDR 0.006, power 0.255
  holm        : FDR 0.006, power 0.260
  bh          : FDR 0.040, power 0.613
  BH controls FDR at 0.05 while recovering ~2x the true movers of Bonferroni.

BH accepts a small, controlled fraction of false discoveries (~5%) in exchange for detecting far more real effects than
the all-or-nothing FWER methods. For a dashboard of many metrics or many variants, BH is the default.

4. Summary

Testing many metrics or variants at once inflates error exactly as peeking does over time:

• The multiplicity problem: with 160 null metrics, naive 0.05 testing produced ~8 false positives per experiment and a family-wise error rate of ~100% — spurious winners guaranteed.
• FWER control (Bonferroni, Holm) bounds the probability of any false positive; it held the FWER at 5% but power fell to ~0.26 (only the strongest effects survive). Holm dominates Bonferroni and should be preferred when FWER is the goal.
• FDR control (Benjamini-Hochberg) bounds the expected proportion of false discoveries; it held the FDR at ~0.04 while recovering ~0.61 of the true movers — roughly twice the power of the FWER methods.

Practical guidance: pre-register a small number of primary metrics (test those with FWER control or no correction) and apply BH to the broader secondary/guardrail dashboard; treat any un-corrected scan of many metrics as hypothesis-generating, not confirmatory. Cross-links: this is the parallel-testing analogue of the sequential peeking problem (1d) — both are multiplicity, across metrics vs across time; the FDR/threshold logic connects to the conformal prediction and model-selection notebooks in the ML arc; and pre-registering primary metrics echoes the pre-analysis discipline behind honest RCTs (1a). The R companion uses base R's p.adjust (bonferroni, holm, BH).