Covariates in Difference-in-Differences: The LaLonde Test in Python

Abstract

Forty years ago Robert LaLonde delivered one of the most influential blows of the credibility revolution: he took a job-training experiment with a known treatment effect, threw away the randomized control group, and showed that a gauntlet of respectable econometric estimators could not recover the answer. This tutorial reproduces — in Python, and in the spirit of a recent essay by Scott Cunningham — a modern version of that test aimed squarely at covariates in difference-in-differences (DiD). Working with the Dehejia-Wahba subsample (185 National Supported Work trainees plus 15,992 Current Population Survey controls, a two-period panel of 1975 and 1978 earnings), the estimand is the average treatment effect on the treated (ATT), and the ground truth is the experimental benchmark of \$1,794. We estimate the ATT eight ways using pyfixest for the regressions and hand-coded inverse-propensity-weighting (Abadie 2005) and doubly-robust (Sant’Anna-Zhao 2020) estimators, cross-checked against the diff-diff package. The result is stark and dollar-accurate. Three specifications that keep covariates out of the counterfactual trend — no covariates, additive covariates, and covariate-by-treatment interactions — all return the naive \$3,621, roughly twice the truth. The instant covariates are allowed to bend the control group’s trend, via X × post (\$1,711) or first-difference saturation (\$1,770, numerically identical to the Heckman-Ichimura-Todd estimator), the estimate snaps to the benchmark; the propensity-based estimators land nearby at \$1,861 and \$1,993. The lesson is that covariates in DiD are not a robustness knob to be twisted for reassurance — they perform a specific job, satisfying conditional parallel trends and relaxing constant treatment effects, and the most common applied specification of all, two-way fixed effects with additive controls, does not do that job.

Overview

Imagine you already know the right answer. A job-training program raised the earnings of disadvantaged workers by about \$1,794 — we know this because the program was evaluated with a randomized controlled trial, and randomization makes the treated and control groups exchangeable. Now throw the experimental control group away and replace it with a large survey of ordinary Americans. Can a difference-in-differences estimator, armed with the usual baseline covariates, still find the \$1,794? And does it matter how you feed those covariates to the model?

This is the LaLonde test, named for LaLonde (1986), who used exactly this setup to show that unbiased-in-principle estimators can be badly biased in application. The exercise here follows a recent post by Scott Cunningham on his Causal Inference: The Mixtape Substack, translating his Stata and R analysis into Python and adding a package cross-check. If you are new to difference-in-differences, start with the companion tutorial, Introduction to Difference-in-Differences in Python; this post is the advanced sequel that asks a sharper question about covariates.

The punchline, which we will earn spec by spec, is a distinction that most applied work blurs. A covariate can enter a DiD in three fundamentally different places: in the level of the outcome (additive controls), in the treatment effect (covariate-by-treatment interactions), or in the trend (covariate-by-time interactions). Only the last one addresses the reason the naive estimate is wrong. Because our control group is a representative slice of America and our treated group is a set of disadvantaged trainees, the two groups are wildly imbalanced; and if groups with different characteristics are on different earnings trajectories, then that imbalance mechanically breaks the parallel-trends assumption. Fixing it requires modeling those trajectories — putting covariates in the trend — not simply adding them to the regression.

Learning objectives:

• Understand why a covariate rescues a DiD estimate only when it enters the control group’s counterfactual trend, not the level or the treatment effect.
• Implement eight covariate specifications in pyfixest — from the naive two-way fixed effects model to a fully saturated first-difference regression — and recover the ATT by g-computation.
• Hand-code the Abadie (2005) inverse-propensity-weighting and Sant’Anna-Zhao (2020) doubly-robust DiD estimators, and cross-check every number against the diff-diff package.
• Judge all eight estimates against the \$1,794 experimental benchmark and diagnose why the ubiquitous additive-controls specification stays stuck at \$3,621.

Key concepts at a glance

The post leans on a small vocabulary repeatedly. Each concept below has three parts. The definition is always visible. The example and analogy sit behind clickable cards: open them when you need them, leave them collapsed for a quick scan. If a later section mentions “conditional parallel trends” or “outcome regression” and the term feels slippery, this is the section to re-read.

1. The LaLonde test. Benchmark an observational estimator against a known experimental ATT by replacing the RCT control group with survey controls and checking whether the method still recovers the truth.

