Causal Inference VI — Difference-in-Differences (R companion)¶
fixest, bacondecomp, and did — the 2×2, the diagnosis, and the modern estimator¶
The staggered-DiD revolution happened largely in R, and this companion uses the packages that implement it:
fixest::feols— the DiD regression (and the twoway-FE / interaction form) with clustered standard errors;bacondecomp::bacon— the Goodman-Bacon decomposition that diagnoses twoway-FE bias by splitting it into its 2×2 comparisons;did(Callaway & Sant'Anna) — group-time ATTs from clean controls, withaggtefor overall and event-study aggregations andggdidfor the plots.
All results match the from-scratch Python notebook: the Card-Krueger DiD of +2.75, the biased twoway-FE estimate on the staggered simulation, and its repair by Callaway-Sant'Anna.
1. Card-Krueger 2×2 — fixest::feols¶
The 2×2 DiD as a regression interaction: feols(y ~ nj*after) on the stacked before/after data. The nj:after coefficient is the difference-in-differences — New Jersey's employment change relative to Pennsylvania's after the 1992 minimum-wage increase — with store-clustered standard errors. It reproduces the famous +2.75 FTE: the minimum wage did not reduce employment.
suppressMessages({library(fixest); library(bacondecomp); library(did)})
ck<-read.csv("card_krueger.csv")
long<-rbind(data.frame(store=ck$sheet,y=ck$fte,nj=ck$state,after=0),
data.frame(store=ck$sheet,y=ck$fte2,nj=ck$state,after=1))
long<-long[!is.na(long$y),]
m<-feols(y~nj*after, data=long, cluster=~store)
cat(sprintf("Card-Krueger DiD (nj:after) = %+.3f (clustered SE %.3f)\n", coef(m)["nj:after"], se(m)["nj:after"]))
njb<-mean(ck$fte[ck$state==1],na.rm=TRUE); nja<-mean(ck$fte2[ck$state==1],na.rm=TRUE)
pab<-mean(ck$fte[ck$state==0],na.rm=TRUE); paa<-mean(ck$fte2[ck$state==0],na.rm=TRUE)
cat(sprintf(" NJ %.2f -> %.2f (%+.2f); PA %.2f -> %.2f (%+.2f); DiD %+.2f\n", njb,nja,nja-njb,pab,paa,paa-pab,(nja-njb)-(paa-pab)))
options(repr.plot.width=7.5, repr.plot.height=4.4)
plot(c(0,1),c(njb,nja),type="b",pch=19,col="#c53030",lwd=2.5,ylim=c(19,24),xaxt="n",xlab="",ylab="FTE employment",main="Card-Krueger: NJ vs PA, before/after the minimum-wage hike")
lines(c(0,1),c(pab,paa),type="b",pch=15,col="#2b6cb0",lwd=2.5); axis(1,at=c(0,1),labels=c("before","after"))
lines(c(0,1),c(njb,njb+(paa-pab)),lty=2,col="grey"); legend("bottomleft",c("NJ (treated)","PA (control)","NJ counterfactual"),col=c("#c53030","#2b6cb0","grey"),lwd=2,pch=c(19,15,NA),lty=c(1,1,2),bty="n")
Warning message: "package 'fixest' was built under R version 4.6.1"
Warning message: "package 'bacondecomp' was built under R version 4.6.1"
Warning message: "package 'did' was built under R version 4.6.1"
Card-Krueger DiD (nj:after) = +2.754 (clustered SE 1.306)
NJ 20.44 -> 21.03 (+0.59); PA 23.33 -> 21.17 (-2.17); DiD +2.75
2. Diagnosing twoway-FE bias — bacondecomp::bacon¶
On staggered-adoption data, feols(y ~ ... | id + t) gives the naive twoway-FE estimate — biased when effects are dynamic. bacondecomp::bacon decomposes that estimate into its underlying 2×2 DiDs and their weights, exposing the culprit: the "Later vs Earlier Treated" comparisons (which use already-treated units as controls) carry real weight and a low estimate, dragging the pooled TWFE below the truth (true ATT ≈ 3.6).
s<-read.csv("staggered_sim.csv"); s$gg<-ifelse(s$g==999,0,s$g)
trueATT<-mean(s$eff[s$treat==1])
twfe<-coef(feols(y~treat | id+t, data=s))["treat"]
bc<-bacon(y~treat, data=s, id_var="id", time_var="t")
cat(sprintf("true ATT = %.3f\n", trueATT))
cat(sprintf("naive TWFE = %.3f (biased)\n\n", twfe))
agg<-aggregate(cbind(weight)~type, bc, sum)
est<-sapply(split(bc,bc$type), function(x) weighted.mean(x$estimate,x$weight))
cat("Goodman-Bacon decomposition:\n"); for(k in names(est)) cat(sprintf(" %-28s weight %.3f avg 2x2 DiD %.3f\n", k, agg$weight[agg$type==k], est[k]))
options(repr.plot.width=8, repr.plot.height=4.4)
cols<-c("Earlier vs Later Treated"="#2b6cb0","Later vs Earlier Treated"="#c53030","Treated vs Untreated"="#2f855a")
plot(bc$weight, bc$estimate, pch=19, col=cols[bc$type], cex=1.4, xlab="Goodman-Bacon weight", ylab="2x2 DiD estimate",
main="Bacon decomposition: 'Later vs Earlier' (red) are the forbidden comparisons")
abline(h=trueATT, lty=2); legend("topright", names(cols), col=cols, pch=19, bty="n", cex=.8)
type weight avg_est 1 Earlier vs Later Treated 0.14072 2.75562 2 Later vs Earlier Treated 0.23453 0.84319 3 Treated vs Untreated 0.62476 3.34762
true ATT = 3.580
naive TWFE = 2.677 (biased)
Goodman-Bacon decomposition:
Earlier vs Later Treated weight 0.141 avg 2x2 DiD 2.756 Later vs Earlier Treated weight 0.235 avg 2x2 DiD 0.843 Treated vs Untreated weight 0.625 avg 2x2 DiD 3.348
3. The fix — did (Callaway & Sant'Anna)¶
att_gt estimates every group-time ATT using not-yet-treated units as clean controls (control_group = "notyettreated"), and aggte collapses them into an overall effect or an event study (type = "dynamic"). The overall CS estimate recovers the true ATT that TWFE missed, and ggdid shows the event study: pre-treatment coefficients near zero (parallel-trends check) and a growing post-treatment effect.
cs<-att_gt(yname="y", tname="t", idname="id", gname="gg", data=s,
control_group="notyettreated", bstrap=FALSE, cband=FALSE)
ov<-aggte(cs, type="group")
cat(sprintf("true ATT = %.3f\n", trueATT))
cat(sprintf("naive TWFE = %.3f (biased)\n", twfe))
cat(sprintf("Callaway-SantAnna (did) = %.3f <- recovers the truth\n", ov$overall.att))
dyn<-aggte(cs, type="dynamic")
options(repr.plot.width=8, repr.plot.height=4.4)
print(ggdid(dyn) + ggplot2::ggtitle("did::ggdid event study — flat pre-trends, growing dynamic effect"))
true ATT = 3.580
naive TWFE = 2.677 (biased)
Callaway-SantAnna (did) = 3.454 <- recovers the truth
4. Summary¶
R's staggered-DiD packages reproduced the full arc of the Python notebook: fixest gave the Card-Krueger +2.75 DiD as a clustered interaction; bacondecomp diagnosed the twoway-FE bias by exposing the "Later vs Earlier Treated" forbidden comparisons (its weighted average matched the biased TWFE estimate exactly); and did (Callaway-Sant'Anna) repaired it, recovering the true ATT and producing a valid event study via ggdid.
This is the current applied-econometrics workflow for policy evaluation with staggered timing: estimate with a modern group-time estimator (did, sunab, did2s), never naive TWFE when effects may be heterogeneous; check pre-trends in the event study; and report the Goodman-Bacon decomposition when a TWFE number must be defended. Cross-links: DiD is the twoway fixed effects of the panel notebook plus a treatment indicator; next, Synthetic Control constructs a bespoke comparison unit when no natural control group exists.