Causal Inference V — Panel Data & Fixed Effects (R companion)¶
plm and fixest on the Stock-Watson traffic-fatalities panel¶
R's plm is the reference package for panel econometrics, and fixest is the fast modern workhorse for high-dimensional fixed effects. This companion reproduces the Python notebook's specification ladder — pooled OLS, within (FE), first-difference, random effects, the Hausman test, and twoway FE with clustered standard errors — with these two packages, on the same beer-tax / traffic-fatalities data. All estimates match the from-scratch Python results: pooled +0.36 (wrong sign), within/twoway ≈ −0.65.
1. Data and the pooled wrong-sign puzzle¶
The Stock-Watson panel (48 states × 7 years, 1982–1988): fatality rate frate (per 10,000), real beertax, and covariates. Pooled OLS ignores state heterogeneity and delivers the notorious positive coefficient — higher beer taxes appear to accompany more traffic deaths — a textbook omitted-variable artifact of comparing intrinsically different states.
suppressMessages({library(plm); library(fixest)})
d<-read.csv("fatalities.csv"); pd<-pdata.frame(d, index=c("state","year"))
cat(sprintf("%d states x %d years = %d obs\n", length(unique(d$state)), length(unique(d$year)), nrow(d)))
pool<-plm(frate~beertax, pd, model="pooling")
cat(sprintf("\nPooled OLS beertax = %+.3f <-- POSITIVE (confounded by state heterogeneity)\n", coef(pool)["beertax"]))
options(repr.plot.width=7.5, repr.plot.height=4.4)
plot(d$beertax, d$frate, pch=19, col=rgb(.6,.6,.6,.5), xlab="beer tax", ylab="fatality rate (per 10,000)",
main="Pooled: high-tax states look deadlier (state heterogeneity, not causation)")
abline(pool, col="#c53030", lwd=2.5)
Warning message: "package 'plm' was built under R version 4.6.1"
Warning message: "package 'fixest' was built under R version 4.6.1"
48 states x 7 years = 336 obs
Pooled OLS beertax = +0.365 <-- POSITIVE (confounded by state heterogeneity)
2. Within (FE) and first-difference — plm and fixest¶
plm(model = "within") is the fixed-effects estimator (demean within states); model = "fd" is first-difference; and fixest::feols(y ~ x | state) fits the same FE with clustered SEs by default. Removing state fixed effects flips the sign to negative — the expected causal effect of beer taxes on traffic deaths — and feols matches the from-scratch and plm estimates exactly.
fe <-plm(frate~beertax, pd, model="within")
fd <-plm(frate~beertax, pd, model="fd")
few<-feols(frate~beertax | state, data=d)
cat(sprintf("pooled OLS beertax = %+.3f\n", coef(plm(frate~beertax,pd,model='pooling'))["beertax"]))
cat(sprintf("within (FE) plm beertax = %+.3f\n", coef(fe)["beertax"]))
cat(sprintf("fixest feols|state beertax= %+.3f (clustered SE %.3f)\n", coef(few)["beertax"], se(few)["beertax"]))
cat(sprintf("first-difference beertax = %+.3f (= FE only at T=2)\n", coef(fd)["beertax"]))
options(repr.plot.width=8, repr.plot.height=3.6)
vals<-c(pooled=coef(plm(frate~beertax,pd,model='pooling'))["beertax"], `first-diff`=coef(fd)["beertax"], `within (FE)`=coef(fe)["beertax"])
bp<-barplot(vals, horiz=TRUE, las=1, col=c("#c53030","#dd6b20","#2f855a"), xlim=c(-0.8,0.5),
xlab="beer-tax coefficient", main="Removing state fixed effects flips the sign")
abline(v=0, col="black"); text(vals+ifelse(vals>=0,0.05,-0.05), bp, sprintf("%+.3f",vals), cex=.85)
pooled OLS beertax = +0.365
within (FE) plm beertax = -0.656
fixest feols|state beertax= -0.656 (clustered SE 0.188)
first-difference beertax = +0.014 (= FE only at T=2)
3. Fixed vs random effects — plm::phtest (Hausman)¶
plm(model = "random") fits random effects (efficient if the heterogeneity is uncorrelated with the regressor), and phtest() runs the Hausman test comparing FE and RE. As in Python, RE is pulled toward the confounded pooled estimate, and Hausman rejects it — so fixed effects are required.
re<-plm(frate~beertax, pd, model="random")
h<-phtest(fe, re)
cat(sprintf("Fixed effects beertax = %+.3f\n", coef(fe)["beertax"]))
cat(sprintf("Random effects beertax = %+.3f (pulled toward pooled)\n", coef(re)["beertax"]))
cat(sprintf("\nHausman test: chisq = %.2f, df = %d, p = %.4f\n", h$statistic, h$parameter, h$p.value))
cat(" p < 0.05 -> RE inconsistent -> USE FIXED EFFECTS.\n")
Fixed effects beertax = -0.656
Random effects beertax = -0.052 (pulled toward pooled)
Hausman test: chisq = 18.35, df = 1, p = 0.0000
p < 0.05 -> RE inconsistent -> USE FIXED EFFECTS.
4. Twoway fixed effects and the DiD bridge — fixest¶
Adding year fixed effects removes national shocks common to all states. feols(y ~ x | state + year) is the twoway FE estimator (clustering by state by default), and plm(effect = "twoways") agrees. This is precisely the difference-in-differences engine of the next notebook — a treatment indicator plus unit and time fixed effects. The effect stays negative after purging year shocks and after adding covariates.
tw <-feols(frate~beertax | state+year, data=d)
twc<-feols(frate~beertax+unemp+income+drinkage+miles | state+year, data=d)
cat(sprintf("twoway FE (state+year) beertax = %+.3f (clustered SE %.3f, p %.3f)\n", coef(tw)["beertax"], se(tw)["beertax"], pvalue(tw)["beertax"]))
cat(sprintf("twoway FE + controls beertax = %+.3f (clustered SE %.3f, p %.3f)\n", coef(twc)["beertax"], se(twc)["beertax"], pvalue(twc)["beertax"]))
options(repr.plot.width=8.5, repr.plot.height=4)
vals<-c(pooled=coef(plm(frate~beertax,pd,model='pooling'))["beertax"], `random eff`=coef(re)["beertax"],
`state FE`=coef(fe)["beertax"], `twoway FE`=coef(tw)["beertax"], `twoway+ctrl`=coef(twc)["beertax"])
bp<-barplot(vals, col=c("#c53030","#dd6b20","#2f855a","#2b6cb0","#6b46c1"), ylab="beer-tax coefficient",
main="Specification ladder: pooled OLS -> twoway FE", ylim=c(-0.8,0.5), names.arg=names(vals), cex.names=.8)
abline(h=0, col="black"); text(bp, vals+ifelse(vals>=0,0.04,-0.05), sprintf("%+.2f",vals), cex=.8)
twoway FE (state+year) beertax = -0.640 (clustered SE 0.197, p 0.001)
twoway FE + controls beertax = -0.517 (clustered SE 0.168, p 0.002)
5. Summary¶
plm and fixest reproduced the panel toolkit end to end and confirmed the from-scratch Python results: pooled OLS +0.36 (wrong sign) → within/twoway FE ≈ −0.65 (higher beer taxes reduce traffic deaths), phtest rejected random effects (mandating FE), and feols delivered twoway FE with clustered standard errors in one line. fixest is the tool to remember for the rest of the arc — its y ~ x | fe1 + fe2 syntax with automatic clustering is exactly what difference-in-differences and event studies require.
The key econometric messages: fixed effects difference away time-invariant confounding, the Hausman test chooses FE vs RE, and twoway FE is difference-in-differences generalized. Cross-links to the BVAR / Multivariate Time Series arc (same panel data, dynamic-forecasting question) and forward to Difference-in-Differences, where the twoway-FE engine meets a treatment that switches on over time — and the modern staggered-adoption critique of naive TWFE.