Multiple Testing and Peeking

26.1 The garden of forking paths (30 min)

Concept. A 5% false positive rate holds for one pre-specified test. Real readouts look at many metrics, many segments and many days. Each extra look is another chance for noise to cross the line. Three habits produce most false discoveries in business experiments:

1. Many comparisons: 4 metrics × 5 regions × 3 devices = 60 chances to find something.
2. Peeking: checking results daily and stopping as soon as they look good.
3. Reporting the best: presenting the segment or metric with the biggest effect, without saying it was chosen because it was biggest.

None of these is dishonest in intent. All of them are fixed by the same tool: decide in the plan which comparisons are confirmatory, and correct for the rest.

26.2 Twenty comparisons in an A/A test (60 min)

Create day26/many_comparisons.py. It runs 300 A/A tests and reads each one 20 ways (4 metrics × 5 regions):

# day26/many_comparisons.py   (python -m day26.many_comparisons)
"""
An A/A test read 20 ways (4 metrics x 5 regions). Nothing works, yet something usually 'wins'.
Then three corrections: Bonferroni, Holm (family-wise error) and Benjamini-Hochberg (false discovery rate).
"""
import numpy as np
import polars as pl
from scipy.stats import ttest_ind
from statsmodels.stats.multitest import multipletests

from causalkit.lab import population, run

pop = population()
region = pop["region"].to_numpy()
METRICS = ["revenue", "orders", "converted", "margin"]
rng = np.random.default_rng(26)

