Session 10: Causal Inference, Built Step by Step

Published course presentation · repository edition

Thierry Warin, PhD

HEC Montréal

Outline · The hour ahead

Part 1 · The question, in pictures
Part 2 · A world where we know the truth
Part 3 · Matching, tested on real data
Part 4 · What matching does not fix
Part 5 · Your recipe

What session 10 is for

The mathematics

State the fundamental problem of causal inference, and write the naive difference as effect plus selection bias
Explain why randomisation solves it, and what it costs
Estimate a propensity score and match treated to control units
Report an ATT, and state the assumption it rests on

The Python

Simulate both potential outcomes and watch a coin flip recover the effect
Fit a propensity score with smf.logit and match with NearestNeighbors
Check overlap and balance (standardised differences) before reading any estimate
Score the naive, regression and matched estimates against an experimental benchmark

The international business

Estimate the effect of a programme that firms chose rather than drew
Say why the firms that joined were already different
Report how far the naive comparison would have been wrong

Your Python so far

Everything you need to type is on one sheet — how to start Python in VS Codium on your platform, how to make a file and run it, the DataFrame basics, and every name from sessions 1 to 10, each with a line showing it in use.

Python Cheatsheet · Session 10 (PDF) · in the repository at 10-…/PYTHON-CHEATSHEET.pdf

One sheet per package, functions explained with pictures: NumPy · matplotlib · statsmodels · scikit-learn · SciPy

Print it, or keep it open beside the editor. Nothing on it is examinable; all of it is assumed.

Part 1 · The question, in pictures

Did the programme work?

A government runs an export-promotion programme: advice, trade missions, a subsidy for the first foreign trade fair. Two years later, the firms that joined have grown their export sales 5.3 points faster than the firms that did not.

The tempting reading: the programme added 5.3 points of growth. Fund it again, and bigger.

The question to ask first: which firms chose to join? If the better firms joined, they would have grown faster anyway, and the 5.3 points is partly the programme and partly the firms.

Today’s question: did the policy do anything, or were the groups different to begin with? By the end of the hour you can separate the two, and say when you cannot.

Every firm has two outcomes; you see one

Six firms, each with sales growth with and without the programme. For each firm only one of the two is observed; the other, and the effect, are question marks.

For firm $i$: $Y_i(1)$ is its sales growth with the programme, $Y_i(0)$ without. The effect for that firm is $\tau_i = Y_i(1) - Y_i(0)$. A firm either joins or does not, so one of the two is always missing, and so is $\tau_i$.

The fundamental problem, and what we estimate instead

No individual effect is ever observed (Holland, 1986). Causal inference is a missing-data problem: we fill in the missing column with other firms who look like this one, and the whole method is in choosing those other firms well.

What we can hope to estimate is an average:

Name	Formula	In words	When you want it
ATE	$\mathrm{E}[Y(1) - Y(0)]$	the effect if every firm joined	rolling the programme out to all
ATT	$\mathrm{E}[Y(1) - Y(0) \mid T = 1]$	the effect on the firms that joined	was it worth it for those who used it?

$T$ is the treatment: 1 if the firm joined, 0 if not. Matching, today’s method, estimates the ATT.

Why the naive comparison fails, in one line

\[ \underbrace{\mathrm{E}[Y \mid T{=}1] - \mathrm{E}[Y \mid T{=}0]}_{\text{what you compute}} = \underbrace{\mathrm{E}[Y(1) - Y(0) \mid T{=}1]}_{\text{ATT: what you want}} + \underbrace{\mathrm{E}[Y(0) \mid T{=}1] - \mathrm{E}[Y(0) \mid T{=}0]}_{\text{selection bias}} \]

A bar of 5.33 split into the effect, 2.00, and selection bias, 3.33.

Selection bias is how differently the two groups would have done without the programme. Here (from the simulation in Part 2) the firms that joined would have grown 3.33 points faster anyway.

The picture behind the bias: a back door