2. Estimand: ATT vs ATE. The ATT is the effect on the treated units; the ATE is the effect on everyone. Under randomization they coincide. Replacing the experimental control with a survey control changes who the comparison represents but leaves the treated group — and therefore the ATT — intact.

3. Conditional parallel trends. DiD identifies the ATT if treated and control outcomes would have moved in parallel absent treatment. When groups are imbalanced, that parallelism is only credible given covariates $X$ — parallel trends conditional on $X$.

4. Covariate imbalance. Systematic differences in $X$ between treated and control. Harmless on its own — but combined with covariate-specific trends, it becomes the bias term that breaks parallel trends.

5. Trend vs level vs effect. The crux. A covariate in the level (additive) or the treatment effect ($X \times D$) leaves the counterfactual trend untouched and is inert. A covariate in the trend ($X \times \text{post}$) bends the control’s counterfactual and corrects the bias.

6. Outcome regression / HIT (1997). First-difference the outcome, fit a flexible model on the controls, and impute each treated unit’s counterfactual change. Because the outcome is now a change, a covariate’s coefficient is its effect on the trend — for free.

7. Doubly robust DiD (Sant’Anna-Zhao 2020). Combine an outcome-regression model of the trend with an inverse-propensity reweighting of the controls. Consistent if either model is correctly specified — two shots at the truth.

The decision that organizes the entire post is a single question — where does the covariate enter? — and its answer sorts every estimator into “inert” or “corrected”:

graph TD
    Q{"Where does covariate X<br/>enter the DiD?"}
    Q -->|"Nowhere (baseline)"| N["Spec 0: Naive TWFE"]
    Q -->|"In the LEVEL<br/>(additive)"| L["Spec A: Additive X"]
    Q -->|"In the EFFECT<br/>(X × treatment)"| E["Spec BT: X × treatment"]
    Q -->|"In the TREND<br/>(X × post)"| T["Spec B: X × post"]
    Q -->|"Trend + effect<br/>(saturated FD)"| S["Spec C = HIT 1997"]
    Q -->|"Propensity<br/>reweighting"| P["IPW / Doubly robust"]
    N --> INERT["Counterfactual trend untouched<br/>ATT stays ~3,621 (inert)"]
    L --> INERT
    E --> INERT
    T --> FIX["Bends the control's<br/>counterfactual trend<br/>ATT snaps to ~1,794 (corrected)"]
    S --> FIX
    P --> FIX
    FIX -.->|recovers| BENCH["RCT benchmark<br/>ATT = 1,794"]
    style Q fill:#d97757,stroke:#141413,color:#fff
    style INERT fill:#9aa0a6,stroke:#141413,color:#fff
    style FIX fill:#6a9bcc,stroke:#141413,color:#fff
    style BENCH fill:#141413,stroke:#141413,color:#fff
    style N fill:#9aa0a6,stroke:#141413,color:#fff
    style L fill:#9aa0a6,stroke:#141413,color:#fff
    style E fill:#9aa0a6,stroke:#141413,color:#fff
    style T fill:#6a9bcc,stroke:#141413,color:#fff
    style S fill:#6a9bcc,stroke:#141413,color:#fff
    style P fill:#6a9bcc,stroke:#141413,color:#fff

Setup and imports

We lean on three packages. pyfixest gives us fast, formula-based OLS with robust standard errors — the workhorse for the six regression specifications. causaldata ships the canonical LaLonde data, so nothing needs to be downloaded by hand. And diff-diff provides a scikit-learn-style DiD toolkit that we use to independently validate every hand-coded number. We fix the seed to 90210 to match the reference exactly, so the bootstrap standard errors are reproducible.

import numpy as np
import pandas as pd
import statsmodels.api as sm
import pyfixest as pf

RANDOM_SEED = 90210          # reproducible cluster bootstrap
N_BOOT = 199                 # bootstrap replications
BENCHMARK = 1794.0           # experimental ATT (Dehejia-Wahba)

# Scott's canonical LaLonde covariate set (age cube; u74 = 1[re74 == 0])
XVARS = ["age", "agesq", "agecube", "educ", "educsq",
         "marr", "nodegree", "black", "hisp", "re74", "u74"]

Every covariate here is a time-invariant baseline characteristic: age and its square and cube, education and its square, indicators for marriage, no high-school degree, race and ethnicity, 1974 earnings, and an unemployment flag for 1974. None of them changes between our two periods. Keep that fact in your pocket — it is the reason the additive specification will turn out to be completely inert.

