Causal Inference I — Randomized Experiments & Randomization Inference¶
The Rubin Causal Model, Fisher's exact test, and Neyman's repeated sampling — where causation is cheap¶
This is the entry point of the Causal Inference arc, and it starts where the whole framework is defined: the randomized experiment. Before any of the observational machinery — matching, instruments, discontinuities, difference-in-differences — a designed experiment is the setting where a causal effect can be identified with almost no assumptions, because the researcher, not nature, decides who is treated. Everything later in the arc is measured against this benchmark.
We introduce three foundations, each on a concrete example:
- The Rubin Causal Model (RCM) — potential outcomes $Y_i(1), Y_i(0)$, the fundamental problem of causal inference, SUTVA, and the reason randomization identifies the average treatment effect. Shown with a simulation where the truth is known.
- Fisher's randomization inference — the sharp null of no effect for anyone, tested exactly by enumerating the assignments that randomization could have produced. Introduced on Fisher's own Lady Tasting Tea and his agricultural Design of Experiments.
- Neyman's repeated-sampling framework — the complementary view: estimate the average effect and put a confidence interval on it, from the sampling distribution induced by re-randomization.
Two of the datasets are the historical originals — the Lady Tasting Tea (Fisher 1935) and Darwin's maize (the cross- vs self-fertilized corn that Fisher used to demonstrate the permutation test). Python-lead; an R companion follows with HistData, coin, ri2, and estimatr.
1. The Rubin Causal Model — potential outcomes and why randomization works¶
Each unit $i$ has two potential outcomes: $Y_i(1)$, the outcome it would show under treatment, and $Y_i(0)$, the outcome under control. The unit-level causal effect is $\tau_i = Y_i(1)-Y_i(0)$. The fundamental problem of causal inference (Holland, 1986) is that we only ever observe one of them — $Y_i = W_i Y_i(1) + (1-W_i)Y_i(0)$, where $W_i\in\{0,1\}$ is the treatment — so $\tau_i$ is never observed for any single unit. The other potential outcome is a missing counterfactual — the same missing-data logic that runs through the Missing Data arc, here made causal.
We can still learn the average treatment effect $\text{ATE}=\mathbb{E}[Y_i(1)-Y_i(0)]$ — if treatment is independent of the potential outcomes. Randomization guarantees exactly that: $W_i \perp (Y_i(1),Y_i(0))$, so the treated and control groups are, in expectation, identical in every pre-treatment characteristic — observed and unobserved. Then the simple difference in means is unbiased for the ATE. We also require SUTVA (no interference between units, and a single version of the treatment).
The simulation below makes the point that no observational method can fully escape: with a confounder $C$ driving both who gets treated and the outcome, the naive difference in means is badly biased; randomizing the same units removes the bias entirely. The known ATE is 4.
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
from scipy import stats
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
rng=np.random.default_rng(7)
N=2000; ATE=4.0
C=rng.normal(0,1,N) # confounder (e.g. baseline motivation/severity)
Y0=50+5*C+rng.normal(0,3,N) # potential outcome under control
Y1=Y0+ATE # constant treatment effect -> ATE known = 4
# (a) CONFOUNDED assignment: sicker/higher-C units more likely treated
p=1/(1+np.exp(-1.5*C)); Wc=(rng.uniform(size=N)<p).astype(int)
Yc=Wc*Y1+(1-Wc)*Y0
naive_conf=Yc[Wc==1].mean()-Yc[Wc==0].mean()
# (b) RANDOMIZED assignment on the SAME units: coin flip, ignores C
Wr=(rng.uniform(size=N)<0.5).astype(int)
Yr=Wr*Y1+(1-Wr)*Y0
naive_rand=Yr[Wr==1].mean()-Yr[Wr==0].mean()
print(f"True ATE = {ATE}")
print(f"(a) confounded assignment: diff-in-means = {naive_conf:.2f} -> BIASED by {naive_conf-ATE:+.2f} (treated had higher C)")
print(f"(b) randomized assignment: diff-in-means = {naive_rand:.2f} -> unbiased (C balanced by design)")
# randomization distribution of the estimator: re-randomize W many times, potential outcomes FIXED
draws=np.array([ (lambda w: (w*Y1+(1-w)*Y0)[w==1].mean()-(w*Y1+(1-w)*Y0)[w==0].mean())((rng.uniform(size=N)<0.5).astype(int)) for _ in range(4000)])
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].hist(C[Wc==1],bins=30,alpha=.55,color=RED,density=True,label="treated (confounded)")
ax[0].hist(C[Wc==0],bins=30,alpha=.55,color=BLUE,density=True,label="control (confounded)")
ax[0].hist(C[Wr==1],bins=30,histtype="step",lw=2,color="black",density=True,label="treated (randomized)")
ax[0].set_xlabel("confounder C"); ax[0].set_ylabel("density"); ax[0].set_title("Confounded assignment imbalances C; randomization does not"); ax[0].legend(fontsize=8)
ax[1].hist(draws,bins=40,color=GREEN,alpha=.75,density=True); ax[1].axvline(ATE,color=RED,lw=2,label=f"true ATE = {ATE}"); ax[1].axvline(draws.mean(),color="black",ls="--",label=f"mean of estimator = {draws.mean():.2f}")
ax[1].set_xlabel("difference-in-means over re-randomizations"); ax[1].set_title("Randomization distribution of the estimator (design-based, unbiased)"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"\nOver 4000 re-randomizations the diff-in-means averages {draws.mean():.2f} = the ATE: randomization makes the SIMPLE")
print("estimator unbiased, with a sampling distribution that comes entirely from the assignment mechanism -- Neyman's view (§4).")
True ATE = 4.0 (a) confounded assignment: diff-in-means = 9.21 -> BIASED by +5.21 (treated had higher C) (b) randomized assignment: diff-in-means = 3.94 -> unbiased (C balanced by design)
Over 4000 re-randomizations the diff-in-means averages 4.00 = the ATE: randomization makes the SIMPLE estimator unbiased, with a sampling distribution that comes entirely from the assignment mechanism -- Neyman's view (§4).
2. Fisher's sharp null and the exact randomization test — the Lady Tasting Tea¶
Fisher's inferential move needs no model and no large sample. His sharp null hypothesis is that treatment changes no unit's outcome at all: $H_0:\ Y_i(1)=Y_i(0)$ for every $i$. Under this null, both potential outcomes are known for everyone (they are equal to the one we observed), so we can compute the test statistic under every treatment assignment the randomization could have generated. That set of values is the exact null (randomization) distribution, and the $p$-value is simply the fraction of assignments giving a statistic as extreme as the one observed — exact, by construction.
The original illustration is the Lady Tasting Tea (Fisher, The Design of Experiments, 1935). Muriel Bristol claimed she could tell whether milk or tea was poured into the cup first. Fisher presented 8 cups — 4 milk-first, 4 tea-first, in random order — and told her there were 4 of each; she had to name the 4 milk-first cups. She got all 4 right. Under the null that she is guessing, the number of correct identifications follows a hypergeometric distribution (drawing 4 from 8 without replacement — the same distribution catalogued in the Statistical Distributions arc), and getting all 4 correct has probability $1/\binom{8}{4}=1/70$. This is Fisher's exact test.
from math import comb
from scipy.stats import hypergeom, fisher_exact
# Under H0 (guessing): X = # milk-first cups correctly named ~ Hypergeom(M=8, n=4 milk-first, N=4 picks)
M,n_milk,picks=8,4,4
xs=np.arange(0,5); pmf=hypergeom.pmf(xs,M,n_milk,picks)
p_exact=hypergeom.sf(3,M,n_milk,picks) # P(X>=4) = P(all 4 correct)
# same thing as Fisher's exact test on the 2x2 outcome table for a perfect score
p_fisher=fisher_exact([[4,0],[0,4]],alternative="greater")[1]
print("Null distribution of correct guesses (hypergeometric):")
for x,pr in zip(xs,pmf): print(f" {x} correct: P = {pr:.4f}"+(" <- observed" if x==4 else ""))
print(f"\nExact one-sided p-value P(all 4 correct | guessing) = {p_exact:.4f} = 1/70 = {1/70:.4f}")
print(f"scipy fisher_exact (package cross-check): {p_fisher:.4f}")
fig,ax=plt.subplots(figsize=(7.5,4))
bars=ax.bar(xs,pmf,color=[RED if x==4 else GREY for x in xs]); ax.set_xlabel("number of milk-first cups correctly identified"); ax.set_ylabel("probability under H0 (guessing)")
ax.set_title("Exact randomization distribution — the Lady Tasting Tea");
for x,pr in zip(xs,pmf): ax.text(x,pr+0.01,f"{pr:.3f}",ha="center",fontsize=8)
ax.text(4,pmf[4]+0.06,"observed\np = 1/70",ha="center",color=RED,fontsize=9)
plt.tight_layout(); plt.show()
print("The p-value is not an approximation -- it is the exact proportion of the 70 equally-likely arrangements in which a")
print("guesser matches the true score. Fisher's genius: the randomization itself supplies the reference distribution.")
Null distribution of correct guesses (hypergeometric): 0 correct: P = 0.0143 1 correct: P = 0.2286 2 correct: P = 0.5143 3 correct: P = 0.2286 4 correct: P = 0.0143 <- observed Exact one-sided p-value P(all 4 correct | guessing) = 0.0143 = 1/70 = 0.0143 scipy fisher_exact (package cross-check): 0.0143
The p-value is not an approximation -- it is the exact proportion of the 70 equally-likely arrangements in which a guesser matches the true score. Fisher's genius: the randomization itself supplies the reference distribution.
3. The real field trial — Darwin's maize (Fisher, Design of Experiments, 1935)¶
Fisher's textbook demonstration of the permutation test used Charles Darwin's data on cross- vs self-fertilized plants (The Effects of Cross and Self Fertilisation in the Vegetable Kingdom, 1876). Darwin grew 15 pairs of Zea mays (corn) plants; within each pair, one plant was cross-fertilized and the other self-fertilized, and the two were grown in the same pot under matched conditions. He measured final plant heights (in inches). The pairing is the design: because the two plants in a pot share soil, light, and water, the difference in height within a pair isolates the fertilization effect from pot-to-pot variation — a matched-pair experiment.
Fisher's paired randomization test treats the sign of each within-pair difference as the randomized quantity: under the sharp null of no fertilization effect, the "cross" and "self" labels within a pair are exchangeable, so each observed difference $d_i$ was equally likely to have come out $+d_i$ or $-d_i$. With 15 pairs there are exactly $2^{15}=32{,}768$ sign patterns — few enough to enumerate the entire null distribution exactly. The test statistic is the mean difference; the $p$-value is the fraction of sign patterns giving a mean at least as large. We compare it to the classical paired $t$-test (Fisher's own benchmark: $t=2.148$).
cross=np.array([23.5,12,21,22,19.125,21.5,22.125,20.375,18.25,21.625,23.25,21,22.125,23,12])
selff=np.array([17.375,20.375,20,20,18.375,18.625,18.625,15.25,16.5,18,16.25,18,12.75,15.5,18])
d=cross-selff; Nn=len(d); Tobs=d.mean()
print(f"Darwin's {Nn} pairs of Zea mays (heights in inches):")
print(pd.DataFrame({"pair":np.arange(1,Nn+1),"cross":cross,"self":selff,"diff":d}).to_string(index=False))
print(f"\nmean height advantage of cross-fertilized = {Tobs:.3f} inches; {int((d>0).sum())}/{Nn} pairs favour cross-fertilization")
# EXACT sign-flip randomization distribution: all 2^15 patterns
signs=((np.arange(2**Nn)[:,None]>>np.arange(Nn))&1)*2-1 # {+1,-1}^15
Tnull=(signs*np.abs(d)).mean(1)
p_one=np.mean(Tnull>=Tobs); p_two=np.mean(np.abs(Tnull)>=abs(Tobs))
tt=stats.ttest_rel(cross,selff)
print(f"\nExact permutation test (enumerated all {2**Nn} sign patterns):")
print(f" one-sided p = {p_one:.5f} two-sided p = {p_two:.5f}")
print(f"Paired t-test (classical benchmark): t = {tt.statistic:.4f}, two-sided p = {tt.pvalue:.5f} (Fisher reported t=2.148)")
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].bar(np.arange(1,Nn+1),d,color=[GREEN if x>0 else RED for x in d]); ax[0].axhline(0,color="k",lw=.7); ax[0].axhline(Tobs,color=BLUE,ls="--",label=f"mean = {Tobs:.2f} in")
ax[0].set_xlabel("pair"); ax[0].set_ylabel("height diff: cross − self (in)"); ax[0].set_title("Darwin's 15 matched pairs (13 favour cross-fertilization)"); ax[0].legend(fontsize=8)
ax[1].hist(Tnull,bins=60,color=GREY,alpha=.8,density=True); ax[1].axvline(Tobs,color=RED,lw=2,label=f"observed {Tobs:.2f}"); ax[1].axvline(-Tobs,color=RED,lw=1,ls=":")
ax[1].set_xlabel("mean difference under H0 (all $2^{15}$ sign flips)"); ax[1].set_ylabel("density"); ax[1].set_title(f"Exact randomization distribution — one-sided p = {p_one:.3f}"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"\nOnly {int(p_one*2**Nn)} of the {2**Nn} equally-likely sign patterns give a mean as large as Darwin's -> p = {p_one:.3f}.")
print("The permutation test needs no normality assumption yet agrees closely with the paired t-test -- Fisher's point exactly.")
Darwin's 15 pairs of Zea mays (heights in inches):
pair cross self diff
1 23.500 17.375 6.125
2 12.000 20.375 -8.375
3 21.000 20.000 1.000
4 22.000 20.000 2.000
5 19.125 18.375 0.750
6 21.500 18.625 2.875
7 22.125 18.625 3.500
8 20.375 15.250 5.125
9 18.250 16.500 1.750
10 21.625 18.000 3.625
11 23.250 16.250 7.000
12 21.000 18.000 3.000
13 22.125 12.750 9.375
14 23.000 15.500 7.500
15 12.000 18.000 -6.000
mean height advantage of cross-fertilized = 2.617 inches; 13/15 pairs favour cross-fertilization
Exact permutation test (enumerated all 32768 sign patterns):
one-sided p = 0.02634 two-sided p = 0.05267
Paired t-test (classical benchmark): t = 2.1480, two-sided p = 0.04970 (Fisher reported t=2.148)
Only 863 of the 32768 equally-likely sign patterns give a mean as large as Darwin's -> p = 0.026. The permutation test needs no normality assumption yet agrees closely with the paired t-test -- Fisher's point exactly.
4. Neyman's repeated-sampling framework — estimating how big, with a confidence interval¶
Fisher asks a yes/no question: is there any effect at all? Neyman (1923) asks the quantitative one: how large is the average effect, and how precisely do we know it? His target is the ATE, his estimator the difference in means (or, for paired data, the mean of the within-pair differences), and his inference comes from the sampling distribution induced by re-randomization — the same distribution we drew in §1. The standard error of the mean paired difference is $\text{SE}=s_d/\sqrt{n}$, giving the familiar $\hat\tau \pm 1.96\,\text{SE}$ interval. (For unpaired two-arm trials Neyman's variance $\frac{s_1^2}{n_1}+\frac{s_0^2}{n_0}$ is conservative — it slightly overstates uncertainty because the covariance of the two potential outcomes is unidentified.)
The two frameworks are complementary, not rival: Fisher tests the sharp null exactly (great for small samples and precise "did anything happen?"), while Neyman estimates the average effect with an interval (the effect size a decision needs). We report both on Darwin's data, and confirm the Neyman interval against the confidence interval a paired $t$-test produces.
tau=d.mean(); sd=d.std(ddof=1); se=sd/np.sqrt(Nn)
ci=(tau-1.96*se, tau+1.96*se)
tci=stats.t.interval(0.95,Nn-1,loc=tau,scale=se) # paired-t CI (package cross-check)
print("NEYMAN (estimation) on Darwin's maize:")
print(f" ATE estimate (mean paired diff) = {tau:.3f} inches")
print(f" SE = s_d/sqrt(n) = {sd:.3f}/sqrt({Nn}) = {se:.3f}")
print(f" 95% CI (normal) = [{ci[0]:.3f}, {ci[1]:.3f}]")
print(f" 95% CI (paired-t) = [{tci[0]:.3f}, {tci[1]:.3f}] (t-based, wider, small n)")
print(f"\nFISHER (testing): sharp-null exact one-sided p = {p_one:.3f} -- reject 'no effect for anyone'.")
print(f"NEYMAN (estimating): cross-fertilization adds ~{tau:.1f} inches on average, 95% CI excludes 0.")
fig,ax=plt.subplots(figsize=(8,3.2))
ax.errorbar([tau],[1],xerr=[[tau-ci[0]],[ci[1]-tau]],fmt="o",color=BLUE,capsize=6,ms=9,label="Neyman 95% CI (normal)")
ax.errorbar([tau],[0.7],xerr=[[tau-tci[0]],[tci[1]-tau]],fmt="s",color=GREEN,capsize=6,ms=8,label="paired-t 95% CI")
ax.axvline(0,color=RED,lw=1.5,ls="--",label="no effect"); ax.set_yticks([]); ax.set_ylim(0.4,1.4); ax.set_xlabel("average height advantage of cross-fertilization (inches)")
ax.set_title("Neyman estimation: the effect size and its uncertainty"); ax.legend(fontsize=8,loc="upper left")
plt.tight_layout(); plt.show()
print("Both frameworks agree the effect is real; Neyman adds the magnitude and a confidence interval, Fisher gives an exact")
print("small-sample p. Together they are the two lenses every later, harder design in this arc is trying to earn back.")
NEYMAN (estimation) on Darwin's maize: ATE estimate (mean paired diff) = 2.617 inches SE = s_d/sqrt(n) = 4.718/sqrt(15) = 1.218 95% CI (normal) = [0.229, 5.004] 95% CI (paired-t) = [0.004, 5.229] (t-based, wider, small n) FISHER (testing): sharp-null exact one-sided p = 0.026 -- reject 'no effect for anyone'. NEYMAN (estimating): cross-fertilization adds ~2.6 inches on average, 95% CI excludes 0.
Both frameworks agree the effect is real; Neyman adds the magnitude and a confidence interval, Fisher gives an exact small-sample p. Together they are the two lenses every later, harder design in this arc is trying to earn back.
5. Summary¶
The randomized experiment is the foundation of the whole arc because it identifies a causal effect by design, not assumption. We built up the three pillars:
- The Rubin Causal Model frames causation as a missing-data problem — each unit has two potential outcomes and we see one. Randomization makes treatment independent of those potential outcomes, so the simple difference in means is unbiased for the ATE; the simulation showed a confounder biasing the naive estimate and randomization erasing that bias.
- Fisher's randomization inference tests the sharp null of no effect for anyone by enumerating the assignments randomization could have produced — exact, model-free, small-sample-valid. The Lady Tasting Tea gave $p=1/70$ from the hypergeometric; Darwin's maize gave an exact permutation $p\approx0.026$ that closely matched the paired $t$-test Fisher used it to justify.
- Neyman's framework estimates the average effect and its sampling variability, yielding a confidence interval (cross-fertilization $\approx 2.6$ inches, CI excluding zero). Fisher tests, Neyman estimates — the two complementary lenses of experimental analysis.
Cross-links. The potential-outcomes-as-missing-data view is the causal face of the Missing Data arc; the hypergeometric null is from the Statistical Distributions catalogue; the permutation logic recurs whenever a reference distribution is built by resampling (it reappears as the purged validation discipline of the ML arc). What comes next: every subsequent subsection relaxes the luxury of randomization. When treatment is not randomly assigned, we must assume something to stand in for it — and the next notebook, Potential Outcomes & Matching, makes the first such assumption (unconfoundedness) explicit and shows how to condition on covariates to approximate the experiment we could not run.