Causal Inference III — Instrumental Variables (R companion)¶

AER/ivreg two-stage least squares with built-in diagnostics, and fixest::feols¶

This companion reproduces the Card returns-to-schooling IV with R's reference tools. R is especially strong here because its IV packages ship the diagnostic tests — weak-instrument F, Wu-Hausman endogeneity, Sargan overidentification — as first-class output, exactly the checks the Python notebook computed by hand:

  • ivreg::ivreg (equivalently AER::ivreg) — the standard 2SLS estimator, with summary(..., diagnostics = TRUE);
  • fixest::feols — fast IV via the y ~ controls | (D ~ Z) syntax, reporting the first-stage F automatically.

All estimates match the from-scratch Python notebook: OLS return ≈ 0.075, 2SLS ≈ 0.132, first-stage F ≈ 13.

1. Data and the first stage (relevance)¶

Card's (1995) sample: 3,010 men, outcome log wage (lwage), endogenous treatment years of education (educ), instrument nearc4 (grew up near a four-year college), with experience, race, region, and urban controls. The first stage confirms relevance: regressing education on the instrument and controls, proximity to a college raises schooling by about 0.3 years (conditional), and the association is highly significant.

In [1]:
suppressMessages({library(ivreg); library(AER); library(fixest)})
d<-read.csv("card.csv")
ctrl<-c("exper","expersq","black","south","smsa","reg661","reg662","reg663","reg664","reg665","reg666","reg667","reg668","smsa66")
cf<-paste(ctrl,collapse="+")
cat(sprintf("n = %d; near a 4-yr college: %.0f%%; mean educ %.1f yrs; mean log wage %.2f\n",
            nrow(d), 100*mean(d$nearc4), mean(d$educ), mean(d$lwage)))
fs<-lm(as.formula(paste("educ~nearc4+",cf)), data=d)
cat(sprintf("\nFIRST STAGE (relevance): nearc4 coefficient = %.3f years of schooling (t = %.2f)\n",
            coef(fs)["nearc4"], summary(fs)$coefficients["nearc4","t value"]))
options(repr.plot.width=7.5, repr.plot.height=4.2)
boxplot(educ~nearc4, data=d, names=c("not near","near 4-yr college"), col=c("#a0aec0","#2f855a"),
        ylab="years of education", main="First stage: the instrument shifts schooling")

Warning message:
"package 'ivreg' was built under R version 4.6.1"

Warning message:
"package 'fixest' was built under R version 4.6.1"

n = 3010; near a 4-yr college: 68%; mean educ 13.3 yrs; mean log wage 6.26

FIRST STAGE (relevance): nearc4 coefficient = 0.320 years of schooling (t = 3.64)
No description has been provided for this image

2. OLS vs 2SLS — ivreg with diagnostics¶