Data: the LaLonde-DW non-experimental panel

The construction is where the fidelity of the whole exercise is won or lost, so it is worth being explicit. The nsw_mixtape dataset is already the Dehejia-Wahba subsample — 185 treated trainees and 260 experimental controls. For the non-experimental analysis we keep the 185 treated and discard the experimental controls, replacing them with the 15,992 CPS controls from cps_mixtape. The experimental controls are set aside for one job only: computing the benchmark. We then reshape the wide earnings columns (re75, re78) into a two-period panel with post = 1 in 1978.

from causaldata import nsw_mixtape, cps_mixtape

nsw = nsw_mixtape.load_pandas().data          # 185 treated + 260 exp. controls
cps = cps_mixtape.load_pandas().data          # 15,992 CPS controls (treat = 0)

# Non-experimental sample = NSW treated only + all CPS controls
wide = pd.concat([nsw[nsw.treat == 1], cps], ignore_index=True)
wide["ever_treated"] = wide["treat"].astype(int)
wide["agesq"]   = wide.age ** 2
wide["agecube"] = wide.age ** 3
wide["educsq"]  = wide.educ ** 2
wide["u74"]     = (wide.re74 == 0).astype(float)
wide["dy"]      = wide.re78 - wide.re75        # first difference (pre 1975, post 1978)

# Long two-period panel: post = 0 -> re75, post = 1 -> re78
pre  = wide.assign(re=wide.re75, post=0.0)
post = wide.assign(re=wide.re78, post=1.0)
panel = pd.concat([pre, post], ignore_index=True)

print(pd.crosstab(panel.ever_treated, panel.post))

post          0.0    1.0
ever_treated
0           15992  15992
1             185    185

The cell counts confirm a clean, balanced two-period panel: 185 treated and 15,992 controls observed in both 1975 and 1978. Because randomization is gone, the identifying assumption is no longer plain parallel trends but conditional parallel trends — trainees and CPS controls with the same covariates are assumed to have moved in parallel absent treatment. Whether that assumption is plausible depends entirely on how different the two groups actually are, which is the next thing to look at.

The covariate imbalance problem

Before running a single regression, we should ask how comparable the groups are. The standard diagnostic is the standardized mean difference (SMD): the gap in a covariate’s mean between treated and control, divided by the pooled standard deviation. A rule of thumb flags anything beyond 0.1 as imbalanced and beyond 0.25 as severe. We compute it for the trainees against the CPS controls, and — as a reference point — against the discarded experimental controls.

def smd(a, b, col):
    s = np.sqrt((a[col].var() + b[col].var()) / 2)
    return 0.0 if s == 0 else (a[col].mean() - b[col].mean()) / s

treated = wide[wide.ever_treated == 1]
cps     = wide[wide.ever_treated == 0]
for col in ["black", "re74", "married", "educ"]:
    print(f"{col:8s}  vs CPS = {smd(treated, cps, col):+.2f}")

The picture is dramatic. Against the CPS controls (orange), the trainees differ by a standardized +2.3 on race, −1.6 on 1974 and 1975 earnings, and −1.3 on marriage — imbalances many times the “severe” threshold. Against the randomized experimental controls (blue), every difference collapses to essentially zero. This is the entire problem in one figure: the CPS is a representative sample of America, so its members are older, better educated, more often married, and far richer than a set of disadvantaged trainees. Imbalance this extreme is exactly the condition under which conditioning on covariates stops being optional.

Raw earnings trends

Imbalance in levels is only half the story. It becomes a bias only if the imbalanced groups are also on different trends. So let us plot mean earnings for each group in 1974, 1975, and 1978.

The CPS controls (black) sit far above everyone else — about \$14,017 in 1974 against the trainees’ \$2,096 — and drift gently upward. The trainees (orange) and the experimental controls (blue) start together near the bottom and rise in tandem into 1978. A DiD that uses the CPS as-is is implicitly assuming the trainees, absent training, would have followed the CPS’s flat high-earning path. That is not credible: low-earning workers and high-earning workers are on different earnings trajectories, and the trainees look nothing like the CPS. This is precisely where covariates must do work — not by shifting levels, but by modeling those divergent trends.

A two-minute DiD refresher