Left: productivity causes both joining the programme and sales growth. Right: a coin decides who joins, and productivity no longer reaches the programme.

Left: productive firms join more and grow more, so the comparison picks up both paths. Right: a coin cannot know which firms are productive, so the back door is shut, and the difference in means measures the programme alone. That is the whole case for randomisation.

Part 2 · A world where we know the truth

Step 1 · Simulate 1,000 firms, with both outcomes

In a simulation we can do what real data never allows: write down both columns.

# Step 1 · 1,000 firms; the programme adds exactly +2 points to every one
import numpy as np
import pandas as pd

rng = np.random.default_rng(10)
n = 1000
productivity = rng.normal(0, 1, n)                        # nobody in the ministry sees this
sales_without = 10 + 3 * productivity + rng.normal(0, 2, n)
sales_with = sales_without + 2                            # the true effect: +2

firms = pd.DataFrame({"productivity": productivity,
                      "sales_without": sales_without, "sales_with": sales_with})
print(firms.head().round(2))
print("true effect, every firm:", (firms["sales_with"] - firms["sales_without"]).mean())

   productivity  sales_without  sales_with
0         -1.10           5.02        7.02
1         -0.73           9.40       11.40
2         -0.78           7.19        9.19
3          0.27          11.02       13.02
4         -0.25          11.34       13.34
true effect, every firm: 2.0

Reading it: each row is a firm with its two potential outcomes. More productive firms grow faster with or without the programme (the 3 * productivity). Every effect is exactly +2: that is the number each method must recover.

Step 2 · Let a coin decide who joins

# Step 2 · Randomise: a coin flip decides; we then see only one outcome per firm
import numpy as np
import pandas as pd

rng = np.random.default_rng(10)
n = 1000
productivity = rng.normal(0, 1, n)
sales_without = 10 + 3 * productivity + rng.normal(0, 2, n)
sales_with = sales_without + 2

coin = rng.integers(0, 2, n)                                    # 1 = joins, 0 = does not
observed = np.where(coin == 1, sales_with, sales_without)       # the column you would get

print("difference in means:", round(observed[coin == 1].mean() - observed[coin == 0].mean(), 2))
print("mean productivity, joined vs not:",
      round(productivity[coin == 1].mean(), 2), round(productivity[coin == 0].mean(), 2))

difference in means: 2.15
mean productivity, joined vs not: -0.04 -0.1

Reading it: the plain difference in means gives 2.15, close to the true +2. The two groups have nearly the same productivity: the coin balanced it without anyone measuring it. The 0.15 is luck, not bias, which the next step shows.

Step 3 · Flip the coin 1,000 times

# Step 3 · Rerun the experiment: is the coin right on average?
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(10)
n = 1000
productivity = rng.normal(0, 1, n)
sales_without = 10 + 3 * productivity + rng.normal(0, 2, n)
sales_with = sales_without + 2

estimates = []
for draw in range(1000):
    coin = rng.integers(0, 2, n)
    observed = np.where(coin == 1, sales_with, sales_without)
    estimates.append(observed[coin == 1].mean() - observed[coin == 0].mean())

print("mean of the estimates:", round(np.mean(estimates), 2))
print("their spread (sd):    ", round(np.std(estimates), 2))
plt.hist(estimates, bins=30)
plt.axvline(2, color="red")
plt.xlabel("difference in means, one coin each")
plt.show()

mean of the estimates: 2.0
their spread (sd):     0.23

Reading it: centred on 2 (the red line): randomisation is unbiased. Any one experiment misses by about 0.23 (the spread). That spread is what a standard error measures: randomisation turns bias into noise, and noise is what confidence intervals are for.

Step 4 · Now let the firms choose

# Step 4 · Self-selection: productive firms are more likely to join
import numpy as np

rng = np.random.default_rng(10)
n = 1000
productivity = rng.normal(0, 1, n)
sales_without = 10 + 3 * productivity + rng.normal(0, 2, n)
sales_with = sales_without + 2

