Randomisation as the Benchmark

13.1 Why a coin flip removes selection bias (40 min)

Concept. Day 11 showed: naive difference = ATT + selection bias, where the bias is E[Y₀ | T=1] − E[Y₀ | T=0]. Selection bias exists because the treated and untreated would have differed anyway.

If treatment is assigned by a coin flip, it can't depend on anything about the customer, including their potential outcomes. Then on average:

E[Y₀ | T=1] = E[Y₀ | T=0]      → selection bias = 0
E[Y₁ − Y₀ | T=1] = E[Y₁ − Y₀]  → ATT = ATE

So the simple difference in means estimates the ATE, with no model and no adjustment. That's why randomised experiments are the benchmark every other method in this course is compared against. In Modules 4–6, the question for every observational method will be: how close does it get to what an experiment would have told us?

"On average" matters. In any single experiment the groups differ a little by chance. Randomisation doesn't make the groups identical; it makes the differences random, so they average out and can be quantified with standard errors.

13.2 Replaying the email as a randomised experiment (50 min)

The truth file lets us do what no real company can: rerun the spring campaign with a coin flip instead of targeting, on the same customers.

Create day13/one_randomization.py:

# day13/one_randomization.py   (python -m day13.one_randomization)
"""What if marketing had flipped a coin instead of targeting? The truth file lets us replay it."""
import numpy as np
import polars as pl

po = pl.read_parquet("data/truth/potential_outcomes.parquet")
rng = np.random.default_rng(2026)

# 1. Assign the email at random to the same number of customers (1,279)
z = np.zeros(po.height, dtype=np.int8)
z[rng.choice(po.height, size=1279, replace=False)] = 1
rand = po.with_columns(z=pl.Series(z)).with_columns(
    y=pl.when(pl.col("z") == 1).then(pl.col("y1")).otherwise(pl.col("y0"))   # what we would observe
)

# 2. The simple difference in means
diff = rand.filter(pl.col("z") == 1)["y"].mean() - rand.filter(pl.col("z") == 0)["y"].mean()
print(f"Randomised difference in means:  {diff:+.3f}")
print(f"True ATE:                        {po['tau'].mean():+.3f}")

# 3. Why it works: without the email, the two random groups look alike
bias = rand.filter(pl.col("z") == 1)["y0"].mean() - rand.filter(pl.col("z") == 0)["y0"].mean()
print(f"Selection bias with randomisation: {bias:+.3f}  (vs +1.268 with targeting)")

python -m day13.one_randomization

Randomised difference in means:  +0.391
True ATE:                        +0.319
Selection bias with randomisation: +0.071  (vs +1.268 with targeting)

The estimate is +0.391 instead of +1.711. It isn't exactly 0.319: in this particular draw, the randomly emailed group happened to contain slightly heavier shoppers (+0.071). That's chance imbalance, not selection bias, and the next section shows it averages out.

13.3 Right on average: the sampling distribution (55 min)

Concept. An estimator's sampling distribution is the spread of estimates you'd get across many repetitions of the same design. It tells you two things: whether the estimator is right on average (unbiased), and how far off a single estimate can be (standard error).

Create day13/many_randomizations.py:

# day13/many_randomizations.py   (python -m day13.many_randomizations)
"""Repeat the coin flip 2,000 times: randomisation is right ON AVERAGE, and we can measure the noise."""
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import polars as pl

po = pl.read_parquet("data/truth/potential_outcomes.parquet")
y0, y1 = po["y0"].to_numpy(), po["y1"].to_numpy()
n, n_treated, reps = po.height, 1279, 2000
rng = np.random.default_rng(13)

estimates = np.empty(reps)
for r in range(reps):
    z = np.zeros(n, dtype=bool)
    z[rng.choice(n, size=n_treated, replace=False)] = True
    estimates[r] = y1[z].mean() - y0[~z].mean()

ate = (y1 - y0).mean()
print(f"True ATE                      {ate:.3f}")
print(f"Average of 2,000 estimates    {estimates.mean():.3f}")
print(f"Standard deviation (noise)    {estimates.std():.3f}")
print(f"Middle 95% of estimates       [{np.percentile(estimates, 2.5):.3f}, {np.percentile(estimates, 97.5):.3f}]")
print("Targeted (real) naive estimate: 1.711")

fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(estimates, bins=50, color="#3987e5", alpha=0.85)
ax.axvline(ate, color="black", linewidth=2, label=f"true ATE = {ate:.2f}")
ax.set_xlabel("Difference in means across 2,000 random assignments (orders per customer)")
ax.set_ylabel("Number of assignments")
ax.legend()
fig.tight_layout()
fig.savefig("figures/day13_randomization_distribution.png", dpi=150)
print("Saved figures/day13_randomization_distribution.png")