If the details of difference-in-differences are hazy, the introductory tutorial covers them properly; here is just enough to fix notation. The two-by-two DiD estimator of the ATT is the difference of two differences,

$$\widehat{\text{ATT}} = (\bar Y_{\text{treated,post}} - \bar Y_{\text{treated,pre}}) - (\bar Y_{\text{control,post}} - \bar Y_{\text{control,pre}}),$$

and, as Cunningham emphasizes, four regressions compute this identical number: a saturated regression with a treatment dummy, a post dummy and their interaction; two-way fixed effects with a post-by-treatment interaction; a first-difference regression on the treatment dummy; and an across-group regression on the post dummy. We use the saturated form because it is the one that lets us add time-invariant covariates in the various ways we want to compare. Throughout, post:ever_treated is the interaction whose coefficient is the DiD estimate of the ATT.

Spec 0 — Naive TWFE (no covariates)

Start with the estimator that ignores covariates entirely. This is the number to beat.

def att(fit):
    td = fit.tidy()
    key = [k for k in td.index if "post" in k.lower() and "ever_treated" in k.lower()][0]
    return td.loc[key, "Estimate"], td.loc[key, "Std. Error"]

s0 = pf.feols("re ~ post * ever_treated", data=panel, vcov="HC1")
print("Spec 0 (naive) ATT = %.0f (SE %.0f)" % att(s0))

Spec 0 (naive) ATT = 3621 (SE 632)

The naive DiD returns \$3,621 — positive, statistically significant, and about twice the true \$1,794. On its own it looks like a perfectly respectable result; nothing about the output warns you that it is off by a factor of two. That is the whole danger LaLonde exposed: a plausible, precise, and wrong estimate. Everything that follows is an attempt to fix it with covariates.

Spec A — Additive X (covariates in the level)

The reflexive fix — and by far the most common specification in applied panel work — is to throw the covariates into the regression additively. This is two-way fixed effects with controls.

XF = " + ".join(XVARS)
sA = pf.feols(f"re ~ post * ever_treated + {XF}", data=panel, vcov="HC1")
print("Spec A (additive X) ATT = %.0f (SE %.0f)" % att(sA))

Spec A (additive X) ATT = 3621 (SE 672)

Nothing moved. Spec A returns \$3,621, identical to the naive estimate to the dollar. This is not a coincidence or a rounding accident — it is mechanical. Because the covariates are time-invariant, the within transformation that defines two-way fixed effects sweeps them out entirely; the DiD coefficient never depended on them in the first place. Adding time-invariant controls to a TWFE DiD is, quite literally, doing nothing to the point estimate. If you have ever added baseline controls to a DiD and reported that “the estimate is robust to covariates,” this is worth sitting with.

Spec BT — X × treatment (covariates in the effect)

Perhaps the additive form was too rigid because it forces the treatment effect to be the same for everyone. So let the effect vary with covariates: interact each covariate with the treatment switch $T = \text{post} \times D$, and recover the ATT by g-computation — predict each unit’s outcome with the switch on and off, and average the difference over the treated-post cell.

d = panel.assign(T=panel.post * panel.ever_treated)
T_ints = " + ".join(f"T:{x}" for x in XVARS)
fit = pf.feols(f"re ~ post + ever_treated + T + {XF} + {T_ints}", data=d, vcov="HC1")
tau = fit.predict(newdata=d.assign(T=1.0)) - fit.predict(newdata=d.assign(T=0.0))
mask = (d.ever_treated == 1) & (d.post == 1)
print("Spec BT (X x treatment) ATT = %.0f" % tau[mask.values].mean())

Spec BT (X x treatment) ATT = 3621

Inert again — \$3,621, to the dollar. Allowing heterogeneous treatment effects saturates the model in levels, which relaxes the constant-effects assumption but does absolutely nothing about conditional parallel trends. The problem was never that the treatment effect varies; it is that the control group’s counterfactual trend is mismodeled. Interacting covariates with treatment is answering a question nobody asked. We have now tried covariates in the level and in the effect, and both left the estimate exactly where it started.

Spec B — X × post (covariates in the trend)

Now put the covariates where the problem actually lives: the trend. Interact each covariate with the post indicator, so that workers with different characteristics are allowed to be on different earnings trajectories over time.