applies = (productivity + rng.normal(0, 1, n) > 0.5).astype(int)
observed = np.where(applies == 1, sales_with, sales_without)

naive = observed[applies == 1].mean() - observed[applies == 0].mean()
bias = sales_without[applies == 1].mean() - sales_without[applies == 0].mean()
print(f"share that joined:     {applies.mean():.2f}")
print(f"naive difference:      {naive:.2f}")
print(f"selection bias:        {bias:.2f}   (needs the column you never see)")
print(f"naive minus bias:      {naive - bias:.2f}")

share that joined:     0.32
naive difference:      5.33
selection bias:        3.33   (needs the column you never see)
naive minus bias:      2.00

Reading it: same firms, same true effect of +2, but the naive comparison now says 5.33. The gap, 3.33, is exactly the selection bias of the formula: firms that joined would have grown 3.33 points faster without the programme. In real data that line cannot be computed.

Step 5 · You can see that the groups differ

# Step 5 · Compare the groups BEFORE treatment
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

rng = np.random.default_rng(10)
n = 1000
productivity = rng.normal(0, 1, n)
sales_without = 10 + 3 * productivity + rng.normal(0, 2, n)
applies = (productivity + rng.normal(0, 1, n) > 0.5).astype(int)

firms = pd.DataFrame({"productivity": productivity, "applies": applies})
print(firms.groupby("applies")["productivity"].mean().round(2))

firms.boxplot(column="productivity", by="applies")
plt.suptitle("")
plt.title("Productivity, before the programme")
plt.show()

applies
0   -0.42
1    0.67
Name: productivity, dtype: float64

Reading it: the firms that joined were about one standard deviation more productive before the programme existed. This is the balance check: compare the groups on what was fixed before treatment. In real data you cannot see the bias, but you can see this.

What randomisation buys, and what it costs

Buys:

balance on everything, measured or not (Step 2)
an unbiased estimate from a difference in means (Step 3)
a standard error that means what it says

Costs:

often impossible: you cannot randomise which countries sign a trade agreement
sometimes unethical, often expensive
tells you about the people in the trial: external validity is argued, not estimated

So the rest of the session asks: without a coin, can we build a control group that looks like the treated group, and does it give the right answer? We test it where the right answer is known.

Part 3 · Matching, tested on real data

LaLonde’s test

Left: the experiment, 185 trained and 260 untrained by lottery. Right: the same 185 trained compared with 15,992 people from a survey.

The National Supported Work programme (US, 1976–77) gave disadvantaged workers a job-training placement by lottery. LaLonde (1986) kept the trainees, threw away the lottery’s controls, and asked whether the methods economists use on survey data would recover the experiment’s answer.

Column	Meaning	Column	Meaning
`treat`	1 = trained	`married`, `nodegree`	1 = married; 1 = no high-school diploma
`age`, `education`	years	`re74`, `re75`	earnings in 1974 and 1975, before training, in dollars
`black`, `hispanic`	1 = yes	`re78`	earnings in 1978, after: the outcome, in dollars

Step 1 · The experiment gives the answer

# Step 1 · The lottery's estimate: our benchmark
import pandas as pd
import statsmodels.formula.api as smf

url = "https://users.nber.org/~rdehejia/data/"
nsw = pd.read_stata(url + "nsw_dw.dta")            # Dehejia & Wahba's version of the experiment

print(nsw["treat"].value_counts())
experiment = smf.ols("re78 ~ treat", data=nsw).fit()
print(experiment.summary().tables[1])

treat
0.0    260
1.0    185
Name: count, dtype: int64
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept   4554.8011    408.046     11.162      0.000    3752.855    5356.747
treat       1794.3424    632.853      2.835      0.005     550.574    3038.110
==============================================================================

Reading it: with a lottery, the regression on treat alone is the difference in means. Training raised 1978 earnings by $1,794 (standard error $633, p = 0.005). This is the number every method on the next slides is trying to hit.

