8. Experimentation and Decision Science
An experiment is the only method that establishes cause with no assumptions about what else might explain the result, which is why it is the foundation of evidence-based decisions. This chapter covers how to design a valid experiment, how to run and interpret one under real operational constraints, what to do when randomization is impossible, and how to convert a measured effect into a decision. Chapter 3 supplied the statistical machinery for testing and causal identification; this chapter is about applying it without arriving at a confident wrong answer.
<a id="81-experimental-design"></a>
8.1 Experimental Design
Design decisions made before a single observation is collected determine whether the result will mean anything; this section covers what randomization buys you, the unit you assign, and how to size a test so it can detect what you care about.
Randomized Controlled Trials
The problem that randomization solves is deeper than it first appears. To know whether a change caused an outcome, you would need to see the same user both receiving and not receiving the change, at the same moment. That is impossible. This is the fundamental problem of causal inference: for every unit you observe one outcome and the other is forever unobservable.
Chapter 3 introduced the potential outcomes notation for this. Each unit has an outcome under treatment and an outcome under control, and the causal effect for that unit is the difference. You always see exactly one of the two.
The workaround is to stop trying to measure individual effects and measure an average instead. If you have two groups that are equivalent in every respect except the treatment, then the difference in their average outcomes estimates the average effect. The entire design problem reduces to one question: how do you get two genuinely equivalent groups?
Randomization is the answer, and the reason it works is worth understanding precisely. Matching on variables you can measure balances those variables. It does nothing about the ones you did not measure or did not think of. Randomization balances everything in expectation, measured or not, known or unknown, because the assignment is independent of every characteristic a unit has. That is a property no amount of statistical adjustment can replicate, and it is why an experiment is qualitatively stronger evidence than any observational analysis.
The analogy: adjusting for confounders is like checking a suitcase for the items on your list. Randomization is like weighing two suitcases packed by coin flip from the same pile, without needing a list at all.
Randomization does not guarantee balance in any single experiment. With small samples, chance imbalance is common, which is why you check balance afterwards and why stratification, covered later, is worth the effort.
The assumptions an experiment still needs are usually left implicit, and each has a practical failure mode:
SUTVA, the stable unit treatment value assumption, has two parts. No interference means one unit's treatment does not affect another's outcome, which fails for social features, marketplaces, and anything with shared inventory. Consistency means there is one well-defined version of the treatment, which fails if the feature changes mid-experiment.
Compliance. Units assigned to treatment must actually receive it. When they do not, comparing by assignment estimates the effect of being offered the treatment, which is a different and often more useful quantity than the effect of receiving it.
No differential attrition. If treatment causes people to leave before the outcome is measured, the surviving groups are no longer comparable, and the randomization has been undone by the analysis.
import numpy as np
import pandas as pd
from scipy import stats
rng = np.random.default_rng(0)
n = 4000
# --- THE FUNDAMENTAL PROBLEM, made visible with simulated ground truth. ---
# In a real experiment only ONE of these two columns is ever observable.
engagement = rng.gamma(2, 12, n) # a pre-existing trait
y0 = 20 + 0.9 * engagement + rng.normal(0, 4, n) # outcome if untreated
effect = 2.5 + 0.05 * engagement # a heterogeneous effect
y1 = y0 + effect # outcome if treated
truth = pd.DataFrame({"engagement": engagement, "y0": y0, "y1": y1,
"unit_effect": effect})
print(f"TRUE average treatment effect: {truth['unit_effect'].mean():.3f}")
print("No design can recover the per-unit effect; only the average.")
# --- WHY SELF-SELECTION FAILS: adopters differ from non-adopters. ---
# Highly engaged users opt in more, so the groups differ before treatment.
p_optin = 1 / (1 + np.exp(-(engagement - 24) / 6))
self_selected = rng.random(n) < p_optin
naive = (y1[self_selected].mean() - y0[~self_selected].mean())
print(f"\nself-selected comparison: {naive:+.3f}"
f" (bias {naive - truth['unit_effect'].mean():+.3f})")
print(f" mean engagement: adopters {engagement[self_selected].mean():.1f}"
f" vs others {engagement[~self_selected].mean():.1f} <- not comparable")
# --- RANDOMIZATION: assignment independent of every characteristic. ---
assigned = rng.random(n) < 0.5
observed = np.where(assigned, y1, y0) # only one outcome per unit
ate_hat = observed[assigned].mean() - observed[~assigned].mean()
se = np.sqrt(observed[assigned].var(ddof=1) / assigned.sum()
+ observed[~assigned].var(ddof=1) / (~assigned).sum())
print(f"\nrandomized estimate: {ate_hat:+.3f}"
f" 95% CI [{ate_hat-1.96*se:+.3f}, {ate_hat+1.96*se:+.3f}]")
print(f" mean engagement: treated {engagement[assigned].mean():.1f}"
f" vs control {engagement[~assigned].mean():.1f} <- balanced")
# --- IT BALANCES WHAT YOU NEVER MEASURED, which is the whole point. ---
unmeasured = rng.normal(50, 15, n) # never recorded anywhere
print(f"\nan UNMEASURED covariate, balanced anyway:"
f" {unmeasured[assigned].mean():.2f} vs {unmeasured[~assigned].mean():.2f}")
print("Matching can only balance columns you have. Randomization balances all.")
# --- BALANCE IS EXPECTED, NOT GUARANTEED: it depends on sample size. ---
print("\nworst standardised covariate imbalance across 500 randomizations:")
for size in (50, 500, 5000):
worst = []
for rep in range(500):
r = np.random.default_rng(rep)
idx = r.choice(n, size, replace=False)
a = r.random(size) < 0.5
if a.sum() < 5 or (~a).sum() < 5:
continue
e = engagement[idx]
# Standardised difference: scale-free, so comparable across variables.
d = ((e[a].mean() - e[~a].mean())
/ np.sqrt((e[a].var(ddof=1) + e[~a].var(ddof=1)) / 2))
worst.append(abs(d))
print(f" n={size:>5} median |std diff| {np.median(worst):.3f}"
f" 95th pct {np.percentile(worst, 95):.3f}"
f" share above 0.1: {np.mean(np.array(worst) > 0.1):.0%}")
# At n=50, a meaningful imbalance is the norm. Randomization is a
# large-sample guarantee, which is why small tests need stratification.
# --- A BALANCE CHECK you run before trusting any result. ---
def balance_table(df: pd.DataFrame, group: str, covariates: list[str]):
rows = []
for c in covariates:
a = df.loc[df[group], c]
b = df.loc[~df[group], c]
pooled = np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
rows.append({
"covariate": c,
"treated": round(a.mean(), 3), "control": round(b.mean(), 3),
"std_diff": round((a.mean() - b.mean()) / pooled, 3),
# A t-test here is a diagnostic, not a hypothesis worth testing:
# under correct randomization it is false-positive by design 5%
# of the time. Judge the standardised difference, not the p-value.
"p": round(stats.ttest_ind(a, b, equal_var=False).pvalue, 3),
"flag": abs((a.mean() - b.mean()) / pooled) > 0.1,
})
return pd.DataFrame(rows)
exp = truth.assign(assigned=assigned, observed=observed, unmeasured=unmeasured,
tenure=rng.integers(1, 60, n).astype(float))
print("\nbalance check:\n",
balance_table(exp, "assigned", ["engagement", "unmeasured", "tenure"])
.to_string(index=False))
# --- DIFFERENTIAL ATTRITION undoes the randomization after the fact. ---
# Treatment makes low-engagement users churn before the outcome is measured.
p_drop = np.where(assigned, 1/(1+np.exp((engagement-14)/3)) * 0.8, 0.05)
stayed = rng.random(n) > p_drop
biased = observed[assigned & stayed].mean() - observed[~assigned & stayed].mean()
print(f"\nafter differential attrition: {biased:+.3f}"
f" vs true {truth['unit_effect'].mean():+.3f}")
print(f" attrition rate: treated {(~stayed & assigned).sum()/assigned.sum():.1%}"
f" vs control {(~stayed & ~assigned).sum()/(~assigned).sum():.1%}")
print("Unequal attrition rates between arms is the alarm. Always report them.")
# --- COMPLIANCE: two different estimands, and only one is unbiased. ---
complied = assigned & (rng.random(n) < 0.7) # 30% never see the feature
partial = np.where(complied, y1, y0)
itt = partial[assigned].mean() - partial[~assigned].mean() # by ASSIGNMENT
per_protocol = partial[complied].mean() - partial[~assigned].mean() # by receipt
print(f"\nintention to treat (by assignment): {itt:+.3f} <- unbiased")
print(f"per protocol (by receipt): {per_protocol:+.3f} <- biased")
print("Per-protocol reintroduces selection, because compliers differ.")
print(f"Scaling ITT by the compliance rate gives the complier effect:"
f" {itt/0.7:+.3f}")
THE FUNDAMENTAL PROBLEM OF CAUSAL INFERENCE
unit outcome if treated outcome if control effect
1 42.1 ???? ????
2 ???? 31.8 ????
3 38.6 ???? ????
^ ^
observed for treated observed for control
Exactly one cell per row is ever knowable. Individual effects are
unidentifiable; the AVERAGE is what an experiment recovers.
WHAT EACH STRATEGY BALANCES
strategy measured covariates UNMEASURED covariates
---------------------------------------------------------------------
compare adopters no no
adjust / regress yes, if modelled NO
match on covariates yes, those matched NO
randomize yes, in expectation YES, in expectation
Only the last row balances what you never thought to record.
BALANCE IS A LARGE-SAMPLE PROPERTY
n=50 median |std diff| 0.19 -> imbalance is the NORM
n=500 median |std diff| 0.06
n=5000 median |std diff| 0.02 -> reliably balanced
Small experiments need stratified assignment; large ones rarely do.
ASSUMPTIONS THAT SURVIVE RANDOMIZATION
assumption fails when symptom
---------------------------------------------------------------------------
no interference social, marketplace, shared effect leaks into control
inventory
consistency the treatment changed mid-test effect is an average of
two different things
compliance users never reach the feature ITT is diluted
no differential treatment causes exits arm attrition rates differ
attrition
Key Takeaways
- Individual causal effects are unobservable, so an experiment estimates an average by constructing two equivalent groups.
- Randomization balances unmeasured characteristics in expectation, which no amount of statistical adjustment can achieve.
- Balance is a large-sample guarantee, so small experiments need stratified assignment and every experiment needs a balance check.
- Analyse by assignment rather than by receipt, because filtering to compliers reintroduces exactly the selection randomization removed.
- Report attrition by arm, since differential dropout undoes randomization after the fact.
🧪 Practice
- Simulate potential outcomes with a known average effect, then show that a self-selected comparison is biased while a randomized one is not.
- Run 500 randomizations at sample sizes of 50 and 5,000 and report how often the standardised covariate difference exceeds 0.1.
- Write a balance table function reporting standardised differences, and explain why a significant p-value in it is not necessarily a problem.
- Simulate a test where 30% of the treatment arm never receives the feature, and compare the intention-to-treat and per-protocol estimates.
- Interview: A colleague says their experiment is valid because the treatment and control groups match on age, region, and tenure. What is your response? (Hint: what did they not match on, and what would have handled it without a list?)
Control and Treatment Groups
The control group is the counterfactual made concrete, and getting it wrong is the most common way an otherwise well-run experiment produces a meaningless number. The principle is simple to state and easy to violate: the control group must differ from treatment in exactly one respect, the thing you are testing.
Several violations are common enough to name.
The control is not actually business as usual. If you rebuild a page to add a feature, and the control gets the old page while treatment gets the new page plus the feature, you are measuring the rebuild and the feature together and cannot separate them. The fix is to ship the rebuild to everyone first, or to run a three-arm test isolating each change.
Different measurement between arms. A new implementation often logs differently. If treatment fires a tracking event that control does not, the observed difference includes a measurement artefact. This is the most insidious variant because it looks exactly like a real effect.
Contamination. Users see both variants because assignment is by device and they switch devices, or because logged-out traffic is randomized per session, or because they share an account. Contamination biases the estimate toward zero, so it makes real effects harder to detect, which is at least a conservative failure.
Then the choices that make up a design.
Split ratio. A 50/50 split maximises statistical power for a fixed total sample, because the variance of the difference depends on the sum of the reciprocals of the group sizes, which is minimised when they are equal. Unequal splits are sometimes justified by risk, since a 99/1 split exposes fewer users to a possibly harmful change, but the cost in power is real and worth quantifying rather than assuming.
Ramping. Production practice is to start small, verify nothing is broken, then increase exposure. The important subtlety is that a ramp changes the population composition over time, so you should analyse each ramp stage separately or restrict the primary analysis to a period of stable allocation. Pooling across a ramp mixes different exposure regimes.
A/A tests. Running an experiment where both arms get the identical experience is the single most valuable piece of infrastructure validation available. Since the true effect is zero, any consistent difference means the system is broken: the randomization is not random, the assignment leaks, the metric pipeline differs by arm, or the variance estimate is wrong. Run them continuously, and check that the distribution of p-values across many A/A tests is uniform, because that is what the theory predicts and a skewed distribution localises the defect.
Holdouts. A long-running control group that never receives any new features measures the cumulative effect of everything shipped over a quarter, which is usually much smaller than the sum of the individually measured wins. That gap is informative rather than embarrassing: it captures interaction between features, novelty decay, and the fact that a metric cannot improve without limit.
import numpy as np
import pandas as pd
from scipy import stats
rng = np.random.default_rng(1)
# --- SPLIT RATIO AND POWER: quantify the cost of an unequal split. ---
# se(difference) ~ sqrt(var/n_a + var/n_b), minimised at n_a == n_b.
N, sigma, effect = 20_000, 10.0, 0.35
print("relative standard error by allocation (total N fixed at 20,000):")
base = None
for share in (0.50, 0.40, 0.25, 0.10, 0.02):
na, nb = int(N * share), N - int(N * share)
se = sigma * np.sqrt(1/na + 1/nb)
base = base or se
# Power for a two-sided test at alpha = 0.05.
z = effect / se
power = stats.norm.cdf(z - 1.959964) + stats.norm.cdf(-z - 1.959964)
print(f" {share:>5.0%}/{1-share:<5.0%} se {se:.4f}"
f" se vs balanced {se/base:5.2f}x power {power:.3f}")
print("\nA 90/10 split needs ~1.7x the users of a 50/50 for the same power.")
print("Use an unequal split for RISK control, and price the power loss.")
# --- A/A TEST: the infrastructure check that finds what nothing else does. ---
def run_aa(n: int, seed: int, broken: bool = False) -> float:
r = np.random.default_rng(seed)
y = r.gamma(2, 20, n)
if broken:
# A realistic defect: the hash is correlated with a user attribute,
# so assignment is not independent of the outcome.
arm = (y > np.median(y)) ^ (r.random(n) < 0.35)
else:
arm = r.random(n) < 0.5
return stats.ttest_ind(y[arm], y[~arm], equal_var=False).pvalue
for label, broken in [("healthy system", False), ("broken assignment", True)]:
pvals = np.array([run_aa(6000, s, broken) for s in range(400)])
# Under a correct system p-values are UNIFORM, so ~5% fall below 0.05.
ks = stats.kstest(pvals, "uniform").pvalue
print(f"\n{label}:")
print(f" share of A/A tests with p < 0.05: {np.mean(pvals < 0.05):.1%}"
f" (expected 5%)")
print(f" uniformity of the p-value distribution: KS p = {ks:.2e}"
f" {'OK' if ks > 0.01 else 'FAIL: assignment is not random'}")
# --- SAMPLE RATIO MISMATCH: the first thing to check on any test. ---
def srm_check(n_a: int, n_b: int, expected_share: float = 0.5) -> dict:
"""A chi-squared test that the observed split matches the intended one.
This is the highest-value single check in experiment analysis."""
total = n_a + n_b
exp_a = total * expected_share
chi2 = ((n_a - exp_a) ** 2 / exp_a
+ (n_b - (total - exp_a)) ** 2 / (total - exp_a))
p = 1 - stats.chi2.cdf(chi2, df=1)
return {"n_a": n_a, "n_b": n_b, "observed_share": round(n_a / total, 4),
"chi2": round(chi2, 2), "p_value": f"{p:.2e}",
"verdict": "INVALID: do not interpret" if p < 0.001 else "ok"}
print("\nsample ratio mismatch checks:")
for a, b in [(50_000, 50_100), (50_000, 51_800), (50_000, 47_000)]:
print(" ", srm_check(a, b))
print("A mismatch means users were lost or misassigned NON-randomly, so the")
print("arms are no longer comparable and the effect estimate is unusable.")
# --- CONTAMINATION biases toward zero: a conservative but real distortion. ---
n = 30_000
truth = 1.8
assigned = rng.random(n) < 0.5
for leak in (0.0, 0.05, 0.20, 0.40):
# A fraction of each arm actually receives the other arm's experience.
crossed = rng.random(n) < leak
received = assigned ^ crossed
y = 40 + truth * received + rng.normal(0, 9, n)
est = y[assigned].mean() - y[~assigned].mean()
print(f"contamination {leak:>4.0%} estimate {est:+.3f}"
f" ({est/truth:.0%} of the true {truth})")
print("Dilution scales with the leak rate. It hides effects; it does not")
print("invent them, which makes it the least dangerous failure mode.")
# --- THE CONFOUNDED CONTROL: two changes, one comparison. ---
n = 24_000
g = rng.integers(0, 3, n) # 0 = old page, 1 = new, 2 = new+feature
redesign_effect, feature_effect = 1.2, 0.8
y = (30 + redesign_effect * (g >= 1) + feature_effect * (g == 2)
+ rng.normal(0, 7, n))
two_arm = y[g == 2].mean() - y[g == 0].mean() # confounded
three_arm = y[g == 2].mean() - y[g == 1].mean() # isolates the feature
print(f"\nold page vs new+feature: {two_arm:+.3f}"
f" <- redesign AND feature together")
print(f"new page vs new+feature: {three_arm:+.3f}"
f" <- the feature alone (true {feature_effect})")
print("A third arm is how you separate two changes that shipped together.")
# --- RAMPING: pooling across stages mixes different allocation regimes. ---
stages = pd.DataFrame([
("day 1-2", 0.01, 4_000, 1.9),
("day 3-5", 0.10, 40_000, 1.8),
("day 6-14", 0.50, 300_000, 0.6), # the real effect, at stable allocation
], columns=["window", "treat_share", "users", "observed_lift_pct"])
pooled = np.average(stages["observed_lift_pct"], weights=stages["users"])
print("\nramp stages:\n", stages.to_string(index=False))
print(f"\nvolume-weighted pooled lift: {pooled:.2f}%")
print(f"stable-allocation lift only: "
f"{stages.iloc[-1]['observed_lift_pct']:.2f}%")
print("Early ramp stages reach an unrepresentative slice and often show")
print("novelty-inflated effects. Analyse the stable window, not the pool.")
# --- HOLDOUT: the sum of wins exceeds the measured total, predictably. ---
shipped = pd.DataFrame({
"feature": [f"f{i}" for i in range(1, 7)],
"measured_lift_pct": [1.4, 0.9, 2.1, 0.6, 1.1, 0.8],
})
print(f"\nsum of individually measured lifts: "
f"{shipped['measured_lift_pct'].sum():.1f}%")
print("quarterly holdout measures: 3.8%")
print("The gap is real: feature interaction, novelty decay, and metric")
print("ceilings. A holdout is the only way to see it.")
THE ONE-DIFFERENCE RULE
VALID INVALID
control: old page control: old page
treatment: old page + feature treatment: NEW page + feature
^ one difference ^ two differences
Fix: ship the rebuild to everyone first, or add a third arm:
A = old page | B = new page | C = new page + feature
B vs C isolates the feature. A vs B isolates the rebuild.
SPLIT RATIO VS POWER (total N fixed)
allocation se vs balanced power at the same effect
50/50 1.00x 0.85
40/60 1.02x 0.84
25/75 1.15x 0.76
10/90 1.67x 0.45
2/98 3.57x 0.13
Unequal splits are a RISK decision that costs power. Price it.
WHAT AN A/A TEST FINDS THAT NOTHING ELSE DOES
true effect is zero, so any consistent difference is a DEFECT:
- assignment hash correlated with a user attribute
- assignment leaking across arms
- metric pipeline differing by arm
- variance underestimated (p-values not uniform)
healthy: 5% of A/A tests show p < 0.05, p-values UNIFORM
broken: far more, and the p-value distribution is skewed
CHECK ORDER ON ANY TEST RESULT
1. sample ratio mismatch -> if it fails, STOP. Nothing else matters.
2. arm attrition rates -> differential dropout?
3. covariate balance -> std diffs within 0.1?
4. metric pipeline parity -> same events logged in both arms?
5. only then look at the effect
Key Takeaways
- Control and treatment must differ in exactly one respect, so a change bundled with a redesign needs a third arm to be separable.
- A 50/50 split maximises power for a fixed sample, and unequal splits are a risk decision whose power cost should be quantified.
- A/A tests validate the experiment system itself, and the p-values across many of them should be uniform.
- Sample ratio mismatch is the first check on any result, and a failure means the effect estimate is unusable rather than merely noisy.
- Contamination dilutes effects toward zero, and a long-running holdout reveals that shipped wins do not sum.
🧪 Practice
- Compute the standard error of the difference for 50/50, 25/75, and 10/90 splits at fixed total sample size, and convert each into power.
- Simulate 400 A/A tests on a healthy system and confirm about 5% reach significance and the p-values are uniform.
- Write a sample ratio mismatch check and run it on splits of 50,000 against 50,100 and 50,000 against 51,800.
- Simulate 20% contamination between arms and report how much of the true effect survives.
- Interview: Your test shows a 3% lift and the treatment arm has 2% more users than control. What do you do? (Hint: what would cause users to be unevenly assigned, and is the surplus likely to be a random subset?)
Randomization Units
The choice of what to randomize is the design decision with the largest consequences and the least attention. Randomizing users, sessions, requests, or groups of users all produce valid experiments, and they answer different questions with different statistical properties.
The two forces pulling in opposite directions are consistency and precision.
Randomizing at a fine unit, such as each page request, gives you enormous sample size and tight confidence intervals. It also gives a single user a different experience on every page load, which is confusing, invalidates any metric that spans more than one request, and makes measuring a change in user behaviour impossible, because the user was never consistently in one condition.
Randomizing at a coarse unit, such as a whole city or company account, gives a consistent experience and handles spillover between people. It also collapses your effective sample size, because the unit of analysis becomes the city, not the person, and twenty cities is twenty observations regardless of how many users they contain.
The rule that resolves most cases: randomize at the level at which the treatment acts and at which you want to measure the outcome, and no finer. If the metric is "did this user convert this week", the unit is the user. If the metric is "did this request return a relevant result", the request may be acceptable.
The second consideration is interference. When one unit's treatment affects another's outcome, individual randomization is invalid because control units are partially treated. The standard cases:
Social features. If treatment users post more, their control friends see more content. The control group is contaminated, and the measured effect understates the true one.
Marketplaces. If treated buyers bid more aggressively, control buyers lose auctions they would have won. Here the effect on control is negative, so the difference overstates the true effect, sometimes by a lot. This is the most dangerous version because the bias inflates rather than deflates.
Shared resources. A treatment that consumes more server capacity slows down control users.
The fix is cluster randomization: assign whole groups so that interference happens mostly within a group rather than across arms. Clusters can be geographic, social communities, or time periods. The cost is precision, and the mechanism is worth knowing by name: the intraclass correlation. When outcomes within a cluster are correlated, each additional member of a cluster contributes less independent information, and the effective sample size shrinks by a factor called the design effect, approximately one plus the product of the intraclass correlation and the cluster size minus one. With clusters of 100 and an intraclass correlation of only 0.05, the design effect is about 6, so your 100,000 users carry the information of roughly 17,000.
Switchback designs, alternating treatment across time periods for everyone, are the standard answer for marketplaces where geographic clusters are too few. They handle interference and control for time trends, at the cost of assuming no carryover between periods.
import numpy as np
import pandas as pd
from scipy import stats
rng = np.random.default_rng(2)
# --- THE SAME EXPERIMENT AT THREE UNITS OF RANDOMIZATION ---
n_users, sessions_per_user = 2_000, 6
user_id = np.repeat(np.arange(n_users), sessions_per_user)
# Each user has a persistent tendency: this is what creates within-user
# correlation, and it is why the unit of analysis matters.
user_effect = np.repeat(rng.normal(0, 8, n_users), sessions_per_user)
true_lift = 1.0
rows = len(user_id)
# 1. RANDOMIZE BY SESSION: huge n, inconsistent experience.
by_session = rng.random(rows) < 0.5
# 2. RANDOMIZE BY USER: consistent experience, correct unit.
user_arm = rng.random(n_users) < 0.5
by_user = np.repeat(user_arm, sessions_per_user)
outcome_session = 30 + true_lift * by_session + user_effect + rng.normal(0, 5, rows)
outcome_user = 30 + true_lift * by_user + user_effect + rng.normal(0, 5, rows)
# --- THE VARIANCE TRAP: analysing user-randomized data at session level. ---
def naive_test(y, arm, label):
a, b = y[arm], y[~arm]
se = np.sqrt(a.var(ddof=1)/len(a) + b.var(ddof=1)/len(b))
diff = a.mean() - b.mean()
return {"analysis": label, "n_rows": len(y), "estimate": round(diff, 3),
"se": round(se, 4), "ci_width": round(2*1.96*se, 4)}
results = [
naive_test(outcome_session, by_session, "session-randomized, session-analysed"),
naive_test(outcome_user, by_user, "user-randomized, SESSION-analysed (WRONG)"),
]
# The correct analysis for user-randomized data: aggregate to the user first.
per_user = pd.DataFrame({"user": user_id, "y": outcome_user,
"arm": by_user}).groupby("user").agg(
y=("y", "mean"), arm=("arm", "first"))
a, b = per_user.loc[per_user["arm"], "y"], per_user.loc[~per_user["arm"], "y"]
se_correct = np.sqrt(a.var(ddof=1)/len(a) + b.var(ddof=1)/len(b))
results.append({"analysis": "user-randomized, user-analysed (RIGHT)",
"n_rows": len(per_user),
"estimate": round(a.mean()-b.mean(), 3),
"se": round(se_correct, 4),
"ci_width": round(2*1.96*se_correct, 4)})
print(pd.DataFrame(results).to_string(index=False))
print(f"\nThe wrong analysis understates the standard error by "
f"{se_correct/results[1]['se']:.1f}x, because 12,000 sessions from")
print("2,000 users are not 12,000 independent observations.")
# --- INTRACLASS CORRELATION AND THE DESIGN EFFECT ---
def icc(values: np.ndarray, cluster: np.ndarray) -> float:
"""Share of total variance that lies BETWEEN clusters."""
df = pd.DataFrame({"y": values, "g": cluster})
grand = df["y"].mean()
means = df.groupby("g")["y"].mean()
sizes = df.groupby("g")["y"].size()
between = (sizes * (means - grand) ** 2).sum() / (len(means) - 1)
within = df.groupby("g")["y"].var(ddof=1).mean()
m_bar = sizes.mean()
return float(max(between - within, 0)
/ max(between + (m_bar - 1) * within, 1e-12))
rho = icc(outcome_user, user_id)
print(f"\nintraclass correlation within users: {rho:.3f}")
print("\ndesign effect = 1 + (m - 1) * rho, and effective n = n / design effect:")
for m in (2, 10, 50, 200):
for r in (0.01, 0.05, 0.20):
de = 1 + (m - 1) * r
print(f" cluster size {m:>3}, rho {r:.2f} -> design effect {de:5.2f}"
f" 100,000 users carry {100_000/de:>8,.0f} units of information")
print()
# --- INTERFERENCE: the two directions of bias. ---
n = 20_000
arm = rng.random(n) < 0.5
direct = 1.5
# SOCIAL SPILLOVER: treatment leaks POSITIVELY into control -> understates.
spill_to_control = 0.6 * arm.mean() * 1.5
y_social = 40 + direct * arm + spill_to_control * ~arm + rng.normal(0, 8, n)
est_social = y_social[arm].mean() - y_social[~arm].mean()
# MARKETPLACE CANNIBALISATION: treatment HURTS control -> overstates.
harm_to_control = -0.9
y_market = 40 + direct * arm + harm_to_control * ~arm + rng.normal(0, 8, n)
est_market = y_market[arm].mean() - y_market[~arm].mean()
print(f"true direct effect: {direct:+.2f}")
print(f"with social spillover: {est_social:+.2f}"
f" ({est_social/direct:.0%} of truth) -- UNDERSTATES")
print(f"with marketplace crowding: {est_market:+.2f}"
f" ({est_market/direct:.0%} of truth) -- OVERSTATES")
print("\nUnderstating is conservative. Overstating ships a change that does")
print("not work at 100% rollout, which is how marketplace tests fail.")
# --- CLUSTER RANDOMIZATION: correct standard errors, far fewer of them. ---
n_clusters, per_cluster = 40, 500
cluster = np.repeat(np.arange(n_clusters), per_cluster)
cluster_effect = np.repeat(rng.normal(0, 6, n_clusters), per_cluster)
cl_arm_c = rng.random(n_clusters) < 0.5
cl_arm = np.repeat(cl_arm_c, per_cluster)
y_cl = 40 + 1.5 * cl_arm + cluster_effect + rng.normal(0, 9, len(cluster))
# WRONG: treat each user as independent.
naive_se = np.sqrt(y_cl[cl_arm].var(ddof=1)/cl_arm.sum()
+ y_cl[~cl_arm].var(ddof=1)/(~cl_arm).sum())
# RIGHT: the cluster is the unit. Aggregate, then test on cluster means.
cm = pd.DataFrame({"c": cluster, "y": y_cl, "arm": cl_arm}).groupby("c").agg(
y=("y", "mean"), arm=("arm", "first"))
ca, cb = cm.loc[cm["arm"], "y"], cm.loc[~cm["arm"], "y"]
cluster_se = np.sqrt(ca.var(ddof=1)/len(ca) + cb.var(ddof=1)/len(cb))
print(f"\n{n_clusters} clusters x {per_cluster} users = {len(cluster):,} users")
print(f" naive per-user se: {naive_se:.4f} (p = "
f"{stats.ttest_ind(y_cl[cl_arm], y_cl[~cl_arm], equal_var=False).pvalue:.2e})")
print(f" cluster-level se: {cluster_se:.4f} (p = "
f"{stats.ttest_ind(ca, cb, equal_var=False).pvalue:.3f})")
print(f" the naive se is {cluster_se/naive_se:.1f}x too small")
print("20,000 users randomized in 40 clusters give you 40 observations.")
# --- SWITCHBACK: alternate in time when clusters are too few. ---
periods = 60
period_arm = np.tile([True, False], periods // 2) # alternate by period
time_trend = np.sin(np.arange(periods) / 7 * 2 * np.pi) * 3
period_y = 50 + 1.5 * period_arm + time_trend + rng.normal(0, 1.5, periods)
sw = pd.DataFrame({"period": range(periods), "arm": period_arm, "y": period_y})
# Pairing adjacent periods differences out the shared time trend.
paired = sw["y"].values.reshape(-1, 2)
diffs = paired[:, 0] - paired[:, 1]
t = stats.ttest_1samp(diffs, 0)
print(f"\nswitchback: {periods} alternating periods")
print(f" paired estimate {diffs.mean():+.3f}"
f" 95% CI [{diffs.mean()-1.96*diffs.std(ddof=1)/np.sqrt(len(diffs)):+.3f},"
f" {diffs.mean()+1.96*diffs.std(ddof=1)/np.sqrt(len(diffs)):+.3f}]"
f" p {t.pvalue:.4f}")
print("Alternating handles interference AND differences out the time trend.")
print("It assumes no carryover between adjacent periods.")
THE UNIT LADDER
unit sample size consistent? interference right when
---------------------------------------------------------------------------
request largest NO unhandled stateless ranking
session large within a unhandled session-scoped
session metrics
USER moderate YES unhandled the default
account smaller yes within account B2B, shared plans
cluster/geo small yes HANDLED marketplaces, social
time period smallest yes handled few clusters exist
RULE: randomize at the level the treatment acts and the outcome is
measured, and no finer.
THE VARIANCE TRAP
2,000 users x 6 sessions = 12,000 rows
analysed as 12,000 independent sessions se = 0.0921 p = 0.0004
analysed correctly at the user level se = 0.2510 p = 0.31
The first is wrong by 2.7x and turns noise into a "significant" result.
Randomize at unit X -> analyse at unit X.
DESIGN EFFECT = 1 + (m - 1) * rho
cluster size m rho=0.01 rho=0.05 rho=0.20
2 1.01 1.05 1.20
10 1.09 1.45 2.80
50 1.49 3.45 10.80
200 2.99 10.95 40.80
effective n = n / design effect. With m=200 and rho=0.05, 100,000
users carry the information of about 9,100.
INTERFERENCE: WHICH DIRECTION IS THE BIAS?
social spillover treatment helps control too
-> control improves -> difference SHRINKS -> understates (safe)
marketplace crowding treatment takes from control
-> control worsens -> difference GROWS -> OVERSTATES (dangerous)
A marketplace test that looks great at 50/50 can be worthless at 100%,
because there is no control left to take from.
Key Takeaways
- Randomize at the level the treatment acts and the outcome is measured, and analyse at the same level you randomized.
- Analysing user-randomized data as independent sessions understates the standard error and manufactures significance.
- The design effect, one plus the intraclass correlation times cluster size minus one, is how much cluster randomization costs in effective sample size.
- Social spillover understates effects, which is conservative, while marketplace crowding overstates them, which ships changes that fail at full rollout.
- Switchback designs alternate treatment over time and difference out the time trend when too few clusters exist.
🧪 Practice
- Simulate user-randomized data with multiple sessions per user, then compare the standard error from a session-level analysis against a user-level one.
- Write an intraclass correlation function and compute the design effect for cluster sizes of 10, 50, and 200.
- Simulate both positive social spillover and negative marketplace crowding, and report the direction and size of the bias in each.
- Build a switchback experiment over 60 alternating periods with a strong time trend, and show that pairing adjacent periods removes it.
- Interview: You are testing a driver-incentive change on a ride-hailing marketplace. What is your randomization unit and why? (Hint: what happens to control drivers when treated drivers accept more rides in the same city?)
Blocking and Stratification
Randomization gives you balance in expectation. Blocking and stratification give you balance by construction, and the payoff is not only fairness but precision: a well-chosen stratification can cut the required sample size substantially at no extra cost.
The intuition is about where variance comes from. If your outcome varies a lot across countries, and your two arms happen to contain slightly different country mixes, some of the observed difference is country mix rather than treatment. That is noise. Remove the country mix as a source of variation and the remaining comparison is cleaner.
Stratified randomization does this by partitioning units into strata, such as country crossed with platform, and randomizing within each stratum separately. Every stratum ends up with the intended split, so the mix is identical by construction. The analysis should then account for the strata, by computing a stratum-weighted average effect or including strata in a regression, because a stratified design analysed as if it were simple randomization leaves the precision gain on the table.
Blocking is the same idea from the experimental-design tradition, where a block is a set of similar units and treatments are assigned within the block. A paired design is the extreme case with blocks of size two, and it is why pairing adjacent time periods in a switchback works so well.
The industrial version of this idea is CUPED, controlled experiment using pre-experiment data, and it deserves particular attention because it is one of the highest-return techniques in practice. The observation is that the best predictor of a user's metric during the experiment is usually the same user's metric before the experiment. Pre-period data cannot have been affected by the treatment, so you can safely subtract off the part of the outcome it explains. The adjusted metric has the same expected difference between arms and a smaller variance, and the variance reduction is the square of the correlation between the pre-period and in-experiment values. A correlation of 0.7 removes about half the variance, which is equivalent to doubling your sample size for free.
Two cautions. The adjustment must use only pre-experiment data, because using anything measured during the experiment can absorb part of the treatment effect and bias the estimate. And stratifying on too many variables creates strata with very few units, which reintroduces imbalance and complicates the analysis; a handful of high-variance dimensions is the sweet spot.
Post-stratification, adjusting after the fact rather than randomizing within strata, recovers most of the same precision gain and is what you use when you did not stratify up front. It is nearly always worth doing.
import numpy as np
import pandas as pd
from scipy import stats
rng = np.random.default_rng(3)
n = 6_000
# A population where an observable dimension drives most of the variance.
country = rng.choice(["US", "DE", "JP", "BR"], n, p=[.45, .25, .2, .1])
base = pd.Series(country).map({"US": 60, "DE": 45, "JP": 80, "BR": 25}).values
true_effect = 2.0
# --- SIMPLE VS STRATIFIED RANDOMIZATION ---
simple_arm = rng.random(n) < 0.5
strat_arm = np.zeros(n, dtype=bool)
for c in np.unique(country): # randomize WITHIN each stratum
idx = np.where(country == c)[0]
take = rng.permutation(idx)[: len(idx) // 2]
strat_arm[take] = True
print("treated share within each country:")
comp = pd.DataFrame({
"simple": pd.Series(simple_arm).groupby(country).mean().round(3),
"stratified": pd.Series(strat_arm).groupby(country).mean().round(3),
})
print(comp.to_string())
print("\nStratified hits 0.5 in every stratum by construction; simple")
print("randomization is close on average and off in any single stratum.")
# --- THE PRECISION GAIN, measured over many replications. ---
def one_run(seed: int, stratify: bool) -> float:
r = np.random.default_rng(seed)
c = r.choice(["US", "DE", "JP", "BR"], n, p=[.45, .25, .2, .1])
b = pd.Series(c).map({"US": 60, "DE": 45, "JP": 80, "BR": 25}).values
if stratify:
arm = np.zeros(n, dtype=bool)
for k in np.unique(c):
idx = np.where(c == k)[0]
arm[r.permutation(idx)[: len(idx)//2]] = True
else:
arm = r.random(n) < 0.5
y = b + true_effect * arm + r.normal(0, 10, n)
return y[arm].mean() - y[~arm].mean()
for label, strat in [("simple", False), ("stratified", True)]:
ests = np.array([one_run(s, strat) for s in range(400)])
print(f"\n{label:11s} mean {ests.mean():+.3f} sd of estimate {ests.std():.4f}")
simple_sd = np.array([one_run(s, False) for s in range(400)]).std()
strat_sd = np.array([one_run(s, True) for s in range(400)]).std()
print(f"\nvariance ratio {(strat_sd/simple_sd)**2:.3f}"
f" -> equivalent to {(simple_sd/strat_sd)**2:.2f}x the sample size")
# --- POST-STRATIFICATION: recover the gain even without stratified assignment. ---
y = base + true_effect * simple_arm + rng.normal(0, 10, n)
df = pd.DataFrame({"country": country, "arm": simple_arm, "y": y})
unadj = y[simple_arm].mean() - y[~simple_arm].mean()
se_unadj = np.sqrt(y[simple_arm].var(ddof=1)/simple_arm.sum()
+ y[~simple_arm].var(ddof=1)/(~simple_arm).sum())
# Stratum-weighted estimate: effect per stratum, weighted by stratum share.
parts = []
for c, g in df.groupby("country"):
a, b_ = g.loc[g["arm"], "y"], g.loc[~g["arm"], "y"]
parts.append({"stratum": c, "w": len(g)/len(df),
"effect": a.mean() - b_.mean(),
"var": a.var(ddof=1)/len(a) + b_.var(ddof=1)/len(b_)})
p = pd.DataFrame(parts)
adj = (p["w"] * p["effect"]).sum()
se_adj = np.sqrt((p["w"]**2 * p["var"]).sum())
print(f"\nunadjusted {unadj:+.3f} se {se_unadj:.4f}")
print(f"post-stratified {adj:+.3f} se {se_adj:.4f}"
f" ({1-(se_adj/se_unadj):.0%} narrower)")
# --- CUPED: the highest-return variance reduction in practice. ---
m = 20_000
# Pre-period metric, and an in-experiment metric strongly related to it.
pre = rng.gamma(3, 15, m)
arm = rng.random(m) < 0.5
lift = 1.5
post = 0.75 * pre + 12 + lift * arm + rng.normal(0, 12, m)
def cuped(post: np.ndarray, pre: np.ndarray, arm: np.ndarray):
"""Subtract off the part of the outcome explained by PRE-period data.
theta is estimated from the pooled sample; pre cannot be affected by
treatment, so the expected difference is unchanged."""
theta = np.cov(post, pre, ddof=1)[0, 1] / pre.var(ddof=1)
adjusted = post - theta * (pre - pre.mean())
a, b = adjusted[arm], adjusted[~arm]
se = np.sqrt(a.var(ddof=1)/len(a) + b.var(ddof=1)/len(b))
return a.mean() - b.mean(), se, theta, adjusted
raw_est = post[arm].mean() - post[~arm].mean()
raw_se = np.sqrt(post[arm].var(ddof=1)/arm.sum()
+ post[~arm].var(ddof=1)/(~arm).sum())
c_est, c_se, theta, adjusted = cuped(post, pre, arm)
corr = np.corrcoef(post, pre)[0, 1]
print(f"\ncorrelation(pre, post) = {corr:.3f}")
print(f"raw estimate {raw_est:+.4f} se {raw_se:.4f}")
print(f"CUPED estimate {c_est:+.4f} se {c_se:.4f} (theta {theta:.3f})")
print(f"variance reduction: {1-(c_se/raw_se)**2:.1%}"
f" (theory predicts rho^2 = {corr**2:.1%})")
print(f"equivalent sample size multiplier: {(raw_se/c_se)**2:.2f}x")
# The whole point, stated in the currency that matters: users needed.
print("\nusers needed for 80% power at a 1% relative lift:")
for label, sd in [("raw metric", post.std(ddof=1)),
("CUPED metric", adjusted.std(ddof=1))]:
mde_abs = 0.01 * post.mean()
need = 2 * (1.959964 + 0.841621)**2 * sd**2 / mde_abs**2
print(f" {label:14s} sd {sd:6.2f} -> {need:>10,.0f} per arm")
# --- THE CUPED MISTAKE: using data from DURING the experiment. ---
during = post * 0.5 + rng.normal(0, 3, m) # contaminated by treatment
theta_bad = np.cov(post, during, ddof=1)[0, 1] / during.var(ddof=1)
bad = post - theta_bad * (during - during.mean())
bad_est = bad[arm].mean() - bad[~arm].mean()
print(f"\nadjusting on an IN-EXPERIMENT covariate: {bad_est:+.4f}"
f" vs true {lift:+.2f}")
print("It absorbed part of the treatment effect. Only pre-period data is safe.")
# --- TOO MANY STRATA: the failure mode of over-stratifying. ---
print("\nstrata count vs units per stratum (n = 6,000):")
dims = [("country", 4), ("platform", 3), ("tenure band", 5), ("plan", 4)]
cells = 1
for name, k in dims:
cells *= k
print(f" + {name:12s} -> {cells:>4} strata, ~{6000/cells:>7.1f} units each"
f" {'ok' if 6000/cells >= 30 else 'TOO FEW: imbalance returns'}")
WHERE THE PRECISION GAIN COMES FROM
total variance = BETWEEN-strata variance + WITHIN-strata variance
|
stratification removes this part from the comparison
outcome by country: JP 80 | US 60 | DE 45 | BR 25
A slight country-mix difference between arms looks like a treatment
effect. Fix the mix and that noise disappears.
THE TOOLKIT, BY WHEN YOU CAN APPLY IT
technique applied needs typical gain
---------------------------------------------------------------------------
stratified assign. before strata known up front 10-40% variance
blocking / pairing before similar units large, if tight
post-stratification after strata recorded most of the above
CUPED after PRE-period metric 30-70% variance
regression adjust. after pre-period covariates similar to CUPED
CUPED IN ONE LINE
adjusted = post - theta * (pre - mean(pre)), theta = cov(post,pre)/var(pre)
variance reduction = rho^2, where rho = corr(pre, post)
rho variance removed equivalent sample multiplier
0.3 9% 1.10x
0.5 25% 1.33x
0.7 49% 1.96x
0.9 81% 5.26x
A metric with rho = 0.7 needs HALF the users. This is free power.
THE RULE THAT KEEPS IT VALID
the covariate must be measured BEFORE assignment
pre-period spend SAFE
pre-period visit count SAFE
in-experiment sessions UNSAFE: absorbs the treatment effect
post-period anything UNSAFE
OVER-STRATIFYING
+ country (4) 4 strata, 1500 each ok
+ platform (3) 12 strata, 500 each ok
+ tenure band (5) 60 strata, 100 each ok
+ plan (4) 240 strata, 25 each TOO FEW
Use a few high-variance dimensions, not every dimension you have.
Key Takeaways
- Stratified randomization balances the mix by construction and removes between-strata variance from the comparison, which raises precision.
- Analyse a stratified design with stratum weights or a regression, or the precision gain is wasted.
- CUPED subtracts the part of the outcome explained by a pre-period metric, reducing variance by the square of their correlation.
- A pre-post correlation of 0.7 halves the required sample size, which makes CUPED one of the highest-return techniques available.
- The adjustment covariate must predate assignment, and stratifying on too many dimensions reintroduces the imbalance it was meant to remove.
🧪 Practice
- Compare simple and stratified randomization over 400 replications and report the ratio of the estimator variances.
- Implement post-stratification as a stratum-weighted effect with its standard error, and compare against the unadjusted difference.
- Implement CUPED and verify empirically that the variance reduction matches the square of the pre-post correlation.
- Show that adjusting on a covariate measured during the experiment biases the estimate toward zero.
- Interview: Your team wants to detect a 0.5% lift but the test would need twelve weeks. What do you propose before asking for more traffic? (Hint: what do you already know about these users from before the test started?)
Factorial Designs
The intuitive way to test several changes is one at a time: run A, then run B, then run C. This feels rigorous and it is both slower and less informative than the alternative, which is a genuinely counterintuitive result worth understanding.
A factorial design varies every factor simultaneously across all combinations. With three binary factors there are eight combinations, and each user is assigned to one. The surprise is that this does not cost you eight times the sample size per effect. To estimate the effect of factor A you compare all four cells where A is on against all four where A is off, which uses every observation in the experiment. The same data simultaneously estimates the effect of B and the effect of C. You get three effects for the price of one test, which is why the technique is called hidden replication.
One-factor-at-a-time testing, by contrast, splits your traffic across three sequential tests, so each effect is estimated on a third of the data, and then it still cannot answer the most important question.
That question is interaction. Two changes may not combine additively. A new recommendation module might help a lot on its own and add nothing once a better search bar is present, because both address the same need. Sequential testing measures each in isolation and then assumes the sum, which is exactly the assumption most likely to be wrong when two features address related problems. A factorial design measures the interaction directly, because it contains cells with both changes present.
The cost is cell count, which grows as two to the power of the number of factors. Five factors is 32 cells, seven is 128, and cells get thin fast. Two responses:
Fractional factorial designs run a carefully chosen subset of cells, giving up the ability to distinguish certain high-order interactions in exchange for a fraction of the size. A half-fraction of five factors runs 16 cells instead of 32. The trade is formalised as resolution: a resolution IV design keeps main effects clean of two-factor interactions but confounds two-factor interactions with each other, which is usually acceptable because high-order interactions are typically small.
Orthogonal arrays and their industrial packaging as Taguchi designs generalise this to factors with more than two levels.
In online experimentation, the common practice is subtly different and worth naming: many companies run concurrent independent experiments, each randomized separately, on overlapping traffic. This is effectively a factorial design where the assignment is independent rather than balanced, which works for main effects and means interactions are not measured at all. Most of the time that is fine. When two experiments touch the same surface, it is not, and the usual mitigation is to make them mutually exclusive by putting them in the same layer.
import numpy as np
import pandas as pd
import itertools
from scipy import stats
import statsmodels.api as sm
rng = np.random.default_rng(4)
# --- A WORLD WITH A REAL INTERACTION ---
# Search and recommendations both help discovery, so together they add
# LESS than the sum of their parts. This is exactly what OFAT cannot see.
TRUE = {"intercept": 40.0, "A_search": 2.0, "B_recs": 1.6,
"C_layout": 0.5, "AB_interaction": -1.1}
def outcome(a, b, c, size, r):
return (TRUE["intercept"] + TRUE["A_search"]*a + TRUE["B_recs"]*b
+ TRUE["C_layout"]*c + TRUE["AB_interaction"]*(a & b)
+ r.normal(0, 10, size))
# --- FULL FACTORIAL: 2^3 = 8 cells, every combination. ---
N = 24_000
cells = list(itertools.product([0, 1], repeat=3))
per_cell = N // len(cells)
rows = []
for a, b, c in cells:
y = outcome(bool(a), bool(b), bool(c), per_cell, rng)
rows.append(pd.DataFrame({"A": a, "B": b, "C": c, "y": y}))
fac = pd.concat(rows, ignore_index=True)
print(f"full factorial: {len(cells)} cells x {per_cell:,} = {len(fac):,} users")
print("\ncell means:")
print(fac.groupby(["A", "B", "C"])["y"].agg(n="size", mean="mean")
.round(2).to_string())
# EVERY observation contributes to EVERY main effect.
print("\nmain effects, each using all 24,000 users:")
for f in ("A", "B", "C"):
on, off = fac.loc[fac[f] == 1, "y"], fac.loc[fac[f] == 0, "y"]
se = np.sqrt(on.var(ddof=1)/len(on) + off.var(ddof=1)/len(off))
print(f" {f}: {on.mean()-off.mean():+.3f} se {se:.3f}"
f" (n_on {len(on):,}, n_off {len(off):,})")
# --- THE INTERACTION, which only a factorial design can measure. ---
m = fac.groupby(["A", "B"])["y"].mean()
interaction = (m[(1, 1)] - m[(1, 0)]) - (m[(0, 1)] - m[(0, 0)])
print(f"\neffect of B when A is OFF: {m[(0,1)] - m[(0,0)]:+.3f}")
print(f"effect of B when A is ON: {m[(1,1)] - m[(1,0)]:+.3f}")
print(f"interaction AB: {interaction:+.3f}"
f" (true {TRUE['AB_interaction']:+.2f})")
# A regression gives all effects with standard errors in one fit.
X = sm.add_constant(fac.assign(AB=fac["A"]*fac["B"])[["A", "B", "C", "AB"]])
fit = sm.OLS(fac["y"], X).fit()
print("\nregression with the interaction term:")
print(pd.DataFrame({"coef": fit.params.round(3), "se": fit.bse.round(3),
"p": fit.pvalues.round(4)}).to_string())
# --- ONE FACTOR AT A TIME: same total traffic, worse on both counts. ---
print("\n--- OFAT with the same 24,000 users split across 3 tests ---")
ofat_n = N // 3
ofat_effects = {}
for f, idx in [("A", 0), ("B", 1), ("C", 2)]:
r = np.random.default_rng(100 + idx)
arm = r.random(ofat_n) < 0.5
# Only factor f varies; the others stay at their baseline (off), which is
# exactly what "one factor at a time" means.
a = arm if f == "A" else np.zeros(ofat_n, bool)
b = arm if f == "B" else np.zeros(ofat_n, bool)
c = arm if f == "C" else np.zeros(ofat_n, bool)
y = outcome(a, b, c, ofat_n, r)
on_y, off_y = y[arm], y[~arm]
se = np.sqrt(on_y.var(ddof=1)/len(on_y) + off_y.var(ddof=1)/len(off_y))
ofat_effects[f] = on_y.mean() - off_y.mean()
print(f" {f}: {ofat_effects[f]:+.3f} se {se:.3f} (n {ofat_n:,})")
print(f"\nOFAT predicts A+B together = "
f"{ofat_effects['A'] + ofat_effects['B']:+.3f}")
print(f"factorial measures A+B together = {m[(1,1)] - m[(0,0)]:+.3f}")
print(f"OFAT overstates the combination by "
f"{ofat_effects['A']+ofat_effects['B'] - (m[(1,1)]-m[(0,0)]):+.3f}")
print("Shipping both based on OFAT would miss the target, and nobody would")
print("know why, because each test was individually correct.")
# --- STANDARD ERROR COMPARISON: factorial wins on precision too. ---
print("\nstandard error per main effect, same total users:")
print(f" factorial (all 24,000 per effect): "
f"{10*np.sqrt(1/12000 + 1/12000):.4f}")
print(f" OFAT (8,000 per effect): "
f"{10*np.sqrt(1/4000 + 1/4000):.4f}")
print("Factorial is more precise AND measures interactions. There is no")
print("trade-off here; OFAT is simply worse when factors can be combined.")
# --- CELL GROWTH AND FRACTIONAL DESIGNS ---
print("\ncell count as factors accumulate (10,000 users total):")
for k in range(2, 8):
full, half = 2**k, 2**(k-1)
print(f" {k} factors -> full {full:>3} cells ({10_000/full:>6.0f}/cell)"
f" half-fraction {half:>3} cells ({10_000/half:>6.0f}/cell)")
# A resolution IV half-fraction of 2^4: choose cells where the product of
# all four coded factors equals +1. Main effects stay clean of 2-factor
# interactions; 2-factor interactions are confounded in pairs.
coded = np.array(list(itertools.product([-1, 1], repeat=4)))
half_frac = coded[coded.prod(axis=1) == 1]
print(f"\n2^4 half-fraction: {len(half_frac)} of {len(coded)} cells")
print(pd.DataFrame(half_frac, columns=list("ABCD")).to_string(index=False))
# Verify the main effects are orthogonal within the fraction.
print("\ncolumn correlations within the fraction (0 = orthogonal, so each")
print("main effect is estimable independently of the others):")
print(pd.DataFrame(half_frac, columns=list("ABCD")).corr().round(3).to_string())
# --- CONCURRENT INDEPENDENT TESTS: what most companies actually run. ---
K = 30_000
ia = rng.random(K) < 0.5 # experiment 1, randomized independently
ib = rng.random(K) < 0.5 # experiment 2, independent of the first
y = outcome(ia, ib, np.zeros(K, bool), K, rng)
print(f"\nconcurrent independent tests, overlap check:"
f" {np.mean(ia & ib):.3f} (expected 0.25)")
for name, flag in [("exp 1 (A)", ia), ("exp 2 (B)", ib)]:
est = y[flag].mean() - y[~flag].mean()
print(f" {name}: {est:+.3f}")
print("Main effects come out fine because the assignments are independent.")
print("The interaction is present in the data and nobody computes it, so two")
print("teams each claim a win and the combined launch underdelivers.")
WHY FACTORIAL IS NOT MORE EXPENSIVE
2^3 factorial, 8 cells x 3,000 users = 24,000 total
cell A B C
1 0 0 0 |
2 0 0 1 | effect of A =
3 0 1 0 | mean(cells 5-8) - mean(cells 1-4)
4 0 1 1 | uses ALL 24,000 users
5 1 0 0 |
6 1 0 1 | and the SAME data gives B and C
7 1 1 0 |
8 1 1 1 |
OFAT with the same traffic: 3 tests x 8,000 = 4,000 per arm per effect
-> se is 1.73x larger, and no interaction is measurable at all.
THE INTERACTION OFAT CANNOT SEE
B off B on effect of B
A off 40.0 41.6 +1.6
A on 42.0 42.5 +0.5
^ B helps much less once A is present
interaction = (+0.5) - (+1.6) = -1.1
OFAT predicts A+B = +2.0 +1.6 = +3.6
reality A+B = +2.5
The shortfall is the interaction, and it is invisible to OFAT.
CELL GROWTH
factors full cells half-fraction users/cell at 10,000 total
2 4 2 2,500 / 5,000
3 8 4 1,250 / 2,500
4 16 8 625 / 1,250
5 32 16 313 / 625
6 64 32 156 / 313
7 128 64 78 / 156
RESOLUTION, FOR FRACTIONAL DESIGNS
III main effects confounded with 2-factor interactions risky
IV main effects clean; 2-factor interactions confounded usual choice
with each other
V main effects and 2-factor interactions all clean needs more cells
ONLINE PRACTICE: CONCURRENT INDEPENDENT TESTS
two experiments randomized independently on the same traffic
-> ~25% of users get both
-> main effects are unbiased
-> the interaction EXISTS in the data and is never estimated
Mitigation: put experiments that touch the same surface in one
mutually exclusive layer.
Key Takeaways
- A factorial design estimates every main effect using the full sample, so testing three factors together is more precise than three sequential tests.
- Only a factorial design measures interaction, and interaction is largest exactly when two features address the same user need.
- One-factor-at-a-time testing predicts the combined effect by addition, which is the assumption most likely to fail.
- Cell count doubles per factor, so use fractional factorial designs and accept confounded high-order interactions.
- Concurrent independent online tests give unbiased main effects and leave interactions unmeasured, so isolate experiments that share a surface.
🧪 Practice
- Simulate a two-by-two factorial with a known negative interaction and recover all three effects from a single regression.
- Compare the standard error of a main effect from a factorial design against a one-factor-at-a-time design with the same total traffic.
- Show that the sum of two individually measured effects overstates the combined effect when a negative interaction exists.
- Construct a half-fraction of a four-factor design and verify the main effect columns are mutually orthogonal.
- Interview: Two teams ship experiments on the same checkout page in the same week, each reporting a 2% lift, and the combined metric moves 2.1%. What happened? (Hint: what quantity did neither test estimate?)
Sample Size and Power Calculation
A test that cannot detect the effect you care about wastes the traffic it consumes and, worse, produces a non-significant result that gets read as evidence of no effect. Power analysis is how you avoid running that test, and it belongs before the experiment rather than after.
Chapter 3 covered power formally. The practical framing here is that four quantities are locked together, and fixing any three determines the fourth:
- The significance level, conventionally 0.05, the false positive rate you accept.
- Power, conventionally 0.80, the probability of detecting a real effect of the size you assumed.
- The effect size you want to detect, called the minimum detectable effect or MDE.
- The sample size.
In practice you fix the first two by convention, and then you have a relationship between the MDE and the sample size. This is the single most useful thing to internalise about power: the MDE shrinks with the square root of the sample size. Halving the effect you can detect requires four times the users. This is why detecting small effects is so expensive, and why variance reduction techniques matter so much.
For a continuous metric the sample size per arm is approximately
n = 2 * (z_alpha/2 + z_beta)^2 * variance / effect^2
and for a binary metric the variance term is the two group proportions times their complements. The constant for the standard 5% and 80% choice is about 7.849, which is worth remembering because it makes back-of-envelope calculation fast.
Three practical points that formulas hide.
The MDE is a business input, not a statistical one. The right question is not "what can we detect" but "what is the smallest effect that would change our decision". If a 0.1% lift is not worth shipping, do not size the test to find it. If your minimum viable effect requires more traffic than you have, the honest conclusion is that this experiment cannot be run, and you should either pick a more sensitive metric, reduce variance, or make the decision another way.
Ratio metrics and skewed metrics need care. Revenue per user is heavily skewed, so its variance is large and the sample requirement is high. Capping outliers, using a binary version such as "purchased at all", or applying CUPED all help. The variance of a ratio where the denominator also varies needs the delta method rather than the naive formula.
Multiple metrics multiply the false positive rate. Testing ten metrics at 5% gives roughly a 40% chance of at least one spurious significant result. Declare one primary metric in advance and treat the rest as secondary, or correct for multiplicity.
Finally, a warning about the aftermath. Post-hoc power calculated from the observed effect is not informative. It is a deterministic function of the p-value and tells you nothing you did not already know. If a test was underpowered, the fix is a wider interval and an honest "we cannot rule out effects between X and Y", not a power calculation run backwards.
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.power import NormalIndPower, TTestIndPower
from statsmodels.stats.proportion import proportion_effectsize
ALPHA, POWER = 0.05, 0.80
z_a = stats.norm.ppf(1 - ALPHA/2) # 1.959964
z_b = stats.norm.ppf(POWER) # 0.841621
K = (z_a + z_b) ** 2
print(f"z_alpha/2 = {z_a:.4f}, z_beta = {z_b:.4f}, (sum)^2 = {K:.4f}")
print("For 5% / 80% the constant is ~7.85. Memorise it for quick estimates.")
# --- BINARY METRIC: sample size per arm from a baseline and a relative lift. ---
def n_binary(p0: float, rel_lift: float, alpha=ALPHA, power=POWER) -> int:
p1 = p0 * (1 + rel_lift)
za, zb = stats.norm.ppf(1-alpha/2), stats.norm.ppf(power)
return int(np.ceil((za+zb)**2 * (p0*(1-p0) + p1*(1-p1)) / (p1-p0)**2))
print("\nusers PER ARM to detect a relative lift on a binary metric:")
base_rates = [0.02, 0.05, 0.20]
lifts = [0.01, 0.02, 0.05, 0.10]
tbl = pd.DataFrame({f"{l:.0%} lift": {f"base {p:.0%}": f"{n_binary(p,l):,}"
for p in base_rates} for l in lifts})
print(tbl.to_string())
# Cross-check against statsmodels, which uses the arcsine effect size.
es = proportion_effectsize(0.05*1.05, 0.05)
sm_n = NormalIndPower().solve_power(effect_size=es, alpha=ALPHA, power=POWER,
ratio=1.0, alternative="two-sided")
print(f"\n5% base, 5% lift: closed form {n_binary(0.05, 0.05):,}"
f" statsmodels {sm_n:,.0f} (different effect-size parameterisations)")
# --- THE SQUARE-ROOT LAW: the most important fact in power analysis. ---
print("\nMDE shrinks as 1/sqrt(n): halving the MDE costs 4x the users.")
p0, sigma_like = 0.05, np.sqrt(0.05*0.95)
for n in (5_000, 20_000, 80_000, 320_000, 1_280_000):
mde_abs = (z_a + z_b) * np.sqrt(2 * p0*(1-p0) / n)
print(f" n = {n:>9,} per arm -> MDE {mde_abs:.5f} absolute"
f" = {mde_abs/p0:6.2%} relative")
# --- CONTINUOUS METRIC, and why skew is so expensive. ---
def n_continuous(sd: float, abs_effect: float, alpha=ALPHA, power=POWER) -> int:
za, zb = stats.norm.ppf(1-alpha/2), stats.norm.ppf(power)
return int(np.ceil(2 * (za+zb)**2 * sd**2 / abs_effect**2))
rng = np.random.default_rng(5)
revenue = rng.lognormal(2.2, 1.5, 200_000) # heavily skewed, realistic
mean_rev = revenue.mean()
target = 0.02 * mean_rev # detect a 2% lift
variants = {
"raw revenue": revenue,
"capped at p99": np.clip(revenue, None, np.quantile(revenue, 0.99)),
"capped at p95": np.clip(revenue, None, np.quantile(revenue, 0.95)),
"log1p revenue": np.log1p(revenue),
}
print(f"\nmean revenue {mean_rev:.2f}, detecting a 2% lift ({target:.3f}):")
for name, v in variants.items():
# For a transformed metric the 2% lift is on that metric's own scale.
eff = 0.02 * v.mean()
print(f" {name:16s} sd {v.std(ddof=1):8.2f} cv {v.std()/v.mean():5.2f}"
f" -> {n_continuous(v.std(ddof=1), eff):>12,} per arm")
# A binary proxy is often far cheaper than the revenue metric itself.
purchased = (revenue > np.quantile(revenue, 0.55)).astype(float)
print(f" {'binary: purchased':16s} p {purchased.mean():.3f}"
f" -> {n_binary(purchased.mean(), 0.02):>12,} per arm")
# --- DURATION: convert users into calendar days, which is the real constraint. ---
def duration_days(n_per_arm: int, daily_users: int, n_arms: int = 2,
min_days: int = 7) -> dict:
days = n_per_arm * n_arms / daily_users
# Always run whole weeks: weekday/weekend behaviour differs systematically.
weeks = max(np.ceil(days / 7), np.ceil(min_days / 7))
return {"raw_days": round(days, 1), "run_weeks": int(weeks),
"run_days": int(weeks * 7)}
print("\nduration at 40,000 eligible users/day:")
for lift in (0.01, 0.02, 0.05, 0.10):
n = n_binary(0.05, lift)
d = duration_days(n, 40_000)
print(f" {lift:>4.0%} lift n {n:>9,}/arm {d['raw_days']:>7.1f} days"
f" -> run {d['run_weeks']} week(s)")
# --- UNEQUAL ALLOCATION costs sample size. ---
print("\ntotal users needed at various allocations (5% base, 5% lift):")
n_bal = n_binary(0.05, 0.05)
for ratio in (1.0, 2.0, 4.0, 9.0):
# n_a * n_b / (n_a + n_b) must match the balanced design's n/2.
n_a = int(np.ceil(n_bal * (1 + 1/ratio) / 2))
n_b = int(np.ceil(n_a * ratio))
print(f" 1:{ratio:<4.0f} arm A {n_a:>9,} arm B {n_b:>9,}"
f" total {n_a+n_b:>10,} ({(n_a+n_b)/(2*n_bal):.2f}x balanced)")
# --- MULTIPLE METRICS inflate the false positive rate. ---
print("\nprobability of at least one false positive at alpha = 0.05:")
for m in (1, 3, 5, 10, 20):
mult = ((stats.norm.ppf(1 - 0.05/(2*m)) + z_b) / (z_a + z_b)) ** 2
print(f" {m:>2} independent metrics: {1-0.95**m:.1%}"
f" Bonferroni alpha {0.05/m:.4f}"
f" -> n multiplier {mult:.2f}x")
# --- POST-HOC POWER IS NOT INFORMATIVE ---
n, p0 = 20_000, 0.05
observed_lift = 0.015
se = np.sqrt(2 * p0*(1-p0)/n)
z_obs = observed_lift * p0 / se
p_val = 2 * (1 - stats.norm.cdf(abs(z_obs)))
posthoc = stats.norm.cdf(abs(z_obs) - z_a) + stats.norm.cdf(-abs(z_obs) - z_a)
print(f"\nobserved lift {observed_lift:.1%} relative, p = {p_val:.3f}")
print(f"post-hoc 'power' from the observed effect: {posthoc:.2f}")
print("This is a function of the p-value alone and adds no information.")
print("Report the interval instead, and say what it rules out:")
obs_abs = observed_lift * p0
print(f" effect {obs_abs:+.5f} absolute, 95% CI "
f"[{obs_abs-1.96*se:+.5f}, {obs_abs+1.96*se:+.5f}]")
print(f" in relative terms: [{(obs_abs-1.96*se)/p0:+.2%},"
f" {(obs_abs+1.96*se)/p0:+.2%}]")
print(" -> we can rule out lifts above ~4%, not 'there is no effect'.")
# --- THE PRE-EXPERIMENT DECISION TABLE you should actually produce. ---
daily = 40_000
plan = []
for lift in (0.005, 0.01, 0.02, 0.03, 0.05):
n = n_binary(0.05, lift)
d = duration_days(n, daily)
plan.append({"target_lift": f"{lift:.1%}", "n_per_arm": f"{n:,}",
"weeks": d["run_weeks"],
"feasible": "yes" if d["run_weeks"] <= 4 else "NO"})
print("\nPRE-EXPERIMENT FEASIBILITY (40,000 users/day, 4-week limit)")
print(pd.DataFrame(plan).to_string(index=False))
print("\nDecide the MDE from the business case, then read off feasibility.")
print("If the smallest worthwhile effect is infeasible, do not run the test.")
THE FOUR LOCKED QUANTITIES
alpha (0.05) ---+
power (0.80) ---+---> fix any three, the fourth follows
MDE ---+
sample size ---+
n per arm (continuous) = 2 * (z_a/2 + z_b)^2 * var / effect^2
n per arm (binary) = (z_a/2 + z_b)^2 * (p0*q0 + p1*q1) / (p1-p0)^2
at 5% / 80%, (z_a/2 + z_b)^2 = 7.849
THE SQUARE-ROOT LAW
n per arm MDE (relative, 5% base)
5,000 10.9%
20,000 5.4% <- 4x the users, half the MDE
80,000 2.7% <- 4x again, half again
320,000 1.4%
1,280,000 0.7%
Detecting HALF the effect costs FOUR TIMES the traffic. This is why
variance reduction (CUPED, capping, better metrics) pays so well.
WHY SKEWED METRICS ARE EXPENSIVE
metric cv n per arm for a 2% lift
raw revenue 2.9 ~3,400,000
capped at p99 1.7 ~1,150,000
capped at p95 1.2 ~580,000
binary: purchased 1.1 ~125,000
Same underlying question, 27x difference in cost. Choose the metric
with the business, not after the power calculation fails.
MDE IS A BUSINESS INPUT
WRONG: "we have 4 weeks of traffic, so our MDE is 2.7%"
RIGHT: "a lift below 1% would not justify the maintenance cost, so
our MDE is 1%; that needs 12 weeks, which we do not have.
Options: a more sensitive metric, CUPED, or decide without
an experiment."
POST-HOC POWER
computed from the OBSERVED effect, it is a one-to-one function of the
p-value and carries no new information.
Instead report: effect +1.5% relative, 95% CI [-0.9%, +3.9%]
-> "rules out lifts above ~4%", NOT "no effect"
Key Takeaways
- Significance level, power, minimum detectable effect, and sample size are locked together, so fixing three determines the fourth.
- The detectable effect shrinks only with the square root of sample size, so halving it costs four times the traffic.
- Set the minimum detectable effect from the smallest effect that would change your decision, not from the traffic you happen to have.
- Skewed metrics such as revenue per user are far more expensive than binary proxies, so choose the metric before the power calculation fails.
- Post-hoc power from the observed effect restates the p-value; report the confidence interval and what it rules out instead.
🧪 Practice
- Write a sample size function for a binary metric and tabulate users per arm across baseline rates of 2%, 5%, and 20% and lifts of 1%, 2%, and 5%.
- Demonstrate the square-root law by computing the minimum detectable effect at four sample sizes each four times the last.
- Compare the sample requirement for raw revenue, revenue capped at the 99th percentile, and a binary purchase indicator.
- Build a feasibility table converting target lifts into run weeks at a given daily traffic, and mark which are achievable in four weeks.
- Interview: A test ran for two weeks, showed a 1.5% lift with p equal to 0.18, and the team concluded the feature does not work. What do you tell them? (Hint: what does the confidence interval include, and what is the difference between no evidence of an effect and evidence of no effect?)
<a id="82-ab-testing-in-practice"></a>
8.2 A/B Testing in Practice
A correctly designed experiment can still produce a wrong decision through how it is monitored, stopped, and read; this section covers the operational discipline that turns a valid design into a trustworthy result.
Metric Selection: Primary and Guardrail
The metric you choose determines what the experiment can tell you, and metric selection is where most experiments are won or lost before any data arrives. The core tension is that the metrics which move fastest are usually the ones that matter least.
Start with the structure. A well-specified experiment has exactly one primary metric, a small number of secondary metrics, and a set of guardrails.
The primary metric is the single number that decides whether you ship. Having one is not a bureaucratic constraint, it is what makes the decision rule unambiguous and the statistics valid. If you have five primary metrics you have five chances to find significance, and you will almost always find one. Declared in advance, in writing, before the data is seen.
A good primary metric has four properties. It is sensitive, meaning it actually moves when the feature works. It is aligned with the outcome you care about, so moving it is genuinely good. It is measurable within the experiment window, which rules out annual retention. And it is hard to game, because a team optimising it will find every loophole.
These properties conflict, which is the whole difficulty. Clicks are sensitive, easy to measure, and weakly aligned: you can raise clicks by making a page confusing. Lifetime value is perfectly aligned and cannot be measured in two weeks. The practical resolution is a surrogate: a short-term metric with a demonstrated relationship to the long-term outcome, where the demonstration is the important part. Validate the surrogate against long-run data before you trust it, because an unvalidated surrogate is just a convenient number.
Guardrail metrics are the ones you are not trying to improve but refuse to damage: page load time, error rate, unsubscribe rate, support contacts, refund rate. Their role is asymmetric, which changes the statistics. For the primary metric you ask "is there evidence of improvement". For a guardrail you ask "can I rule out meaningful harm", which is a non-inferiority test against a tolerance, not a test against zero. A guardrail that is merely non-significant is not reassuring if the confidence interval extends well into harmful territory.
Two further distinctions worth having.
Ratio metrics with a user-level denominator need care. Clicks per session computed over all sessions has a denominator that treatment can change. If treatment increases sessions, the ratio can fall while total clicks rise. Decide whether you are measuring per-user or per-session behaviour and be consistent.
Overall Evaluation Criterion. When several metrics matter, some organisations combine them into one weighted composite so there is a single decision number. This forces the trade-offs to be made explicit in advance rather than argued after the results arrive, which is its real value.
import numpy as np
import pandas as pd
from scipy import stats
rng = np.random.default_rng(6)
# --- SENSITIVITY VS ALIGNMENT: the central trade-off, quantified. ---
n = 40_000
metrics = pd.DataFrame([
# (name, baseline, sd/baseline, correlation with long-run value)
("clicks per user", 8.0, 0.55, 0.15),
("sessions per user", 4.2, 0.70, 0.35),
("conversion rate", 0.05, 4.36, 0.62),
("7-day retention", 0.42, 1.18, 0.78),
("90-day revenue", 42.00, 2.90, 0.95),
], columns=["metric", "baseline", "cv", "corr_with_ltv"])
z = stats.norm.ppf(0.975) + stats.norm.ppf(0.80)
metrics["mde_rel"] = z * np.sqrt(2 / n) * metrics["cv"]
metrics["n_for_2pct"] = (2 * z**2 * metrics["cv"]**2 / 0.02**2).astype(int)
print(metrics.round(3).to_string(index=False))
print("\nThe metrics most aligned with long-run value are the least sensitive.")
print("Clicks detect a 2% move in 40k users; 90-day revenue needs ~165k.")
print("Choosing clicks because it is cheap means optimising the wrong thing.")
# --- VALIDATING A SURROGATE: the step that makes a proxy legitimate. ---
# Historical experiments where BOTH the short-term proxy and the long-term
# outcome were eventually measured. This is the evidence a surrogate needs.
past = pd.DataFrame({
"experiment": [f"e{i}" for i in range(1, 13)],
"proxy_lift_pct": [2.1, -0.8, 4.0, 1.2, 0.3, 3.1, -1.5, 2.8, 0.9, 5.2, 1.7, -0.4],
"ltv_lift_pct": [1.4, -0.6, 2.2, 0.9, 0.1, 1.9, -1.1, 1.6, 0.4, 2.6, 1.0, -0.2],
})
slope, intercept, r, p, se = stats.linregress(past["proxy_lift_pct"],
past["ltv_lift_pct"])
print(f"\nsurrogate validation over {len(past)} past experiments:")
print(f" correlation r = {r:.3f} (r^2 = {r**2:.2f}), p = {p:.2e}")
print(f" a 1% proxy lift historically delivered {slope:.2f}% long-run lift")
sign_match = np.mean(np.sign(past["proxy_lift_pct"])
== np.sign(past["ltv_lift_pct"]))
print(f" sign agreement: "
f"{sign_match:.0%}")
print("Without this table, a surrogate is an assumption wearing a number.")
# --- GUARDRAILS ARE NON-INFERIORITY TESTS, not tests against zero. ---
def guardrail_verdict(est: float, se: float, tolerance: float,
higher_is_worse: bool = True) -> dict:
"""Ask 'can I rule out harm beyond the tolerance?', not 'is it zero?'."""
lo, hi = est - 1.96*se, est + 1.96*se
harm_bound = hi if higher_is_worse else lo
ruled_out = ((harm_bound < tolerance) if higher_is_worse
else (harm_bound > -tolerance))
return {"estimate": round(est, 4), "ci": (round(lo, 4), round(hi, 4)),
"tolerance": tolerance,
"verdict": "PASS: harm beyond tolerance ruled out" if ruled_out
else "FAIL: cannot rule out meaningful harm"}
print("\nguardrail: page load time, tolerance +20 ms")
for label, est, se in [("clearly fine", 3.0, 4.0),
("not significant but UNDERPOWERED", 8.0, 30.0),
("real regression", 45.0, 8.0)]:
v = guardrail_verdict(est, se, tolerance=20.0)
print(f" {label:34s} {v['estimate']:>7.1f} ms CI {v['ci']}"
f"\n {'':34s} {v['verdict']}")
print("\nThe middle case has p > 0.05 on a test against zero, and its interval")
print("reaches +67 ms. 'Not significant' is not the same as 'no harm'.")
# --- THE RATIO-METRIC DENOMINATOR TRAP ---
m = 30_000
arm = rng.random(m) < 0.5
# Treatment brings users back more often (more sessions) but each new session
# is a lower-intent visit, so clicks PER SESSION falls while TOTAL clicks rise.
sessions = rng.poisson(np.where(arm, 5.0, 4.0), m) + 1
clicks = rng.poisson(sessions * np.where(arm, 1.75, 1.90))
df = pd.DataFrame({"arm": arm, "sessions": sessions, "clicks": clicks})
per_session_ratio = (df.groupby("arm").apply(
lambda g: g["clicks"].sum() / g["sessions"].sum(), include_groups=False))
per_user_clicks = df.groupby("arm")["clicks"].mean()
print("\nthe same experiment, two denominators:")
print(f" clicks per SESSION: control {per_session_ratio[False]:.3f}"
f" -> treated {per_session_ratio[True]:.3f}"
f" ({per_session_ratio[True]/per_session_ratio[False]-1:+.1%}) DOWN")
print(f" clicks per USER: control {per_user_clicks[False]:.3f}"
f" -> treated {per_user_clicks[True]:.3f}"
f" ({per_user_clicks[True]/per_user_clicks[False]-1:+.1%}) UP")
print("Both are correct. Treatment increased the denominator. Decide which")
print("question you are asking BEFORE you see this.")
# --- A METRIC SPEC, declared before the test runs. ---
SPEC = {
"experiment": "checkout_autofill_v1",
"hypothesis": "Address autofill reduces form abandonment, raising "
"checkout completion.",
"primary": {
"name": "checkout_completion_rate",
"definition": "completed_orders / checkout_starts, user-level, "
"first checkout only",
"baseline": 0.412, "mde_rel": 0.02, "direction": "increase",
"alpha": 0.05, "power": 0.80,
},
"secondary": [
{"name": "time_to_complete_s", "direction": "decrease"},
{"name": "orders_per_user", "direction": "increase"},
],
"guardrails": [
{"name": "page_load_p95_ms", "tolerance_abs": 20, "worse": "higher"},
{"name": "refund_rate", "tolerance_rel": 0.05, "worse": "higher"},
{"name": "support_contacts_per_1k", "tolerance_rel": 0.10, "worse": "higher"},
{"name": "payment_error_rate", "tolerance_rel": 0.02, "worse": "higher"},
],
"decision_rule": "Ship if primary is significantly positive AND every "
"guardrail passes non-inferiority. Do not ship on a "
"secondary metric alone.",
}
print(f"\nmetric spec for {SPEC['experiment']}:")
print(f" primary: 1 ({SPEC['primary']['name']})")
print(f" secondary: {len(SPEC['secondary'])} (not decision-making)")
print(f" guardrails: {len(SPEC['guardrails'])}")
print(f" rule: {SPEC['decision_rule']}")
# --- WHY ONE PRIMARY: the false positive rate with several. ---
print("\nchance of at least one significant 'primary' under no true effect:")
for k in (1, 2, 3, 5, 8):
print(f" {k} primary metrics: {1-0.95**k:.1%}")
print("Five 'primaries' means a 23% chance of shipping nothing on evidence.")
# --- OVERALL EVALUATION CRITERION: force the trade-off in advance. ---
WEIGHTS = {"conversion_rate": 0.55, "orders_per_user": 0.30,
"session_length_s": 0.15}
observed = {"conversion_rate": 0.021, "orders_per_user": 0.014,
"session_length_s": -0.035} # relative lifts
oec = sum(WEIGHTS[k] * observed[k] for k in WEIGHTS)
print(f"\nOEC = {oec:+.4f}")
for k, w in WEIGHTS.items():
print(f" {k:20s} weight {w:.2f} lift {observed[k]:+.3f}"
f" contribution {w*observed[k]:+.4f}")
print("Session length fell and the composite is still positive. The weights")
print("encode that trade-off, and they were agreed BEFORE the results.")
THE METRIC HIERARCHY
PRIMARY (exactly one) decides ship / no-ship
| declared in writing before data
| sized for in the power calculation
SECONDARY (a few) explain the mechanism, never decide alone
|
GUARDRAILS (a set) must not be harmed
tested for NON-INFERIORITY, not against zero
THE SENSITIVITY / ALIGNMENT TRADE-OFF
metric sensitivity alignment with value n for a 2% lift
---------------------------------------------------------------------------
clicks per user highest weakest (r=0.15) 12,000
sessions per user high weak (r=0.35) 19,000
conversion rate moderate good (r=0.62) 746,000
7-day retention low strong (r=0.78) 55,000
90-day revenue lowest strongest (r=0.95) 330,000
Picking the cheap metric optimises the thing you do not care about.
The fix is a VALIDATED surrogate, not a convenient one.
GUARDRAIL LOGIC IS ASYMMETRIC
PRIMARY: is there evidence of a gain?
0 |-----o-----|
|----------------------------------------^ CI excludes 0 -> ship
GUARDRAIL: can I rule out harm beyond the tolerance?
tolerance
0 |
|----o----| | PASS: whole CI inside tolerance
|----------o----------| FAIL: CI crosses the tolerance,
| even though it also crosses 0
"Not significant" on an underpowered guardrail is not reassurance.
THE RATIO DENOMINATOR TRAP
clicks per SESSION: 1.90 -> 1.75 DOWN 8%
clicks per USER: 7.60 -> 8.75 UP 15%
Treatment raised sessions. Both numbers are right; they answer
different questions. Fix the denominator in the spec.
Key Takeaways
- Declare exactly one primary metric in writing before the data arrives, or multiple chances at significance guarantee a false win.
- Sensitivity and alignment with long-term value conflict, so short-term surrogates must be validated against historical long-run outcomes.
- Guardrails are non-inferiority tests against a tolerance, so a non-significant but wide interval is a failure rather than a pass.
- A ratio metric whose denominator the treatment can change will disagree with its per-user counterpart, so fix the denominator in advance.
- A weighted composite criterion forces trade-offs to be agreed before results rather than argued afterwards.
🧪 Practice
- For five candidate metrics with different coefficients of variation, compute the sample size needed to detect a 2% relative lift in each.
- Build a surrogate validation table from historical experiments and report the correlation and sign agreement between proxy and long-run lift.
- Write a guardrail verdict function taking an estimate, standard error, and tolerance, and demonstrate a case that is non-significant yet fails.
- Construct a simulation where clicks per session falls while clicks per user rises, and explain the mechanism.
- Interview: Product wants to declare a win on any of conversion, engagement, or retention, whichever moves. What is wrong and what do you propose? (Hint: what is the false positive rate across three tests, and what makes a decision rule valid?)
Test Duration and Peeking Problems
The most common statistical error in industrial experimentation is checking results repeatedly and stopping when they look significant. It feels harmless and even responsible, and it inflates the false positive rate severely.
The reason is worth building up carefully, because the intuition that "more data cannot hurt" is what makes the error so durable.
A p-value below 0.05 means: if there were no effect, data this extreme would occur 5% of the time. That guarantee applies to one test at one prespecified moment. When you look twenty times over a two-week test, you get twenty chances for the running estimate to wander below the threshold. Under no true effect, the running difference is essentially a random walk, and a random walk crosses any fixed boundary far more often than it ends beyond it. Checking daily for two weeks pushes the false positive rate from 5% to roughly 20-30%. Checking continuously, with a rule to stop as soon as significance appears, drives it toward 100% given enough time.
So the classical fix is procedural: fix the duration in advance and analyse once. Compute the sample size, run to it, then look.
That prescription collides with reality. Businesses want early reads on harmful changes, and waiting two weeks to discover a payment bug is unacceptable. Three resolutions exist, and choosing among them is a design decision.
Separate monitoring from inference. Watch guardrails and error rates continuously for operational safety, and analyse the primary metric once at the planned end. Stopping early for harm is legitimate; stopping early for a win is not, because the asymmetry is deliberate.
Use a method designed for repeated looks. Sequential and Bayesian approaches, covered in the next topic, provide valid inference under continuous monitoring.
Pre-register interim analyses with spending functions. Group sequential designs allocate the total false positive budget across a few planned checkpoints using stricter thresholds early, so the overall rate stays at 5%.
Separately from peeking, duration has substantive requirements that have nothing to do with statistics.
Run whole weeks. Behaviour differs systematically by day of week, so a test ending on a Wednesday over-weights weekdays. Always run in multiples of seven days.
Cover at least one full business cycle, including any weekly or monthly pattern relevant to your product.
Allow for novelty and primacy effects. Users react to change as change, and that reaction decays. A test that is too short measures the reaction rather than the steady state.
Beware very long tests. Cookie churn, users switching devices, and accumulating contamination all degrade a long-running experiment, and the population composition drifts.
import numpy as np
import pandas as pd
from scipy import stats
# --- PEEKING: simulate the actual error rate under no true effect. ---
def peeking_simulation(n_total=20_000, n_looks=1, n_sims=4000,
true_effect=0.0, seed=0):
"""Stop as soon as p < 0.05 at any scheduled look. Measure how often
we wrongly declare a winner when the truth is no effect."""
r = np.random.default_rng(seed)
checkpoints = np.linspace(n_total/n_looks, n_total, n_looks).astype(int)
false_positives = 0
for _ in range(n_sims):
a = r.normal(0, 1, n_total)
b = r.normal(true_effect, 1, n_total)
for cp in checkpoints:
half = cp // 2
if half < 30:
continue
t = stats.ttest_ind(a[:half], b[:half], equal_var=False)
if t.pvalue < 0.05:
false_positives += 1
break # stop-on-significance, as people do
return false_positives / n_sims
print("FALSE POSITIVE RATE under NO true effect, stop-on-significance:")
for looks in (1, 2, 5, 10, 20, 50):
rate = peeking_simulation(n_looks=looks, n_sims=2000, seed=looks)
print(f" {looks:>2} look(s): {rate:.1%}"
f" {'<- the nominal 5%' if looks == 1 else ''}")
print("\nOne look gives the 5% you asked for. Twenty looks nearly quintuples it.")
# --- WHY: the running estimate is a random walk that crosses boundaries. ---
rng = np.random.default_rng(11)
n = 8_000
print("\nfive A/A tests, tracking the running p-value:")
for sim in range(5):
r = np.random.default_rng(200 + sim)
a, b = r.normal(0, 1, n), r.normal(0, 1, n) # NO effect
crossed_at, min_p = None, 1.0
for cp in range(200, n, 200):
p = stats.ttest_ind(a[:cp], b[:cp], equal_var=False).pvalue
min_p = min(min_p, p)
if p < 0.05 and crossed_at is None:
crossed_at = cp
final = stats.ttest_ind(a, b, equal_var=False).pvalue
print(f" sim {sim}: min p over the run {min_p:.4f}"
f" first crossed at n={crossed_at}"
f" FINAL p {final:.3f}")
print("\nSeveral runs dip below 0.05 mid-flight and end nowhere near it.")
print("Stopping at the dip converts a null result into a 'win'.")
# --- AN ASCII VIEW OF THE RANDOM WALK ---
r = np.random.default_rng(7)
a, b = r.normal(0, 1, 6000), r.normal(0, 1, 6000)
track = []
for cp in range(300, 6001, 300):
d = b[:cp].mean() - a[:cp].mean()
se = np.sqrt(a[:cp].var(ddof=1)/cp + b[:cp].var(ddof=1)/cp)
track.append(d / se) # the z statistic over time
print("\nz statistic over time (A/A, true effect = 0):")
print(" n: " + "".join(f"{i*300:>5}" for i in range(1, 21)))
for level, label in [(2, " +2 "), (0, " 0 "), (-2, " -2 ")]:
row = ""
for z in track:
if abs(z - level) < 0.5:
row += " * "
else:
row += " . " if level == 0 else " "
print(f" {label}|" + row)
print(f" crossings of |z| > 1.96: "
f"{sum(1 for z in track if abs(z) > 1.96)} of {len(track)} checkpoints")
# --- ALPHA SPENDING: a valid way to look more than once. ---
def obrien_fleming_bounds(n_looks: int, alpha: float = 0.05) -> pd.DataFrame:
"""Approximate O'Brien-Fleming boundaries: strict early, near-nominal
at the end, so the OVERALL false positive rate stays at alpha."""
t = np.arange(1, n_looks + 1) / n_looks # information fraction
z_final = stats.norm.ppf(1 - alpha/2)
z_bounds = z_final / np.sqrt(t) # the defining shape
return pd.DataFrame({
"look": range(1, n_looks + 1),
"info_fraction": t.round(2),
"z_boundary": z_bounds.round(3),
"p_boundary": (2 * (1 - stats.norm.cdf(z_bounds))).round(5),
})
print("\nO'Brien-Fleming boundaries for 5 planned looks:")
print(obrien_fleming_bounds(5).to_string(index=False))
print("\nAt the first look you need p < 0.00001; at the last, p < 0.05.")
print("Total false positive spend stays at 5%. Pocock boundaries instead")
print("use one constant, less strict threshold at every look.")
# --- POWER COST of interim analyses: you pay a little for the option. ---
print("\nthe price of planned interim looks:")
for looks, penalty in [(1, 1.00), (2, 1.01), (3, 1.02), (5, 1.04), (10, 1.07)]:
print(f" {looks:>2} looks -> about {penalty:.2f}x the sample size"
f" for the same power")
print("A few percent more traffic buys the ability to stop early. Cheap.")
# --- DAY-OF-WEEK: why partial weeks distort the estimate. ---
rng = np.random.default_rng(9)
days = pd.date_range("2024-04-01", periods=21, freq="D")
dow_multiplier = {0: 1.05, 1: 1.08, 2: 1.06, 3: 1.02, 4: 0.95, 5: 0.78, 6: 0.74}
daily = pd.DataFrame({"date": days})
daily["dow"] = daily["date"].dt.dayofweek
daily["conv"] = [0.05 * dow_multiplier[d] * (1 + rng.normal(0, .02))
for d in daily["dow"]]
print("\nconversion by day of week (weekends are much lower):")
print(daily.groupby("dow")["conv"].mean().round(4).to_string())
for stop in (11, 14, 17, 21):
sub = daily.head(stop)
print(f" stopping after {stop:>2} days ({stop/7:.1f} weeks):"
f" mean conv {sub['conv'].mean():.5f}"
f" {'(whole weeks)' if stop % 7 == 0 else '(PARTIAL WEEK: biased)'}")
print("An 11-day test over-weights weekdays and overstates the baseline.")
# --- A DURATION PLAN that respects both statistics and behaviour. ---
def plan_duration(n_needed_per_arm: int, daily_eligible: int,
treat_share: float = 0.5, novelty_days: int = 7) -> dict:
stat_days = n_needed_per_arm / (daily_eligible * treat_share)
# Take the binding constraint, then round UP to whole weeks.
required = max(stat_days, novelty_days, 7)
weeks = int(np.ceil(required / 7))
return {"days_for_power": round(stat_days, 1),
"novelty_floor_days": novelty_days,
"binding_constraint": ("power" if stat_days >= novelty_days
else "novelty decay"),
"run_weeks": weeks, "run_days": weeks * 7,
"analyse": "once, at the end (guardrails monitored daily)"}
for n_need, daily in [(120_000, 50_000), (12_000, 50_000)]:
print(f"\nplan for n={n_need:,}/arm at {daily:,} users/day:")
for k, v in plan_duration(n_need, daily).items():
print(f" {k:22s} {v}")
THE PEEKING PENALTY
looks false positive rate (true effect = 0)
1 5.0% <- what you asked for
2 8.3%
5 14.2%
10 19.7%
20 25.6%
50 33.8%
continuous -> approaches 100% given unlimited time
WHY: THE RUNNING z STATISTIC IS A RANDOM WALK
+2 | * *
| * * * *
0 |. . *. . . *.*. . . .*. . . *. .
| * * * *
-2 |* *
+--------------------------------> n
A walk CROSSES a boundary far more often than it ENDS beyond it.
Stop-on-crossing turns a null into a win.
THE ASYMMETRY THAT IS ALLOWED
stop early for HARM legitimate: you are protecting users, and a
false alarm costs you a shipped feature
stop early for a WIN invalid: this is the peeking error, and it costs
you a wrong belief about the world
O'BRIEN-FLEMING BOUNDARIES (5 looks, total alpha 5%)
look info z boundary p boundary
1 0.20 4.383 0.00001 <- nearly impossible to cross
2 0.40 3.099 0.00194
3 0.60 2.530 0.01141
4 0.80 2.191 0.02845
5 1.00 1.960 0.05000 <- the usual threshold at the end
Strict early, nominal at the end. Cost: ~2-4% more sample size.
DURATION CONSTRAINTS THAT ARE NOT STATISTICAL
whole weeks day-of-week effects are large and systematic
>= 1 business cycle weekly and monthly patterns must be covered
>= novelty window or you measure the reaction, not the steady state
not TOO long cookie churn, device switching, population drift
Key Takeaways
- A p-value's guarantee holds for one prespecified analysis, so repeated looks with stop-on-significance inflate the false positive rate several times over.
- The running test statistic is a random walk, and a walk crosses a boundary far more often than it ends beyond it.
- Monitor guardrails continuously for harm, but analyse the primary metric once at the planned end; the asymmetry is deliberate.
- Alpha-spending boundaries such as O'Brien-Fleming permit planned interim looks for a few percent more sample size.
- Run whole weeks covering at least one business cycle and the novelty window, and avoid tests so long that the population drifts.
🧪 Practice
- Simulate stop-on-significance under no true effect at 1, 5, and 20 looks, and report the false positive rate at each.
- Track the running p-value of a single A/A test and find a case where it dips below 0.05 mid-flight and ends above it.
- Compute O'Brien-Fleming boundaries for five looks and confirm the final boundary is close to the nominal threshold.
- Show that stopping an experiment after 11 days biases a conversion estimate relative to stopping after 14.
- Interview: Your manager checks the experiment dashboard every morning and asks to ship as soon as it goes green. How do you handle this? (Hint: what can they legitimately watch daily, and what would let you give them a valid early read?)
Sequential and Bayesian Testing
The peeking problem has two principled solutions rather than one, and they come from different statistical traditions. Both let you look at results continuously and act on what you see, which is what operational reality demands.
Sequential testing keeps the frequentist framework and changes the boundary. The insight is that if you know in advance you will look many times, you can compute thresholds that account for it. Three variants matter in practice.
Group sequential designs schedule a fixed number of looks and spend the false positive budget across them, as covered in the previous topic.
Always-valid inference, built on sequential probability ratio tests and mixture martingale methods, produces confidence sequences that are valid at every moment simultaneously. You can look continuously, stop whenever you like, and the coverage guarantee holds. This is what most modern experimentation platforms use for continuous monitoring, and the cost is that the intervals are wider than fixed-horizon ones, typically by 20-40%. You are paying for optionality in the currency of precision.
Sequential probability ratio tests compare the likelihood of the data under two hypotheses and stop when the ratio crosses a boundary, which is optimal in expected sample size when the two hypotheses are well specified.
Bayesian testing changes the question. Instead of "how surprising is this data if there is no effect", it asks "given the data and my prior beliefs, what is the probability distribution over the true effect". The output is a posterior, from which you read directly interpretable quantities: the probability the treatment is better, the probability the effect exceeds a threshold worth shipping for, and the expected loss from making the wrong choice.
Two things about Bayesian testing are commonly misunderstood.
First, it does not make peeking free in the way people assume. The posterior is a valid summary of the evidence at any point, and in that sense continuous monitoring is coherent. But a decision rule such as "stop when the posterior probability of being better exceeds 95%" still has a frequentist error rate, and that rate can exceed what you would accept, particularly with a weak prior and small samples. If you care about long-run error rates, you must simulate the decision rule you actually intend to use.
Second, the prior matters most when data is scarce, which is exactly when you are tempted to stop early. An informative prior built from historical experiments is genuinely valuable: most tested changes have small effects, and a prior reflecting that shrinks extreme early estimates toward zero, which is the correct behaviour.
The practical comparison is less about correctness than about what you want to say. A frequentist result controls error rates and cannot tell you the probability that the treatment is better. A Bayesian result gives you that probability and depends on a prior. Expected loss is the Bayesian output most useful for decisions, because it converts the posterior into the currency of the choice being made.
import numpy as np
import pandas as pd
from scipy import stats
rng = np.random.default_rng(13)
# --- BAYESIAN A/B TEST for a conversion rate, with a conjugate prior. ---
# Beta prior + binomial likelihood -> Beta posterior. Chapter 3.4 has the
# derivation; here it is the decision machinery that matters.
def beta_binomial_test(conv_a, n_a, conv_b, n_b,
prior_alpha=1.0, prior_beta=1.0, draws=200_000,
min_effect=0.0, seed=0):
r = np.random.default_rng(seed)
post_a = r.beta(prior_alpha + conv_a, prior_beta + n_a - conv_a, draws)
post_b = r.beta(prior_alpha + conv_b, prior_beta + n_b - conv_b, draws)
lift = post_b - post_a
rel = lift / post_a
return {
"p_b_better": float((lift > 0).mean()),
"p_beats_threshold": float((rel > min_effect).mean()),
"median_rel_lift": float(np.median(rel)),
"cred_int_95": (float(np.quantile(rel, .025)), float(np.quantile(rel, .975))),
# Expected loss: if we ship B and B is worse, how much do we lose?
"expected_loss_ship_b": float(np.maximum(-lift, 0).mean()),
"expected_loss_keep_a": float(np.maximum(lift, 0).mean()),
}
res = beta_binomial_test(conv_a=1_020, n_a=20_000, conv_b=1_115, n_b=20_000,
min_effect=0.02)
print("BAYESIAN RESULT (control 5.10%, treatment 5.58%)")
for k, v in res.items():
print(f" {k:24s} {v if isinstance(v, tuple) else round(v, 5)}")
print("\nRead it directly: a 93% chance the true lift exceeds the 2% we")
print("require to ship, and the expected loss from shipping is tiny.")
# The frequentist answer to the same data, for comparison.
from statsmodels.stats.proportion import proportions_ztest
stat, pval = proportions_ztest([1_115, 1_020], [20_000, 20_000])
print(f"\nfrequentist: z {stat:.3f}, p {pval:.4f}")
print("It controls the error rate and cannot state P(B is better).")
# --- THE PRIOR MATTERS MOST WHEN DATA IS SCARCE ---
print("\nsame small sample, three priors (control 10/200, treatment 18/200):")
priors = {"uninformative Beta(1,1)": (1, 1),
"weak, centred on 6%: Beta(3,47)": (3, 47),
"strong, from history: Beta(30,470)": (30, 470)}
for label, (pa, pb) in priors.items():
r = beta_binomial_test(10, 200, 18, 200, pa, pb, min_effect=0.0, seed=1)
print(f" {label:34s} P(better) {r['p_b_better']:.3f}"
f" median lift {r['median_rel_lift']:+.1%}")
print("\nA historically grounded prior pulls an extreme early estimate toward")
print("zero, which is correct: most changes have small effects.")
# --- DOES A BAYESIAN STOPPING RULE CONTROL ERROR? SIMULATE IT. ---
def bayes_stopping_error(true_rate=0.05, true_lift=0.0, n_max=20_000,
check_every=1_000, threshold=0.95,
prior=(1.0, 1.0), sims=1_500, seed=0):
"""Stop as soon as P(B better) > threshold. Measure how often we stop
for B when the truth is no difference."""
r = np.random.default_rng(seed)
wrong = 0
for _ in range(sims):
a_conv = b_conv = 0
for n in range(check_every, n_max + 1, check_every):
a_conv += r.binomial(check_every, true_rate)
b_conv += r.binomial(check_every, true_rate * (1 + true_lift))
pa = r.beta(prior[0] + a_conv, prior[1] + n - a_conv, 4_000)
pb = r.beta(prior[0] + b_conv, prior[1] + n - b_conv, 4_000)
if (pb > pa).mean() > threshold:
wrong += 1
break
return wrong / sims
print("\nfalse 'B wins' rate for a Bayesian stop-when-P>0.95 rule,")
print("under NO true effect:")
for prior, label in [((1.0, 1.0), "flat Beta(1,1)"),
((30.0, 570.0), "informative Beta(30,570)")]:
rate = bayes_stopping_error(prior=prior, sims=400, seed=3)
print(f" {label:26s} {rate:.1%}")
print("\nThe flat prior's rule fires well above 5%. Bayesian inference is")
print("coherent under peeking; a Bayesian DECISION RULE still has an error")
print("rate you must simulate if you care about it.")
# --- ALWAYS-VALID CONFIDENCE SEQUENCES: valid at every moment. ---
def confidence_sequence(x: np.ndarray, y: np.ndarray, alpha=0.05, rho=None):
"""A mixture-martingale confidence sequence for a difference of means.
Width shrinks like sqrt(log(n)/n) rather than sqrt(1/n), which is the
price of being valid at EVERY n simultaneously."""
n = len(x)
diff = y.mean() - x.mean()
var = (x.var(ddof=1) + y.var(ddof=1)) / 2
rho = rho or max(n / 10, 1) # tuning: where it is tightest
# Radius of the always-valid interval at sample size n.
radius = np.sqrt(2 * var * (n * rho + 1) / (n ** 2 * rho)
* np.log(np.sqrt(n * rho + 1) / alpha))
return diff, diff - radius, diff + radius
print("\nfixed-horizon vs always-valid interval width (true effect 0.05):")
x_all = rng.normal(0, 1, 40_000)
y_all = rng.normal(0.05, 1, 40_000)
for n in (2_000, 10_000, 40_000):
x, y = x_all[:n], y_all[:n]
se = np.sqrt(x.var(ddof=1)/n + y.var(ddof=1)/n)
fixed_w = 2 * 1.96 * se
_, lo, hi = confidence_sequence(x, y)
print(f" n={n:>6} fixed width {fixed_w:.4f}"
f" always-valid width {hi-lo:.4f}"
f" ratio {(hi-lo)/fixed_w:.2f}x")
print("\n20-40% wider, and you may stop at any moment without penalty.")
print("That is the trade: precision for optionality.")
# --- EXPECTED LOSS: the output that maps onto a decision. ---
print("\ndecision by expected loss, threshold = 0.0005 absolute (0.05 pts):")
CARE_THRESHOLD = 0.0005
for label, (ca, na, cb, nb) in {
"clear winner": (1_000, 20_000, 1_160, 20_000),
"too close": (1_000, 20_000, 1_015, 20_000),
"clear loser": (1_000, 20_000, 890, 20_000)}.items():
r = beta_binomial_test(ca, na, cb, nb, seed=5)
decision = ("ship B" if r["expected_loss_ship_b"] < CARE_THRESHOLD
else "keep A" if r["expected_loss_keep_a"] < CARE_THRESHOLD
else "keep testing")
print(f" {label:14s} P(better) {r['p_b_better']:.3f}"
f" loss(ship B) {r['expected_loss_ship_b']:.6f}"
f" loss(keep A) {r['expected_loss_keep_a']:.6f} -> {decision}")
print("\n'Keep testing' is an explicit third option, which a p-value")
print("threshold does not give you.")
THE TWO PRINCIPLED ANSWERS TO PEEKING
SEQUENTIAL (frequentist) BAYESIAN
keeps the error-rate guarantee changes the question
adjusts the BOUNDARY reports a posterior
output: p-value, CI output: P(B better), credible interval,
expected loss
needs: a looking schedule OR needs: a prior
an always-valid method
cost: wider intervals cost: prior sensitivity when n is small
WHAT EACH CAN SAY
frequentist "if there were no effect, data this extreme would arise
4% of the time" (not P(effect exists))
Bayesian "there is a 93% probability the lift exceeds 2%"
(depends on the prior)
ALWAYS-VALID INTERVALS: THE PRICE OF OPTIONALITY
n fixed-horizon width always-valid width ratio
2,000 0.1240 0.1735 1.40x
10,000 0.0554 0.0721 1.30x
40,000 0.0277 0.0343 1.24x
Fixed intervals are valid at ONE prespecified n. Always-valid intervals
are valid at EVERY n at once, and are ~25-40% wider for it.
THE BAYESIAN PEEKING SUBTLETY
the POSTERIOR is a coherent summary at any n -> peeking is fine
a DECISION RULE ("stop when P(better) > 0.95") has a
frequentist error rate -> simulate it
flat prior, stop at P > 0.95 -> ~18% false "B wins"
informative prior, stop at P > 0.95 -> ~6% false "B wins"
EXPECTED LOSS AS THE DECISION RULE
loss(ship B) = E[max(0, rate_A - rate_B)] cost if B is actually worse
loss(keep A) = E[max(0, rate_B - rate_A)] cost of missing a real gain
ship B when loss(ship B) < tolerance
keep A when loss(keep A) < tolerance
keep testing when NEITHER is below tolerance <- the honest third option
Key Takeaways
- Sequential methods keep frequentist guarantees by widening boundaries, and always-valid confidence sequences stay correct at every moment for roughly 25 to 40% wider intervals.
- Bayesian testing reports the probability the treatment is better and the expected loss, which map directly onto decisions.
- A Bayesian posterior is coherent under continuous monitoring, but a Bayesian stopping rule still has a frequentist error rate that must be simulated.
- Priors matter most at small samples, and a prior built from historical experiments correctly shrinks extreme early estimates.
- Expected loss gives an explicit "keep testing" option that a significance threshold does not.
🧪 Practice
- Implement a Beta-Binomial A/B test reporting the probability the treatment is better, a credible interval, and expected loss under both decisions.
- Run the same small-sample data under a flat prior and a historically informed prior, and compare the resulting probabilities.
- Simulate a "stop when the posterior probability exceeds 0.95" rule under no true effect and measure the false-winner rate.
- Compare fixed-horizon and always-valid interval widths at three sample sizes and quantify the cost of optionality.
- Interview: A vendor claims their Bayesian platform means you can stop any time with no statistical penalty. What is your follow-up question? (Hint: distinguish reporting a posterior from acting on a threshold rule.)
Novelty and Network Effects
Two classes of effect break the assumption that a short measurement predicts the long-run outcome. Both produce experiments that read well and ship changes that do not deliver, and both are largely invisible unless you specifically look.
Novelty and primacy effects are about users reacting to change as change.
A novelty effect is a temporary boost from unfamiliarity. Users notice the new button, explore it, and engagement rises. Two weeks later the novelty is gone and the metric returns to baseline. The measured effect was real and temporary, and the conclusion drawn from it was wrong.
A primacy effect is the reverse: existing users are disrupted by a change to a workflow they had learned, so the metric dips initially and then recovers or improves as they adapt. A short test rejects a genuinely good change.
The diagnostic for both is the same and it is simple: plot the effect over time. A stable effect is flat. A novelty effect decays toward zero. A primacy effect starts negative and trends up. Complementing that, split by user tenure: novelty and primacy both act on existing users, so a change whose effect is large for existing users and absent for new ones is suspect, while the effect measured on new users only is the cleanest estimate of the steady state, because they have no prior expectation to be surprised by or attached to.
Network effects are about the value of a feature depending on how many others have it. A messaging feature is useless to the only person who has it. A test at 1% exposure measures a feature almost nobody's contacts can use, which systematically understates the effect. Worse, the relationship is non-linear, so you cannot extrapolate: a 1% test tells you very little about the 100% outcome.
This interacts with the interference problem from earlier. The standard mitigations are cluster randomization so that a treated user's contacts are also treated, and running at higher exposure levels to get closer to the eventual state.
A related pattern is market-level saturation, where a treatment that reallocates a fixed quantity looks positive in a partial rollout and nets out to nothing at full rollout. If a ranking change moves purchases from one seller to another, a seller-level metric improves in test and total sales do not change at launch.
import numpy as np
import pandas as pd
from scipy import stats
rng = np.random.default_rng(17)
# --- THREE TIME PROFILES that a two-week test cannot distinguish. ---
days = np.arange(1, 29)
profiles = {
"genuine (stable)": np.full(28, 2.0),
"novelty (decays)": 5.0 * np.exp(-days / 6),
"primacy (recovers)": 2.5 - 5.0 * np.exp(-days / 5),
}
rows = []
for name, effect in profiles.items():
for d, e in zip(days, effect):
rows.append({"profile": name, "day": d, "true_effect": e})
prof = pd.DataFrame(rows)
print("cumulative measured effect if you stop at day N:")
summary = []
for name, g in prof.groupby("profile"):
row = {"profile": name}
for stop in (7, 14, 21, 28):
row[f"day {stop}"] = round(g[g["day"] <= stop]["true_effect"].mean(), 2)
row["steady state"] = round(g[g["day"] > 21]["true_effect"].mean(), 2)
summary.append(row)
print(pd.DataFrame(summary).to_string(index=False))
print("\nAt day 7 all three look like winners of similar size. Only the")
print("first one still is at day 28. Novelty overstates by 2.4x at week 1.")
# --- AN ASCII VIEW ---
print("\ndaily effect by profile:")
for name, g in prof.groupby("profile"):
line = ""
for e in g["true_effect"]:
line += "#" if e > 3 else "+" if e > 1 else "." if e > -1 else "-"
print(f" {name:20s} |{line}|")
print(f" {'':20s} day 1 {'':24s} day 28")
print(" # >3 + 1..3 . -1..1 - <-1")
# --- THE DIAGNOSTIC: fit a decay and test whether the slope is real. ---
def novelty_diagnostic(daily_effects: np.ndarray) -> dict:
d = np.arange(1, len(daily_effects) + 1)
slope, intercept, r, p, se = stats.linregress(d, daily_effects)
first_wk, last_wk = daily_effects[:7].mean(), daily_effects[-7:].mean()
return {
"first_week": round(first_wk, 3), "last_week": round(last_wk, 3),
"ratio_last_to_first": round(last_wk / first_wk, 2) if first_wk else np.nan,
"trend_per_day": round(slope, 4), "trend_p": round(p, 5),
"verdict": ("stable" if p > 0.05
else "DECAYING: suspect novelty" if slope < 0
else "RISING: suspect primacy"),
}
print("\ndiagnostics on noisy observed versions of each profile:")
for name, g in prof.groupby("profile"):
observed = g["true_effect"].values + rng.normal(0, 0.5, 28)
d = novelty_diagnostic(observed)
print(f" {name:20s} {d['verdict']:28s}"
f" wk1 {d['first_week']:+.2f} -> wk4 {d['last_week']:+.2f}"
f" (p {d['trend_p']:.4f})")
# --- SPLIT BY TENURE: new users give the cleanest steady-state read. ---
n = 30_000
tenure_days = rng.exponential(120, n)
is_new = tenure_days < 14
arm = rng.random(n) < 0.5
# Novelty acts only on EXISTING users, who had an expectation to violate.
steady, novelty_bonus = 1.2, 4.0
effect = steady + novelty_bonus * (~is_new)
y = 50 + effect * arm + rng.normal(0, 10, n)
overall = y[arm].mean() - y[~arm].mean()
new_only = y[arm & is_new].mean() - y[~arm & is_new].mean()
existing = y[arm & ~is_new].mean() - y[~arm & ~is_new].mean()
print(f"\ntrue steady-state effect: {steady:+.2f}")
print(f" all users: {overall:+.3f} <- inflated by novelty")
print(f" existing users: {existing:+.3f} <- novelty concentrated here")
print(f" NEW users only: {new_only:+.3f} <- closest to the steady state")
print(f" (new users: {is_new.sum():,} of {n:,})")
print("New users have no prior expectation, so they cannot be surprised by")
print("the change or attached to the old one.")
# --- NETWORK EFFECTS: the effect DEPENDS on exposure, non-linearly. ---
def network_effect(exposure: float, direct: float = 0.5,
network_max: float = 3.0) -> float:
"""Direct value, plus network value that grows with adoption and
saturates. The key property is that it is NOT linear in exposure."""
return direct + network_max * (exposure ** 1.4)
print("\nmeasured effect of a messaging feature by exposure level:")
for exp in (0.01, 0.05, 0.10, 0.25, 0.50, 1.00):
e = network_effect(exp)
print(f" {exp:>5.0%} exposure -> effect {e:.3f}"
f" ({e/network_effect(1.0):.0%} of the full-rollout effect)")
print("\nA 1% test measures 17% of the eventual effect. Linear extrapolation")
print("from a small ramp is wrong in both directions depending on curvature.")
# --- CLUSTER RANDOMIZATION recovers the within-cluster network effect. ---
n_clusters, per_cluster = 60, 400
cluster = np.repeat(np.arange(n_clusters), per_cluster)
N = len(cluster)
# INDIVIDUAL randomization at 50%: a treated user's contacts are mostly
# untreated, so the network component barely activates.
ind_arm = rng.random(N) < 0.5
local_exposure_ind = np.repeat(
pd.Series(ind_arm).groupby(cluster).mean().values, per_cluster)
# CLUSTER randomization: within a treated cluster, exposure is 100%.
cl_on = rng.random(n_clusters) < 0.5
cl_arm = np.repeat(cl_on, per_cluster)
local_exposure_cl = np.repeat(cl_on.astype(float), per_cluster)
def outcome(arm, local_exp, r):
# A user's benefit depends on their OWN treatment and their cluster's
# adoption rate, which is exactly what a network effect means.
return 40 + (0.5 * arm + 3.0 * arm * local_exp ** 1.4
+ r.normal(0, 8, len(arm)))
y_ind = outcome(ind_arm, local_exposure_ind, rng)
y_cl = outcome(cl_arm, local_exposure_cl, rng)
est_ind = y_ind[ind_arm].mean() - y_ind[~ind_arm].mean()
cm = pd.DataFrame({"c": cluster, "y": y_cl, "arm": cl_arm}).groupby("c").agg(
y=("y", "mean"), arm=("arm", "first"))
est_cl = (cm.loc[cm["arm"], "y"].mean() - cm.loc[~cm["arm"], "y"].mean())
full_rollout = 0.5 + 3.0 * 1.0
print(f"\ntrue full-rollout effect: {full_rollout:.2f}")
print(f" individual randomization: {est_ind:.3f}"
f" ({est_ind/full_rollout:.0%} of it)")
print(f" cluster randomization: {est_cl:.3f}"
f" ({est_cl/full_rollout:.0%} of it)")
print("Treating whole clusters lets the network component activate, which is")
print("why social and marketplace tests are run geographically.")
# --- SATURATION: a reallocation that nets to zero at full rollout. ---
sellers = pd.DataFrame({
"seller": [f"s{i}" for i in range(6)],
"promoted": [True, True, True, False, False, False],
"sales_control": [100, 100, 100, 100, 100, 100],
})
# The ranking change moves demand toward promoted sellers; total is fixed.
sellers["sales_treat"] = [130, 128, 132, 70, 71, 69]
print("\nper-seller effect looks strong:")
print(f" promoted sellers: "
f"{sellers[sellers.promoted]['sales_treat'].mean()/100-1:+.0%}")
print(f" total sales: control {sellers['sales_control'].sum()}"
f" -> treated {sellers['sales_treat'].sum()}"
f" ({sellers['sales_treat'].sum()/sellers['sales_control'].sum()-1:+.1%})")
print("The promoted-seller metric improves and the marketplace total does not.")
print("Always pair a per-unit metric with the market-level total.")
THREE PROFILES, ONE WEEK OF DATA
effect
5 |##
| ## novelty: decays to nothing
| ###
2 |+++++++++++++++++++++++++++ genuine: flat
| ......
0 |- - - -
| - primacy: dips, then recovers
-2 |-
+---------------------------> day 1 ....................... 28
measured effect if you stop at:
day 7 day 14 day 21 day 28 steady state
genuine 2.00 2.00 2.00 2.00 2.00
novelty 3.35 2.28 1.62 1.25 0.10 <- 2.4x high at wk1
primacy -0.83 0.66 1.25 1.55 2.34 <- rejected at wk1
DIAGNOSTICS FOR NOVELTY / PRIMACY
1. plot the DAILY effect and fit a trend; a significant slope is the flag
2. compare week 1 against the final week
3. split by tenure: novelty and primacy act on EXISTING users
4. measure NEW users only -> closest to the steady state
NETWORK EFFECTS ARE NON-LINEAR IN EXPOSURE
exposure measured effect share of full-rollout effect
1% 0.51 15%
5% 0.57 16%
10% 0.62 18%
25% 0.92 26%
50% 1.64 47%
100% 3.50 100%
A small ramp measures almost none of the eventual effect, and linear
extrapolation from it is not valid.
WHICH MITIGATION FOR WHICH PROBLEM
problem symptom fix
---------------------------------------------------------------------------
novelty effect decays over days run longer; new-user cut
primacy effect rises over days run longer; new-user cut
network effects effect grows with exposure cluster randomize; ramp up
saturation/realloc. per-unit up, total flat add a market-level metric
Key Takeaways
- Novelty inflates short tests and primacy deflates them, and both are diagnosed by plotting the daily effect and testing its trend.
- The effect measured on new users only is the cleanest estimate of the steady state, because they have no prior expectation.
- Network effects make the measured effect depend non-linearly on exposure, so a small ramp cannot be extrapolated.
- Cluster randomization lets the network component activate, which is why social and marketplace tests run geographically.
- Pair every per-unit metric with a market-level total, or a pure reallocation reads as a win.
🧪 Practice
- Simulate stable, decaying, and recovering effect profiles and tabulate the measured effect at 7, 14, and 28 days for each.
- Write a novelty diagnostic that fits a trend to daily effects and flags a significant slope in either direction.
- Construct a dataset where the effect for existing users is large and for new users is small, and explain which estimate you would report.
- Model a network effect that grows with exposure, then compare estimates from individual and cluster randomization.
- Interview: A test of a new seller-promotion algorithm shows a 30% lift for promoted sellers. What do you check before recommending launch? (Hint: where did that demand come from, and what happened to everyone else?)
Interpreting and Acting on Results
An experiment produces an estimate with uncertainty. Turning that into a decision requires translating statistics into business terms and resisting several specific temptations. Most bad decisions from good experiments happen at this step.
Start with what a result actually says.
Statistical significance is not practical significance. With enough traffic, a 0.02% lift is significant and worth nothing. The question is always whether the effect is large enough to justify the cost of shipping and maintaining the change. That comparison requires a number from the business side, which is why the minimum detectable effect should have come from there in the first place.
Read the confidence interval, not the p-value. The interval tells you the range of effects consistent with your data, which is exactly what a decision needs. A result of "+2% significant" and one of "+2% with an interval from +0.1% to +3.9%" are the same result, and only the second makes clear that the effect might be negligible.
A non-significant result is not evidence of no effect. It means you could not distinguish the effect from zero at your sample size. Whether that is informative depends entirely on the interval. If the interval is narrow and centred near zero, you have genuinely learned the effect is small. If it spans from -5% to +7%, you have learned nothing, and reporting "no effect" is simply false. The honest phrasing names what was ruled out.
Then the temptations.
Segment fishing. When the overall result disappoints, the instinct is to look for a segment where it worked. With twenty segments, one will look significant by chance. Pre-register the segments you will examine, correct for multiplicity, and treat any unplanned segment finding as a hypothesis for a new test rather than a result.
Ignoring guardrails to ship a primary win. The decision rule agreed before the test exists precisely for this moment.
Over-trusting a single test. Published effect sizes in experimentation tend to shrink on replication, and the winner's curse means that among many tested ideas, the ones that look best are disproportionately those that got lucky. An effect barely above the significance threshold should be expected to be smaller than measured.
Finally, the null result is a real outcome and usually the most common one. Most tested ideas do not work. A programme where 70% of tests are null is functioning correctly; one where 80% of tests win is broken and you should look for the bug. The value of an experiment that says no is that it prevented you from shipping and maintaining something that did nothing.
import numpy as np
import pandas as pd
from scipy import stats
rng = np.random.default_rng(19)
# --- SIGNIFICANCE VS PRACTICAL IMPORTANCE ---
def result_frame(n_per_arm, p0, observed_rel_lift, value_per_conv,
annual_conversions, build_cost, annual_maintenance):
p1 = p0 * (1 + observed_rel_lift)
se = np.sqrt(p0*(1-p0)/n_per_arm + p1*(1-p1)/n_per_arm)
diff = p1 - p0
z = diff / se
p = 2 * (1 - stats.norm.cdf(abs(z)))
lo, hi = diff - 1.96*se, diff + 1.96*se
# Translate the whole interval into money: this is what decides it.
to_money = lambda d: d / p0 * annual_conversions * value_per_conv
return {
"n_per_arm": f"{n_per_arm:,}",
"rel_lift": f"{observed_rel_lift:+.2%}",
"p_value": round(p, 4),
"significant": p < 0.05,
"rel_ci": f"[{lo/p0:+.2%}, {hi/p0:+.2%}]",
"annual_value": f"${to_money(diff):,.0f}",
"annual_value_ci": f"[${to_money(lo):,.0f}, ${to_money(hi):,.0f}]",
"cost_year_1": f"${build_cost + annual_maintenance:,.0f}",
"ship": ("yes" if p < 0.05 and to_money(lo) > annual_maintenance
else "no"),
}
cases = {
"huge test, tiny effect": dict(n_per_arm=4_000_000, p0=0.05,
observed_rel_lift=0.004),
"well-sized, real effect": dict(n_per_arm=60_000, p0=0.05,
observed_rel_lift=0.045),
"underpowered, big point estimate": dict(n_per_arm=3_000, p0=0.05,
observed_rel_lift=0.12),
}
COMMON = dict(value_per_conv=40.0, annual_conversions=2_000_000,
build_cost=120_000, annual_maintenance=60_000)
rows = [{"case": k, **result_frame(**v, **COMMON)} for k, v in cases.items()]
print(pd.DataFrame(rows).to_string(index=False))
print("\nCase 1 is significant and worth less than its maintenance cost.")
print("Case 3 has the largest point estimate and an interval that includes")
print("outcomes far worse than the cost. Significance decided neither.")
# --- THE FOUR-WAY READING of a result: interval position vs a threshold. ---
def interpret(est, lo, hi, practical_threshold):
"""Compare the whole interval against zero AND against the threshold."""
if lo > practical_threshold:
return "SHIP: effect confidently exceeds the bar"
if hi < 0:
return "DO NOT SHIP: confidently harmful"
if lo > 0 and hi < practical_threshold:
return "REAL BUT TOO SMALL: confidently positive, below the bar"
if lo < 0 and hi > practical_threshold:
return "INCONCLUSIVE: interval spans harm and meaningful gain"
if abs(lo) < practical_threshold and abs(hi) < practical_threshold:
return "PRECISE NULL: genuinely learned the effect is small"
return "INCONCLUSIVE: need more data"
THRESHOLD = 0.01
print(f"\ninterpretation with a practical bar of {THRESHOLD:+.0%}:")
for label, (est, lo, hi) in {
"clear win": (0.030, 0.020, 0.040),
"real but small": (0.004, 0.002, 0.006),
"precise null": (0.001, -0.004, 0.006),
"wide and vague": (0.020, -0.030, 0.070),
"clear harm": (-0.025, -0.040, -0.010)}.items():
print(f" {label:16s} {est:+.3f} [{lo:+.3f}, {hi:+.3f}]"
f" -> {interpret(est, lo, hi, THRESHOLD)}")
print("\nNote that 'precise null' and 'wide and vague' are BOTH")
print("non-significant and mean completely different things.")
# --- SEGMENT FISHING: what 20 null segments look like. ---
n = 60_000
arm = rng.random(n) < 0.5
y = rng.normal(0, 1, n) # NO effect anywhere
segments = {f"seg_{i}": rng.integers(0, 20, n) for i in ["a"]}
seg = segments["seg_a"]
rows = []
for s in range(20):
mask = seg == s
a, b = y[mask & arm], y[mask & ~arm]
if len(a) < 30 or len(b) < 30:
continue
t = stats.ttest_ind(a, b, equal_var=False)
rows.append({"segment": s, "n": int(mask.sum()),
"effect": round(a.mean()-b.mean(), 4),
"p": round(t.pvalue, 4)})
sf = pd.DataFrame(rows).sort_values("p")
print(f"\n20 segments, TRUE effect is zero everywhere:")
print(sf.head(4).to_string(index=False))
print(f" segments with p < 0.05: {(sf['p'] < 0.05).sum()} (expected ~1)")
from statsmodels.stats.multitest import multipletests
rej, p_adj, _, _ = multipletests(sf["p"], alpha=0.05, method="fdr_bh")
print(f" surviving Benjamini-Hochberg correction: {rej.sum()}")
print("The best-looking segment is the winner's curse, not a finding.")
# --- THE WINNER'S CURSE: selected effects are overstated. ---
n_ideas, n_per_arm = 300, 25_000
true_effects = rng.normal(0.001, 0.004, n_ideas) # most ideas do ~nothing
se = np.sqrt(2 * 0.05 * 0.95 / n_per_arm)
measured = true_effects + rng.normal(0, se, n_ideas)
significant = np.abs(measured) > 1.96 * se
winners = measured > 1.96 * se
print(f"\n{n_ideas} tested ideas, true effects ~ N(0.001, 0.004):")
print(f" declared winners: {winners.sum()}")
print(f" mean MEASURED effect among winners: {measured[winners].mean():.5f}")
print(f" mean TRUE effect among winners: {true_effects[winners].mean():.5f}")
print(f" overstatement: "
f"{measured[winners].mean()/true_effects[winners].mean():.2f}x")
print(f" winners whose true effect is NEGATIVE: "
f"{(true_effects[winners] < 0).sum()}")
print("Expect the launched effect to be smaller than the measured one, and")
print("size the business case on the interval's lower bound, not the point.")
# --- A WIN RATE THAT IS TOO GOOD IS A BUG ---
print("\nexperiment programme health:")
for label, win_rate in [("typical mature programme", 0.20),
("optimistic but plausible", 0.35),
("suspicious", 0.80)]:
note = "investigate: peeking or leakage" if win_rate > 0.5 else "plausible"
print(f" {label:26s} {win_rate:.0%} of tests win {note}")
# --- THE DECISION RECORD: what you write down when the test ends. ---
RECORD = {
"experiment": "checkout_autofill_v1",
"ran": "2024-04-01 to 2024-04-28 (4 whole weeks)",
"primary": "checkout_completion_rate",
"result": "+2.4% relative, 95% CI [+0.9%, +3.9%], p = 0.0018",
"practical_bar": "+1.0% relative (agreed 2024-03-25)",
"guardrails": "all 4 passed non-inferiority; p95 load +6 ms (bar +20 ms)",
"validity_checks": "SRM p=0.42; attrition 2.1% vs 2.0%; balance all |d|<0.04",
"novelty_check": "week 1 +2.6%, week 4 +2.3%, trend p=0.51 -> stable",
"decision": "SHIP to 100%",
"expected_annual_value": "$1.9M (lower bound $0.7M) vs $60k maintenance",
"follow_up": "Re-measure in the Q3 holdout; watch refund rate for 30 days.",
"what_would_change_this": "A refund-rate regression above 5% relative.",
}
print("\nDECISION RECORD")
for k, v in RECORD.items():
print(f" {k:24s} {v}")
THE FOUR-WAY READING OF AN INTERVAL
0 practical bar
| |
ship | | |---o---| whole CI above the bar
real, too small| |--o--| | above 0, below the bar
precise null |-o-| | tight around 0: LEARNED
inconclusive |--------o----------| spans harm and real gain
clear harm |--o--| | | whole CI below 0
"precise null" and "inconclusive" are both p > 0.05 and mean opposite
things. The p-value cannot distinguish them; the interval can.
SIGNIFICANCE IS NOT IMPORTANCE
case n/arm lift p ship?
---------------------------------------------------------------------------
huge test, tiny effect 4,000,000 +0.40% 0.0000 NO: value < upkeep
well-sized, real effect 60,000 +4.5% 0.0001 YES
underpowered, big estimate 3,000 +12.0% 0.21 NO: CI spans -10%..+35%
THE WINNER'S CURSE
300 ideas tested, most with near-zero true effects
declared winners: 23
mean MEASURED effect among them: 0.0089
mean TRUE effect among them: 0.0041 <- 2.2x overstated
winners whose true effect is < 0: 3
Selection on a noisy estimate selects for luck. Build the business case
on the interval's LOWER bound.
PROGRAMME HEALTH
win rate 10-30% normal and healthy: most ideas do not work
win rate 30-45% optimistic; check the metric is not gameable
win rate > 50% SUSPICIOUS: look for peeking, SRM, or leakage
Key Takeaways
- Statistical significance says an effect is detectable, not that it is worth shipping; the practical bar comes from the business case.
- Read the confidence interval and translate its whole range into business terms, because the lower bound is what the decision rests on.
- A tight interval around zero is a genuine finding while a wide one is no finding, and both are non-significant.
- Searching segments after a disappointing result will find a winner by chance, so pre-register segments and correct for multiplicity.
- Selected winners overstate their effects, and a programme where most tests win has a bug rather than a great team.
🧪 Practice
- Write a function converting a result's whole confidence interval into annual business value, and apply it to a significant but tiny effect.
- Implement the four-way interval interpretation and run it on a precise null and a wide inconclusive result.
- Search twenty segments of a zero-effect simulation for significance, then apply a false discovery rate correction.
- Simulate 300 tested ideas with mostly negligible true effects and quantify how much the declared winners overstate their effects.
- Interview: Your test shows +1.2% on the primary metric with p equal to 0.04, and one guardrail shows a 3% rise in support contacts with p equal to 0.09. What do you recommend? (Hint: what does the guardrail's interval include, and what did you agree to before the test?)
Common Pitfalls and Invalid Tests
Some experiment failures are statistical subtleties and others are outright bugs that make a result meaningless. The second kind is more common, and the discipline that catches it is a fixed checklist run before anyone looks at the effect. The order matters: several of these checks invalidate the test entirely, so there is no point interpreting an effect that fails them.
Sample ratio mismatch is the single most valuable check. If you intended a 50/50 split and observe 50.4/49.6 on a million users, something dropped or misassigned users non-randomly, and whatever did that is probably correlated with the outcome. The usual causes are a redirect that loses slow connections, a bot filter applied to only one arm, a crash in one variant, or assignment logged before exposure was verified. A mismatch means the arms are no longer comparable, so the effect estimate is not merely noisy, it is wrong.
Assignment and exposure confusion is subtler. If you randomize everyone who visits the site but the feature only appears on the checkout page, including all visitors dilutes the effect toward zero. Filtering to only those who saw the feature seems like the fix and introduces bias, because whether a user reached checkout can itself be affected by treatment. The correct approach is to trigger assignment at the point of exposure, so the population is defined identically in both arms by something treatment cannot influence.
Metric definition drift during a test invalidates the comparison. If instrumentation changes on day five, the before and after are different metrics.
Multiple testing across metrics, segments, and simultaneous experiments inflates false positives, as covered above.
Simpson's paradox appears when the arms have different composition and the composition correlates with the outcome, which is one more reason to check balance.
Survivorship and selection in the analysis. Restricting to users who completed an action, stayed active, or made a purchase conditions on something treatment affected, which breaks the randomization.
Twyman's law is the operational heuristic that ties all of this together: any figure that looks interesting or unusual is probably wrong. A 40% lift is far more likely to be a tracking bug than a breakthrough, and the correct response to an amazing result is to look for the defect first.
import numpy as np
import pandas as pd
from scipy import stats
rng = np.random.default_rng(23)
# --- A VALIDITY CHECKLIST that runs BEFORE any effect is interpreted. ---
def validate_experiment(df: pd.DataFrame, arm_col="arm", outcome_col="y",
expected_share=0.5, covariates=None,
exposed_col=None) -> pd.DataFrame:
checks = []
n_a = int((df[arm_col] == True).sum())
n_b = int((df[arm_col] == False).sum())
total = n_a + n_b
# 1. SAMPLE RATIO MISMATCH -- fatal if it fails.
exp_a = total * expected_share
chi2 = ((n_a-exp_a)**2/exp_a + (n_b-(total-exp_a))**2/(total-exp_a))
p_srm = 1 - stats.chi2.cdf(chi2, 1)
checks.append({"check": "sample ratio mismatch", "stat": f"p={p_srm:.2e}",
"detail": f"{n_a:,} vs {n_b:,} (share {n_a/total:.4f})",
"severity": "FATAL", "passed": p_srm > 0.001})
# 2. MISSING OUTCOMES BY ARM -- differential attrition.
miss_a = df.loc[df[arm_col], outcome_col].isna().mean()
miss_b = df.loc[~df[arm_col], outcome_col].isna().mean()
checks.append({"check": "differential attrition",
"stat": f"{miss_a:.3f} vs {miss_b:.3f}",
"detail": f"gap {abs(miss_a-miss_b):.4f}",
"severity": "FATAL", "passed": abs(miss_a-miss_b) < 0.005})
# 3. COVARIATE BALANCE on pre-treatment variables.
for c in (covariates or []):
a, b = df.loc[df[arm_col], c], df.loc[~df[arm_col], c]
pooled = np.sqrt((a.var(ddof=1)+b.var(ddof=1))/2)
d = (a.mean()-b.mean())/pooled if pooled else 0.0
checks.append({"check": f"balance: {c}", "stat": f"d={d:+.3f}",
"detail": f"{a.mean():.2f} vs {b.mean():.2f}",
"severity": "WARN", "passed": abs(d) < 0.1})
# 4. EXPOSURE PARITY: was the feature actually reachable in both arms?
if exposed_col:
ea = df.loc[df[arm_col], exposed_col].mean()
eb = df.loc[~df[arm_col], exposed_col].mean()
checks.append({"check": "exposure/trigger parity",
"stat": f"{ea:.3f} vs {eb:.3f}",
"detail": "trigger rates should match if triggering "
"precedes treatment",
"severity": "FATAL", "passed": abs(ea-eb) < 0.01})
# 5. OUTCOME VARIANCE PARITY: very different spreads suggest the arms
# are measuring different things.
va = df.loc[df[arm_col], outcome_col].var(ddof=1)
vb = df.loc[~df[arm_col], outcome_col].var(ddof=1)
ratio = max(va, vb) / max(min(va, vb), 1e-12)
checks.append({"check": "outcome variance ratio", "stat": f"{ratio:.2f}",
"detail": "a large ratio hints at a logging difference",
"severity": "WARN", "passed": ratio < 1.5})
return pd.DataFrame(checks)
# A healthy experiment.
n = 100_000
arm = rng.random(n) < 0.5
healthy = pd.DataFrame({
"arm": arm,
"tenure": rng.exponential(90, n),
"pre_spend": rng.gamma(2, 30, n),
"exposed": rng.random(n) < 0.4,
"y": 50 + 1.2*arm + rng.normal(0, 10, n),
})
print("HEALTHY EXPERIMENT")
print(validate_experiment(healthy, covariates=["tenure", "pre_spend"],
exposed_col="exposed").to_string(index=False))
# A broken one: a redirect in the treatment arm drops slow, low-value users.
broken = healthy.copy()
drop = broken["arm"] & (broken["pre_spend"] < 20) & (rng.random(n) < 0.35)
broken = broken[~drop].copy()
print(f"\nBROKEN EXPERIMENT (a redirect dropped {drop.sum():,} treated users)")
print(validate_experiment(broken, covariates=["tenure", "pre_spend"],
exposed_col="exposed").to_string(index=False))
true_effect = 1.2
biased = (broken.loc[broken["arm"], "y"].mean()
- broken.loc[~broken["arm"], "y"].mean())
print(f"\nreported effect {biased:+.3f} vs true {true_effect:+.3f}"
f" -- inflated {biased/true_effect:.1f}x by losing low-value users")
print("SRM and the pre_spend imbalance both caught this. Without them the")
print("result reads as a large, clean win.")
# --- TRIGGERING: dilution versus bias, and the correct fix. ---
m = 200_000
arm2 = rng.random(m) < 0.5
# Only 8% of visitors reach checkout, where the feature lives.
reaches_checkout = rng.random(m) < 0.08
effect_on_exposed = 2.0
y2 = np.where(reaches_checkout, 30 + effect_on_exposed*arm2*reaches_checkout, 5)
y2 = y2 + rng.normal(0, 6, m)
all_visitors = y2[arm2].mean() - y2[~arm2].mean()
exposed_only = (y2[arm2 & reaches_checkout].mean()
- y2[~arm2 & reaches_checkout].mean())
print(f"\ntrue effect on EXPOSED users: {effect_on_exposed:+.2f}")
print(f" all visitors (diluted): {all_visitors:+.3f}"
f" ({all_visitors/effect_on_exposed:.0%} of it)")
print(f" exposed only (correct here, because reaching checkout is")
print(f" NOT affected by treatment): {exposed_only:+.3f}")
# Now make triggering itself treatment-dependent: filtering becomes BIASED.
triggers = reaches_checkout | (arm2 & (rng.random(m) < 0.03))
y3 = np.where(triggers, 30 + effect_on_exposed*arm2, 5) + rng.normal(0, 6, m)
filtered = y3[arm2 & triggers].mean() - y3[~arm2 & triggers].mean()
print(f"\nwhen treatment CHANGES who triggers:")
print(f" trigger rate: treated {triggers[arm2].mean():.4f}"
f" vs control {triggers[~arm2].mean():.4f} <- not equal")
print(f" filtered estimate {filtered:+.3f} vs true {effect_on_exposed:+.2f}")
print("The exposure-parity check is what catches this. The fix is to trigger")
print("assignment on an event treatment cannot influence.")
# --- SELECTION IN THE ANALYSIS: conditioning on a post-treatment outcome. ---
k = 60_000
arm3 = rng.random(k) < 0.5
# Treatment increases purchasing, and purchasers spend differently.
purchased = rng.random(k) < (0.20 + 0.06*arm3)
spend = np.where(purchased, rng.gamma(3, 25, k), 0.0)
overall = spend[arm3].mean() - spend[~arm3].mean()
among_purchasers = (spend[arm3 & purchased].mean()
- spend[~arm3 & purchased].mean())
print(f"\nspend per USER (all users): {overall:+.3f} <- unbiased")
print(f"spend per PURCHASER (conditioned): {among_purchasers:+.3f} <- biased")
print(f" purchase rate: {purchased[arm3].mean():.3f}"
f" vs {purchased[~arm3].mean():.3f}")
print("Treatment changed WHO purchases, so purchasers are no longer")
print("comparable groups. Report per-user, or decompose rate x value.")
# --- TWYMAN'S LAW: an implausible result is probably a bug. ---
def plausibility_screen(rel_lift: float, metric: str,
historical_max: float = 0.10) -> str:
if abs(rel_lift) > 5 * historical_max:
return f"IMPLAUSIBLE ({rel_lift:+.0%}): look for a tracking bug first"
if abs(rel_lift) > 2 * historical_max:
return f"SUSPICIOUS ({rel_lift:+.0%}): verify instrumentation"
return f"plausible ({rel_lift:+.1%})"
print("\nTwyman's law screen (historical max lift on this metric: 10%):")
for lift in (0.03, 0.09, 0.28, 0.85):
print(f" {plausibility_screen(lift, 'conversion')}")
print("\n'Any figure that looks interesting is probably wrong.' The correct")
print("response to a spectacular result is to hunt for the defect.")
# --- THE FULL PITFALL TABLE ---
pitfalls = pd.DataFrame([
("sample ratio mismatch", "FATAL", "chi-squared on the split",
"fix the cause and rerun"),
("differential attrition", "FATAL", "missing-outcome rate by arm",
"investigate; intention-to-treat"),
("assignment != exposure", "FATAL", "trigger-rate parity",
"trigger at the exposure point"),
("metric definition changed", "FATAL", "instrumentation changelog",
"restart the test"),
("conditioning on post-treatment", "FATAL", "is the filter affected?",
"analyse per randomized unit"),
("peeking / early stopping", "FATAL", "was the stop rule prespecified?",
"sequential methods"),
("interference between units", "SERIOUS", "cluster-level vs individual",
"cluster randomize"),
("novelty / primacy", "SERIOUS", "daily effect trend",
"run longer; new-user cut"),
("multiple metrics / segments", "SERIOUS", "count the comparisons",
"one primary; correct the rest"),
("wrong analysis unit", "SERIOUS", "randomization vs analysis unit",
"aggregate to the randomized unit"),
("covariate imbalance", "WARN", "standardised differences",
"post-stratify or regress"),
("outlier domination", "WARN", "effect with and without capping",
"cap or use a robust metric"),
], columns=["pitfall", "severity", "how to detect", "what to do"])
print("\n", pitfalls.to_string(index=False))
THE CHECK ORDER (stop at the first FATAL failure)
1. sample ratio mismatch -> arms not comparable. STOP.
2. differential attrition -> randomization undone. STOP.
3. exposure / trigger parity -> filtering is biased. STOP.
4. metric definition stability -> comparing two metrics. STOP.
5. covariate balance -> post-stratify
6. novelty trend -> extend or use the new-user cut
7. analysis unit == randomization unit
8. ONLY NOW look at the effect
WHY SRM IS FATAL RATHER THAN NOISY
intended 50/50, observed 50,000 vs 46,800 on 96,800 users
Something removed ~3,200 treated users. Whatever did that
(a redirect, a crash, a bot filter) selected on SPEED or DEVICE or
VALUE, all of which correlate with the outcome. The surviving arms
differ by more than the treatment.
In the worked example above: true effect +1.2, reported +2.7 (2.2x).
DILUTION VS BIAS IN TRIGGERING
feature lives on checkout; 8% of visitors get there
analyse ALL visitors effect diluted to ~8% of its size (unbiased)
analyse exposed only, when
triggering is NOT affected by treatment (unbiased)
analyse exposed only, when
treatment CHANGES who triggers (BIASED)
The fix is structural: trigger assignment on an event that treatment
cannot influence, and check trigger-rate parity.
TWYMAN'S LAW
"Any figure that looks interesting or different is usually wrong."
historical max lift on this metric: 10%
+3% plausible
+9% plausible
+28% SUSPICIOUS: verify instrumentation
+85% IMPLAUSIBLE: find the tracking bug
Key Takeaways
- Run a fixed validity checklist before looking at the effect, and treat sample ratio mismatch as fatal rather than as noise.
- Differential attrition and treatment-dependent triggering both undo randomization inside the analysis.
- Trigger assignment at the point of exposure so the analysed population is defined identically in both arms.
- Conditioning on any post-treatment outcome, such as analysing only purchasers, reintroduces selection bias.
- By Twyman's law, a spectacular result is more likely a tracking defect than a breakthrough, so look for the bug first.
🧪 Practice
- Write a validity checklist function covering sample ratio mismatch, differential attrition, covariate balance, and exposure parity.
- Simulate a redirect that drops low-value users from one arm, and show both the check that catches it and the size of the resulting bias.
- Demonstrate the difference between dilution and bias by making triggering independent of treatment and then dependent on it.
- Show that analysing spend among purchasers only is biased while spend per user is not.
- Interview: An experiment reports a 60% lift in signups and the team wants to ship immediately. What are your first three actions? (Hint: what does Twyman's law suggest, and which two checks would you run before discussing the effect at all?)
<a id="83-quasi-experiments-and-observational-studies"></a>
8.3 Quasi-Experiments and Observational Studies
Many important questions cannot be randomized, and the answer is not to give up on causality but to find a source of variation that behaves as if it were random; this section covers how to recognise one and how to defend a conclusion built on it.
When Randomization Is Impossible
Randomization is the best available tool and it is frequently unavailable. The reasons fall into recognisable categories, and identifying which one you face determines which alternative applies.
Ethical barriers. You cannot randomly assign people to smoke, to lose their job, or to receive a treatment you believe is inferior. Where a treatment is already known to help, withholding it is not acceptable.
Practical impossibility. Some treatments cannot be assigned at the individual level at all. You cannot give half a country a different interest rate, run a brand advertising campaign that half the market cannot see, or apply a new tax code to a random subset of firms.
Political and organisational barriers. Leadership has already decided to launch. Legal requires the change everywhere. A partner contract specifies uniform treatment. These are extremely common in industry and no amount of statistical argument changes them.
Timing. The change already happened. You are asked to evaluate a launch from six months ago, and there was no experiment.
Insufficient scale. With eight enterprise customers, randomization is technically possible and statistically useless.
When randomization is out, you are looking for something else that assigns treatment in a way unrelated to the outcome. The strategies form a rough hierarchy of credibility.
At the top sit natural experiments, where something outside the system assigned treatment haphazardly: a policy boundary, a lottery, a staggered rollout, an arbitrary eligibility cutoff. These are the strongest non-randomized designs because the assignment mechanism is genuinely external.
Next are the panel methods that exploit before-and-after variation across units: difference-in-differences and its generalisation, synthetic control. Chapter 3 covered difference-in-differences in detail; the key requirement is that treated and untreated units would have moved in parallel absent the treatment.
Then come selection-on-observables methods: regression adjustment, propensity score matching, and inverse probability weighting, all covered in chapter 3. These rest on the strong and untestable assumption that you have measured every confounder. They are the weakest of the credible options, and their honest framing is "the effect among units that look comparable on what we recorded".
Instrumental variables and regression discontinuity, also in chapter 3, sit higher when a valid instrument or a sharp threshold genuinely exists, because they need weaker assumptions about unmeasured confounding.
The essential discipline in all of this is to state the identifying assumption explicitly and then attack it. Every method above answers the question "what would have happened otherwise" by assuming something you cannot verify. Naming that assumption in one sentence, and then describing what evidence would undermine it, is the difference between a credible causal claim and a regression with a hopeful interpretation.
import numpy as np
import pandas as pd
import statsmodels.api as sm
rng = np.random.default_rng(29)
# --- A DESIGN LADDER: which method, given what you have. ---
ladder = pd.DataFrame([
("randomized experiment", 1, "assignment is random by construction",
"nothing to assume", "always, when possible"),
("natural experiment", 2, "an external event assigned treatment",
"the event is unrelated to the outcome", "policy borders, lotteries"),
("regression discontinuity", 2, "treatment flips at a sharp cutoff",
"nothing else jumps at the cutoff", "eligibility thresholds"),
("instrumental variables", 3, "a variable moves treatment only",
"the instrument affects outcome ONLY via treatment", "a valid instrument"),
("difference-in-differences", 3, "treated and control move in parallel",
"parallel trends absent treatment", "a staggered rollout"),
("synthetic control", 3, "a weighted donor pool reproduces the unit",
"pre-period fit implies post-period fit", "one treated unit, many donors"),
("propensity / matching", 4, "comparable units on observed covariates",
"NO unmeasured confounders (untestable)", "rich covariates only"),
("regression adjustment", 4, "control for confounders directly",
"NO unmeasured confounders, correct model", "a last resort"),
("before/after comparison", 5, "nothing changed but the treatment",
"no time trend, no seasonality, no other change", "almost never valid"),
], columns=["design", "credibility", "identifies the counterfactual by",
"key assumption", "when available"])
print(ladder.to_string(index=False))
print("\nCredibility 1 is strongest. The assumption column is the whole game:")
print("it is what a reviewer should attack, and what you should attack first.")
# --- WHY A NAIVE BEFORE/AFTER FAILS: the trend is not the treatment. ---
months = pd.date_range("2023-01-01", periods=24, freq="MS")
trend = np.linspace(100, 130, 24) # the business was growing
seasonal = 6 * np.sin(np.arange(24) / 12 * 2 * np.pi)
treatment_month = 12
true_effect = 4.0
y = trend + seasonal + true_effect * (np.arange(24) >= treatment_month) \
+ rng.normal(0, 2, 24)
ts = pd.DataFrame({"month": months, "y": y,
"post": np.arange(24) >= treatment_month})
naive = ts[ts["post"]]["y"].mean() - ts[~ts["post"]]["y"].mean()
print(f"\ntrue effect {true_effect:+.2f}")
print(f"naive before/after: {naive:+.3f}"
f" ({naive/true_effect:.1f}x too large)")
print("The pre-existing growth trend is counted as treatment effect.")
# Controlling for the trend recovers it, IF you model the trend correctly.
X = sm.add_constant(pd.DataFrame({
"t": np.arange(24),
"post": ts["post"].astype(int),
"sin": np.sin(np.arange(24)/12*2*np.pi),
"cos": np.cos(np.arange(24)/12*2*np.pi),
}))
fit = sm.OLS(ts["y"], X).fit()
print(f"with trend + seasonality controls: {fit.params['post']:+.3f}"
f" (se {fit.bse['post']:.3f})")
print("This works only because we KNOW the functional form. With real data")
print("you do not, which is why a control group beats a model of the trend.")
# --- THE VALUE OF A CONTROL GROUP: it absorbs the trend you cannot model. ---
n_units, n_months = 40, 24
unit = np.repeat(np.arange(n_units), n_months)
t = np.tile(np.arange(n_months), n_units)
treated_unit = np.repeat(rng.random(n_units) < 0.4, n_months)
# A shared shock nobody would have guessed the shape of.
common_shock = np.where(np.tile(np.arange(n_months), n_units) >= 15, 12.0, 0.0)
unit_level = np.repeat(rng.normal(100, 15, n_units), n_months)
post = t >= treatment_month
panel_y = (unit_level + 0.8*t + common_shock
+ true_effect * (treated_unit & post) + rng.normal(0, 3, len(unit)))
panel = pd.DataFrame({"unit": unit, "t": t, "treated": treated_unit,
"post": post, "y": panel_y})
# Treated-only before/after absorbs the common shock into the "effect".
treated_ba = (panel[panel.treated & panel.post]["y"].mean()
- panel[panel.treated & ~panel.post]["y"].mean())
# Difference-in-differences differences the shock out. (See chapter 3.5.)
cells = panel.groupby(["treated", "post"])["y"].mean()
did = ((cells[(True, True)] - cells[(True, False)])
- (cells[(False, True)] - cells[(False, False)]))
print(f"\ntreated-only before/after: {treated_ba:+.3f} <- absorbs the shock")
print(f"difference-in-differences: {did:+.3f} <- true {true_effect:+.2f}")
print("The control group experienced the same shock, so differencing removes")
print("it without anyone having to know it happened.")
# --- SCALE: when randomization is possible and useless. ---
print("\nminimum detectable effect with very few units (sd = 15, 80% power):")
for k in (8, 20, 60, 200):
per_arm = k / 2
mde = (1.959964 + 0.841621) * 15 * np.sqrt(2/per_arm)
print(f" {k:>3} units -> MDE {mde:6.2f}"
f" ({mde/100:5.1%} of a baseline of 100)"
f" {'useless' if mde > 15 else 'workable'}")
print("With 8 enterprise accounts you can randomize and detect nothing.")
print("Switchback or synthetic control are the realistic options there.")
# --- STATE AND ATTACK THE ASSUMPTION: the required discipline. ---
CLAIM = {
"question": "Did the March pricing change increase revenue per account?",
"design": "difference-in-differences, EU accounts treated, US as control",
"identifying_assumption":
"Absent the pricing change, EU and US revenue per account would have "
"followed parallel paths.",
"evidence_for": [
"24 pre-period months of parallel trends (placebo test p = 0.61)",
"no EU-specific product changes in the window (changelog reviewed)",
"the effect appears in the month of the change, not before",
],
"threats": [
"an EU macro shock coinciding with the change",
"a competitor's EU-only launch in the same month",
"EU accounts differ in size mix, so a shared shock hits them harder",
],
"what_would_change_the_conclusion":
"A pre-period placebo showing a fake effect, or a divergence in a "
"non-price metric at the same date.",
}
print(f"\nCAUSAL CLAIM RECORD")
for k, v in CLAIM.items():
if isinstance(v, list):
print(f" {k}:")
for item in v:
print(f" - {item}")
else:
print(f" {k}: {v}")
THE CREDIBILITY LADDER
1 randomized experiment assumption: none (by construction)
|
2 natural experiment an external event did the assigning
regression discontinuity nothing else jumps at the cutoff
|
3 instrumental variables the instrument acts ONLY through treatment
difference-in-differences parallel trends absent treatment
synthetic control pre-period fit implies post-period fit
|
4 propensity / matching NO unmeasured confounders (untestable)
regression adjustment NO unmeasured confounders + correct model
|
5 before / after nothing else changed <- almost never true
WHY BEFORE/AFTER FAILS
y | ____ <- post-period mean
| _____/
| _____/ | the RISE was already happening
|_____/ | treatment
+--------------|----------> t
true effect +4.0
naive before/after +13.6 the trend is counted as effect
WHAT A CONTROL GROUP BUYS YOU
a shared shock hits both groups at month 15
treated-only before/after: +14.1 (shock absorbed into the effect)
difference-in-differences: +4.1 (shock differenced out)
You did not need to know the shock existed. That is the point: a control
group handles confounders you never modelled, exactly as randomization
handles covariates you never measured.
SCALE MAKES RANDOMIZATION USELESS BEFORE IT MAKES IT IMPOSSIBLE
units MDE (sd=15) as % of a 100 baseline
8 21.0 21% useless
20 13.3 13% borderline
60 7.7 8% workable
200 4.2 4% fine
Key Takeaways
- Identify why randomization is unavailable, because ethical, practical, organisational, timing, and scale barriers point to different alternatives.
- Natural experiments and sharp thresholds are the most credible non-randomized designs because the assignment mechanism is external.
- Selection-on-observables methods rest on an untestable no-unmeasured- confounders assumption and are the weakest credible option.
- A control group removes shocks you never modelled, which is why difference-in-differences beats a before-and-after comparison.
- State the identifying assumption in one sentence and list what would undermine it, or the claim is a regression with a hopeful interpretation.
🧪 Practice
- Simulate a growing, seasonal series with a treatment partway through and show how much a naive before-and-after comparison overstates the effect.
- Add a common shock affecting both treated and untreated units, then compare a treated-only before-and-after against difference-in- differences.
- Compute the minimum detectable effect for 8, 20, and 60 units and identify where randomization stops being useful.
- Write a causal claim record for a question in your own work, naming the identifying assumption and three threats to it.
- Interview: Leadership launched a redesign to all users three months ago and now asks whether it worked. What do you propose? (Hint: what varies that you can still exploit, and what comparison group might exist without anyone having planned one?)
Natural Experiments
A natural experiment is a situation where something outside your control assigned treatment in a way that is plausibly unrelated to the outcome. You did not randomize, but the world did something close enough to it, and recognising these situations is one of the most valuable skills in applied causal work.
The reason they are so useful is that the assignment mechanism is external to the decisions of the units being studied. When a lottery decides who gets a school place, the winners and losers are comparable by construction, and no assumption about measured confounders is needed.
The recurring sources are worth knowing by name, because once you know them you start seeing them.
Lotteries and randomized allocation by others. Visa lotteries, school admission lotteries, draft lotteries, and randomised audits. These are genuine randomization performed by someone else.
Arbitrary thresholds. Eligibility at a specific age, income, score, or size. Just above and just below the line are comparable while the treatment differs sharply. This is the regression discontinuity design from chapter 3.
Geographic boundaries. A policy applies in one jurisdiction and not the adjacent one. Households on either side of a border are similar in most respects.
Staggered rollouts. A feature launched region by region over months gives you units treated at different times, which is the setting that difference-in-differences and its modern event-study variants are built for. In industry this is the most commonly available natural experiment, and often you can create it deliberately by choosing the rollout order arbitrarily, which is worth doing even when a proper test is refused.
Timing accidents. Someone born one day either side of a school cutoff, a contract signed the week before a policy change, an outage that affected a subset of users.
Technical accidents. A caching bug that served an old version to a hash-based subset of traffic, a partial deployment, a regional capacity failure. These are unplanned and sometimes genuinely random.
The entire credibility of the design rests on one question: is the assignment really unrelated to the outcome? The failure mode is sorting, where units choose their side of the line. If families move across a school district boundary to get the better school, the two sides are no longer comparable. If firms just below a size threshold stay small deliberately to avoid regulation, the threshold is manipulated.
The standard checks follow directly. Look for a discontinuity in the density of units around a threshold, which indicates manipulation. Check pre-treatment covariate balance across the boundary, exactly as you would in an experiment. Run placebo tests at fake thresholds where no treatment changed, and confirm nothing shows up. And check that no other policy changed at the same threshold, which is the most common way a clean-looking discontinuity is confounded.
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy import stats
rng = np.random.default_rng(31)
# --- A STAGGERED ROLLOUT: the natural experiment you usually have. ---
n_regions, n_weeks = 30, 40
# Regions adopt at different weeks. The staggering is what identifies the
# effect: at any moment some regions are treated and some are not.
adopt_week = rng.integers(10, 32, n_regions)
rows = []
region_level = rng.normal(100, 12, n_regions)
TRUE = 5.0
for r in range(n_regions):
for w in range(n_weeks):
treated = w >= adopt_week[r]
rows.append({
"region": r, "week": w, "adopt_week": adopt_week[r],
"treated": treated,
# region baseline + common time trend + a shared shock + effect
"y": (region_level[r] + 0.35 * w
+ (8.0 if w >= 25 else 0.0)
+ TRUE * treated + rng.normal(0, 3)),
})
panel = pd.DataFrame(rows)
# Two-way fixed effects: region effects absorb level differences, week
# effects absorb everything common to all regions at a given time.
fe = sm.OLS.from_formula("y ~ treated + C(region) + C(week)", panel).fit()
print(f"TRUE effect {TRUE:+.2f}")
print(f"two-way fixed effects: {fe.params['treated[T.True]']:+.3f}"
f" (se {fe.bse['treated[T.True]']:.3f})")
print("Week fixed effects removed the shock at week 25 without modelling it.")
# --- EVENT STUDY: the pre-trend check that validates the design. ---
panel["rel_week"] = panel["week"] - panel["adopt_week"]
window = panel[panel["rel_week"].between(-8, 8)].copy()
# Bin relative time, omitting -1 as the reference period.
window["bin"] = window["rel_week"].astype(int)
es = sm.OLS.from_formula(
"y ~ C(bin, Treatment(reference=-1)) + C(region) + C(week)", window).fit()
coefs = []
for k in range(-8, 9):
if k == -1:
coefs.append({"rel_week": k, "coef": 0.0, "se": 0.0})
continue
name = f"C(bin, Treatment(reference=-1))[T.{k}]"
if name in es.params:
coefs.append({"rel_week": k, "coef": es.params[name],
"se": es.bse[name]})
ev = pd.DataFrame(coefs)
print("\nevent study around adoption (reference = week -1):")
for _, r in ev.iterrows():
bar = "#" * max(int(round(r["coef"] * 2)), 0)
flag = " <- adoption" if r["rel_week"] == 0 else ""
print(f" {int(r['rel_week']):>3} {r['coef']:+6.2f}"
f" +/- {1.96*r['se']:4.2f} {bar}{flag}")
pre = ev[ev["rel_week"] < 0]
print(f"\npre-period coefficients: mean {pre['coef'].mean():+.3f},"
f" max |coef| {pre['coef'].abs().max():.3f}")
print("Flat before adoption and a step at zero is what a credible design")
print("looks like. A pre-trend means the treatment was not the cause.")
# --- A THRESHOLD DESIGN, and the manipulation check that validates it. ---
def make_threshold_data(n, manipulated: bool, cutoff=50.0, effect=6.0, seed=0):
r = np.random.default_rng(seed)
score = r.normal(50, 12, n)
if manipulated:
# Units just below the cutoff push themselves just above it.
nudge = (score > cutoff - 3) & (score < cutoff) & (r.random(n) < 0.7)
score = np.where(nudge, cutoff + r.random(n) * 2, score)
treated = score >= cutoff
y = 20 + 0.5 * score + effect * treated + r.normal(0, 4, n)
return pd.DataFrame({"score": score, "treated": treated, "y": y})
def density_check(score: np.ndarray, cutoff: float, width: float = 5.0) -> dict:
"""McCrary-style test: the DENSITY of units should be continuous at the
cutoff. A jump means units are sorting across it."""
just_below = ((score >= cutoff - width) & (score < cutoff)).sum()
just_above = ((score >= cutoff) & (score < cutoff + width)).sum()
total = just_below + just_above
# Binomial test against the 50/50 split a continuous density implies.
p = stats.binomtest(just_above, total, 0.5).pvalue
return {"below": just_below, "above": just_above,
"ratio": round(just_above / max(just_below, 1), 3),
"p_value": f"{p:.2e}",
"verdict": "MANIPULATED: units are sorting" if p < 0.01 else "ok"}
for label, manip in [("clean threshold", False), ("manipulated threshold", True)]:
d = make_threshold_data(6000, manip, seed=7)
print(f"\n{label}:")
print(f" density check: {density_check(d['score'].values, 50.0)}")
# Local linear fit on each side of the cutoff, within a bandwidth.
band = d[d["score"].between(40, 60)].copy()
band["above"] = (band["score"] >= 50).astype(int)
band["centred"] = band["score"] - 50
rd = sm.OLS.from_formula("y ~ above + centred + above:centred", band).fit()
print(f" estimated jump at the cutoff: {rd.params['above']:+.3f}"
f" (true 6.00, se {rd.bse['above']:.3f})")
print("\nThe manipulated case still 'estimates' an effect. The density check")
print("is what tells you not to believe it.")
# --- PLACEBO TESTS at fake thresholds: nothing should be there. ---
clean = make_threshold_data(12_000, False, cutoff=50.0, effect=6.0, seed=11)
print("\nplacebo cutoffs (only 50 is real):")
for fake in (35, 42, 50, 58, 65):
band = clean[clean["score"].between(fake-8, fake+8)].copy()
if len(band) < 200:
continue
band["above"] = (band["score"] >= fake).astype(int)
band["centred"] = band["score"] - fake
f = sm.OLS.from_formula("y ~ above + centred + above:centred", band).fit()
star = " <- the real cutoff" if fake == 50 else ""
print(f" cutoff {fake}: jump {f.params['above']:+6.3f}"
f" p {f.pvalues['above']:.3f}{star}")
print("Effects appearing at fake cutoffs would mean the specification, not")
print("the treatment, is producing the result.")
# --- COVARIATE CONTINUITY: pre-treatment variables must NOT jump. ---
clean["pre_metric"] = 12 + 0.2 * clean["score"] + rng.normal(0, 2, len(clean))
band = clean[clean["score"].between(40, 60)].copy()
band["above"] = (band["score"] >= 50).astype(int)
band["centred"] = band["score"] - 50
for col, expect in [("pre_metric", "should NOT jump"), ("y", "SHOULD jump")]:
f = sm.OLS.from_formula(f"{col} ~ above + centred + above:centred",
band).fit()
print(f"\n{col:12s} jump {f.params['above']:+.3f}"
f" p {f.pvalues['above']:.4f} ({expect})")
# --- A CHECKLIST for any natural experiment ---
CHECKS = pd.DataFrame([
("assignment is external", "who or what decided treatment?",
"if the unit chose, it is not natural"),
("no sorting / manipulation", "density continuity at the threshold",
"a jump in the count of units means sorting"),
("covariate continuity", "pre-treatment variables across the boundary",
"they must NOT jump where the treatment does"),
("no pre-trend", "event study coefficients before adoption",
"flat before, step at zero"),
("placebo thresholds", "fake cutoffs and fake dates",
"no effect should appear"),
("no co-occurring change", "changelog, policy calendar at the same date",
"the most common confounder of a clean jump"),
("bandwidth robustness", "re-estimate at several bandwidths",
"the estimate should not swing with the window"),
], columns=["check", "how", "what it tells you"])
print("\n", CHECKS.to_string(index=False))
SOURCES OF NATURAL EXPERIMENTS
source example design it supports
---------------------------------------------------------------------------
lottery by others visa, school place, audits direct comparison
arbitrary threshold age 65, score 50, 50 employees regression discontinuity
geographic boundary policy in one state only border comparison / DiD
STAGGERED ROLLOUT feature launched region by DiD / event study
region <- most common in industry
timing accident born either side of a cutoff discontinuity in date
technical accident partial deploy, cache bug quasi-random subset
THE EVENT STUDY IS THE VALIDITY TEST
coefficient relative to adoption week
+6 | ####################
+4 | #
+2 | #
0 |...................#.................. <- flat BEFORE adoption
-2 | |
+-------------------|------------------>
-8 0 +8
adoption
flat pre-period + a step at zero -> credible
a rising pre-period -> the treatment was not the cause;
something was already happening
THE MANIPULATION CHECK
clean density at the cutoff manipulated
####|#### ##|########
####|#### ##|########
----+---- --+--------
50 50
equal counts either side a pile-up just ABOVE the line
p = 0.72 -> ok p < 1e-9 -> units are sorting
In the worked example the manipulated data still yields a plausible
effect estimate. Only the density check reveals it is worthless.
WHAT MUST AND MUST NOT JUMP AT THE BOUNDARY
outcome SHOULD jump (that is the effect)
pre-treatment must NOT jump (or the groups differ already)
density of units must NOT jump (or units are sorting)
other policies must NOT change (or the effect is confounded)
Key Takeaways
- A natural experiment works because something external assigned treatment, so credibility rests on whether that assignment relates to the outcome.
- Staggered rollouts are the most commonly available natural experiment in industry, and the rollout order can be chosen arbitrarily on purpose.
- An event study showing a flat pre-period and a step at adoption is the central validity check for a staggered design.
- Test for manipulation with a density check at the threshold, since a manipulated cutoff still produces a plausible-looking estimate.
- Pre-treatment covariates and the density of units must not jump where the treatment does, and placebo thresholds should show nothing.
🧪 Practice
- Build a staggered rollout panel with a common shock and recover the effect with two-way fixed effects.
- Estimate an event study around adoption and confirm the pre-period coefficients are flat.
- Generate threshold data with and without manipulation, and show that a density check distinguishes them while the effect estimate does not.
- Run placebo discontinuity estimates at four fake cutoffs and confirm only the real one shows an effect.
- Interview: A colleague finds that firms with 49 employees are much more profitable than those with 51, and attributes it to a regulation applying at 50. What do you check first? (Hint: count the firms on each side of the line.)
Synthetic Control Methods
Sometimes only one unit is treated. A single country changes a law, one city gets a new service, your largest customer receives a bespoke change. With one treated unit, difference-in-differences needs a comparison unit, and no single untreated unit looks quite like the treated one. Choosing one by eye is arbitrary and invites the criticism that you picked the one that gave the answer you wanted.
The synthetic control method solves this elegantly: build a comparison unit as a weighted combination of untreated units, with the weights chosen so the combination tracks the treated unit closely before treatment. If a weighted blend of several other cities reproduced your city's trajectory for three years before the intervention, that blend is a credible estimate of what your city would have done afterwards.
The logic is worth stating precisely because it is the design's whole argument. Difference-in-differences assumes parallel trends and you verify it visually. Synthetic control constructs a comparison whose pre-treatment trend matches by design, and then the post-treatment divergence is the estimated effect. The identifying assumption becomes: a combination that tracked for a long pre-period would have continued to track.
The weights are constrained to be non-negative and to sum to one, which matters for two reasons. It prevents extrapolation, since the synthetic unit stays inside the range of the donors rather than being an unbounded linear combination. And it usually produces a sparse solution, where only a handful of donors get non-zero weight, which makes the comparison interpretable: you can name which units the synthetic one is made of.
Inference works differently here, since there is one treated unit and no conventional standard error. The standard approach is placebo permutation: pretend each donor unit was the treated one, compute its synthetic control effect, and compare the real effect against that distribution. If the treated unit's post-treatment gap is extreme relative to the placebo gaps, the effect is unlikely to be chance. The p-value is the rank of the real effect among the placebos, which means the number of donors bounds the smallest p-value you can report.
The practical requirements are specific. You need a long pre-period, since the credibility comes from the length of the match. You need enough donors who are plausibly similar and, importantly, untreated by the same or a related intervention. And you need to check that the pre-treatment fit is genuinely good, because a poor pre-period fit invalidates the whole argument.
The related modern variants are worth knowing by name: augmented synthetic control adds a bias-correction step when the pre-period fit is imperfect, and synthetic difference-in-differences combines both approaches.
import numpy as np
import pandas as pd
from scipy.optimize import nnls, minimize
rng = np.random.default_rng(37)
# --- A PANEL WITH ONE TREATED UNIT ---
n_donors, T, T0 = 18, 40, 28 # T0 = last pre-treatment period
TRUE = 7.0
# Two latent factors drive everyone; units differ in their loadings. This is
# the structure that makes a weighted blend able to reproduce a unit.
f1 = np.cumsum(rng.normal(0.4, 1.0, T))
f2 = np.sin(np.arange(T) / 6) * 4
loadings = rng.uniform(0.3, 1.4, (n_donors + 1, 2))
levels = rng.uniform(80, 120, n_donors + 1)
Y = np.zeros((n_donors + 1, T))
for i in range(n_donors + 1):
Y[i] = (levels[i] + loadings[i, 0] * f1 + loadings[i, 1] * f2
+ rng.normal(0, 1.2, T))
treated_idx = 0
Y[treated_idx, T0:] += TRUE # the intervention, post-T0 only
names = ["TREATED"] + [f"donor_{i:02d}" for i in range(n_donors)]
panel = pd.DataFrame(Y, index=names,
columns=[f"t{t}" for t in range(T)])
# --- FIT THE WEIGHTS on the PRE-period only. ---
def synthetic_control(Y: np.ndarray, treated: int, T0: int):
"""Non-negative weights summing to 1 that best match the treated unit's
pre-period path. The simplex constraint prevents extrapolation and
usually yields a sparse, interpretable solution."""
donors = np.delete(np.arange(Y.shape[0]), treated)
X = Y[donors, :T0].T # pre-period donor paths
y = Y[treated, :T0]
def loss(w):
return np.mean((y - X @ w) ** 2)
w0 = np.full(len(donors), 1 / len(donors))
res = minimize(loss, w0, method="SLSQP",
bounds=[(0, 1)] * len(donors),
constraints=[{"type": "eq", "fun": lambda w: w.sum() - 1}],
options={"maxiter": 800, "ftol": 1e-12})
w = np.clip(res.x, 0, None)
w = w / w.sum()
synth = Y[donors].T @ w # the synthetic unit, all periods
return w, donors, synth
w, donors, synth = synthetic_control(Y, treated_idx, T0)
pre_rmse = np.sqrt(np.mean((Y[treated_idx, :T0] - synth[:T0]) ** 2))
post_gap = (Y[treated_idx, T0:] - synth[T0:]).mean()
print(f"pre-period RMSE: {pre_rmse:.3f}"
f" (sd of the treated series: {Y[treated_idx,:T0].std():.2f})")
print(f"estimated effect (mean post gap): {post_gap:+.3f} (true {TRUE:+.2f})")
# --- THE WEIGHTS ARE SPARSE AND NAMEABLE ---
wt = pd.Series(w, index=[names[d] for d in donors]).sort_values(ascending=False)
print(f"\ndonors with weight > 0.01 ({(wt > 0.01).sum()} of {len(wt)}):")
print(wt[wt > 0.01].round(3).to_string())
print("Being able to say WHICH units the comparison is made of is what makes")
print("this defensible to a stakeholder.")
# --- THE PICTURE: tracking before, diverging after. ---
print("\ntreated vs synthetic (scaled to the plot width):")
lo = min(Y[treated_idx].min(), synth.min())
hi = max(Y[treated_idx].max(), synth.max())
for label, series in [("treated ", Y[treated_idx]), ("synthetic", synth)]:
row = ""
for t in range(T):
pos = int((series[t] - lo) / (hi - lo) * 30)
row += "#" if t >= T0 else "."
print(f" {label} |{row}|")
print(f" {'':10s} {'.'*T0}{'^'*(T-T0)}")
print(f" {'':10s} pre-period (matched) post (divergence = effect)")
print("\ngap (treated - synthetic) by period:")
gap = Y[treated_idx] - synth
for block, label in [(slice(0, T0), "pre "), (slice(T0, T), "post")]:
vals = gap[block]
print(f" {label}: mean {vals.mean():+.3f} sd {vals.std():.3f}"
f" range [{vals.min():+.2f}, {vals.max():+.2f}]")
# --- INFERENCE BY PLACEBO PERMUTATION: pretend each donor was treated. ---
def placebo_distribution(Y, T0, treated_idx):
effects = []
for cand in range(Y.shape[0]):
w_c, donors_c, synth_c = synthetic_control(Y, cand, T0)
pre = np.sqrt(np.mean((Y[cand, :T0] - synth_c[:T0]) ** 2))
post = (Y[cand, T0:] - synth_c[T0:]).mean()
effects.append({"unit": cand, "pre_rmse": pre, "post_effect": post,
# A ratio makes units with poor pre-fit comparable.
"ratio": abs(post) / max(pre, 1e-9)})
return pd.DataFrame(effects)
pl = placebo_distribution(Y, T0, treated_idx)
real = pl[pl["unit"] == treated_idx].iloc[0]
# Exclude donors whose pre-period fit is poor: their gaps are noise, not
# evidence, and including them makes the test too easy to pass.
usable = pl[pl["pre_rmse"] <= 2 * real["pre_rmse"]]
rank = (usable["ratio"] >= real["ratio"]).sum()
print(f"\nplacebo permutation over {len(usable)} units with acceptable pre-fit:")
print(f" treated ratio |effect|/pre_rmse: {real['ratio']:.2f}")
others = usable[usable.unit != treated_idx]["ratio"]
print(f" placebo ratios: median {others.median():.2f},"
f" max {others.max():.2f}")
print(f" rank of the treated unit: {rank} of {len(usable)}")
print(f" permutation p-value: {rank/len(usable):.3f}")
print(f"\nWith {len(usable)} units the smallest reportable p-value is"
f" {1/len(usable):.3f}. Donor count bounds your precision.")
# --- IN-TIME PLACEBO: pretend the treatment happened earlier. ---
fake_T0 = 18
w_f, d_f, synth_f = synthetic_control(Y[:, :T0], treated_idx, fake_T0)
fake_gap = (Y[treated_idx, fake_T0:T0] - synth_f[fake_T0:T0]).mean()
print(f"\nin-time placebo at t{fake_T0} (before the real treatment):")
print(f" fake 'effect' {fake_gap:+.3f} vs real {post_gap:+.3f}")
print("A large fake effect would mean the method finds jumps that are not there.")
# --- WHEN IT FAILS: a poor pre-period fit invalidates the argument. ---
# A treated unit outside the donors' range cannot be reproduced by a convex
# combination, which is exactly what the simplex constraint prevents.
Y_bad = Y.copy()
Y_bad[treated_idx, :] += 60 # now far above every donor
w_b, d_b, synth_b = synthetic_control(Y_bad, treated_idx, T0)
bad_rmse = np.sqrt(np.mean((Y_bad[treated_idx, :T0] - synth_b[:T0]) ** 2))
print(f"\ntreated unit outside the donor range:")
print(f" pre-period RMSE {bad_rmse:.2f} (vs {pre_rmse:.2f} in the good case)")
print(f" 'effect' {(Y_bad[treated_idx, T0:] - synth_b[T0:]).mean():+.2f}"
f" -- meaningless, because the pre-period never matched")
print("Report the pre-period fit alongside every synthetic control estimate.")
print("A poor fit is a failed design, not a result with wide error bars.")
THE IDEA
one treated unit, many untreated donors
treated ......../\../\..___________ <- what happened
synthetic ......../\../\..
^^^^^^^^^^^^^^^^ |
weights fitted | treatment
HERE, pre-period |
+--------- <- what would have happened
the GAP after treatment is the effect
WEIGHTS ON THE SIMPLEX (non-negative, summing to 1)
donor_11 0.312 |###############
donor_04 0.268 |#############
donor_15 0.191 |#########
donor_02 0.140 |#######
donor_09 0.089 |####
others 0.000
-> 5 of 18 donors carry the comparison, and you can NAME them.
Why the simplex: no extrapolation (the synthetic unit stays inside the
donors' range) and a sparse, interpretable solution.
VS DIFFERENCE-IN-DIFFERENCES
DiD synthetic control
comparison one or more untreated a WEIGHTED blend of donors
units, equally weighted
parallel ASSUMED, checked by eye CONSTRUCTED by fitting
trends
works with many treated units ONE treated unit
inference standard errors placebo permutation
fails when trends were not parallel pre-period fit is poor
INFERENCE BY PERMUTATION
pretend each donor was the treated unit, compute its gap:
placebo ratios |effect|/pre_rmse: 0.3 0.4 0.6 0.7 0.9 1.1 1.4 ...
TREATED ratio: 5.8 <- extreme
p = rank / n_units. With 18 usable donors the smallest possible
p-value is 1/18 = 0.056, so donor count bounds your precision.
THE TWO PLACEBO TESTS
in-space each donor as the fake treated unit -> is our gap unusual?
in-time a fake treatment date before the real one -> do we find
effects that are not there?
REQUIREMENTS
long pre-period credibility scales with the length of the match
enough donors bounds the p-value; 15+ is typical
donors untreated by this or any related intervention
GOOD pre-period fit report the RMSE; a poor fit is a FAILED design
Key Takeaways
- Synthetic control builds a comparison unit as a weighted blend of donors chosen to reproduce the treated unit's pre-treatment path.
- It constructs parallel trends rather than assuming them, which is why it works with a single treated unit.
- Weights constrained to be non-negative and sum to one prevent extrapolation and give a sparse, nameable comparison.
- Inference comes from placebo permutation over donors, so the number of donors bounds the smallest reportable p-value.
- Always report the pre-treatment fit, because a poor fit means the design failed rather than that the estimate is imprecise.
🧪 Practice
- Implement synthetic control with weights on the simplex and recover a known effect from a simulated panel with one treated unit.
- Report the fitted weights and confirm that only a handful of donors receive meaningful weight.
- Build the placebo permutation distribution by treating each donor as the treated unit, and compute the resulting p-value.
- Run an in-time placebo at an earlier fake treatment date and confirm no effect appears.
- Interview: Your synthetic control shows a large effect and a pre-period root mean squared error nearly as large as the effect itself. What do you conclude? (Hint: what does the whole argument rest on, and what does a poor match imply about the post-period prediction?)
Uplift Modeling
Standard predictive modelling answers "who is likely to do X". Uplift modelling answers a different and usually more useful question: "whose behaviour would change if we intervened". The distinction sounds subtle and it is the difference between a campaign that works and one that wastes its budget.
Consider a retention discount. The obvious approach is to predict churn risk and send discounts to the highest-risk customers. This is almost always wrong, because among high-risk customers are many who will leave regardless of the discount. Sending them money changes nothing and costs the full amount. Meanwhile, some moderate-risk customers would stay if nudged, and they are the only ones where the spend does anything.
The framework is to divide the population into four types by what they do under each condition:
Persuadables convert only if treated. These are the entire value of the campaign.
Sure things convert either way. Treating them is pure cost, and if the treatment is a discount you are giving margin away for a sale you already had.
Lost causes convert under neither condition. Treating them is wasted spend.
Sleeping dogs, or do-not-disturb, convert only if not treated. These are real and frequently overlooked: a retention email can remind a dormant customer that they are paying for something, and a discount offer can signal that the product is not worth full price. Treating them actively causes harm.
The quantity you want is the uplift, the difference between a unit's probability of converting under treatment and under control. The fundamental problem from earlier applies with full force: you never observe both, so uplift cannot be measured per individual and must be modelled from experimental data. That last part is essential. Uplift modelling needs a randomized experiment as its training data, because it is estimating a causal quantity.
The standard approaches:
Two-model, sometimes called T-learner: fit one model on the treated and one on the control, then predict the difference. Simple, and it can be unreliable because it takes a difference of two separately estimated surfaces, and the errors do not cancel.
Single model with interactions, or S-learner: one model with treatment as a feature plus its interactions, predicting with the flag set both ways.
Class transformation: a clever reframing for the balanced case that turns uplift into a standard classification problem on a redefined target.
Direct uplift trees: tree methods whose splitting criterion maximises the difference in treatment effect between children rather than outcome purity. Causal forests are the modern version and estimate heterogeneous effects directly.
Evaluation cannot use accuracy, because the target is unobservable. The standard tools are the Qini curve and the uplift curve, which show incremental conversions as a function of how much of the population you treat, ranked by predicted uplift. A good model's curve rises steeply then flattens, and the optimal treatment fraction is where it peaks. The area between the curve and the random-targeting diagonal is the summary statistic.
import numpy as np
import pandas as pd
from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier
from sklearn.model_selection import train_test_split
rng = np.random.default_rng(41)
n = 40_000
# --- BUILD A POPULATION WITH ALL FOUR TYPES ---
X = pd.DataFrame({
"tenure": rng.exponential(200, n),
"engagement": rng.beta(2, 5, n),
"price_sensitivity": rng.beta(2, 2, n),
"support_tickets": rng.poisson(1.2, n),
})
# Baseline retention probability without any treatment.
p_control = 1 / (1 + np.exp(-(-0.4 + 0.004*X["tenure"] + 2.2*X["engagement"]
- 0.15*X["support_tickets"])))
# The discount's effect VARIES: it helps price-sensitive mid-engagement users,
# does nothing for the already-loyal, and HURTS some low-engagement users
# (a sleeping-dog effect: the offer reminds them they are paying).
true_uplift = (0.28 * X["price_sensitivity"] * (X["engagement"].between(.15, .55))
- 0.12 * (X["engagement"] < 0.12)
+ 0.01)
p_treat = np.clip(p_control + true_uplift, 0.01, 0.99)
treated = rng.random(n) < 0.5 # RANDOMIZED: required
y = np.where(treated, rng.random(n) < p_treat, rng.random(n) < p_control).astype(int)
data = X.assign(treated=treated, y=y, true_uplift=true_uplift,
p_control=p_control)
# --- THE FOUR TYPES ---
def classify(uplift, p_c, thresh=0.03):
if uplift > thresh:
return "persuadable"
if uplift < -thresh:
return "sleeping dog"
return "sure thing" if p_c > 0.6 else "lost cause"
data["type"] = [classify(u, p) for u, p in zip(true_uplift, p_control)]
summary = data.groupby("type").agg(
n=("y", "size"), share=("y", lambda s: len(s)/n),
mean_uplift=("true_uplift", "mean"),
mean_p_control=("p_control", "mean"))
print("POPULATION COMPOSITION")
print(summary.round(3).to_string())
print(f"\ntrue average uplift: {true_uplift.mean():+.4f}")
print("Only persuadables justify the spend. Sleeping dogs are actively harmed.")
# --- WHY TARGETING BY RISK FAILS ---
tr, te = train_test_split(data, test_size=0.4, random_state=0,
stratify=data["treated"])
FEATS = ["tenure", "engagement", "price_sensitivity", "support_tickets"]
# The conventional approach: predict churn risk from control-group behaviour.
risk_model = GradientBoostingClassifier(random_state=0, max_depth=3).fit(
tr.loc[~tr["treated"], FEATS], 1 - tr.loc[~tr["treated"], "y"])
te = te.copy()
te["churn_risk"] = risk_model.predict_proba(te[FEATS])[:, 1]
# --- UPLIFT MODEL 1: two models (T-learner) ---
m_t = GradientBoostingClassifier(random_state=0, max_depth=3).fit(
tr.loc[tr["treated"], FEATS], tr.loc[tr["treated"], "y"])
m_c = GradientBoostingClassifier(random_state=0, max_depth=3).fit(
tr.loc[~tr["treated"], FEATS], tr.loc[~tr["treated"], "y"])
te["uplift_two_model"] = (m_t.predict_proba(te[FEATS])[:, 1]
- m_c.predict_proba(te[FEATS])[:, 1])
# --- UPLIFT MODEL 2: single model with interactions (S-learner) ---
def add_interactions(df, flag):
out = df[FEATS].copy()
out["treated"] = flag.astype(int)
for f in FEATS: # treatment x feature interactions
out[f"{f}_x_t"] = out[f] * out["treated"]
return out
s_model = GradientBoostingClassifier(random_state=0, max_depth=3).fit(
add_interactions(tr, tr["treated"]), tr["y"])
te["uplift_s_learner"] = (
s_model.predict_proba(add_interactions(te, np.ones(len(te), bool)))[:, 1]
- s_model.predict_proba(add_interactions(te, np.zeros(len(te), bool)))[:, 1])
# --- UPLIFT MODEL 3: class transformation ---
# With 50/50 assignment, define z = 1 when (treated AND converted) or
# (control AND not converted). Then 2*P(z=1) - 1 estimates the uplift.
tr_ct = tr.copy()
tr_ct["z"] = ((tr_ct["treated"] & (tr_ct["y"] == 1))
| (~tr_ct["treated"] & (tr_ct["y"] == 0))).astype(int)
ct_model = RandomForestClassifier(n_estimators=250, min_samples_leaf=40,
random_state=0).fit(tr_ct[FEATS], tr_ct["z"])
te["uplift_class_transform"] = 2 * ct_model.predict_proba(te[FEATS])[:, 1] - 1
print("\ncorrelation of each score with the TRUE uplift:")
for col in ["churn_risk", "uplift_two_model", "uplift_s_learner",
"uplift_class_transform"]:
print(f" {col:24s} {te[col].corr(te['true_uplift']):+.3f}")
print("\nChurn risk is NEGATIVELY or weakly related to uplift: the people")
print("most likely to leave are largely those a discount cannot save.")
# --- QINI / UPLIFT CURVE: the correct evaluation. ---
def uplift_curve(df, score_col, n_bins=20):
"""Incremental conversions gained by treating the top-k% ranked by score.
Computed from the randomized holdout, which is the only valid source."""
d = df.sort_values(score_col, ascending=False).reset_index(drop=True)
out = []
for k in range(1, n_bins + 1):
cut = int(len(d) * k / n_bins)
top = d.iloc[:cut]
n_t, n_c = top["treated"].sum(), (~top["treated"]).sum()
if n_t == 0 or n_c == 0:
continue
rate_t = top.loc[top["treated"], "y"].mean()
rate_c = top.loc[~top["treated"], "y"].mean()
# Incremental conversions if this whole top-k were treated.
incremental = (rate_t - rate_c) * cut
out.append({"treated_frac": k / n_bins, "n": cut,
"uplift_rate": rate_t - rate_c,
"incremental": incremental})
return pd.DataFrame(out)
print("\nincremental conversions by fraction treated (top-ranked first):")
curves = {}
for col in ["churn_risk", "uplift_two_model", "uplift_class_transform"]:
c = uplift_curve(te, col)
curves[col] = c
peak = c.loc[c["incremental"].idxmax()]
print(f"\n {col}:")
for frac in (0.1, 0.3, 0.5, 0.8, 1.0):
row = c[np.isclose(c["treated_frac"], frac)]
if len(row):
r = row.iloc[0]
print(f" treat top {frac:>4.0%}: uplift rate {r['uplift_rate']:+.4f}"
f" incremental {r['incremental']:+8.1f}")
print(f" OPTIMAL: treat top {peak['treated_frac']:.0%}"
f" for {peak['incremental']:+.1f} incremental conversions")
# Treating everyone versus treating the top slice, in money.
full = curves["uplift_two_model"].iloc[-1]
two_m = curves["uplift_two_model"]
best = two_m.loc[two_m["incremental"].idxmax()]
COST_PER_TREATMENT, VALUE_PER_RETENTION = 12.0, 90.0
for label, row in [("treat everyone", full), ("treat the top slice", best)]:
cost = row["n"] * COST_PER_TREATMENT
value = row["incremental"] * VALUE_PER_RETENTION
print(f"\n{label}: treat {row['n']:,} users"
f" cost ${cost:,.0f} value ${value:,.0f}"
f" NET ${value-cost:+,.0f}")
print("\nTargeting by uplift beats treating everyone on net value, and the")
print("curve tells you exactly where to stop.")
# --- SLEEPING DOGS: verify the bottom of the ranking is genuinely harmed. ---
bottom = te.nsmallest(int(len(te) * 0.10), "uplift_two_model")
rate_t = bottom.loc[bottom["treated"], "y"].mean()
rate_c = bottom.loc[~bottom["treated"], "y"].mean()
print(f"\nbottom 10% by predicted uplift:")
print(f" retention: treated {rate_t:.3f} vs control {rate_c:.3f}"
f" ({rate_t-rate_c:+.4f})")
print("A negative measured uplift in the bottom decile confirms the model")
print("found real sleeping dogs. Excluding them is free value.")
THE FOUR TYPES
converts if TREATED
no yes
+-------------+-------------+
converts yes | SLEEPING | SURE |
if CONTROL | DOG | THING |
| harm: treat | waste: they |
| and lose | convert |
| them | anyway |
+-------------+-------------+
no | LOST CAUSE | PERSUADABLE |
| waste: they | THE ONLY |
| never | ONES WORTH |
| convert | TREATING |
+-------------+-------------+
uplift = P(convert | treated) - P(convert | control)
persuadable > 0 sure thing ~ 0 lost cause ~ 0
sleeping dog < 0
WHY RISK TARGETING FAILS
high churn risk = "likely to leave"
high uplift = "leaving would CHANGE if we acted"
correlation of churn risk with true uplift: -0.08
A discount cannot save someone who is leaving for reasons it does not
address, and those are exactly the highest-risk users.
THE APPROACHES
method fits weakness
---------------------------------------------------------------------------
two-model (T) one model per arm differences two noisy
surfaces
single model (S) treatment + interactions treatment flag can be
swamped by other features
class transformation one classifier on a needs ~50/50 assignment
redefined target
uplift trees / splits to maximise the more machinery
causal forest treatment-effect difference
ALL of them REQUIRE randomized training data: uplift is causal.
THE UPLIFT CURVE TELLS YOU WHERE TO STOP
incremental
conversions | ____
| _/ \____
| _/ \____
| _/ \___ <- treating more starts
| / \ DESTROYING value
|/ . . . . . . . . . . . . . (dotted = random targeting)
+-------------------------->
0% optimal ~30% 100% treated
Peak = the optimal treatment fraction. Past it you are paying to treat
sure things and waking sleeping dogs.
Key Takeaways
- Uplift modelling estimates whose behaviour would change if treated, which is a different question from who is likely to act.
- Targeting the highest-risk users wastes budget on those who leave regardless, and the two rankings can be nearly uncorrelated.
- Sleeping dogs, who respond negatively to treatment, are real, and excluding them is free value.
- Uplift is a causal quantity, so the training data must come from a randomized experiment.
- Evaluate with a Qini or uplift curve, whose peak gives the optimal fraction of the population to treat.
🧪 Practice
- Simulate a population containing persuadables, sure things, lost causes, and sleeping dogs, and report each group's share.
- Show that a churn-risk score is weakly or negatively correlated with true uplift on that population.
- Implement two-model and class-transformation uplift estimators and compare their correlation with the true uplift.
- Build an uplift curve and find the treatment fraction maximising net value given a cost per treatment and a value per retention.
- Interview: Marketing wants to send a win-back offer to every customer who has not purchased in 60 days. What do you propose instead and what do you need? (Hint: which of the four types is that segment full of, and what data would let you separate them?)
Sensitivity Analysis
Every non-randomized causal claim rests on an assumption that cannot be verified. Sensitivity analysis asks the only useful follow-up question: how wrong would that assumption have to be to overturn the conclusion? It converts an unverifiable assumption from a weakness you hope nobody mentions into a quantified statement you can defend.
The framing matters because the alternative is worse. Saying "we assume no unmeasured confounding" invites a reviewer to say "but what about X", and the conversation has nowhere to go. Saying "a confounder would need to be more strongly associated with both treatment and outcome than any variable we measured, and we measured income and prior usage" ends the discussion with evidence.
Several specific tools exist.
Rosenbaum bounds ask how much hidden bias in the assignment odds would be needed to make a significant result insignificant. The answer is a number, usually written as gamma, and it is interpretable: gamma of 1.2 means a weak confounder would break the finding, while gamma of 3 means it would take an implausibly strong one.
The E-value is the more modern and more interpretable version. It reports the minimum strength of association, on the risk-ratio scale, that an unmeasured confounder would need to have with both treatment and outcome to explain the observed effect away. Crucially, you then compare the E-value against the associations of the confounders you did measure, which grounds it in something concrete.
Placebo and negative control tests check the specification rather than the confounding. Run the analysis on an outcome the treatment cannot plausibly affect. If an effect appears there, your method finds effects that are not there. Similarly, run it on a period before the treatment existed.
Specification curve analysis addresses a different threat: the many defensible analytical choices you made. Rather than reporting one specification, run all reasonable combinations of outlier rules, control sets, functional forms, and sample restrictions, and show the distribution of estimates. If the conclusion holds across the whole curve, it is robust. If it depends on three specific choices, the reader deserves to know.
Leave-one-out and subsample stability checks whether the result depends on a single unit, period, or region.
Bounding the estimate under extreme assumptions, such as assuming all missing outcomes went the worst possible way, gives a worst-case interval that is often surprisingly informative.
The overarching principle: report the fragility, not just the estimate. A result that survives a hostile sensitivity analysis is far more persuasive than one presented with a single confident number, and a result that does not survive is something you need to know before someone else finds out.
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy import stats
import itertools
rng = np.random.default_rng(43)
# --- A CONFOUNDED OBSERVATIONAL STUDY ---
n = 8_000
# The confounder drives BOTH treatment uptake and the outcome.
confounder = rng.normal(0, 1, n)
p_treat = 1 / (1 + np.exp(-(0.8 * confounder)))
treated = rng.random(n) < p_treat
TRUE_EFFECT = 0.0 # there is NO real effect
y = 10 + 2.0 * confounder + TRUE_EFFECT * treated + rng.normal(0, 2, n)
obs = pd.DataFrame({"treated": treated, "y": y, "confounder": confounder})
naive = sm.OLS(obs["y"], sm.add_constant(obs[["treated"]].astype(float))).fit()
adjusted = sm.OLS(obs["y"], sm.add_constant(
obs[["treated", "confounder"]].astype(float))).fit()
print(f"TRUE effect: {TRUE_EFFECT:+.2f}")
print(f"naive (confounder omitted): {naive.params['treated']:+.4f}"
f" p {naive.pvalues['treated']:.2e}")
print(f"adjusted (confounder in): {adjusted.params['treated']:+.4f}"
f" p {adjusted.pvalues['treated']:.3f}")
print("\nThe naive estimate is strongly 'significant' and entirely spurious.")
print("In real work the confounder column does not exist, which is why we")
print("must ask how strong an unmeasured one would have to be.")
# --- THE E-VALUE: how strong must a confounder be? ---
def e_value(risk_ratio: float) -> float:
"""Minimum association (on the risk-ratio scale) that an unmeasured
confounder must have with BOTH treatment and outcome to explain away
an observed risk ratio. VanderWeele & Ding (2017)."""
rr = risk_ratio if risk_ratio >= 1 else 1 / risk_ratio
return rr + np.sqrt(rr * (rr - 1))
print("\nE-values for observed risk ratios:")
for rr in (1.1, 1.25, 1.5, 2.0, 3.0, 5.0):
print(f" RR {rr:>4.2f} -> E-value {e_value(rr):.2f}"
f" (a confounder needs RR >= {e_value(rr):.2f} with BOTH)")
print("\nInterpretation: an observed RR of 1.25 is explained away by a")
print("confounder only moderately associated with both. An RR of 3 needs a")
print("very strong one. Compare against the RRs of what you DID measure.")
# Ground the E-value in measured associations, which is the crucial step.
measured_rrs = {"prior usage": 2.4, "income band": 1.6, "region": 1.3,
"device type": 1.15}
observed_rr = 1.5
ev = e_value(observed_rr)
print(f"\nobserved RR {observed_rr}, E-value {ev:.2f}")
print("strongest measured confounder associations:")
for k, v in sorted(measured_rrs.items(), key=lambda kv: -kv[1]):
verdict = "COULD explain it away" if v >= ev else "too weak to explain it"
print(f" {k:14s} RR {v:.2f} -> {verdict}")
print("\nBecause prior usage (RR 2.4) exceeds the E-value, a similar")
print("unmeasured variable is entirely plausible. That is a fragile finding.")
# --- ROSENBAUM-STYLE BOUNDS: hidden bias in the assignment odds. ---
def rosenbaum_bound(estimate: float, se: float, gamma: float) -> dict:
"""How the inference degrades if the assignment odds could differ by a
factor of gamma between apparently identical units."""
inflation = np.log(gamma)
z_lower = (abs(estimate) - inflation * se * 2) / se
p_worst = 2 * (1 - stats.norm.cdf(max(z_lower, 0)))
return {"gamma": gamma, "worst_case_p": round(p_worst, 4),
"still_significant": p_worst < 0.05}
est, se = 0.42, 0.11
print(f"\nRosenbaum sensitivity for an estimate of {est} (se {se}):")
for g in (1.0, 1.1, 1.25, 1.5, 2.0, 3.0):
r = rosenbaum_bound(est, se, g)
print(f" gamma {g:>4.2f}: worst-case p {r['worst_case_p']:.4f}"
f" {'significant' if r['still_significant'] else 'NOT significant'}")
print("The gamma at which significance disappears is the headline number.")
# --- NEGATIVE CONTROL OUTCOME: does the method find fake effects? ---
# An outcome the treatment cannot possibly affect, but which shares the
# confounder. If an effect shows up here, the specification is the problem.
obs["negative_control"] = 5 + 1.5 * confounder + rng.normal(0, 2, n)
nc = sm.OLS(obs["negative_control"],
sm.add_constant(obs[["treated"]].astype(float))).fit()
print(f"\nnegative control outcome: effect {nc.params['treated']:+.4f}"
f" p {nc.pvalues['treated']:.2e}")
print("A large effect on an outcome the treatment cannot touch is proof that")
print("confounding, not causation, is producing the main result.")
# --- SPECIFICATION CURVE: report the distribution, not one number. ---
panel = pd.DataFrame({
"y": y, "treated": treated.astype(float), "confounder": confounder,
"covariate_a": rng.normal(0, 1, n), "covariate_b": rng.normal(0, 1, n),
"region": rng.integers(0, 4, n),
})
specs = []
for controls in [[], ["covariate_a"], ["confounder"],
["confounder", "covariate_a"],
["confounder", "covariate_a", "covariate_b"]]:
for winsorize in (False, True):
for drop_region in (None, 0, 1):
d = (panel if drop_region is None
else panel[panel["region"] != drop_region])
yy = d["y"]
if winsorize:
lo, hi = yy.quantile([.01, .99])
yy = yy.clip(lo, hi)
Xs = sm.add_constant(d[["treated"] + controls])
f = sm.OLS(yy, Xs).fit()
specs.append({
"n_controls": len(controls),
"controls": ",".join(controls) or "none",
"winsorized": winsorize,
"dropped_region": drop_region,
"estimate": f.params["treated"],
"p": f.pvalues["treated"],
"sig": f.pvalues["treated"] < 0.05,
})
sc = pd.DataFrame(specs)
print(f"\nSPECIFICATION CURVE over {len(sc)} defensible specifications:")
print(f" estimate: median {sc['estimate'].median():+.4f},"
f" range [{sc['estimate'].min():+.4f}, {sc['estimate'].max():+.4f}]")
print(f" significant in {sc['sig'].mean():.0%} of specifications")
print("\nby whether the true confounder is controlled:")
print(sc.groupby(sc["controls"].str.contains("confounder"))
.agg(n=("estimate", "size"), median_est=("estimate", "median"),
share_sig=("sig", "mean")).round(4).to_string())
print("\nThe result flips entirely on ONE choice. A specification curve makes")
print("that visible instead of letting the author pick a column.")
# --- LEAVE-ONE-OUT: does one group carry the whole result? ---
print("\nleave-one-region-out stability:")
full_est = sm.OLS(panel["y"], sm.add_constant(
panel[["treated", "confounder"]])).fit().params["treated"]
for r in sorted(panel["region"].unique()):
d = panel[panel["region"] != r]
e = sm.OLS(d["y"], sm.add_constant(
d[["treated", "confounder"]])).fit().params["treated"]
print(f" without region {r}: {e:+.4f}"
f" (change {e-full_est:+.4f})")
# --- WORST-CASE BOUNDS on missing outcomes ---
print("\nworst-case bounds with 8% of outcomes missing:")
missing_rate, obs_effect = 0.08, 0.42
spread = panel["y"].std()
print(f" observed effect {obs_effect:+.3f}")
print(f" if ALL missing treated outcomes were the worst and all missing")
print(f" control outcomes the best: {obs_effect - missing_rate*spread*2:+.3f}")
print(f" the reverse: {obs_effect + missing_rate*spread*2:+.3f}")
print("If the conclusion survives both extremes, missingness is not a threat.")
THE QUESTION SENSITIVITY ANALYSIS ANSWERS
NOT "is our assumption true?" (unanswerable)
BUT "how false would it need to be to change our conclusion?"
THE E-VALUE
E-value = RR + sqrt(RR * (RR - 1))
observed RR E-value what a confounder needs with BOTH
1.10 1.43 weak association suffices -> FRAGILE
1.25 1.75 moderate
1.50 2.22 moderate-to-strong
2.00 3.41 strong
3.00 5.45 very strong -> ROBUST
5.00 9.47 implausible -> very robust
THE STEP PEOPLE SKIP: GROUND IT
observed RR 1.50 -> E-value 2.22
measured: prior usage RR 2.40 <- EXCEEDS the E-value
income band RR 1.60
region RR 1.30
A variable as predictive as prior usage would explain the whole result,
and such variables plainly exist. Verdict: fragile.
NEGATIVE CONTROL LOGIC
treatment ---> outcome the claim
treatment -X-> negative control cannot be affected
confounder --> BOTH
an effect on the negative control proves confounding is active
SPECIFICATION CURVE
estimate
+0.5 |#################### significant
| ####
0.0 |........................########## not significant
-0.1 |
+-------------------------------->
specifications, sorted by estimate
controls include the confounder? median est share significant
no +0.41 100%
yes +0.01 6%
The entire conclusion rests on one modelling choice. Show the curve.
THE SENSITIVITY TOOLKIT
threat tool
---------------------------------------------------------------------------
unmeasured confounding E-value, Rosenbaum bounds
specification is wrong negative control outcome, placebo period
analyst choices specification curve
one unit drives it leave-one-out
missing outcomes worst-case bounds
Key Takeaways
- Sensitivity analysis asks how wrong an unverifiable assumption would need to be to overturn the conclusion, which is the only answerable version.
- The E-value gives the minimum confounder strength needed to explain an effect away, and it must be compared against your measured associations to mean anything.
- A negative control outcome the treatment cannot affect exposes a specification that finds effects that are not there.
- A specification curve over all defensible analytical choices reveals whether a conclusion rests on one of them.
- Report fragility alongside the estimate, because a result that survives hostile checks is far more persuasive than a single confident number.
🧪 Practice
- Simulate a confounded study with no true effect and show that the naive estimate is strongly significant.
- Implement the E-value and tabulate it for observed risk ratios from 1.1 to 5, then compare against a set of measured confounder associations.
- Build a negative control outcome that the treatment cannot affect and demonstrate that a spurious effect appears on it.
- Run a specification curve over at least twenty defensible specifications and report the share that reach significance.
- Interview: You present an observational finding and a reviewer asks how you know there is no unmeasured confounding. What is your answer? (Hint: you cannot know, so what quantity can you give them instead, and what do you compare it against?)
<a id="84-decision-analysis"></a>
8.4 Decision Analysis
An estimate is not a decision; converting one into the other requires costs, preferences, and a rule; this section covers how to frame that conversion, how to set thresholds from real consequences, and how to learn while acting.
Cost-Benefit Framing
Before any modelling, the most valuable thing you can do is write down what the decision actually costs and gains. This sounds like finance rather than data science, and skipping it is why many technically excellent analyses change nothing.
The framing has a specific structure worth following.
Enumerate the options, including doing nothing. "Do nothing" is a real option with real consequences and it is the baseline every other option must beat. Its omission is the most common error in decision framing, because analysis is usually commissioned by someone who has already decided to act.
Identify who bears each cost and benefit. A change that saves the company money by increasing customer effort has moved cost rather than eliminated it, and whether that is acceptable is a judgement the analysis should surface rather than hide.
Distinguish fixed from marginal. A one-time build cost and a per-transaction cost behave completely differently as volume changes, and conflating them produces a break-even calculation that is wrong at every volume except the one assumed.
Separate recoverable from sunk. Money already spent is irrelevant to what you should do next. This is worth naming explicitly because organisations reliably fall for it: "we have already invested six months" is not an argument for continuing.
Account for the cost of the decision process itself. Running an experiment costs traffic, engineering time, and delay. Sometimes the cost of deciding well exceeds the value of deciding correctly, and then you should just pick one. A useful test: if two options differ by less than the cost of telling them apart, and the choice is reversible, choosing quickly is the optimal decision.
Value the option, not just the outcome. A reversible decision is worth more than an irreversible one at the same expected value, because you can stop if it goes badly. This asymmetry is the quantitative justification for staged rollouts.
Two further concepts pay for themselves.
The value of information is the maximum you should pay to learn something. It is the improvement in expected outcome from making a better-informed decision, and it is zero when the information would not change your action. Run this test before commissioning any analysis: if you would ship regardless, do not spend three weeks measuring.
Asymmetric loss. Most real decisions have unequal costs of error. Missing a fraudulent transaction costs far more than reviewing a legitimate one. When errors are asymmetric, the optimal decision is not the one that minimises the error count.
import numpy as np
import pandas as pd
# --- A DECISION FRAME: options, costs, and who bears them. ---
options = pd.DataFrame([
("do nothing", 0, 0, 0.00, 0, "none"),
("manual review", 15_000, 240_000, 0.55, 0, "support team"),
("rules engine", 90_000, 40_000, 0.68, 0, "users (false blocks)"),
("ML model", 260_000, 110_000, 0.84, 30_000, "users (false blocks)"),
("ML + human", 310_000, 280_000, 0.89, 30_000, "support + users"),
], columns=["option", "build_cost", "annual_run_cost", "fraud_caught_rate",
"annual_retrain", "cost_borne_by"])
ANNUAL_FRAUD_LOSS = 4_200_000
options["fraud_prevented"] = options["fraud_caught_rate"] * ANNUAL_FRAUD_LOSS
options["annual_cost"] = options["annual_run_cost"] + options["annual_retrain"]
options["year1_net"] = (options["fraud_prevented"] - options["annual_cost"]
- options["build_cost"])
options["steady_state_net"] = options["fraud_prevented"] - options["annual_cost"]
print(options[["option", "build_cost", "annual_cost", "fraud_prevented",
"year1_net", "steady_state_net"]].to_string(index=False))
print(f"\nbest in year 1: "
f"{options.loc[options['year1_net'].idxmax(), 'option']}")
print(f"best in steady state: "
f"{options.loc[options['steady_state_net'].idxmax(), 'option']}")
print("The ranking can differ by horizon, so state the horizon first.")
# --- MARGINAL ANALYSIS: compare against the next-best, not against zero. ---
opts = options.sort_values("fraud_caught_rate").reset_index(drop=True)
opts["marginal_benefit"] = opts["fraud_prevented"].diff()
opts["marginal_cost"] = opts["annual_cost"].diff() + opts["build_cost"].diff()
opts["step_worth_it"] = opts["marginal_benefit"] > opts["marginal_cost"]
print("\nMARGINAL steps up the ladder:")
print(opts[["option", "fraud_caught_rate", "marginal_benefit",
"marginal_cost", "step_worth_it"]].to_string(index=False))
print("\nEvery option has a positive total net. Only some STEPS pay for")
print("themselves, and the last step is where the money is usually wasted.")
# --- SUNK COST: what is already spent does not affect what to do next. ---
already_spent, to_finish = 180_000, 140_000
value_if_finished, alternative = 200_000, 90_000
print(f"\nsunk cost check:")
print(f" already spent (IRRELEVANT): ${already_spent:,}")
print(f" net of FINISHING: ${value_if_finished - to_finish:,}")
print(f" net of SWITCHING: ${alternative:,}")
print(f" -> {'finish' if value_if_finished - to_finish > alternative else 'switch'}")
print(" The project loses money either way and that loss is fixed. Only the")
print(" forward-looking comparison is a decision.")
# --- VALUE OF INFORMATION: zero if it would not change the action. ---
def value_of_information(p_success, value_success, cost_failure,
study_cost, reliability=1.0):
ev_act = p_success * value_success - (1 - p_success) * cost_failure
ev_without = max(ev_act, 0.0) # 0 = the do-nothing option
# With perfect information you act only in the states where acting pays.
ev_with_perfect = p_success * value_success
voi = (ev_with_perfect - ev_without) * reliability
return {"ev_act_now": round(ev_act),
"decision_without_study": "act" if ev_act > 0 else "do nothing",
"voi": round(voi), "study_cost": study_cost,
"worth_studying": voi > study_cost}
print("\nVALUE OF INFORMATION (study reliability 0.8, cost $40,000)")
for label, kw in {
"genuinely uncertain": dict(p_success=0.50, value_success=800_000,
cost_failure=600_000),
"obviously good": dict(p_success=0.95, value_success=800_000,
cost_failure=50_000),
"obviously bad": dict(p_success=0.03, value_success=800_000,
cost_failure=600_000)}.items():
v = value_of_information(**kw, study_cost=40_000, reliability=0.8)
print(f" {label:20s} decision without study: "
f"{v['decision_without_study']:11s} VoI ${v['voi']:>8,}"
f" study? {v['worth_studying']}")
print("\nWhen the decision is already clear, measuring it is a pure cost.")
# --- ASYMMETRIC LOSS: fewest errors is not lowest cost. ---
COST_FN, COST_FP = 800.0, 15.0
policies = pd.DataFrame([
("aggressive", 120, 3_000),
("balanced", 340, 900),
("conservative", 800, 150),
], columns=["policy", "false_negatives", "false_positives"])
policies["total_errors"] = (policies["false_negatives"]
+ policies["false_positives"])
policies["total_cost"] = (policies["false_negatives"] * COST_FN
+ policies["false_positives"] * COST_FP)
print(f"\nasymmetric costs: missed fraud ${COST_FN:.0f},"
f" wasted review ${COST_FP:.0f} (ratio {COST_FN/COST_FP:.0f}:1)")
print(policies.to_string(index=False))
print(f"\nfewest errors: "
f"{policies.loc[policies['total_errors'].idxmin(), 'policy']}"
f" lowest cost: "
f"{policies.loc[policies['total_cost'].idxmin(), 'policy']}")
print("They disagree, so optimising accuracy optimises the wrong objective.")
# --- REVERSIBILITY has value independent of the expected outcome. ---
def staged_rollout_value(p_good, upside, downside,
first_stage=0.05, second_stage=0.25, detect=0.85):
"""A staged rollout caps the downside: a bad outcome is detected early,
so you only pay it on the fraction actually exposed."""
all_at_once = p_good * upside - (1 - p_good) * downside
exposure_if_bad = first_stage + (1 - detect) * (second_stage - first_stage)
staged = p_good * upside - (1 - p_good) * downside * exposure_if_bad
return {"all_at_once_ev": round(all_at_once), "staged_ev": round(staged),
"option_value": round(staged - all_at_once),
"exposure_if_bad": round(exposure_if_bad, 3)}
print("\nREVERSIBILITY AS OPTION VALUE (p(good)=0.6, +$500k / -$900k)")
for k, v in staged_rollout_value(0.6, 500_000, 900_000).items():
print(f" {k:20s} {v}")
print("\nAll-at-once has a negative expected value and staged is positive.")
print("Nothing about the change improved; the option to stop did that.")
THE DECISION FRAME
1. OPTIONS list them all, including DO NOTHING
2. COSTS fixed vs marginal; recoverable vs SUNK
3. WHO BEARS the company, the user, the support team?
4. HORIZON year 1 or steady state? the ranking can change
5. REVERSIBILITY worth more at equal expected value
6. DECIDING COST traffic, engineering time, delay
7. VoI would the answer change the action? if not, VoI = 0
MARGINAL, NOT TOTAL
option catch rate total net is the STEP worth it?
do nothing 0% $0 -
manual review 55% $2.07M yes
rules engine 68% $2.82M yes
ML model 84% $3.42M yes
ML + human 89% $3.43M NO: +$210k benefit, +$220k cost
Every option beats zero. Only some steps beat the previous option.
VALUE OF INFORMATION IS ZERO WHEN THE ACTION IS FIXED
scenario decision without study VoI study $40k?
genuinely uncertain act (barely) $240,000 YES
obviously good act $12,000 no
obviously bad do nothing $19,200 no
ASYMMETRIC LOSS
missed fraud $800, wasted review $15 -> 53:1
policy FN FP errors COST
aggressive 120 3,000 3,120 $141,000 <- lowest cost
balanced 340 900 1,240 $285,500
conservative 800 150 950 $642,250 <- fewest errors
Accuracy and cost point in OPPOSITE directions.
REVERSIBILITY AS OPTION VALUE
all at once: 0.6*500k - 0.4*900k = -$60,000 do not do it
staged 5/25/100, 85% detection: = +$253,000 do it
The downside is only paid on the exposed fraction. That is why you ramp.
Key Takeaways
- Always include "do nothing" as an option, and compare each option against the next-best rather than against zero.
- Separate fixed from marginal costs and ignore sunk costs, since only the forward-looking comparison is a decision.
- The value of information is zero when the answer would not change the action, which is the test to run before commissioning analysis.
- When error costs are asymmetric, minimising the error count and minimising total cost point in different directions.
- Reversibility carries quantifiable option value, which is the formal argument for staged rollouts.
🧪 Practice
- Build an options table with build and run costs, and show that the best option can differ between a one-year and a steady-state horizon.
- Compute marginal benefit and marginal cost for each step up an options ladder, and identify the last step worth taking.
- Implement a value-of-information function and find a case where a study is not worth its cost because the decision would not change.
- Construct three policies where the one with the fewest total errors is not the one with the lowest total cost.
- Interview: A team has spent $180,000 over six months and needs $140,000 more to finish, for an expected $200,000 of value. What do you advise? (Hint: which of those three numbers is irrelevant?)
Expected Value and Utility
Expected value is the probability-weighted average of outcomes, and it is the default rule for choosing under uncertainty. It is also not quite what anyone actually wants to maximise, and understanding the gap is what separates a mechanical calculation from good advice.
The rule is straightforward: multiply each outcome by its probability, sum, and pick the largest. Its justification is the law of large numbers. Over many repetitions the average outcome converges to the expected value, so maximising it is optimal for a decision you will face repeatedly.
That clause is the crux. For a one-off decision whose outcomes are large relative to your capacity to absorb them, expected value gives bad advice. A bet with a positive expected value that could bankrupt you is a bad bet: take it a thousand times and you are rich, take it once and you may be finished, with no thousand times available.
The reason is that utility is not linear in money. Losing your last hundred thousand hurts far more than gaining a hundred thousand helps, because the consequences are not symmetric. Ruin ends the sequence of future decisions. Utility theory handles this by applying a concave function to outcomes before taking the expectation, so losses weigh more heavily. Logarithmic utility is the standard choice, and it has the useful property of refusing any bet with a possibility of total ruin regardless of the upside.
The practical forms this takes:
Risk aversion is rational for one-shot decisions with ruinous downsides. Declining a positive-expected-value bet that could end the company is correct rather than timid.
Repeated small decisions should use expected value. Pricing, bidding, routing, and any high-frequency choice is exactly what expected value was built for, and adding risk aversion there costs money.
Variance matters when it affects your ability to keep operating. A strategy with lower expected value and much lower variance can win if it keeps you solvent long enough to compound.
The Kelly criterion formalises the repeated-bet case: it gives the fraction of capital to stake that maximises the long-run growth rate, and that fraction is far smaller than the one maximising expected value. Staking to maximise expected value each round loses everything with probability one.
Certainty equivalent is the guaranteed amount you would accept instead of a gamble, and the gap between it and the expected value is the risk premium. This turns risk preference into a number you can discuss rather than a disposition.
import numpy as np
import pandas as pd
rng = np.random.default_rng(47)
# --- EXPECTED VALUE over discrete branches ---
def expected_value(outcomes):
"""outcomes = [(probability, payoff), ...]"""
assert abs(sum(p for p, _ in outcomes) - 1) < 1e-9, "probabilities must sum to 1"
return sum(p * v for p, v in outcomes)
decisions = {
"ship now": [(0.60, 400_000), (0.40, -150_000)],
"test then ship": [(0.60, 340_000), (0.40, -20_000)],
"do nothing": [(1.00, 0)],
}
print("EXPECTED VALUE BY OPTION")
for name, outs in decisions.items():
print(f" {name:16s} EV ${expected_value(outs):>10,.0f}"
f" worst case ${min(v for _, v in outs):>10,.0f}")
print("\nTesting lowers the upside and shrinks the downside far more.")
# --- WHERE EXPECTED VALUE GIVES BAD ADVICE ---
CAPITAL = 1_000_000
bet = [(0.55, 1_200_000), (0.45, -1_000_000)]
print(f"\na bet with EV ${expected_value(bet):,.0f} on capital ${CAPITAL:,},")
print(f"where losing means losing EVERYTHING:")
def simulate_all_in(n_rounds=20, n_paths=20_000):
finals = np.empty(n_paths)
for i in range(n_paths):
c = CAPITAL
for _ in range(n_rounds):
c += 1_200_000 if rng.random() < 0.55 else -1_000_000
if c <= 0:
c = 0.0
break # ruin is absorbing: no recovery
finals[i] = c
return finals
finals = simulate_all_in()
print(f" after 20 all-in rounds: {np.mean(finals > 0):.1%} survive")
print(f" median final ${np.median(finals):,.0f}"
f" mean final ${finals.mean():,.0f}")
print(f" {np.mean(finals == 0):.0%} of paths end at zero, and the MEAN hides it.")
print("Maximising expected value round by round leads to ruin almost surely.")
# --- UTILITY: a concave function weights losses more heavily. ---
def log_utility(w, floor=1.0):
return np.log(np.maximum(w, floor))
def crra_utility(w, gamma=2.0, floor=1.0):
"""Constant relative risk aversion. gamma=0 is risk-neutral (linear),
gamma=1 is logarithmic, and larger gamma is more risk-averse."""
w = np.maximum(w, floor)
return np.log(w) if abs(gamma - 1) < 1e-9 else (w ** (1 - gamma)) / (1 - gamma)
print("\nthe same bet under linear and concave utility (capital $1,000,000):")
rows = []
for gamma, label in [(0.0, "risk neutral"), (1.0, "logarithmic"),
(2.0, "moderately averse"), (4.0, "strongly averse")]:
eu = sum(p * crra_utility(CAPITAL + v, gamma) for p, v in bet)
eu_decline = crra_utility(CAPITAL, gamma)
rows.append({"utility": label, "gamma": gamma,
"EU(take the bet)": round(eu, 6),
"EU(decline)": round(eu_decline, 6),
"take it?": eu > eu_decline})
print(pd.DataFrame(rows).to_string(index=False))
print("\nRisk-neutral takes a bet that a logarithmic decision maker refuses,")
print("and the bet's expected value never changed.")
# --- CERTAINTY EQUIVALENT AND RISK PREMIUM ---
def certainty_equivalent(outcomes, wealth, gamma=2.0):
"""The guaranteed amount that leaves you as well off as the gamble."""
eu = sum(p * crra_utility(wealth + v, gamma) for p, v in outcomes)
if abs(gamma - 1) < 1e-9:
return float(np.exp(eu) - wealth)
return float((eu * (1 - gamma)) ** (1 / (1 - gamma)) - wealth)
gamble = [(0.5, 200_000), (0.5, -100_000)]
ev = expected_value(gamble)
print(f"\ngamble: 50% +$200,000 / 50% -$100,000, EV ${ev:,.0f}")
print("certainty equivalent by wealth and risk aversion:")
for wealth in (300_000, 1_000_000, 10_000_000):
row = f" wealth ${wealth:>11,}: "
for g in (1.0, 2.0, 4.0):
ce = certainty_equivalent(gamble, wealth, g)
row += f" gamma {g}: CE ${ce:>9,.0f} (premium ${ev-ce:>8,.0f})"
print(row)
print("\nThe same gamble is worth far less to a smaller balance sheet. The")
print("risk premium is the number to put in front of a stakeholder.")
# --- KELLY: stake for long-run GROWTH, not for expected value. ---
def kelly_fraction(p_win: float, win_mult: float, loss_mult: float) -> float:
"""Fraction of capital to stake that maximises log growth.
win_mult and loss_mult are gains/losses per unit staked."""
return (p_win * win_mult - (1 - p_win) * loss_mult) / (win_mult * loss_mult)
p, b = 0.55, 1.0 # even-money bet, 55% edge
k = kelly_fraction(p, b, 1.0)
print(f"\neven-money bet won {p:.0%} of the time: Kelly stake {k:.1%} of capital")
print("(expected-value maximisation would say 100%.)")
def growth_simulation(fraction, n_rounds=250, n_paths=4_000, p_win=0.55):
"""Median terminal wealth is the right summary for a compounding process;
the mean is dominated by a few extreme paths."""
finals = np.empty(n_paths)
for i in range(n_paths):
w = 1.0
for _ in range(n_rounds):
stake = w * fraction
w += stake if rng.random() < p_win else -stake
if w <= 1e-12:
w = 0.0
break
finals[i] = w
return finals
print("\nterminal wealth after 250 rounds, starting at 1.0:")
for frac in (0.02, 0.05, k, 0.20, 0.50, 1.00):
f = growth_simulation(frac)
label = " <- Kelly" if abs(frac - k) < 1e-9 else ""
print(f" stake {frac:>5.1%} median {np.median(f):>12,.2f}"
f" ruined {np.mean(f < 0.01):>5.1%}{label}")
print("\nStaking above Kelly lowers median growth AND raises ruin risk. The")
print("all-in strategy maximises expected value and is ruined every time.")
# --- WHEN TO USE WHICH RULE ---
guide = pd.DataFrame([
("ad bidding, per impression", "millions/day", "tiny", "EXPECTED VALUE"),
("pricing a single SKU", "continuous", "small", "expected value"),
("a feature launch", "dozens/year", "moderate", "EV + staged rollout"),
("entering a new market", "once", "company-scale", "utility, risk-averse"),
("a bet that could end the firm", "once", "ruinous", "REFUSE regardless of EV"),
], columns=["decision", "repetitions", "stake relative to capacity", "rule"])
print("\n", guide.to_string(index=False))
WHEN EXPECTED VALUE IS THE RIGHT RULE
repeated + small stakes -> expected value (the law of large numbers
actually applies)
one-shot + large stakes -> utility (there is no "long run" to average
over)
any possibility of RUIN -> refuse, whatever the expected value
WHY RUIN BREAKS EXPECTED VALUE
bet: 55% win $1.2M / 45% lose $1.0M on $1.0M of capital
EV = +$210,000 per round, which sounds compelling
after 20 all-in rounds:
survive 22%
median final $0
MEAN final $2,900,000 <- driven entirely by the surviving paths
The mean is real and unreachable: 78% of paths are at zero, and zero is
absorbing.
UTILITY MAKES LOSSES WEIGH MORE
utility
| ____---------- logarithmic: the slope falls as wealth rises,
| ___--- so a loss costs more than an equal gain gains
| _--
| / linear (risk-neutral): a dollar is a dollar
| / everywhere
+---------------------> wealth
CERTAINTY EQUIVALENT: RISK PREFERENCE AS A NUMBER
gamble: 50% +$200,000 / 50% -$100,000, EV = +$50,000
wealth gamma=1 gamma=2 gamma=4
$300,000 CE +$16,000 CE -$9,000 CE -$46,000
$1,000,000 CE +$39,000 CE +$28,000 CE +$6,000
$10,000,000 CE +$49,000 CE +$48,000 CE +$45,000
The same gamble; a small balance sheet should decline it.
KELLY: STAKE FOR GROWTH, NOT FOR EXPECTED VALUE
stake median terminal wealth ruined
2% 9.4 0.0%
5% 58.1 0.0%
10% 196.3 0.0% <- Kelly for a 55% edge
20% 31.7 0.2%
50% 0.00 71%
100% 0.00 100% <- maximises EV every round
Overbetting lowers median growth and raises ruin at the same time.
Key Takeaways
- Expected value is optimal for decisions faced repeatedly at stakes you can absorb, because that is when averages apply.
- For one-shot decisions with ruinous downsides, risk aversion is rational and a concave utility formalises it.
- Ruin is absorbing, so maximising expected value each round leads to losing everything with probability one.
- The certainty equivalent turns risk preference into a number stakeholders can discuss, and the gap to expected value is the risk premium.
- The Kelly criterion stakes for long-run growth, and overbetting lowers median growth while raising ruin risk simultaneously.
🧪 Practice
- Simulate repeatedly taking a positive-expected-value bet that stakes all capital, and report the survival rate and median outcome.
- Evaluate the same bet under risk-neutral and logarithmic utility and show they disagree.
- Compute the certainty equivalent of a gamble at three wealth levels and report the risk premium for each.
- Simulate long-run growth across several staking fractions and confirm the Kelly fraction maximises median terminal wealth.
- Interview: An investment has a 60% chance of tripling the company's value and a 40% chance of bankrupting it. Expected value is strongly positive. What do you advise? (Hint: how many times can the company play this game?)
Decision Thresholds from Business Costs
A classifier outputs a score and a decision needs an action, so something must convert one into the other. That something is a threshold, and the default of 0.5 is almost never right. Choosing it from costs rather than convention is one of the highest-value and most frequently skipped steps in applied modelling.
The reason 0.5 is arbitrary is that it implicitly assumes the two errors cost the same. They rarely do. If missing a fraudulent transaction costs 500 and reviewing a legitimate one costs 10, you should accept about fifty false alarms to prevent one miss, and that means flagging anything above roughly a 2% probability rather than 50%.
That number is not a coincidence. For a well-calibrated probability, the cost-minimising threshold has a closed form: it is the cost of a false positive divided by the sum of the two error costs. Only the ratio matters, which is liberating in practice, because stakeholders who cannot tell you absolute costs can usually tell you a relative one. "A miss is about fifty times worse than a false alarm" is a question people can answer.
Two conditions qualify this.
The probabilities must be calibrated. The formula operates on a probability that means what it says. A model producing scores that rank well but are not calibrated needs recalibration first, or the threshold you compute will not be the threshold you get. Tree ensembles and support vector machines are commonly miscalibrated, and boosting tends to push probabilities toward the extremes.
Capacity is often the real constraint. If a review team can process 400 cases a day, the threshold is whatever produces 400 cases, and the cost calculation becomes a question of which 400 rather than how many. This is extremely common and usually more binding than the cost ratio.
Beyond the binary threshold, several extensions matter.
Three-way decisions with an abstain region. Automatically approve above one threshold, automatically decline below another, and route the middle band to humans. This is often better than any single threshold because it puts human effort exactly where the model is uncertain.
Per-segment thresholds. When the cost of an error varies by case, for example with transaction size, the optimal threshold varies too, and a single global threshold is leaving value on the table. The caution is fairness: a threshold that differs by a protected attribute is a decision that needs explicit justification.
Threshold drift. The optimal threshold depends on the base rate and on the cost structure, both of which change. Monitor and revisit it rather than setting it once.
import numpy as np
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.calibration import CalibratedClassifierCV, calibration_curve
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import brier_score_loss
rng = np.random.default_rng(53)
X, y = make_classification(n_samples=40_000, n_features=18, n_informative=7,
weights=[0.985, 0.015], flip_y=0.01, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.4, random_state=0,
stratify=y)
model = GradientBoostingClassifier(random_state=0, max_depth=3).fit(Xtr, ytr)
p = model.predict_proba(Xte)[:, 1]
print(f"base rate {yte.mean():.4f}, {yte.sum()} positives in {len(yte):,} cases")
# --- THE CLOSED FORM for a calibrated probability. ---
def optimal_threshold(cost_fp: float, cost_fn: float) -> float:
"""Act when p * cost_fn > (1-p) * cost_fp, i.e. p > cfp / (cfp + cfn).
Only the RATIO of the two costs matters."""
return cost_fp / (cost_fp + cost_fn)
print("\nthe cost-optimal threshold depends only on the cost RATIO:")
for cfn, cfp in [(10, 10), (50, 10), (100, 10), (500, 10), (2_000, 10)]:
t = optimal_threshold(cfp, cfn)
print(f" miss ${cfn:>5} / false alarm ${cfp:>3}"
f" (ratio {cfn/cfp:>5.0f}:1) -> threshold {t:.4f}")
print("\nStakeholders cannot state absolute costs and CAN state a ratio.")
print("'A miss is about 50x worse' is an answerable question.")
# --- VERIFY EMPIRICALLY, by sweeping the threshold. ---
COST_FN, COST_FP = 500.0, 10.0
grid = np.unique(np.quantile(p, np.linspace(0.5, 0.9999, 600)))
sweep = []
for t in grid:
flagged = p >= t
fn = int(((yte == 1) & ~flagged).sum())
fp = int(((yte == 0) & flagged).sum())
tp = int(((yte == 1) & flagged).sum())
sweep.append({"threshold": t, "flagged": int(flagged.sum()),
"tp": tp, "fp": fp, "fn": fn,
"precision": tp / max(flagged.sum(), 1),
"recall": tp / max(yte.sum(), 1),
"total_cost": fn * COST_FN + fp * COST_FP})
sw = pd.DataFrame(sweep)
best = sw.loc[sw["total_cost"].idxmin()]
theory = optimal_threshold(COST_FP, COST_FN)
print(f"\ntheoretical optimum {theory:.4f},"
f" empirical optimum {best['threshold']:.4f}")
print(f"cost at the empirical optimum: ${best['total_cost']:,.0f}")
at_half = sw.iloc[(sw["threshold"] - 0.5).abs().argmin()]
at_base = sw.iloc[(sw["threshold"] - yte.mean()).abs().argmin()]
print("\ncost comparison across threshold choices:")
for label, row in [("default 0.5", at_half), ("the base rate", at_base),
("cost-optimal", best)]:
print(f" {label:16s} t={row['threshold']:.4f}"
f" flagged {row['flagged']:>5}"
f" recall {row['recall']:.3f} precision {row['precision']:.3f}"
f" cost ${row['total_cost']:>10,.0f}")
saving = at_half["total_cost"] - best["total_cost"]
print(f"\nmoving off the default saves ${saving:,.0f}"
f" ({saving/at_half['total_cost']:.0%}) with no model change at all.")
# --- CALIBRATION IS A PRECONDITION for the closed form. ---
print("\nis the model calibrated? (predicted vs observed by decile of score)")
frac_pos, mean_pred = calibration_curve(yte, p, n_bins=8, strategy="quantile")
for mp, fp_ in zip(mean_pred, frac_pos):
flag = "" if abs(mp - fp_) < 0.3 * max(fp_, 1e-6) else " <- off"
print(f" predicted {mp:.4f} observed {fp_:.4f}{flag}")
calibrated = CalibratedClassifierCV(
GradientBoostingClassifier(random_state=0, max_depth=3),
method="isotonic", cv=3).fit(Xtr, ytr)
p_cal = calibrated.predict_proba(Xte)[:, 1]
print(f"\nBrier score: raw {brier_score_loss(yte, p):.6f}"
f" calibrated {brier_score_loss(yte, p_cal):.6f}")
cost_raw = (((yte == 1) & (p < theory)).sum() * COST_FN
+ ((yte == 0) & (p >= theory)).sum() * COST_FP)
cost_cal = (((yte == 1) & (p_cal < theory)).sum() * COST_FN
+ ((yte == 0) & (p_cal >= theory)).sum() * COST_FP)
print(f"cost at the theoretical threshold: raw ${cost_raw:,.0f}"
f" calibrated ${cost_cal:,.0f}")
print("The closed form only lands on the true optimum for calibrated scores.")
# --- CAPACITY: usually the binding constraint. ---
DAILY_CAPACITY, DAILY_CASES = 400, 20_000
target_rate = DAILY_CAPACITY / DAILY_CASES
cap_threshold = float(np.quantile(p_cal, 1 - target_rate))
flagged = p_cal >= cap_threshold
print(f"\ncapacity constraint: {DAILY_CAPACITY} reviews of {DAILY_CASES:,} cases")
print(f" required flag rate {target_rate:.3%} -> threshold {cap_threshold:.4f}")
print(f" recall achieved {((yte==1) & flagged).sum()/yte.sum():.3f}"
f" precision {((yte==1) & flagged).sum()/max(flagged.sum(),1):.3f}")
print(f" cost-optimal threshold would flag "
f"{(p_cal >= theory).mean():.2%} of cases"
f" = {(p_cal >= theory).mean()*DAILY_CASES:,.0f}/day"
f" ({(p_cal >= theory).mean()*DAILY_CASES/DAILY_CAPACITY:.0f}x capacity)")
print("When capacity binds, the question is WHICH cases, not how many.")
# --- THREE-WAY DECISIONS: auto-approve, review, auto-decline. ---
def three_way(p, y, t_low, t_high, cost_fn, cost_fp, review_cost, review_acc=0.95):
auto_approve = p < t_low
auto_decline = p >= t_high
review = ~auto_approve & ~auto_decline
# Missed positives among auto-approved, plus review misses.
fn = ((y == 1) & auto_approve).sum() + ((y == 1) & review).sum() * (1-review_acc)
fp = ((y == 0) & auto_decline).sum() + ((y == 0) & review).sum() * (1-review_acc)
return {"auto_approve_pct": auto_approve.mean(),
"review_pct": review.mean(), "auto_decline_pct": auto_decline.mean(),
"reviews": int(review.sum()),
"total_cost": fn*cost_fn + fp*cost_fp + review.sum()*review_cost}
print("\nthree-way policies (review costs $8 and is 95% accurate):")
single = three_way(p_cal, yte, theory, theory, COST_FN, COST_FP, 8.0)
print(f" single threshold review {single['review_pct']:.2%}"
f" cost ${single['total_cost']:>10,.0f}")
for lo, hi in [(0.005, 0.08), (0.01, 0.15), (0.005, 0.30)]:
r = three_way(p_cal, yte, lo, hi, COST_FN, COST_FP, 8.0)
print(f" approve<{lo:.3f} decline>{hi:.2f}"
f" review {r['review_pct']:.2%} ({r['reviews']:,})"
f" cost ${r['total_cost']:>10,.0f}")
print("An abstain band puts human effort where the model is uncertain, which")
print("beats any single cut when review is cheaper than the errors it prevents.")
# --- PER-CASE COSTS: the threshold should vary with the amount at risk. ---
amounts = rng.lognormal(4.2, 1.3, len(yte))
print("\nwhen the cost of a miss scales with the transaction amount:")
global_t = theory
total_global = (((yte == 1) & (p_cal < global_t)) * amounts).sum() \
+ ((yte == 0) & (p_cal >= global_t)).sum() * COST_FP
# Per-case threshold: cfp / (cfp + amount_at_risk).
per_case_t = COST_FP / (COST_FP + amounts)
total_percase = (((yte == 1) & (p_cal < per_case_t)) * amounts).sum() \
+ ((yte == 0) & (p_cal >= per_case_t)).sum() * COST_FP
print(f" one global threshold: ${total_global:,.0f}")
print(f" per-case thresholds: ${total_percase:,.0f}"
f" ({1-total_percase/total_global:.0%} lower)")
print(f" threshold range: {per_case_t.min():.5f} to {per_case_t.max():.3f}")
print("A $12 transaction and a $9,000 one do not deserve the same scrutiny.")
print("Caution: a threshold that varies by a protected attribute needs an")
print("explicit fairness justification, which chapter 16 covers.")
THE CLOSED FORM
act when p * cost_FN > (1 - p) * cost_FP
cost_FP
threshold = -----------------
cost_FP + cost_FN
Only the RATIO matters:
cost ratio (miss : false alarm) optimal threshold
1 : 1 0.5000 <- the only case 0.5 is right
5 : 1 0.1667
10 : 1 0.0909
50 : 1 0.0196
200 : 1 0.0050
THE COST OF USING 0.5
threshold flagged recall precision total cost
default 0.5 142 0.243 0.394 $1,097,000
the base rate 903 0.612 0.198 $585,000
cost-optimal 2,418 0.847 0.104 $412,000
A 62% cost reduction, no model change, no extra data. This is usually
the cheapest available improvement to a deployed classifier.
PRECONDITIONS AND CONSTRAINTS
CALIBRATION the formula acts on a probability that means what it says;
recalibrate boosted trees and SVMs first
CAPACITY if reviewers can handle 400/day and the optimum flags
2,400, the threshold IS the capacity quantile, and the
question becomes WHICH 400
DRIFT the optimum moves with the base rate and the cost
structure; monitor it
THE THREE-WAY DECISION
p: 0 ---------- t_low ------------ t_high ---------- 1
| AUTO | SEND TO | AUTO |
| APPROVE | HUMAN | DECLINE |
| (cheap) | (costly, but | (cheap) |
| | accurate) | |
Human effort goes exactly where the model is uncertain. Beats any
single cut whenever review costs less than the errors it prevents.
Key Takeaways
- The default threshold of 0.5 assumes equal error costs, and the cost- optimal threshold is the false positive cost divided by the sum of both error costs.
- Only the cost ratio matters, which is fortunate because stakeholders can state a ratio even when they cannot state absolute costs.
- The closed form requires calibrated probabilities, so recalibrate before applying it.
- Review capacity is often more binding than the cost ratio, turning the question from how many to flag into which ones.
- An abstain band routing uncertain cases to humans beats any single threshold when review costs less than the errors it prevents.
🧪 Practice
- Implement the closed-form optimal threshold and tabulate it for cost ratios from 1:1 to 200:1.
- Sweep the threshold on a real classifier's scores and confirm the empirical cost minimum matches the closed form.
- Compare total cost at 0.5, at the base rate, and at the cost-optimal threshold, and express the saving as a percentage.
- Implement a three-way policy with an abstain band and find the band that minimises total cost including review cost.
- Interview: Your fraud model is deployed at a 0.5 threshold and the review team says they see almost no fraud in their queue. What is happening and what do you change? (Hint: at 0.5, how many of the actual positives are being flagged at all?)
Multi-Armed Bandits
An A/B test separates learning from acting. You run the test, accept that half your traffic gets the worse variant for the whole duration, then act on the result. Bandit algorithms interleave the two: they shift traffic toward better-performing options while still exploring, so the cost of learning is paid down as it is incurred.
The name comes from a gambler facing several slot machines with unknown payouts. Each pull gives information and a payoff, and the tension is the explore-exploit trade-off: play the best-so-far arm to earn, or play an uncertain arm to learn. Pure exploitation locks onto an arm that got lucky early; pure exploration wastes traffic on arms already known to be bad.
The quantity being minimised is regret, the cumulative difference between what you earned and what you would have earned always playing the best arm. Good algorithms achieve regret growing logarithmically in the number of trials rather than linearly, which means the per-decision cost of learning tends to zero.
The main algorithms, in increasing sophistication:
Epsilon-greedy plays the best-known arm with probability one minus epsilon and a random arm otherwise. Trivial to implement and never stops wasting a fixed fraction on known-bad arms unless you decay epsilon over time.
Upper confidence bound plays the arm with the highest optimistic estimate, adding a bonus that grows with uncertainty. The principle is "optimism in the face of uncertainty", and it explores automatically: an arm tried rarely has a large bonus, so it gets tried. It has strong theoretical guarantees and no randomness.
Thompson sampling maintains a posterior over each arm's value, draws one sample from each, and plays the arm with the highest draw. Exploration emerges naturally from posterior uncertainty, it performs excellently in practice, it handles delayed feedback gracefully, and it is about ten lines of code. This is usually the right default.
Contextual bandits condition the choice on features of the current request, which turns the problem into personalisation: the best arm for one user differs from another's.
The essential thing is knowing when not to use a bandit. They are not a strict improvement over A/B tests.
Use a bandit when you have many alternatives, feedback is fast, the metric is short-horizon, and you care mainly about cumulative performance. Content selection, headline choice, ad creative, and layout ranking all fit.
Use an A/B test when you need an unbiased effect estimate with a confidence interval, when feedback is delayed relative to the decision, when effects vary over time, or when the result must be defensible to a regulator or a sceptical reviewer. A bandit's adaptive allocation makes conventional inference invalid, because the sample size in each arm depends on the outcomes observed, so a naive t-test on bandit data is biased.
There is also a subtler trap: a bandit optimising a short-term metric will happily converge to something that damages a long-term one, and because it stops exploring you may never learn that.
import numpy as np
import pandas as pd
# --- THE SETUP: five arms, unknown true rates. ---
TRUE_RATES = np.array([0.050, 0.055, 0.062, 0.048, 0.058])
BEST = TRUE_RATES.max()
N_ARMS, HORIZON = len(TRUE_RATES), 40_000
print(f"true rates: {TRUE_RATES}")
print(f"best arm: {TRUE_RATES.argmax()} at {BEST:.3f}")
def pull(arm: int, r) -> int:
return int(r.random() < TRUE_RATES[arm])
# --- STRATEGY 1: A/B TEST (equal allocation throughout, then commit) ---
def ab_test(horizon, r, test_frac=0.5):
n_test = int(horizon * test_frac)
successes, trials = np.zeros(N_ARMS), np.zeros(N_ARMS)
reward = 0
for t in range(n_test): # equal split while testing
a = t % N_ARMS
y = pull(a, r); reward += y
successes[a] += y; trials[a] += 1
winner = int((successes / np.maximum(trials, 1)).argmax())
for _ in range(horizon - n_test): # then exploit only
reward += pull(winner, r)
return reward, winner, trials
# --- STRATEGY 2: EPSILON-GREEDY ---
def epsilon_greedy(horizon, r, eps=0.1, decay=False):
successes, trials = np.zeros(N_ARMS), np.zeros(N_ARMS)
reward = 0
for t in range(1, horizon + 1):
e = eps / np.sqrt(t / 100) if decay else eps
if r.random() < e or trials.min() == 0:
a = r.integers(N_ARMS) # explore
else:
a = int((successes / np.maximum(trials, 1)).argmax()) # exploit
y = pull(a, r); reward += y
successes[a] += y; trials[a] += 1
return reward, int(trials.argmax()), trials
# --- STRATEGY 3: UPPER CONFIDENCE BOUND ---
def ucb1(horizon, r, c=2.0):
successes, trials = np.zeros(N_ARMS), np.zeros(N_ARMS)
reward = 0
for t in range(1, horizon + 1):
if trials.min() == 0:
a = int(trials.argmin()) # each arm once first
else:
mean = successes / trials
# The bonus grows with total trials and shrinks with this arm's
# trials: optimism in the face of uncertainty.
bonus = np.sqrt(c * np.log(t) / trials)
a = int((mean + bonus).argmax())
y = pull(a, r); reward += y
successes[a] += y; trials[a] += 1
return reward, int((successes/np.maximum(trials,1)).argmax()), trials
# --- STRATEGY 4: THOMPSON SAMPLING (the usual default) ---
def thompson(horizon, r, prior=(1.0, 1.0)):
alpha = np.full(N_ARMS, prior[0])
beta = np.full(N_ARMS, prior[1])
reward = 0
for _ in range(horizon):
# Draw one sample per arm from its posterior and play the best draw.
draws = r.beta(alpha, beta)
a = int(draws.argmax())
y = pull(a, r); reward += y
alpha[a] += y
beta[a] += 1 - y
trials = alpha + beta - prior[0] - prior[1]
return reward, int((alpha/(alpha+beta)).argmax()), trials
print(f"\nover {HORIZON:,} pulls, averaged across 20 runs:")
oracle = BEST * HORIZON
results = []
for name, fn, kw in [("A/B test (50% test)", ab_test, {}),
("epsilon-greedy 0.1", epsilon_greedy, {}),
("epsilon-greedy, decayed", epsilon_greedy, {"decay": True}),
("UCB1", ucb1, {}),
("Thompson sampling", thompson, {})]:
rewards, correct, alloc = [], 0, np.zeros(N_ARMS)
for seed in range(20):
r = np.random.default_rng(seed)
rew, win, trials = fn(HORIZON, r, **kw)
rewards.append(rew)
correct += int(win == TRUE_RATES.argmax())
alloc += trials
results.append({
"strategy": name,
"mean_reward": int(np.mean(rewards)),
"regret": int(oracle - np.mean(rewards)),
"picked_best": f"{correct}/20",
"best_arm_share": round(alloc[TRUE_RATES.argmax()] / alloc.sum(), 3),
"worst_arm_share": round(alloc[TRUE_RATES.argmin()] / alloc.sum(), 3),
})
print(pd.DataFrame(results).to_string(index=False))
print(f"\noracle (always the best arm): {oracle:.0f} conversions")
print("Regret is the cost of not knowing. Thompson and UCB pay far less of it")
print("than a fixed 50% test, because they stop sending traffic to bad arms.")
# --- REGRET GROWTH: logarithmic, not linear. ---
print("\ncumulative regret as the horizon grows (Thompson, 10 runs):")
prev = None
for h in (2_000, 8_000, 32_000):
regs = []
for seed in range(10):
r = np.random.default_rng(100 + seed)
rew, _, _ = thompson(h, r)
regs.append(BEST * h - rew)
m = np.mean(regs)
growth = f" ({m/prev:.2f}x for 4x the horizon)" if prev else ""
print(f" horizon {h:>6,}: regret {m:>7.1f}"
f" per-pull {m/h:.5f}{growth}")
prev = m
print("Regret grows far slower than the horizon, so the per-decision cost of")
print("learning falls toward zero. That is the whole appeal.")
# --- WHY NAIVE INFERENCE ON BANDIT DATA IS BIASED ---
print("\nestimating an arm's rate from bandit data vs from a fixed split:")
for label, fn in [("Thompson", thompson), ("A/B (equal split)", ab_test)]:
ests = []
for seed in range(40):
r = np.random.default_rng(500 + seed)
if label == "Thompson":
alpha = np.ones(N_ARMS); beta = np.ones(N_ARMS)
for _ in range(6_000):
a = int(r.beta(alpha, beta).argmax())
y = pull(a, r); alpha[a] += y; beta[a] += 1 - y
worst = TRUE_RATES.argmin()
n = alpha[worst] + beta[worst] - 2
ests.append((alpha[worst] - 1) / max(n, 1))
else:
s = t = 0
worst = TRUE_RATES.argmin()
for _ in range(6_000 // N_ARMS):
y = pull(worst, r); s += y; t += 1
ests.append(s / t)
print(f" {label:20s} worst-arm estimate {np.mean(ests):.4f}"
f" (true {TRUE_RATES.min():.4f})"
f" sd {np.std(ests):.4f}")
print("The bandit's estimate of a bad arm rests on very few pulls, so it is")
print("noisy and its sample size is OUTCOME-DEPENDENT. A naive t-test on")
print("bandit data does not have the coverage it claims.")
# --- WHEN TO USE WHICH ---
guide = pd.DataFrame([
("many variants (5+)", "bandit", "regret scales with arms in a fixed test"),
("need a CI on the effect", "A/B TEST", "adaptive allocation breaks inference"),
("fast feedback (minutes)", "bandit", "needs quick updates to adapt"),
("delayed feedback (weeks)", "A/B test", "cannot adapt before committing"),
("short-lived content", "bandit", "no time for a full test"),
("a permanent product change", "A/B test", "you need a defensible estimate"),
("regulated / auditable", "A/B TEST", "adaptive designs are hard to defend"),
("maximise cumulative reward", "bandit", "that IS the bandit objective"),
("long-horizon metric", "A/B test", "a bandit optimises the short one"),
("effects change over time", "A/B test, or a non-stationary bandit",
"a standard bandit assumes stationarity"),
], columns=["situation", "choose", "why"])
print("\n", guide.to_string(index=False))
THE EXPLORE-EXPLOIT TRADE-OFF
exploit only lock onto whichever arm got lucky first
explore only keep paying to learn what you already know
good algorithm explore in proportion to UNCERTAINTY
ALLOCATION OVER TIME
A/B test (50% test, then commit) Thompson sampling
arm arm
0 |#########| 0 |####|
1 |#########| 1 |#####|
2 |#########|################# 2 |###|##################
3 |#########| 3 |##|
4 |#########| 4 |######|
+---------+-----------------> +-------------------->
equal while testing, then all traffic shifts CONTINUOUSLY
to the winner toward the better arms
RESULTS OVER 40,000 PULLS (true best rate 6.2%)
strategy reward regret best-arm share
A/B test, 50% test 2,401 79 0.60
epsilon-greedy 0.1 2,389 91 0.68
epsilon-greedy, decayed 2,436 44 0.85
UCB1 2,443 37 0.87
Thompson sampling 2,452 28 0.91
REGRET GROWS LOGARITHMICALLY
horizon regret per-pull
2,000 12.4 0.00620
8,000 21.8 0.00272 (1.76x for 4x the horizon)
32,000 33.1 0.00103 (1.52x for 4x the horizon)
The per-decision cost of learning falls toward zero.
WHY NAIVE INFERENCE FAILS ON BANDIT DATA
a bandit deliberately UNDER-SAMPLES bad arms, so:
- each arm's sample size depends on the outcomes observed
- the estimate for a rarely-pulled arm is noisy and biased
A t-test on bandit data does not have its stated coverage. If you need
a confidence interval, run an experiment.
CHOOSE A BANDIT WHEN CHOOSE AN A/B TEST WHEN
many arms you need an unbiased CI
fast feedback feedback is delayed
short-lived content it is a permanent change
cumulative reward is the goal it must be auditable
the metric is long-horizon
Key Takeaways
- Bandits interleave learning and acting, shifting traffic toward better arms so the cost of exploration is paid down as it is incurred.
- Regret grows logarithmically rather than linearly, so the per-decision cost of learning tends toward zero.
- Thompson sampling is the usual default: it is short, handles delayed feedback, and performs excellently.
- Adaptive allocation makes conventional inference invalid, because each arm's sample size depends on observed outcomes.
- Choose an A/B test when you need a defensible effect estimate, when feedback is delayed, or when the metric is long-horizon.
🧪 Practice
- Implement epsilon-greedy, upper confidence bound, and Thompson sampling on five arms, and compare cumulative regret over 40,000 pulls.
- Show that regret grows sublinearly by measuring it at three horizons each four times the last.
- Compare the estimated rate of the worst arm from bandit data against a fixed equal split, and explain the difference in precision.
- Add a decay schedule to epsilon-greedy and quantify the improvement over a fixed epsilon.
- Interview: A product manager wants to replace all A/B tests with bandits to stop wasting traffic. Where do you push back? (Hint: what does the company lose when it can no longer state a confidence interval, and what happens when the outcome takes three weeks to observe?)
Simulation and Monte Carlo Methods
Many decisions involve enough interacting uncertainties that no formula gives the answer. Simulation replaces the formula: you specify the uncertainty in the inputs, sample from it many times, propagate each sample through the logic, and read the distribution of outcomes. It is the most generally useful quantitative tool in decision analysis and the least demanding mathematically.
The reason it matters so much is the flaw of averages. Plugging the average value of each input into a model does not give the average output whenever the logic is non-linear, and almost all real logic is non-linear. Capacity limits, minimums and maximums, thresholds, options, and any if-statement all break it. A project where every task takes its average duration finishes far earlier than the average project, because delays on parallel paths compound while early finishes do not help. Building a plan around average inputs systematically understates cost and duration and overstates profit.
The method is mechanical:
- Specify a distribution for each uncertain input, not a point estimate.
- Encode the dependencies, because inputs are rarely independent, and correlation between them changes the tails dramatically.
- Sample and propagate many times through the actual decision logic.
- Read the distribution, including the tails, not just the mean.
Choosing distributions is where judgement enters. Use a triangular or PERT distribution when you have a plausible minimum, mode, and maximum from an expert, which is the common case. Use a lognormal for quantities that cannot be negative and have a long right tail, such as durations and costs. Use empirical resampling from historical data when you have it, which avoids assuming a shape at all.
The outputs worth reading are not the mean. The percentiles tell you the downside, and a plan should usually be made against the 80th or 90th percentile of cost and duration rather than the median. The probability of a specific bad outcome is directly decision-relevant. And sensitivity analysis on the simulation, correlating each input's sampled value with the output, tells you which uncertainty is actually driving your risk and therefore where to spend effort reducing it.
Two cautions. Simulation inherits your assumptions: a precise-looking distribution built on invented inputs is false confidence dressed up as analysis, and a sensible practice is to run the simulation under several assumption sets. And check convergence by confirming your reported statistic is stable as you add samples, particularly for tail quantiles, which need far more samples than the mean.
import numpy as np
import pandas as pd
rng = np.random.default_rng(59)
N_SIMS = 100_000
# --- THE FLAW OF AVERAGES, demonstrated on a capacity-limited system. ---
# Revenue = price * min(demand, capacity). The min() is the non-linearity.
PRICE, CAPACITY = 40.0, 1_000
mean_demand = 1_000
point_estimate = PRICE * min(mean_demand, CAPACITY)
demand = rng.lognormal(np.log(mean_demand) - 0.5*0.35**2, 0.35, N_SIMS)
revenue = PRICE * np.minimum(demand, CAPACITY)
print("THE FLAW OF AVERAGES")
print(f" plug in average demand ({mean_demand}): revenue ${point_estimate:,.0f}")
print(f" simulated MEAN revenue: ${revenue.mean():,.0f}")
print(f" gap: ${point_estimate - revenue.mean():,.0f}"
f" ({1 - revenue.mean()/point_estimate:.1%} overstatement)")
print(f" mean simulated demand was {demand.mean():,.0f}, so the inputs were")
print(f" right and the OUTPUT is still wrong: capacity truncates the upside")
print(f" while the downside is unbounded.")
# --- A PROJECT PLAN: parallel paths make averages worse still. ---
def triangular(low, mode, high, size, r):
return r.triangular(low, mode, high, size)
tasks = pd.DataFrame([
("design", 8, 12, 25, "start"),
("backend", 15, 22, 45, "design"),
("frontend", 12, 18, 40, "design"),
("data model", 6, 10, 28, "design"),
("integration", 5, 8, 20, "backend+frontend+data model"),
("QA", 4, 7, 18, "integration"),
], columns=["task", "low", "mode", "high", "depends_on"])
print("\nPROJECT DURATION")
print(tasks.to_string(index=False))
deterministic = (tasks.loc[tasks.task == "design", "mode"].iloc[0]
+ max(tasks.loc[tasks.task.isin(["backend", "frontend",
"data model"]), "mode"])
+ tasks.loc[tasks.task == "integration", "mode"].iloc[0]
+ tasks.loc[tasks.task == "QA", "mode"].iloc[0])
sims = {t.task: triangular(t.low, t.mode, t.high, N_SIMS, rng)
for t in tasks.itertuples()}
# The critical path is the MAX of the parallel branch, sampled per run.
parallel_max = np.maximum.reduce([sims["backend"], sims["frontend"],
sims["data model"]])
total = sims["design"] + parallel_max + sims["integration"] + sims["QA"]
print(f"\n all-modes estimate: {deterministic:.0f} days")
print(f" simulated mean: {total.mean():.1f} days")
print(f" simulated median: {np.median(total):.1f} days")
for q in (50, 80, 90, 95):
print(f" p{q}: {np.percentile(total, q):.1f} days")
print(f"\n P(finish by the all-modes estimate of {deterministic:.0f} days):"
f" {np.mean(total <= deterministic):.1%}")
print(" Committing to the sum of most-likely durations gives a plan with")
print(" roughly a 1-in-10 chance of being met. Plan to p80.")
# Which parallel task is the bottleneck? Not always the longest on average.
for name in ("backend", "frontend", "data model"):
share = np.mean(sims[name] == parallel_max)
print(f" {name:12s} is the binding path {share:.1%} of the time")
# --- CORRELATION CHANGES THE TAILS ---
def correlated_normals(rho, size, r):
"""Two standard normals with correlation rho, via a Cholesky step."""
z1 = r.normal(size=size)
z2 = rho * z1 + np.sqrt(1 - rho**2) * r.normal(size=size)
return z1, z2
print("\nCORRELATION BETWEEN RISKS")
for rho in (0.0, 0.5, 0.9):
z1, z2 = correlated_normals(rho, N_SIMS, rng)
cost_a = np.exp(np.log(500_000) + 0.30 * z1)
cost_b = np.exp(np.log(300_000) + 0.35 * z2)
tot = cost_a + cost_b
print(f" rho {rho:.1f}: mean ${tot.mean():>12,.0f}"
f" p95 ${np.percentile(tot, 95):>12,.0f}"
f" P(>$1.2M) {np.mean(tot > 1_200_000):.2%}")
print(" The mean barely moves and the TAIL moves a lot. Assuming")
print(" independence between risks that share a cause understates disaster.")
# --- DECISION UNDER SIMULATION: compare options on their distributions. ---
def simulate_option(fixed_cost, unit_margin_dist, demand_dist, capacity, n, r):
margin = unit_margin_dist(n, r)
dem = demand_dist(n, r)
units = np.minimum(dem, capacity)
return units * margin - fixed_cost
opts = {
"small plant": dict(fixed_cost=1_200_000, capacity=30_000),
"large plant": dict(fixed_cost=2_600_000, capacity=70_000),
"outsource": dict(fixed_cost=200_000, capacity=90_000),
}
margin_dists = {"small plant": lambda n, r: r.normal(85, 12, n),
"large plant": lambda n, r: r.normal(92, 12, n),
"outsource": lambda n, r: r.normal(52, 9, n)}
demand_dist = lambda n, r: r.lognormal(np.log(40_000) - .5*.45**2, .45, n)
print("\nCAPACITY DECISION")
rows = []
profits = {}
for name, kw in opts.items():
pr = simulate_option(kw["fixed_cost"], margin_dists[name], demand_dist,
kw["capacity"], N_SIMS, rng)
profits[name] = pr
rows.append({"option": name, "mean": int(pr.mean()),
"median": int(np.median(pr)),
"p5": int(np.percentile(pr, 5)),
"p95": int(np.percentile(pr, 95)),
"P(loss)": round(float((pr < 0).mean()), 4),
"P(>$2M)": round(float((pr > 2_000_000).mean()), 3)})
print(pd.DataFrame(rows).to_string(index=False))
best_mean = max(profits, key=lambda k: profits[k].mean())
safest = min(profits, key=lambda k: (profits[k] < 0).mean())
print(f"\nhighest mean profit: {best_mean}")
print(f"lowest chance of a loss: {safest}")
print("If these differ, the choice depends on risk tolerance, and the")
print("distributions make that trade-off explicit rather than implicit.")
# --- WHICH UNCERTAINTY DRIVES THE RISK? ---
print("\nsensitivity: correlation of each sampled input with the outcome")
d = demand_dist(N_SIMS, rng)
m = margin_dists["large plant"](N_SIMS, rng)
u = np.minimum(d, opts["large plant"]["capacity"])
pr = u * m - opts["large plant"]["fixed_cost"]
for label, series in [("demand", d), ("unit margin", m)]:
print(f" {label:12s} correlation with profit"
f" {np.corrcoef(series, pr)[0,1]:+.3f}")
print(" Spend your effort reducing the uncertainty that actually moves the")
print(" outcome. Here that means a better demand forecast, not a better")
print(" margin estimate.")
# --- CONVERGENCE: tail quantiles need far more samples than the mean. ---
print("\nconvergence check (large plant profit):")
big = simulate_option(opts["large plant"]["fixed_cost"],
margin_dists["large plant"], demand_dist,
opts["large plant"]["capacity"], 400_000, rng)
for n in (1_000, 10_000, 100_000, 400_000):
sub = big[:n]
print(f" n={n:>7,} mean ${sub.mean():>10,.0f}"
f" p5 ${np.percentile(sub, 5):>10,.0f}"
f" p99 ${np.percentile(sub, 99):>10,.0f}")
print(" The mean stabilises quickly; the tail percentiles keep moving.")
print(" Always report the sample count alongside a tail statistic.")
THE FLAW OF AVERAGES
revenue = price * min(demand, capacity)
average demand in -> $40,000 (the plan)
simulated mean -> $35,100 (reality)
demand distribution revenue
#### | capacity cap
######## | |
############ |____________|________ upside TRUNCATED
############## | /
----|-------------- | /
capacity | / downside UNBOUNDED
+---------------------> demand
Any min, max, threshold, or if-statement breaks the average.
PROJECT DURATION: PARALLEL PATHS COMPOUND DELAY
sum of most-likely durations: 47 days
simulated median: 56 days
simulated p80: 68 days
P(finishing by 47 days): 11%
Three parallel tasks: the branch takes the MAX, so one slipping delays
everything while the others finishing early helps nobody. Plan to p80.
CORRELATION MOVES THE TAIL, NOT THE MEAN
rho mean cost p95 cost P(> $1.2M)
0.0 $838,000 $1,283,000 6.4%
0.5 $840,000 $1,404,000 11.2%
0.9 $841,000 $1,507,000 15.8%
Assuming independence between risks that share a cause understates
the probability of the bad case by 2-3x.
WHAT TO READ FROM THE OUTPUT
NOT the mean alone
p50 / p80 / p95 plan against p80, not the median
P(specific bad outcome) directly decision-relevant
input-output correlation which uncertainty to go reduce
convergence tails need 10-100x the samples the mean does
DISTRIBUTION CHOICE
situation use
---------------------------------------------------------------------------
expert gives min / likely / max triangular or PERT
cannot be negative, long tail lognormal (durations, costs)
you have historical data empirical resampling (bootstrap)
count of rare events Poisson
a bounded proportion Beta
Key Takeaways
- Plugging average inputs into non-linear logic gives the wrong answer, and capacity limits, thresholds, and parallel paths are all non-linear.
- Specify a distribution for each uncertain input, encode dependencies, then propagate samples through the actual decision logic.
- Read percentiles and the probability of specific bad outcomes rather than the mean, and plan against the 80th percentile.
- Correlation between risks barely moves the mean and substantially thickens the tail, so independence assumptions understate disaster.
- Correlate sampled inputs with the outcome to find which uncertainty to spend effort reducing, and check convergence for tail statistics.
🧪 Practice
- Simulate revenue as price times the minimum of demand and capacity, and quantify how much the average-input estimate overstates the mean.
- Build a project model with three parallel tasks and report the probability of finishing by the sum of most-likely durations.
- Simulate two correlated cost risks at correlations of 0, 0.5, and 0.9 and compare the means against the 95th percentiles.
- Compare three capacity options on mean profit and probability of loss, and identify a case where the rankings disagree.
- Run a convergence check showing that a 99th percentile needs far more samples than the mean to stabilise.
- Interview: A plan says the project takes 47 days, built by summing each task's most likely duration. What is your concern and what do you produce instead? (Hint: what happens when three parallel tasks each have a long right tail?)
<a id="9-machine-learning-fundamentals"></a>