True ATE                      0.319
Average of 2,000 estimates    0.317
Standard deviation (noise)    0.073
Middle 95% of estimates       [0.180, 0.461]
Targeted (real) naive estimate: 1.711
Saved figures/day13_randomization_distribution.png

Open the figure. The 2,000 estimates form a bell curve centred on the true ATE. Its standard deviation, 0.073, is the standard error of the experiment. The real, targeted comparison (1.711) is about 19 standard errors away: no amount of noise explains it.

Try it. How does the noise change with the number of emailed customers? Rerun the simulation with 200, 1,279 and 2,000 emailed customers (1,000 repetitions each).

import numpy as np
import polars as pl

po = pl.read_parquet("data/truth/potential_outcomes.parquet")
y0, y1 = po["y0"].to_numpy(), po["y1"].to_numpy()
rng = np.random.default_rng(13)

for n_treated in [200, 1279, 2000]:
    est = []
    for _ in range(1000):
        z = np.zeros(po.height, dtype=bool)
        z[rng.choice(po.height, size=n_treated, replace=False)] = True
        est.append(y1[z].mean() - y0[~z].mean())
    est = np.array(est)
    print(f"{n_treated:>5} emailed: mean {est.mean():.3f}, SD {est.std():.3f}, "
          f"share of estimates below 0.1: {(est < 0.1).mean():.1%}")

  200 emailed: mean 0.320, SD 0.163, share of estimates below 0.1: 7.7%
 1279 emailed: mean 0.313, SD 0.073, share of estimates below 0.1: 0.3%
 2000 emailed: mean 0.314, SD 0.065, share of estimates below 0.1: 0.1%

13.4 One experiment, analysed properly (40 min)

Concept. In practice you run one experiment. Two tools give its uncertainty:

1. Regression of the outcome on the treatment reproduces the difference in means, and gives a standard error and confidence interval. Use robust (HC2) standard errors: they don't assume equal variance in both groups.
2. Adding pre-treatment covariates doesn't change what is estimated (randomisation already removed bias) but absorbs noise, giving a tighter interval. Module 3 develops this as CUPED.