Step 2 · Swap in survey controls: the naive estimate

# Step 2 · The 185 trainees against 15,992 people from the Current Population Survey
import pandas as pd

url = "https://users.nber.org/~rdehejia/data/"
nsw = pd.read_stata(url + "nsw_dw.dta")
cps = pd.read_stata(url + "cps_controls.dta")
data = pd.concat([nsw[nsw["treat"] == 1], cps], ignore_index=True)

print(data["treat"].value_counts())
print(data.groupby("treat")["re78"].mean().round(0))
naive = data.loc[data["treat"] == 1, "re78"].mean() - data.loc[data["treat"] == 0, "re78"].mean()
print(f"naive estimate: {naive:,.0f} dollars")

treat
0.0    15992
1.0      185
Name: count, dtype: int64
treat
0.0    14847.0
1.0     6349.0
Name: re78, dtype: float32
naive estimate: -8,498 dollars

Reading it: the naive comparison says training cut earnings by $8,498: the wrong sign and five times too large. Nothing is wrong with the arithmetic. The survey’s average worker earned far more than these trainees ever would have; that is selection bias, $10,000 of it.

Step 3 · Check balance: how different are the groups?

# Step 3 · Means and standardised differences, before matching
import numpy as np
import pandas as pd

url = "https://users.nber.org/~rdehejia/data/"
nsw = pd.read_stata(url + "nsw_dw.dta")
cps = pd.read_stata(url + "cps_controls.dta")
data = pd.concat([nsw[nsw["treat"] == 1], cps], ignore_index=True)
X = ["age", "education", "black", "hispanic", "married", "nodegree", "re74", "re75"]

treated, control = data[data["treat"] == 1], data[data["treat"] == 0]
sd = np.sqrt((treated[X].var() + control[X].var()) / 2)
balance = pd.DataFrame({"treated": treated[X].mean(), "control": control[X].mean()})
balance["std. diff."] = (treated[X].mean() - control[X].mean()) / sd
print(balance.astype(float).round(2))

           treated   control  std. diff.
age          25.82     33.23       -0.80
education    10.35     12.03       -0.68
black         0.84      0.07        2.43
hispanic      0.06      0.07       -0.05
married       0.19      0.71       -1.23
nodegree      0.71      0.30        0.90
re74       2095.57  14016.80       -1.57
re75       1532.06  13650.80       -1.75

Reading it: the standardised difference is the gap in means in standard deviations, so years and dollars compare. Below 0.1 is balanced; most here are out by 0.7 to 2.4. Trainees are younger, mostly Black, unmarried, without a diploma, and earned under $2,100 a year before, against $14,000.

The idea of matching: one number to compare on

Treated people on top, controls below, on a propensity-score axis. Each treated person is joined by a line to the control with the closest score; controls far from every treated person are grey and unused.

Propensity score $e(x) = P(T = 1 \mid X = x)$: the probability of being trained, given the eight characteristics. Fit it with a logit (session 7).

Rosenbaum and Rubin (1983): if the eight $X$ explain who got trained, comparing people with the same score is as good as comparing people with the same eight $X$. Eight dimensions become one.

Matching: pair each trainee with the survey person whose score is closest: the nearest neighbour of session 8, on one number. The matched controls are the “other firms” that fill in the missing column.

Step 4 · Estimate the propensity score

# Step 4 · A logit of treatment on the pre-training characteristics
import pandas as pd
import statsmodels.formula.api as smf

url = "https://users.nber.org/~rdehejia/data/"
nsw = pd.read_stata(url + "nsw_dw.dta")
cps = pd.read_stata(url + "cps_controls.dta")
data = pd.concat([nsw[nsw["treat"] == 1], cps], ignore_index=True)

score = smf.logit("treat ~ age + education + black + hispanic + married + nodegree + re74 + re75",
                  data=data).fit(disp=0)