post_ints = " + ".join(f"post:{x}" for x in XVARS)
sB = pf.feols(f"re ~ post * ever_treated + {XF} + {post_ints}", data=panel, vcov="HC1")
print("Spec B (X x post) ATT = %.0f (SE %.0f)" % att(sB))

Spec B (X x post) ATT = 1711 (SE 704)

The estimate collapses from \$3,621 to \$1,711 — a \$1,910 move that lands within \$83 of the \$1,794 benchmark. Same covariates as Spec A; the only change is that they now multiply post instead of sitting in the level. That single change lets the model say “high-earning, well-educated workers were on a steeper path than low-earning ones,” which is what the counterfactual for the trainees actually requires. This is the pivot of the entire post: covariates rescue the estimate, but only from the trend.

Spec C — Saturated first differences = HIT (1997)

Spec B corrected the trend but still imposes a constant treatment effect. We can do both jobs at once. First-difference the outcome, then regress the change $\Delta y = y_{1978} - y_{1975}$ on the covariates, the treatment indicator, and their full interaction, and recover the ATT by g-computation over the treated.

D_ints = " + ".join(f"ever_treated:{x}" for x in XVARS)
sC = pf.feols(f"dy ~ {XF} + {D_ints} + ever_treated", data=wide, vcov="HC1")
tau = (sC.predict(newdata=wide.assign(ever_treated=1.0))
       - sC.predict(newdata=wide.assign(ever_treated=0.0)))
print("Spec C (saturated FD) ATT = %.0f" % tau[wide.ever_treated.values == 1].mean())

Spec C (saturated FD) ATT = 1770

\$1,770 — within \$24 of the truth. The subtlety worth savoring: because the outcome is now a change, a covariate’s own coefficient is no longer a level effect but its effect on the trend. So the first-difference regression does two things simultaneously — the covariate main effects bend the control’s counterfactual trend (what Spec B did), and the treatment interactions relax the constant-effect assumption (what Spec BT tried to do). One regression, both corrections. The pure-levels saturation of Spec BT never touched the trend, which is why it sat at \$3,621 while this lands at the benchmark.

HIT (1997) by hand

Here is the result that Cunningham rightly calls “not terribly intuitive.” The fully saturated first-difference regression above is numerically identical to a multi-step outcome-regression procedure — Heckman, Ichimura and Todd (1997) — that only ever fits a model on the controls. Fit $\Delta y$ on covariates using the control group alone, impute the counterfactual change for each treated unit, and average the treated units’ actual-minus-imputed gains.

import statsmodels.formula.api as smf
mH = smf.ols(f"dy ~ {XF}", data=wide[wide.ever_treated == 0]).fit()   # controls only
hit = (wide.dy - mH.predict(wide))[wide.ever_treated == 1].mean()      # impute to treated
print("HIT by hand ATT = %.0f" % hit)

HIT by hand ATT = 1770

Identical to Spec C — \$1,770. A regression that touches only the controls and then imputes to the treated gives exactly what a saturated regression on the whole sample gives. This is the Oaxaca-Blinder logic applied to a change, and it clarifies what “outcome regression” means in the DiD context: learn how ordinary workers’ earnings evolved, then ask how much more each trainee gained than a statistical twin would have. The numerical equivalence is a small piece of econometric magic, and it is reassuring that two conceptually different recipes agree to the dollar.

IPW — Abadie (2005) by hand

Regression adjustment models the outcome. The alternative is to model treatment — estimate each unit’s probability of being a trainee and reweight the controls to resemble the treated. Abadie (2005) gives the propensity-weighted DiD estimator of the ATT with weights

$$w_i = \frac{D_i - \hat p(X_i)}{1 - \hat p(X_i)} \cdot \frac{1}{\Pr(D = 1)},$$

applied to the first-differenced outcome.

Xc = sm.add_constant(wide[XVARS])
wide["phat"] = sm.Logit(wide.ever_treated, Xc).fit(disp=0).predict(Xc)
p = wide.ever_treated.mean()
w = (wide.ever_treated - wide.phat) / (1 - wide.phat) / p
print("IPW (Abadie 2005) ATT = %.0f" % (w * wide.dy).mean())

IPW (Abadie 2005) ATT = 1861

The inverse-propensity-weighted estimate is \$1,861, another near-miss on the benchmark from a completely different modeling philosophy. Instead of specifying how earnings trend with covariates, it specifies how treatment depends on them and lets the reweighting handle the rest. That two independent strategies — outcome regression and propensity weighting — both land near \$1,800 is the kind of convergence that builds confidence in the answer.