any_hit = {"none": 0, "bonferroni": 0, "holm": 0, "fdr_bh": 0}
REPS = 300
for r in range(REPS):
    z = rng.permutation(np.repeat([0, 1], pop.height // 2))
    out = run(pop, z, "aa_test", seed=3000 + r)
    pvals = []
    for m in METRICS:
        y = out[m].to_numpy()
        for g in np.unique(region):
            sel = region == g
            pvals.append(ttest_ind(y[sel & (z == 1)], y[sel & (z == 0)], equal_var=False).pvalue)
    pvals = np.array(pvals)
    any_hit["none"] += (pvals < 0.05).any()
    for method in ["bonferroni", "holm", "fdr_bh"]:
        any_hit[method] += multipletests(pvals, alpha=0.05, method=method)[0].any()

print(f"{len(METRICS) * 5} comparisons per A/A test, {REPS} A/A tests")
with pl.Config(fmt_str_lengths=70, tbl_width_chars=110):
    print(pl.DataFrame({"correction": list(any_hit),
                        "share of A/A tests with at least one significant result":
                            [round(v / REPS, 3) for v in any_hit.values()]}))
print("Without correction, 1 - 0.95^20 = 0.64 would be expected if the 20 tests were independent.")

python -m day26.many_comparisons

20 comparisons per A/A test, 300 A/A tests
│ correction ┆ share of A/A tests with at least one significant result │
│ none       ┆ 0.467                                                   │
│ bonferroni ┆ 0.033                                                   │
│ holm       ┆ 0.033                                                   │
│ fdr_bh     ┆ 0.063                                                   │
Without correction, 1 - 0.95^20 = 0.64 would be expected if the 20 tests were independent.

Nearly half of all A/A tests (where nothing happened) produce at least one "significant" regional result. It's below the 64% of 20 independent tests because the four metrics are correlated (a region's revenue and margin move together), but it's still a coin toss.

26.3 Corrections: Bonferroni, Holm, Benjamini–Hochberg (45 min)

In an A/A test every discovery is false, so FDR equals FWER and all three keep the rate near 5%. The difference shows when there are real effects.

Try it. Run the same 20 comparisons on email_v2, which really works in every region. On average, how many of the 20 real effects does each method find?

# day26/real_effect_corrections.py   (python -m day26.real_effect_corrections)
import numpy as np
from scipy.stats import ttest_ind
from statsmodels.stats.multitest import multipletests

from causalkit.lab import population, run

pop = population()
region = pop["region"].to_numpy()
rng = np.random.default_rng(262)
found = {"none": 0, "holm": 0, "fdr_bh": 0}
for r in range(100):
    z = rng.permutation(np.repeat([0, 1], pop.height // 2))
    out = run(pop, z, "email_v2", seed=8000 + r)
    p = []
    for m in ["revenue", "orders", "converted", "margin"]:
        y = out[m].to_numpy()
        for g in np.unique(region):
            s = region == g
            p.append(ttest_ind(y[s & (z == 1)], y[s & (z == 0)], equal_var=False).pvalue)
    p = np.array(p)
    found["none"] += (p < 0.05).sum()
    for method in ["holm", "fdr_bh"]:
        found[method] += multipletests(p, method=method)[0].sum()
print({k: round(float(v) / 100, 2) for k, v in found.items()}, "true effects found per experiment (of 20)")

{'none': 7.34, 'holm': 2.36, 'fdr_bh': 4.76} true effects found per experiment (of 20)

26.4 Peeking: checking every day (45 min)

Concept. A p-value is valid at the sample size you planned. If you compute it every day and stop the first time it dips below 0.05, you give noise 28 chances instead of one.

Create day26/peeking.py. It simulates 2,000 A/A tests where 700 new customers per group arrive each day, with revenue drawn from last month's real distribution:

# day26/peeking.py   (python -m day26.peeking)
"""
Peeking: an A/A test checked every day for 28 days, stopped the first time p < 0.05.
Customers' 4-week revenue is resampled from last month's real distribution.
"""
import numpy as np
from scipy.stats import norm

from causalkit.lab import population, run

pop = population()
revenue = run(pop, np.zeros(pop.height), "aa_test", seed=99)["revenue"].to_numpy()
rng = np.random.default_rng(260)

SIMS, DAYS, PER_DAY = 2000, 28, 700            # 700 new customers per group per day (~40k in total)
stopped_by = np.zeros(SIMS, dtype=bool)
final_sig = np.zeros(SIMS, dtype=bool)
first_day = np.full(SIMS, -1)
for s in range(SIMS):
    a = rng.choice(revenue, size=(DAYS, PER_DAY))           # group A, by day
    b = rng.choice(revenue, size=(DAYS, PER_DAY))           # group B: same distribution (A/A)
    n = PER_DAY * np.arange(1, DAYS + 1)
    sa, sb = a.sum(1).cumsum(), b.sum(1).cumsum()
    qa, qb = (a**2).sum(1).cumsum(), (b**2).sum(1).cumsum()
    va, vb = qa / n - (sa / n) ** 2, qb / n - (sb / n) ** 2
    zstat = (sa / n - sb / n) / np.sqrt(va / n + vb / n)
    p = 2 * norm.sf(np.abs(zstat))
    hits = np.flatnonzero(p < 0.05)
    if hits.size:
        stopped_by[s] = True
        first_day[s] = hits[0] + 1
    final_sig[s] = p[-1] < 0.05

print(f"Analyse once, on day 28:                false positive rate {final_sig.mean():.3f}")
print(f"Peek daily, stop at the first p < 0.05: false positive rate {stopped_by.mean():.3f}")
print(f"Of the peeking 'wins', share declared in the first week: {(first_day[stopped_by] <= 7).mean():.0%}")

python -m day26.peeking

Analyse once, on day 28:                false positive rate 0.057
Peek daily, stop at the first p < 0.05: false positive rate 0.292
Of the peeking 'wins', share declared in the first week: 58%

Daily peeking turns a 5% error rate into 29%. Most false wins appear in the first week, when samples are small and noisy, which is exactly when an eager team is most tempted to call the result. Day 27 shows how to look early legitimately.

26.5 Day 26 mini-project: the winner's curse (75 min)

Goal. email_v2 lifts orders by exactly 5% in every region. A dashboard shows the regional results and the team reports the best one: "+X% in our best region!" How big does that winner look on average, and which region tends to win?

Create day26/winners_curse.py:

# day26/winners_curse.py   (python -m day26.winners_curse)
"""
Mini-project: the 'best region' in an experiment. email_v2 lifts orders by 5% in EVERY region.
Report the region with the biggest estimated lift, 300 times: how big does the winner look?
"""
import numpy as np
import polars as pl

from causalkit.lab import population, run

pop = population()
region = pop["region"].to_numpy()
regions = np.unique(region)
rng = np.random.default_rng(261)

winners, winner_names, all_lifts = [], [], []
for r in range(300):
    z = rng.permutation(np.repeat([0, 1], pop.height // 2))
    y = run(pop, z, "email_v2", seed=4000 + r)["orders"].to_numpy()
    lifts = {g: y[(region == g) & (z == 1)].mean() / y[(region == g) & (z == 0)].mean() - 1 for g in regions}
    best = max(lifts, key=lifts.get)
    winners.append(lifts[best]); winner_names.append(best); all_lifts.extend(lifts.values())

print(f"True lift in every region:            +5.0%")
print(f"Average estimated lift, all regions:  {np.mean(all_lifts):+.1%}")
print(f"Average lift of the 'winning' region: {np.mean(winners):+.1%}")
print(pl.Series("winning region", winner_names).value_counts().sort("count", descending=True))

python -m day26.winners_curse

True lift in every region:            +5.0%
Average estimated lift, all regions:  +5.2%
Average lift of the 'winning' region: +8.7%
│ winning region ┆ count │
│ ATL            ┆ 81    │
│ AB             ┆ 61    │
│ BC             ┆ 61    │
│ QC             ┆ 60    │
│ ON             ┆ 37    │

• The average region is estimated correctly (+5.2%); the selected region overstates the effect by 74% (+8.7%).
• The smallest regions win most often (ATL and AB are 12% of customers each; ON is 38%). Small samples are noisier, so their extremes are more extreme. "Our best results come from small markets" is often just this.

Add a Reporting policy section to docs/EXPERIMENT_PLAN_email_v2.md: one primary metric (confirmatory, no correction), secondary metrics with Holm, segments as exploratory with Benjamini–Hochberg and a sentence warning that the largest segment effect is biased upward.

26.6 Save your work with Git (15 min)

git add .
git commit -m "Day 26: multiple comparisons, corrections, peeking, winner's curse"
git push

Day 26 checklist

• I can explain why many comparisons inflate false positives, and estimate by how much
• I can choose between Bonferroni, Holm and Benjamini–Hochberg
• I can show by simulation what daily peeking does to the error rate
• I can explain the winner's curse and why small segments win most often
• I have a written reporting policy for metrics and segments