data["pscore"] = score.predict(data)
print(data.groupby("treat")["pscore"].describe().round(3))

         count   mean    std    min   25%    50%    75%    max
treat                                                         
0.0    15992.0  0.008  0.043  0.000  0.00  0.000  0.002  0.488
1.0      185.0  0.271  0.168  0.001  0.13  0.257  0.455  0.487

Reading it: a trainee’s median score is 0.26; a survey person’s is 0.0001. The logit’s coefficients are not the result and are never reported: the score is a tool for building the comparison, judged only by the balance it produces.

Step 5 · Check overlap: is there anyone to compare with?

# Step 5 · The two distributions of the score
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf

url = "https://users.nber.org/~rdehejia/data/"
nsw = pd.read_stata(url + "nsw_dw.dta")
cps = pd.read_stata(url + "cps_controls.dta")
data = pd.concat([nsw[nsw["treat"] == 1], cps], ignore_index=True)
score = smf.logit("treat ~ age + education + black + hispanic + married + nodegree + re74 + re75",
                  data=data).fit(disp=0)
data["pscore"] = score.predict(data)

bins = np.linspace(0, 0.5, 26)
plt.hist(data.loc[data["treat"] == 0, "pscore"], bins=bins, alpha=0.6, label="survey (control)")
plt.hist(data.loc[data["treat"] == 1, "pscore"], bins=bins, alpha=0.6, label="trained")
plt.yscale("log")
plt.xlabel("propensity score")
plt.legend()
plt.show()

print("controls with a score above 0.05:", (data.loc[data["treat"] == 0, "pscore"] > 0.05).sum())

controls with a score above 0.05: 534

Reading it: the vertical axis is logarithmic. Almost all 15,992 survey people sit in the first bar, at a score near 0: nobody like them was trained. But every bar where trainees are also has survey people, a few hundred in all: there is overlap, so a comparison is possible.

Step 6 · Match each trainee to the nearest control

# Step 6 · 1:1 nearest-neighbour matching on the score, with replacement
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.neighbors import NearestNeighbors

url = "https://users.nber.org/~rdehejia/data/"
nsw = pd.read_stata(url + "nsw_dw.dta")
cps = pd.read_stata(url + "cps_controls.dta")
data = pd.concat([nsw[nsw["treat"] == 1], cps], ignore_index=True)
score = smf.logit("treat ~ age + education + black + hispanic + married + nodegree + re74 + re75",
                  data=data).fit(disp=0)
data["pscore"] = score.predict(data)

treated, control = data[data["treat"] == 1], data[data["treat"] == 0]
nn = NearestNeighbors(n_neighbors=1).fit(control[["pscore"]])
distance, index = nn.kneighbors(treated[["pscore"]])
matched = control.iloc[index.ravel()]

print(f"{len(treated)} trainees matched to {matched.index.nunique()} different controls")
print(f"largest gap in score between a pair: {distance.max():.4f}")

185 trainees matched to 124 different controls
largest gap in score between a pair: 0.0038

Reading it: the same NearestNeighbors as session 8, now with one column. 185 trainees use 124 distinct controls: a good control can serve twice (with replacement). No pair is more than 0.004 apart in score, so every match is close. 15,868 survey people are never used.

Step 7 · Did matching balance the groups?

# Step 7 · Standardised differences, before and after matching
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.neighbors import NearestNeighbors

url = "https://users.nber.org/~rdehejia/data/"
nsw = pd.read_stata(url + "nsw_dw.dta")
cps = pd.read_stata(url + "cps_controls.dta")
data = pd.concat([nsw[nsw["treat"] == 1], cps], ignore_index=True)
X = ["age", "education", "black", "hispanic", "married", "nodegree", "re74", "re75"]
score = smf.logit("treat ~ age + education + black + hispanic + married + nodegree + re74 + re75",
                  data=data).fit(disp=0)
data["pscore"] = score.predict(data)