ivreg fits 2SLS with the two-part formula outcome ~ treatment + controls | instrument + controls. Its summary(diagnostics = TRUE) reports the tests that decide whether IV is trustworthy here. As in Python, 2SLS (~0.132) exceeds OLS (~0.075); the Wu-Hausman test is the formal check of whether OLS and IV differ significantly (here only weakly, given IV's larger standard error — an honest caveat), and weak-instruments reports the first-stage F.

In [2]:
ols<-lm(as.formula(paste("lwage~educ+",cf)), data=d)
iv <-ivreg(as.formula(paste("lwage~educ+",cf,"|nearc4+",cf)), data=d)
s<-summary(iv, diagnostics=TRUE)
cat(sprintf("OLS  return to schooling = %.4f  (SE %.4f)\n", coef(ols)["educ"], summary(ols)$coefficients["educ","Std. Error"]))
cat(sprintf("2SLS return to schooling = %.4f  (SE %.4f)\n\n", coef(iv)["educ"], s$coefficients["educ","Std. Error"]))
cat("ivreg diagnostics:\n"); print(round(s$diagnostics,3))
options(repr.plot.width=7.5, repr.plot.height=3.2)
est<-c(OLS=coef(ols)["educ"], IV=coef(iv)["educ"]); se<-c(summary(ols)$coefficients["educ","Std. Error"], s$coefficients["educ","Std. Error"])
plot(est, c(1,0.7), xlim=c(0,0.26), ylim=c(0.4,1.3), pch=c(19,15), col=c("#2b6cb0","#c53030"), cex=1.6, yaxt="n", ylab="",
     xlab="return to a year of schooling (log points)", main="OLS vs 2SLS — Card returns to schooling")
arrows(est-1.96*se, c(1,0.7), est+1.96*se, c(1,0.7), angle=90, code=3, length=.08, lwd=2, col=c("#2b6cb0","#c53030"))
legend("topright", c(sprintf("OLS %.3f",est[1]), sprintf("2SLS %.3f",est[2])), pch=c(19,15), col=c("#2b6cb0","#c53030"), bty="n")

OLS  return to schooling = 0.0747  (SE 0.0035)

2SLS return to schooling = 0.1315  (SE 0.0550)


ivreg diagnostics:

                 df1  df2 statistic p-value
Weak instruments   1 2994    13.256    0.00
Wu-Hausman         1 2993     1.168    0.28
Sargan             0   NA        NA      NA
No description has been provided for this image

3. LATE and the complier share¶

With heterogeneous effects, 2SLS estimates the Local Average Treatment Effect — the return for compliers, those induced into more schooling by the instrument (Imbens & Angrist, 1994), under relevance, exclusion, and monotonicity. The complier share is the first-stage jump in the probability of crossing into college; the rest of the sample are always-takers and never-takers whom the instrument does not move. This is why the Card LATE can exceed the OLS average: proximity moves credit-constrained students, whose return is high.

In [3]:
d$col<-as.integer(d$educ>=13)
comp<-mean(d$col[d$nearc4==1])-mean(d$col[d$nearc4==0])
always<-mean(d$col[d$nearc4==0]); never<-1-mean(d$col[d$nearc4==1])
cat(sprintf("At the 'entered college (educ>=13)' margin the instrument moves:\n"))
cat(sprintf("  compliers %.3f | always-takers %.3f | never-takers %.3f\n", comp, always, never))
options(repr.plot.width=7.5, repr.plot.height=4)
barplot(c(compliers=comp, `always-takers`=always, `never-takers`=never), col=c("#2f855a","#2b6cb0","#a0aec0"),
        ylab="share of sample", main="Who does the instrument move? (LATE = complier effect)")

At the 'entered college (educ>=13)' margin the instrument moves:

  compliers 0.122 | always-takers 0.422 | never-takers 0.456
No description has been provided for this image

4. Weak instruments — the diagnostic, and the Angrist-Krueger cautionary tale¶

A weak first stage biases 2SLS toward OLS and wrecks inference. The Staiger-Stock rule wants first-stage $F>10$; Card's nearc4 gives $F\approx13$, adequate. The famous failure is Angrist & Krueger's (1991) quarter-of-birth instrument: quarter of birth shifts schooling only trivially through compulsory-attendance laws, and Bound, Jaeger & Baker (1995) showed that with hundreds of weak quarter×year×state interactions the 2SLS estimates were driven by bias — even randomly generated instruments reproduced similar numbers. fixest::feols reports the first-stage F automatically; a compact simulation with a known effect reproduces the pathology — as instrument strength falls, the estimate drifts from the truth toward the confounded OLS value.

In [4]:
fe<-feols(as.formula(paste("lwage~",cf,"|educ~nearc4")), data=d)
cat(sprintf("fixest feols 2SLS return = %.4f ; first-stage F = %.2f\n\n", coef(fe)["fit_educ"], fitstat(fe,"ivf1")$ivf1$stat))
set.seed(0); beta<-0.5; reps<-300; strengths<-c(.02,.05,.1,.15,.25,.4,.6); ivm<-Fm<-numeric(length(strengths))
tsls<-function(y,D,Z){z1<-cbind(1,Z); x1<-cbind(1,D); ph<-z1%*%solve(crossprod(z1),crossprod(z1,x1)); solve(crossprod(ph),crossprod(ph,y))[2]}
for(k in seq_along(strengths)){pi<-strengths[k]; iv_<-Fs<-numeric(reps)
  for(r in 1:reps){m<-800; U<-rnorm(m); Z<-rnorm(m); D<-pi*Z+U+rnorm(m); Y<-beta*D+U+rnorm(m)
    iv_[r]<-tsls(Y,D,Z); ff<-lm(D~Z); Fs[r]<-summary(ff)$coefficients["Z","t value"]^2}
  ivm[k]<-mean(iv_); Fm[k]<-mean(Fs)}
olm<-{U<-rnorm(20000);D<-0.3*rnorm(20000)+U+rnorm(20000);Y<-beta*D+U+rnorm(20000);coef(lm(Y~D))[2]}
options(repr.plot.width=8, repr.plot.height=4.4)
plot(Fm, ivm, type="b", pch=19, col="#c53030", lwd=2, log="x", xlab="first-stage F (log scale)", ylab="mean 2SLS estimate",
     main="Weak instruments bias 2SLS toward OLS", ylim=range(c(ivm,beta,olm)))
abline(h=beta, col="#2f855a", lty=2, lwd=2); abline(h=olm, col="#2b6cb0", lty=3, lwd=2); abline(v=10, col="black", lty=3)
legend("bottomright", c("2SLS estimate","true 0.5","OLS (confounded)","F=10"), col=c("#c53030","#2f855a","#2b6cb0","black"),
       lty=c(1,2,3,3), pch=c(19,NA,NA,NA), lwd=2, bty="n", cex=.85)

fixest feols 2SLS return = 0.1315 ; first-stage F = 13.26
No description has been provided for this image

5. Summary¶

R's IV packages reproduced the Card analysis and, crucially, supplied the diagnostics as standard output: ivreg gave OLS 0.075 vs 2SLS 0.132 with the weak-instrument F (≈13) and Wu-Hausman endogeneity test in one call, and fixest::feols confirmed the estimate and first-stage F. The complier-share decomposition made the LATE interpretation concrete, and a compact simulation reproduced the weak-instrument bias behind the Bound-Jaeger-Baker critique of the Angrist-Krueger quarter-of-birth instrument.

The IV workflow to carry forward: argue the exclusion restriction, test relevance (first-stage F), report the Wu-Hausman and — when over-identified — Sargan tests, and interpret the estimate as a LATE. Next, Regression Discontinuity, where identification comes from a sharp threshold rule rather than an excluded instrument, and the fuzzy-RD design is estimated by exactly this 2SLS machinery.