DR — Sant’Anna-Zhao (2020) by hand

Why choose between the two? The doubly-robust estimator of Sant’Anna and Zhao (2020) combines the outcome regression on the controls with the propensity reweighting, and is consistent if either model is correctly specified.

wide["dyhat"] = mH.predict(wide)            # outcome regression from the HIT step
dr_t = (wide.ever_treated * (wide.dy - wide.dyhat) / p).mean()
dr_c = ((1 - wide.ever_treated) * (wide.phat / (1 - wide.phat)) * (wide.dy - wide.dyhat) / p).mean()
print("DR (Sant'Anna-Zhao 2020) ATT = %.0f" % (dr_t - dr_c))

DR (Sant'Anna-Zhao 2020) ATT = 1993

The doubly-robust estimate is \$1,993. It sits a little further from the benchmark than the regression-adjustment estimators, but it comes with the strongest theoretical guarantee: it would still be consistent even if we had gotten the trend model or the propensity model wrong (just not both). In a real application, where you never know which model is right, that insurance is exactly the point.

The payoff: the covariate arc

Eight estimates, one picture. Ordered as an argument, the specifications tell a story that no single number could.

The three inert specifications cluster at \$3,621; the four that touch the trend, plus the two propensity estimators, cluster around the \$1,794 line. The grouping is the thesis. And the ladder view makes the “snap” impossible to miss:

The estimate is flat at \$3,621 across the first three specifications, then drops off a cliff to \$1,711 the instant covariates enter the trend, and stays near the benchmark thereafter. That cliff — between “X × treatment” and “X × post” — is the single most important feature of the whole analysis.

Cross-check with the diff-diff package

Hand-coded estimators are pedagogically transparent, but a reader is entitled to ask whether we coded them correctly. So we run the same designs through the independent diff-diff package and compare.

from diff_diff import DifferenceInDifferences, CallawaySantAnna

naive = DifferenceInDifferences(cluster="id", seed=90210).fit(
    panel, outcome="re", treatment="ever_treated", time="post", unit="id")

cs_df = panel.assign(first_treat=np.where(panel.ever_treated == 1, 1, 0))
cs = CallawaySantAnna(estimation_method="dr", seed=90210).fit(
    cs_df, outcome="re", unit="id", time="post", first_treat="first_treat", covariates=XVARS)

naive 2x2      by pyfixest $3,621  vs diff-diff $3,621
additive X     by pyfixest $3,621  vs diff-diff $3,621
doubly robust  by hand   $1,993    vs diff-diff (Callaway-Sant'Anna) $1,979

The package agrees exactly on the naive and additive designs and lands within \$14 of the by-hand doubly-robust estimate (\$1,993 vs \$1,979). The small remaining gap is expected: diff-diff’s Callaway-Sant’Anna implementation makes its own choices about propensity estimation and weight normalization, exactly the kind of default-level difference the reference flags between hand-coded and packaged estimators. For our purposes the takeaway is that the transparent code and the battle-tested package tell the same story.

Robustness: how much should we trust \$1,770 over \$1,711?

One caution deserves to be front and center, and it comes from a comment on Cunningham’s original post by Alexis Diamond. It is tempting to rank the corrected estimators — is \$1,770 “better” than \$1,711 because it is closer to the benchmark? Look again at the confidence intervals in the table: every corrected estimate spans roughly \$400 to \$3,100. The benchmark itself is estimated with a standard error of about \$671. Against that much noise, the \$59 gap between Spec B and Spec C is meaningless. What is trustworthy is not any single dollar figure but the pattern: three specifications that ignore the trend are all wrong in the same direction, and every specification that models the trend moves decisively toward the truth. Stability across sensible specifications — and, as Diamond argues, across many datasets, not one lucky one — is the signal. This is also why we bootstrap the standard errors (199 id-clustered resamples, seed 90210): the point estimates are only as informative as their uncertainty allows.

Discussion

The exercise settles a question that applied researchers wave away too often: are covariates in a difference-in-differences a robustness check, or do they do real work? If they were a robustness check, we would hope they leave the estimate unchanged — a moving estimate would be a warning sign. But here the covariates are not decoration; they are load-bearing. Under covariate imbalance combined with covariate-specific trends, they are what makes conditional parallel trends hold, and leaving them out — or putting them in the wrong place — produces an estimate that is off by a factor of two while looking perfectly precise.

The sharper lesson is that not all ways of including covariates are equal. The most common specification in the entire applied panel literature — two-way fixed effects with additive controls — is exactly the one that fails here, because time-invariant controls vanish under the within transformation and never touch the counterfactual trend. The specifications that succeed are the ones that let covariates bend the control group’s trajectory: X × post, first-difference saturation, and their outcome-regression and doubly-robust cousins. When you next see a DiD table where “we control for baseline covariates” is offered as reassurance, the right question is not whether covariates are included but where — in the level, the effect, or the trend.

Summary and next steps

1. The reconstruction is exact. Rebuilt from the causaldata package, the naive DiD returns \$3,621 and every one of the eight estimators matches Cunningham’s reported figures to the dollar — confirming both the data and the code.
2. Covariate placement, not inclusion, is what matters. Additive covariates (\$3,621) and covariate-by-treatment interactions (\$3,621) are inert; covariate-by-time interactions (\$1,711) and first-difference saturation (\$1,770) recover the \$1,794 benchmark. IPW (\$1,861) and doubly-robust (\$1,993) land nearby from a different modeling angle.
3. Trust the pattern, not the decimal. With 95% intervals spanning thousands of dollars and a benchmark that is itself noisy, the credible finding is the clean split between trend-ignoring and trend-modeling specifications — not a ranking among the corrected estimates.
4. Limitation and where to go next. This is one dataset with a small treated group and a famously thin covariate set; as Diamond notes, reliability comes from stability across many RCT-versus-observational benchmarks, which is the mission of the emerging rctvsobs.org repository. A natural next step is to extend the analysis to staggered-timing settings with the Callaway-Sant’Anna estimator covered in the introductory DiD tutorial, where doubly-robust covariate adjustment becomes the default rather than an afterthought.

Exercises

1. Swap the covariate set. Drop re74 and u74 from XVARS and re-run Spec B. How much of the correction survives when the model can no longer condition on prior earnings? What does that tell you about which covariate is doing the work?
2. Verify the HIT equivalence yourself. Confirm numerically that Spec C (saturated first differences on the full sample) and the control-only imputation both return \$1,770. Then break the equivalence by fitting the outcome regression on the whole sample instead of the controls — does the number change, and why?
3. Stress-test the propensity model. Trim units with estimated propensity scores above 0.9 and re-estimate the IPW and doubly-robust ATTs. How sensitive is each estimator to extreme weights, and which one would you trust more in a real application?

References

1. Cunningham, S. (2026). Covariates, diff in diff and LaLonde test. Scott’s Mixtape Substack.
2. LaLonde, R. J. (1986). Evaluating the Econometric Evaluations of Training Programs with Experimental Data. American Economic Review, 76(4), 604–620.
3. Dehejia, R. H. & Wahba, S. (2002). Propensity Score-Matching Methods for Nonexperimental Causal Studies. Review of Economics and Statistics, 84(1), 151–161.
4. Heckman, J. J., Ichimura, H. & Todd, P. E. (1997). Matching as an Econometric Evaluation Estimator: Evidence from Evaluating a Job Training Programme. Review of Economic Studies, 64(4), 605–654.
5. Abadie, A. (2005). Semiparametric Difference-in-Differences Estimators. Review of Economic Studies, 72(1), 1–19.
6. Sant’Anna, P. H. C. & Zhao, J. (2020). Doubly Robust Difference-in-Differences Estimators. Journal of Econometrics, 219(1), 101–122.
7. Gerber, I. (2026). diff-diff: Difference-in-Differences Causal Inference for Python. GitHub repository.
8. pyfixest: Fast High-Dimensional Fixed Effects Estimation in Python.
9. causaldata: Example Data Sets for Causal Inference Textbooks.
10. Cunningham, S. (2021). Causal Inference: The Mixtape. Yale University Press.

Acknowledgements

AI tools (Claude Code, Gemini, NotebookLM) were used to make the contents of this post more accessible to students. The analysis reproduces and builds on Scott Cunningham’s essay “Covariates, diff in diff and LaLonde test.” Nevertheless, the content in this post may still have errors. Caution is needed when applying the contents of this post to true research projects.