treated, control = data[data["treat"] == 1], data[data["treat"] == 0]
nn = NearestNeighbors(n_neighbors=1).fit(control[["pscore"]])
distance, index = nn.kneighbors(treated[["pscore"]])
matched = control.iloc[index.ravel()]

sd = np.sqrt((treated[X].var() + control[X].var()) / 2)
balance = pd.DataFrame({"before": (treated[X].mean() - control[X].mean()) / sd,
                        "after": (treated[X].mean() - matched[X].mean()) / sd})
print(balance.round(2))

           before  after
age         -0.80   0.31
education   -0.68   0.00
black        2.43   0.05
hispanic    -0.05   0.07
married     -1.23   0.09
nodegree     0.90   0.00
re74        -1.57   0.01
re75        -1.75   0.03

Reading it: a huge improvement: earnings, race and diploma are now balanced. But age is still off by 0.31: trainees are older than their matches. The score is not good enough yet, and this is decided before anyone looks at earnings in 1978.

Step 8 · Improve the score, then check again

# Step 8 · Add age squared and "earned nothing" indicators to the score
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
from sklearn.neighbors import NearestNeighbors

url = "https://users.nber.org/~rdehejia/data/"
nsw = pd.read_stata(url + "nsw_dw.dta")
cps = pd.read_stata(url + "cps_controls.dta")
data = pd.concat([nsw[nsw["treat"] == 1], cps], ignore_index=True)
X = ["age", "education", "black", "hispanic", "married", "nodegree", "re74", "re75"]
score = smf.logit("treat ~ age + education + black + hispanic + married + nodegree"
                  " + re74 + re75 + I(age**2) + I(re74 == 0) + I(re75 == 0)",
                  data=data).fit(disp=0)
data["pscore"] = score.predict(data)

treated, control = data[data["treat"] == 1], data[data["treat"] == 0]
nn = NearestNeighbors(n_neighbors=1).fit(control[["pscore"]])
distance, index = nn.kneighbors(treated[["pscore"]])
matched = control.iloc[index.ravel()]

sd = np.sqrt((treated[X].var() + control[X].var()) / 2)
before = (treated[X].mean() - control[X].mean()) / sd
after = (treated[X].mean() - matched[X].mean()) / sd
print(after.round(2))

plt.scatter(before, X, label="before")
plt.scatter(after, X, label="after")
plt.axvline(-0.1, color="grey", linestyle="--")
plt.axvline(0.1, color="grey", linestyle="--")
plt.xlabel("standardised difference")
plt.legend()
plt.show()

age         -0.09
education   -0.08
black       -0.03
hispanic     0.07
married      0.06
nodegree     0.08
re74         0.01
re75         0.01
dtype: float32

Reading it: I(age**2) lets the probability of training rise and then fall with age; I(re74 == 0) marks people who earned nothing. Now every difference is inside ±0.1 (the dashed lines). This picture, a love plot, is the evidence that the comparison is fair.

Step 9 · Only now, look at the outcome

# Step 9 · The ATT: trainees' 1978 earnings minus their matches'
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.neighbors import NearestNeighbors

url = "https://users.nber.org/~rdehejia/data/"
nsw = pd.read_stata(url + "nsw_dw.dta")
cps = pd.read_stata(url + "cps_controls.dta")
data = pd.concat([nsw[nsw["treat"] == 1], cps], ignore_index=True)
score = smf.logit("treat ~ age + education + black + hispanic + married + nodegree"
                  " + re74 + re75 + I(age**2) + I(re74 == 0) + I(re75 == 0)",
                  data=data).fit(disp=0)
data["pscore"] = score.predict(data)

treated, control = data[data["treat"] == 1], data[data["treat"] == 0]
nn = NearestNeighbors(n_neighbors=1).fit(control[["pscore"]])
distance, index = nn.kneighbors(treated[["pscore"]])
matched = control.iloc[index.ravel()]

