import matplotlib
matplotlib.use("Agg") # headless: render to buffers, never to a screen
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
rng = np.random.default_rng(42) # one generator, pinned, drives the whole notebook
plt.rcParams["figure.figsize"] = (9, 5)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3
print("numpy", np.__version__)
print("pandas", pd.__version__)nb3.1 — Potential Outcomes & Selection Bias: Playing God with the Science Table
Week 3, Chapter 3.1. Companion notebook.
In Chapter 3.1 we proved an exact identity with pencil and paper:
\underbrace{\mathbb{E}[Y_i\mid D_i=1]-\mathbb{E}[Y_i\mid D_i=0]}_{\text{naive difference }\Delta} \;=\; \underbrace{\text{ATT}}_{\text{what we want}} \;+\; \underbrace{\mathbb{E}[Y_i(0)\mid D_i=1]-\mathbb{E}[Y_i(0)\mid D_i=0]}_{\text{selection bias}}.
In the real world we can never check this, because for each person we see only one of the two potential outcomes — the Fundamental Problem of Causal Inference. The other is the missing counterfactual, and the selection-bias term is built entirely out of missing pieces.
But this is a simulation, and a simulation lets us play God. We will write down the full science table — both Y_i(0) and Y_i(1) for every single unit — so that the true ATE, the true ATT, and the selection-bias term are all things we can compute exactly. Then we hide one column the way reality does, run the naive estimator, and watch it lie by precisely the amount the decomposition predicts. The plan:
- Build Maya’s world. Simulate applicants with a hidden “savviness” trait that drives the untreated approval Y_i(0); give everyone a known treatment effect; reveal both potential outcomes.
- Induce selection. Let savvier applicants self-select into the financial-literacy program, so D_i is correlated with Y_i(0). Verify \Delta = \text{ATT} + \text{selection bias} on the dollar using the hidden counterfactual column, and reproduce the chapter’s 0.07 + 0.14 = 0.21 illustration.
- Watch the naive estimator lie, and dial selection strength up and down to grow and shrink the bias.
- Turn on the coin. Re-assign D_i by a fair random draw and watch the selection-bias term collapse to sampling noise around zero, so the naive difference recovers the ATE.
- Your Turn: change the selection rule, predict the sign and size of the bias before running.
1. Build the full science table: both parallel worlds at once
The trick we are about to pull. For every applicant we generate both Y_i(0) (approval if she does not enroll) and Y_i(1) (approval if she does). Real data shows only the one the treatment switch selects; we keep both, so we can compute the truths reality hides.
Each applicant i has a hidden, standardized savviness trait s_i — a blend of organization, credit awareness, and follow-through. We make s_i drive the untreated potential outcome: savvier people are more approvable even with no program at all. We model approval on a continuous “approvability” scale first (think a latent credit-readiness score), which keeps the algebra transparent; later we will also show the same story for a binary approve/deny outcome.
We hand every unit a known individual treatment effect \tau_i. To keep the chapter’s ATT story clean we start with a constant effect \tau_i = \tau = 0.07 for everyone — so the program lifts approvability by exactly 7 points regardless of who you are. With a constant effect the ATT equals the ATE by construction, which means any gap between the naive number and 0.07 is pure selection bias, nowhere to hide.
N = 200_000 # a big population so Monte Carlo noise is tiny
TAU_TRUE = 0.07 # known, constant individual treatment effect (the chapter's ATT)
def build_science_table(n, tau, rng):
'''Generate the FULL science table: savviness s, and BOTH potential outcomes Y0, Y1.
Y0 = baseline approvability driven by savviness (this is what differs across people).
Y1 = Y0 + tau (constant treatment effect, so ATT == ATE by construction).
Returns a DataFrame with the god's-eye columns we could never see in real data.
'''
s = rng.normal(0.0, 1.0, size=n) # hidden savviness trait, standardized
# Untreated potential outcome: centered near the chapter's control mean 0.41,
# rising with savviness, plus idiosyncratic noise. (Continuous "approvability".)
y0 = 0.41 + 0.20 * s + rng.normal(0.0, 0.10, size=n)
y1 = y0 + tau # constant effect: everyone gains tau
return pd.DataFrame({"s": s, "Y0": y0, "Y1": y1})
sci = build_science_table(N, TAU_TRUE, rng)
sci.head()With both columns in hand, the estimands are not “targets we hope to recover” — they are just column statistics. The ATE is the mean of \tau_i = Y_i(1)-Y_i(0) over everyone; the ATT is the same mean restricted to the treated (we have not assigned treatment yet, so for now ATE is all we can read). Because the effect is constant, \tau_i is the same number for every row.
tau_i = sci["Y1"] - sci["Y0"] # individual treatment effects (the hidden middle column)
ATE_true = tau_i.mean() # average over the WHOLE population
print(f"True individual effect tau_i (constant): {tau_i.iloc[0]:+.4f}")
print(f"std of tau_i across population : {tau_i.std():.2e} (zero up to float noise)")
print(f"True ATE = E[Y1 - Y0] : {ATE_true:+.4f}")
assert abs(ATE_true - TAU_TRUE) < 1e-9, "constant-effect ATE must equal tau exactly"
print("\nCHECK PASSED: with a constant effect, the ATE equals tau on the nose.")2. Induce selection: who signs up for a free six-week evening course?
Now the move that breaks the naive comparison. We do not assign the program randomly. We let people self-select, exactly as Maya’s credit union did: the savvier you are, the more likely you enroll. Concretely, enrollment probability rises with savviness s_i through a logistic rule,
\Pr(D_i = 1 \mid s_i) = \frac{1}{1 + e^{-(\alpha + \gamma\, s_i)}},
and then we flip a biased coin for each person. The slope \gamma is the selection strength: bigger \gamma means savviness matters more for who enrolls, which means the treated group is more positively selected on Y_i(0) — and that is precisely the engine of the bias. Crucially, s_i drives both enrollment and Y_i(0): it is the confounder the chapter warned about. The intercept \alpha just sets roughly what fraction of the population enrolls.
def assign_self_selected(sci, gamma, alpha, rng):
'''Treatment self-selected on savviness via a logistic rule. Returns a NEW frame (a .copy())
with D, the observed outcome Yobs, and the propensity p attached.'''
df = sci.copy() # never mutate the science table in place
logit = alpha + gamma * df["s"].to_numpy()
p = 1.0 / (1.0 + np.exp(-logit)) # P(D=1 | s)
D = (rng.random(len(df)) < p).astype(int) # biased-coin draw per person
df["p"] = p
df["D"] = D
# Observation rule: Yobs = D*Y1 + (1-D)*Y0 -- reality shows only the selected column.
df["Yobs"] = np.where(D == 1, df["Y1"], df["Y0"])
return df
GAMMA = 1.2 # selection strength: savvier -> much likelier to enroll
ALPHA = -0.2 # intercept -> roughly 45-48% enroll
maya = assign_self_selected(sci, GAMMA, ALPHA, rng)
share_treated = maya["D"].mean()
print(f"Share enrolled (D=1): {share_treated:.3f}")
print(f"Mean savviness, enrollees (D=1) : {maya.loc[maya.D == 1, 's'].mean():+.3f}")
print(f"Mean savviness, non-enrollees (D=0): {maya.loc[maya.D == 0, 's'].mean():+.3f}")
print("Enrollees are the savvier crowd -- selection is now baked in.")2.2 Reproduce the chapter’s Maya numbers (0.62 / 0.55 / 0.41 → 0.07 + 0.14 = 0.21)
The chapter tells Maya’s story with a specific set of approval rates: enrollees approve at 0.62 with the program and would have approved at 0.55 without it; non-enrollees approve at 0.41. That gives \text{ATT}=0.62-0.55=0.07, \text{selection bias}=0.55-0.41=0.14, and a naive gap of 0.21. Those are binary approval rates, so let us reproduce them in a binary world to show the decomposition is not an artifact of our continuous scale. We target the three numbers directly: build a population whose treated group has \mathbb{E}[Y(0)\mid D=1]=0.55 and whose controls have \mathbb{E}[Y(0)\mid D=0]=0.41, with a constant +0.07 treatment effect on the approval probability.
def build_binary_maya(n_each, p0_treated, p0_control, tau, rng):
'''Two groups of binary potential outcomes engineered to hit the chapter's rates.
p0_treated = P(approve | no program) among those who WILL enroll (-> 0.55)
p0_control = P(approve | no program) among those who will NOT (-> 0.41)
tau = additive lift in approval probability from the program (-> 0.07).
'''
rows = []
for D, p0, k in [(1, p0_treated, n_each), (0, p0_control, n_each)]:
y0 = (rng.random(k) < p0).astype(int) # Y(0): approve w/o program
p1 = np.clip(p0 + tau, 0.0, 1.0) # Y(1) prob = Y(0) prob + tau
y1 = (rng.random(k) < p1).astype(int) # Y(1): approve w/ program
sub = pd.DataFrame({"D": D, "Y0": y0, "Y1": y1})
rows.append(sub)
df = pd.concat(rows, ignore_index=True)
df["Yobs"] = np.where(df["D"] == 1, df["Y1"], df["Y0"]) # observation rule
return df
bin_maya = build_binary_maya(300_000, p0_treated=0.55, p0_control=0.41, tau=0.07, rng=rng)
bt = bin_maya[bin_maya["D"] == 1]
bc = bin_maya[bin_maya["D"] == 0]
b_E_Y1_D1 = bt["Y1"].mean()
b_E_Y0_D1 = bt["Y0"].mean() # hidden counterfactual
b_E_Y0_D0 = bc["Y0"].mean()
b_ATT = b_E_Y1_D1 - b_E_Y0_D1
b_bias = b_E_Y0_D1 - b_E_Y0_D0
b_naive = bt["Yobs"].mean() - bc["Yobs"].mean()
print(f"E[Y(1)|D=1] (observed enrollee rate) : {b_E_Y1_D1:.3f} (chapter: 0.62)")
print(f"E[Y(0)|D=1] (enrollee counterfactual) : {b_E_Y0_D1:.3f} (chapter: 0.55)")
print(f"E[Y(0)|D=0] (non-enrollee rate) : {b_E_Y0_D0:.3f} (chapter: 0.41)")
print(f"ATT : {b_ATT:.3f} (chapter: 0.07)")
print(f"selection bias : {b_bias:.3f} (chapter: 0.14)")
print(f"naive gap : {b_naive:.3f} (chapter: 0.21)")# Reproduce the chapter's three numbers to within Monte Carlo noise (large N, so tight).
assert abs(b_E_Y1_D1 - 0.62) < 0.01, "enrollee observed rate should be ~0.62"
assert abs(b_ATT - 0.07) < 0.01, "ATT should be ~0.07"
assert abs(b_bias - 0.14) < 0.01, "selection bias should be ~0.14"
assert abs(b_naive - 0.21) < 0.01, "naive gap should be ~0.21"
assert abs(b_naive - (b_ATT + b_bias)) < 1e-9, "decomposition holds exactly"
print("CHECK PASSED: reproduced Maya's 0.62 / 0.55 / 0.41 -> ATT 0.07 + bias 0.14 = naive 0.21.")3. Watch the naive estimator lie, and dial the selection strength
Back to the continuous world. The naive difference overstates the truth because the treated were already more approvable. The size of that overstatement is controlled by \gamma, the selection strength. Sweep \gamma from 0 (enrollment unrelated to savviness) up to strong positive selection, and for each value compute the true selection bias (using the hidden column) and the naive estimate. We expect: at \gamma=0 the bias is ~zero and the naive number nails \tau=0.07; as \gamma grows, the bias grows and the naive number inflates by exactly that amount.
gammas = [0.0, 0.3, 0.6, 1.2, 2.0, 3.0]
rows = []
for g in gammas:
d = assign_self_selected(sci, gamma=g, alpha=ALPHA, rng=rng)
t = d[d["D"] == 1]
c = d[d["D"] == 0]
bias_g = t["Y0"].mean() - c["Y0"].mean() # true selection bias (hidden column)
att_g = (t["Y1"] - t["Y0"]).mean() # true ATT (= tau, constant effect)
naive_g = t["Yobs"].mean() - c["Yobs"].mean() # what Maya would report
rows.append({
"gamma": g,
"share_treated": d["D"].mean(),
"selection_bias": bias_g,
"ATT": att_g,
"naive_diff": naive_g,
"naive - (ATT+bias)": naive_g - (att_g + bias_g),
})
sweep = pd.DataFrame(rows)
print(sweep.to_string(index=False, float_format=lambda v: f"{v:.4f}"))
assert sweep["naive - (ATT+bias)"].abs().max() < 1e-9, "decomposition holds at every gamma"
assert abs(sweep.loc[sweep.gamma == 0.0, "selection_bias"].iloc[0]) < 0.01, "gamma=0 -> ~no bias"
assert sweep.loc[sweep.gamma >= 0.3, "selection_bias"].is_monotonic_increasing, \
"bias should grow with selection strength"
print("\nCHECK PASSED: bias is ~0 at gamma=0, grows with gamma, and naive == ATT+bias throughout.")fig, ax = plt.subplots()
ax.plot(sweep["gamma"], sweep["naive_diff"], "o-", color="crimson",
label=r"naive difference $\hat\Delta$")
ax.plot(sweep["gamma"], sweep["selection_bias"], "s--", color="darkorange",
label="selection-bias term (hidden)")
ax.axhline(TAU_TRUE, color="green", ls=":", label=r"true ATT = ATE = $0.07$")
ax.set_xlabel(r"selection strength $\gamma$ (how strongly savviness drives enrollment)")
ax.set_ylabel("percentage-point effect on approvability")
ax.set_title("The naive number = true effect + a selection bias that grows with self-selection")
ax.legend()
fig.tight_layout()
print("At gamma=0 the red curve sits on 0.07. As gamma rises, naive = 0.07 + the orange bias.")4. Turn on the coin: randomization zeroes the selection-bias term
Now the gold standard. Keep the exact same population — the same science table, the same savviness, the same potential outcomes — but throw away self-selection and assign the program by a fair coin. The coin does not know or care about savviness, so the treated and control groups become, in expectation, the same kind of people: \big(Y_i(0),Y_i(1)\big)\perp\!\!\!\perp D_i. That independence forces \mathbb{E}[Y(0)\mid D=1]=\mathbb{E}[Y(0)\mid D=0], so the selection-bias term collapses to sampling noise around zero, ATT collapses into ATE, and the naive difference reads the true effect straight off.
def assign_randomized(sci, p_treat, rng):
'''Treatment by a fair (or fixed-probability) coin, independent of everything. Returns a .copy().'''
df = sci.copy()
df["D"] = (rng.random(len(df)) < p_treat).astype(int) # coin ignores savviness entirely
df["Yobs"] = np.where(df["D"] == 1, df["Y1"], df["Y0"])
return df
rct = assign_randomized(sci, p_treat=0.5, rng=rng)
rt = rct[rct["D"] == 1]
rc = rct[rct["D"] == 0]
rct_bias = rt["Y0"].mean() - rc["Y0"].mean() # should be ~0 now
rct_att = (rt["Y1"] - rt["Y0"]).mean()
rct_naive = rt["Yobs"].mean() - rc["Yobs"].mean()
print(f"Mean savviness, treated : {rt['s'].mean():+.4f}")
print(f"Mean savviness, control : {rc['s'].mean():+.4f} (coin balances the confounder)")
print("-" * 52)
print(f"selection bias (randomized) : {rct_bias:+.4f} <- collapsed to sampling noise")
print(f"true ATE : {ATE_true:+.4f}")
print(f"naive difference : {rct_naive:+.4f} <- now recovers the ATE")
print(f"|naive - ATE| : {abs(rct_naive - ATE_true):.4f}")# Under randomization the bias is sampling noise (shrinks with N); with N=200k it should be tiny.
assert abs(rct_bias) < 0.01, "randomization should zero the selection bias (up to MC noise)"
assert abs(rct_naive - ATE_true) < 0.01, "naive difference should recover the ATE under randomization"
print("CHECK PASSED: random assignment -> selection bias ~ 0 and naive difference ~ ATE.")4.1 Side by side: self-selection vs. the coin, on the identical population
The cleanest way to see what randomization buys is to put the two assignment mechanisms next to each other on the same science table. Same people, same potential outcomes, same true effect of 0.07 — only the rule for who gets treated changes. Self-selection inflates the naive number; the coin does not.
compare = pd.DataFrame([
{"assignment": "self-selected (gamma=1.2)",
"selection_bias": selection_bias, "ATT_or_ATE": ATT_true, "naive_diff": naive_diff},
{"assignment": "randomized (fair coin)",
"selection_bias": rct_bias, "ATT_or_ATE": rct_att, "naive_diff": rct_naive},
])
compare["truth"] = TAU_TRUE
print(compare.to_string(index=False, float_format=lambda v: f"{v:+.4f}"))
fig, ax = plt.subplots()
xpos = np.arange(len(compare))
ax.bar(xpos - 0.2, compare["naive_diff"], width=0.4, color="crimson", label="naive difference")
ax.bar(xpos + 0.2, compare["selection_bias"], width=0.4, color="darkorange", label="selection bias")
ax.axhline(TAU_TRUE, color="green", ls=":", label=r"truth $=0.07$")
ax.set_xticks(xpos)
ax.set_xticklabels(compare["assignment"], rotation=8)
ax.set_ylabel("percentage-point effect")
ax.set_title("Same population, two assignment rules: the coin kills the bias the choice created")
ax.legend()
fig.tight_layout()
print("Left bars: self-selection inflates naive far above 0.07. Right bars: the coin lands on 0.07.")Your Turn — change the selection rule, predict the bias before you run
Everything above used positive selection: savvier (and so more approvable) people enrolled more, which made the selection bias positive and the naive number an overstatement. But the sign of the bias is the sign of the correlation between Y_i(0) and D_i — and you control that correlation through the selection rule. Flip the rule and you flip the bias.
The selection bias is \mathbb{E}[Y_i(0)\mid D_i=1]-\mathbb{E}[Y_i(0)\mid D_i=0]: the treated group’s average untreated outcome minus the control group’s. Since Y_i(0) rises with savviness s_i in our world, the bias is positive when the treated are savvier on average, and negative when the treated are less savvy on average.
Predict first. Suppose the program instead recruits the people who need it most — a caseworker steers low-savviness applicants in, so enrollment falls with s_i. Use \gamma = -1.2 (negative selection strength). Before running the cell:
- Will the treated group be more or less savvy than the controls? So is \mathbb{E}[Y_i(0)\mid D_i=1] above or below \mathbb{E}[Y_i(0)\mid D_i=0]?
- Will the selection bias be positive or negative?
- Will the naive difference overstate or understate the true +0.07 effect? By roughly how much — should |\text{bias}| be about the same size as the \gamma=+1.2 case (it was about +0.19), since we only flipped the sign?
Write your three predictions down, then run.
# Negative selection: a caseworker steers LOW-savviness applicants in. We flip gamma's sign.
neg = assign_self_selected(sci, gamma=-1.2, alpha=ALPHA, rng=rng)
nt = neg[neg["D"] == 1]
nc = neg[neg["D"] == 0]
neg_E_Y0_D1 = nt["Y0"].mean()
neg_E_Y0_D0 = nc["Y0"].mean()
neg_bias = neg_E_Y0_D1 - neg_E_Y0_D0
neg_att = (nt["Y1"] - nt["Y0"]).mean()
neg_naive = nt["Yobs"].mean() - nc["Yobs"].mean()
print(f"Mean savviness, treated : {nt['s'].mean():+.3f}")
print(f"Mean savviness, control : {nc['s'].mean():+.3f} (now the treated are LESS savvy)")
print("-" * 52)
print(f"E[Y(0)|D=1] : {neg_E_Y0_D1:+.4f}")
print(f"E[Y(0)|D=0] : {neg_E_Y0_D0:+.4f}")
print(f"selection bias : {neg_bias:+.4f} <- now NEGATIVE")
print(f"true ATT : {neg_att:+.4f}")
print(f"naive difference : {neg_naive:+.4f} <- now UNDERstates 0.07")
print(f"naive == ATT + bias ? : {abs(neg_naive - (neg_att + neg_bias)) < 1e-9}")
# The decomposition still holds exactly; only the SIGN of the bias flipped.
assert abs(neg_naive - (neg_att + neg_bias)) < 1e-9, "decomposition holds for negative selection too"
assert neg_bias < 0, "negative selection should make the bias negative"
assert neg_naive < TAU_TRUE, "negative selection should make the naive number understate the truth"
print("\nCHECK PASSED: flipping the selection rule flips the bias sign; naive now UNDERstates the effect.")
# TRY NEXT:
# (a) Make the effect HETEROGENEOUS: set Y1 = Y0 + tau + 0.15*s in build_science_table, so savvier
# people gain MORE. Now ATT != ATE -- re-derive both from the columns and see which the naive
# number is closer to under positive selection.
# (b) SUTVA-breaker: after randomizing, add a spillover so each control's Yobs rises a little when
# its treated "neighbors" are many (e.g., add 0.05 * (#treated in a random pair)). Re-run the
# naive difference and watch a CLEAN randomized estimate go wrong anyway -- interference, no fix.