Create day13/regression_ci.py (it uses pre_orders from Module 1's published customer_features):

# day13/regression_ci.py   (python -m day13.regression_ci)
"""One randomised experiment, analysed properly: estimate, standard error, confidence interval."""
import numpy as np
import polars as pl
import statsmodels.formula.api as smf

po = pl.read_parquet("data/truth/potential_outcomes.parquet")
features = pl.read_parquet("data/published/customer_features.parquet").select("customer_id", "pre_orders")

rng = np.random.default_rng(2026)
z = np.zeros(po.height, dtype=np.int8)
z[rng.choice(po.height, size=1279, replace=False)] = 1
df = (po.with_columns(z=pl.Series(z))
        .with_columns(y=pl.when(pl.col("z") == 1).then(pl.col("y1")).otherwise(pl.col("y0")))
        .join(features, on="customer_id")
        .to_pandas())

# 1. Difference in means as a regression: the coefficient on z IS the difference
m1 = smf.ols("y ~ z", data=df).fit(cov_type="HC2")      # HC2 = robust standard errors
# 2. Add a pre-treatment covariate: same target, less noise
m2 = smf.ols("y ~ z + pre_orders", data=df).fit(cov_type="HC2")

for name, m in [("y ~ z", m1), ("y ~ z + pre_orders", m2)]:
    lo, hi = m.conf_int().loc["z"]
    print(f"{name:<22} effect {m.params['z']:+.3f}  SE {m.bse['z']:.3f}  95% CI [{lo:+.3f}, {hi:+.3f}]")
print(f"True ATE {po['tau'].mean():+.3f}")

y ~ z                  effect +0.391  SE 0.073  95% CI [+0.249, +0.534]
y ~ z + pre_orders     effect +0.355  SE 0.064  95% CI [+0.230, +0.481]
True ATE +0.319

• The robust standard error (0.073) matches the spread across 2,000 replays in 13.3: one experiment can measure its own uncertainty.
• Adding pre_orders shrinks the standard error by 12% and moves the estimate closer to the truth, because it corrects part of this draw's chance imbalance.
• Both intervals contain the true ATE.

Covariates in experiments: only ever add variables measured before treatment. Day 17 shows what happens otherwise.

13.5 Day 13 mini-project: balance and a permutation test (70 min)

Goal. Two checks every experiment report should contain:

1. A balance table: are treated and control similar on pre-treatment characteristics? Measured with standardised mean differences (SMD): the gap in means divided by the pooled standard deviation, so different variables are comparable. |SMD| > 0.1 is the usual warning sign.
2. A permutation test (randomisation inference): if the email did nothing, the treatment labels are arbitrary. Shuffle them thousands of times; how often is the shuffled difference as large as the real one? That share is a p-value that needs no distributional assumptions.

Create day13/balance_and_permutation.py:

# day13/balance_and_permutation.py   (python -m day13.balance_and_permutation)
"""Mini-project: (1) balance under targeting vs randomisation, (2) a permutation test."""
import numpy as np
import polars as pl

features = pl.read_parquet("data/published/customer_features.parquet")
po = pl.read_parquet("data/truth/potential_outcomes.parquet").select("customer_id", "y0", "y1")
df = features.join(po, on="customer_id")

rng = np.random.default_rng(2026)
z = np.zeros(df.height, dtype=np.int8)
z[rng.choice(df.height, size=1279, replace=False)] = 1
df = df.with_columns(z=pl.Series(z))

COVARIATES = ["pre_orders", "pre_revenue", "is_loyalty", "age", "tenure_days"]


def smd(data: pl.DataFrame, group: str, col: str) -> float:
    """Standardised mean difference: gap in means, in units of the pooled standard deviation."""
    a = data.filter(pl.col(group) == 1)[col]
    b = data.filter(pl.col(group) == 0)[col]
    return (a.mean() - b.mean()) / np.sqrt((a.var() + b.var()) / 2)


# 1. Balance table
table = pl.DataFrame({
    "covariate": COVARIATES,
    "smd_targeted": [round(smd(df, "treated", c), 3) for c in COVARIATES],
    "smd_randomised": [round(smd(df, "z", c), 3) for c in COVARIATES],
})
print(table)
print("Rule of thumb: |SMD| above 0.1 signals imbalance.\n")

# 2. Permutation test for the randomised experiment
y = np.where(df["z"].to_numpy() == 1, df["y1"].to_numpy(), df["y0"].to_numpy())
zz = df["z"].to_numpy().astype(bool)
observed = y[zz].mean() - y[~zz].mean()

reps = 5000
null = np.empty(reps)
for r in range(reps):
    shuffled = rng.permutation(zz)             # pretend the email did nothing: labels are arbitrary
    null[r] = y[shuffled].mean() - y[~shuffled].mean()

p_value = (np.abs(null) >= abs(observed)).mean()
print(f"Observed difference      {observed:+.3f}")
shown = f"{p_value:.4f}" if p_value > 0 else f"< {1 / reps:.4f}"
print(f"Permutation p-value      {shown}  (none of {reps:,} shuffles was as extreme)" if p_value == 0
      else f"Permutation p-value      {shown}  ({reps:,} shuffles)")

│ covariate   ┆ smd_targeted ┆ smd_randomised │
│ pre_orders  ┆ 0.627        ┆ 0.037          │
│ pre_revenue ┆ 0.461        ┆ 0.005          │
│ is_loyalty  ┆ 0.621        ┆ 0.031          │
│ age         ┆ 0.038        ┆ 0.025          │
│ tenure_days ┆ -0.04        ┆ 0.01           │
Rule of thumb: |SMD| above 0.1 signals imbalance.

Observed difference      +0.391
Permutation p-value      < 0.0002  (none of 5,000 shuffles was as extreme)

Under targeting, three covariates are badly imbalanced (SMD above 0.6); under randomisation, every SMD is below 0.04. The randomised experiment's effect is far too large to be chance.

Finish by adding a short section to docs/ESTIMANDS.md: "What if marketing had randomised?" In four sentences: what the experiment would have estimated, how precisely, what it would have cost (emails sent to customers marketing considered less promising), and why that cost may be worth paying.

13.6 Save your work with Git (15 min)

git add .
git commit -m "Day 13: randomisation benchmark, sampling distribution, permutation test"
git push

Coming next (Days 14–17). Randomisation is the gold standard, but NorthMart's email wasn't randomised. Days 14–17 introduce DAGs, the tool for deciding what to adjust for when randomisation isn't available, and for spotting the adjustments that make things worse.

Day 13 checklist

• I can explain why randomisation makes selection bias zero on average, and ATT equal ATE
• I can simulate a sampling distribution and read bias and standard error from it
• I can analyse one experiment with robust standard errors and pre-treatment covariates
• I can build a balance table with standardised mean differences
• I can run a permutation test and interpret its p-value