att = treated["re78"].mean() - matched["re78"].mean()
print(f"trained, 1978:          {treated['re78'].mean():8,.0f}")
print(f"their matches, 1978:    {matched['re78'].mean():8,.0f}")
print(f"ATT (matching):         {att:+8,.0f}   experiment: +1,794")

trained, 1978:             6,349
their matches, 1978:       4,348
ATT (matching):           +2,001   experiment: +1,794

Reading it: the matched controls earned what the trainees would have earned without training: the missing column, filled in. The difference, +$2,001, is within a third of a standard error of the experiment’s +$1,794. From −$8,498 to the right answer, using only survey data.

Step 9 · Every estimate, against the truth

Five estimates on one axis: naive −8,498, regression +699, first matching +1,831, balanced matching +2,001, and the experiment +1,794 with its 95% interval.

Both matching estimates land inside the experiment’s 95% interval. Balance moved the estimate by $170; it is the balance check, not the number, that makes it defensible.

Regression with the same eight controls gives +$699. A linear model compares trainees with survey people earning $14,000, by extrapolating along a straight line. Matching only compares like with like.

Part 4 · What matching does not fix

Matching sees only what you measured

# Back to the firms: suppose the ministry sees only firm size, a noisy sign of productivity
import numpy as np
import pandas as pd
from sklearn.neighbors import NearestNeighbors

rng = np.random.default_rng(10)
n = 1000
productivity = rng.normal(0, 1, n)
sales_without = 10 + 3 * productivity + rng.normal(0, 2, n)
sales_with = sales_without + 2
applies = (productivity + rng.normal(0, 1, n) > 0.5).astype(int)
size = productivity + rng.normal(0, 1, n)                     # what the data set contains

firms = pd.DataFrame({"productivity": productivity, "size": size, "applies": applies,
                      "sales": np.where(applies == 1, sales_with, sales_without)})
treated, control = firms[firms["applies"] == 1], firms[firms["applies"] == 0]
for matched_on in ["size", "productivity"]:
    nn = NearestNeighbors(n_neighbors=1).fit(control[[matched_on]])
    distance, index = nn.kneighbors(treated[[matched_on]])
    att = treated["sales"].mean() - control.iloc[index.ravel()]["sales"].mean()
    print(f"matched on {matched_on:13} ATT = {att:.2f}   (truth: 2)")

matched on size          ATT = 3.94   (truth: 2)
matched on productivity  ATT = 2.21   (truth: 2)

Reading it: matching on the variable that drove selection recovers about 2 (2.21). Matching on a noisy proxy for it leaves 3.94: more than half of the bias survives, and size would look perfectly balanced. Balance can only be checked on what you measured. The assumption that nothing else matters, conditional independence, can never be tested.

Not every variable belongs in the score

Left: firm size before the programme causes both joining and growth; control for it. Right: new export contracts are caused by the programme and cause growth; do not control for them.

The rule: match and control only on variables fixed before treatment. A variable the programme itself changes (new contracts, the trainee’s 1978 job) is part of the effect: controlling for it removes some of what you are trying to measure. In LaLonde’s data, re74 and re75 qualify; re78 never could.

What matching does not fix

Unmeasured confounders. Matching is regression’s assumption in a different form: everything that drives both treatment and outcome is in $X$. If it is not, both are biased, in the same direction.
No overlap. Treated units with no comparable control cannot be matched honestly. Drop them, and say that the estimate now applies only to the rest.
A lucky data set. Smith and Todd (2005) showed the LaLonde result moves with the sample and the score. One success is evidence the method can work, not that it did in your project.
Standard errors. The OLS standard error on a matched sample ignores that the score was estimated and the controls reused; bootstrap it, or report it as approximate.
The order of work. Choose the score by balance before looking at the outcome. Trying scores until the ATT looks good is p-hacking with extra steps.

Part 5 · Your recipe

Seven steps, and the one rule

Seven steps from the question to the ATT, with an arrow from the balance check back to the score.

Step	What you do	In Python
1–2	name the estimand; argue which pre-treatment variables drive treatment	on paper
3	the naive difference, as the size of the problem	`.groupby("treat")["y"].mean()`
4	the propensity score	`smf.logit(...).fit()`, `.predict()`
5	match, and check overlap	`NearestNeighbors(n_neighbors=1)`, `.kneighbors()`, `plt.hist`
6	balance: every standardised difference inside ±0.1	`.mean()`, `.var()`, `plt.scatter`
7	the ATT, last	`treated["y"].mean() - matched["y"].mean()`

The whole recipe, in one cell to copy

# Propensity-score matching, start to finish: change the data, X and the formula for your project
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.neighbors import NearestNeighbors

url = "https://users.nber.org/~rdehejia/data/"                    # 1 · data: treated + candidate controls
data = pd.concat([pd.read_stata(url + "nsw_dw.dta").query("treat == 1"),
                  pd.read_stata(url + "cps_controls.dta")], ignore_index=True)
X = ["age", "education", "black", "hispanic", "married", "nodegree", "re74", "re75"]
y = "re78"

naive = data.loc[data["treat"] == 1, y].mean() - data.loc[data["treat"] == 0, y].mean()

score = smf.logit("treat ~ " + " + ".join(X) +                    # 2 · the score
                  " + I(age**2) + I(re74 == 0) + I(re75 == 0)", data=data).fit(disp=0)
data["pscore"] = score.predict(data)

treated, control = data[data["treat"] == 1], data[data["treat"] == 0]
nn = NearestNeighbors(n_neighbors=1).fit(control[["pscore"]])   # 3 · match
distance, index = nn.kneighbors(treated[["pscore"]])
matched = control.iloc[index.ravel()]

sd = np.sqrt((treated[X].var() + control[X].var()) / 2)          # 4 · balance, before the outcome
print(((treated[X].mean() - matched[X].mean()) / sd).round(2))

att = treated[y].mean() - matched[y].mean()                       # 5 · the ATT, last
print(f"naive {naive:,.0f}   ATT {att:,.0f}")

age         -0.09
education   -0.08
black       -0.03
hispanic     0.07
married      0.06
nodegree     0.08
re74         0.01
re75         0.01
dtype: float32
naive -8,498   ATT 2,001

What to report

The estimand, named: ATT or ATE.
The identifying assumption, in words: which variables you claim explain who was treated, and why nothing else does. An argument, not a table.
The naive estimate beside the matched one: the distance between them is the size of the problem.
Overlap: the score’s range in both groups, and what you did with units outside it.
Balance before and after, every covariate, as standardised differences (the love plot).
That every control is pre-treatment, and the score was fixed before the outcome was looked at.

A balanced match makes the comparison fair on what you measured. It says nothing about what you did not.

Your turn

Apply the recipe to LaLonde’s other survey: replace cps_controls.dta by psid_controls.dta (2,490 people from the Panel Study of Income Dynamics) in the cell above.

What is the naive estimate? Why is it even further from $1,794 than with the CPS?
With the first score (no squares), is every standardised difference inside ±0.1? Which is not?
Improve the score for balance without printing the ATT. If some difference stays above 0.1, say which, and what it means. Then print the ATT once.
How far is your ATT from the experiment? How many distinct controls did you use?

Work in pairs; one shares their screen at the end. The point is the order of work, not the number.

Before the next session

Required: the next session’s article, its Dataverse package, and the self-check. That is the pre-session deck.

Optional: the QMIB Lab for this session — a knowledge check on the hour you have just sat through: warin.ca/qmib-labs

Never a prerequisite, never marked for correctness. A completed report counts toward participation because attempting it is the engagement being measured.

Name	Formula	In words	When you want it
ATE	\(\mathrm{E}[Y(1) - Y(0)]\)	the effect if every firm joined	rolling the programme out to all
ATT	\(\mathrm{E}[Y(1) - Y(0) \mid T = 1]\)	the effect on the firms that joined	was it worth it for those who used it?