Causal Inference I(c) — Noncompliance and Cluster Designs (R companion)¶
estimatr — iv_robust for CACE, cluster-robust difference_in_means, on two real experiments¶
estimatr gives design-based estimators for exactly these problems: difference_in_means for the intention-to-treat effect, iv_robust for the complier average causal effect (assignment instrumenting receipt), and difference_in_means(clusters=) / lm_robust(clusters=) for cluster-robust inference. This companion reproduces the Python notebook on the Sommer–Zeger vitamin A trial (real one-sided noncompliance) and Project STAR (real within-school clustering), plus the known-truth coverage simulation.
1. Noncompliance and the intention-to-treat effect — the vitamin A trial¶
Sommer & Zeger (1991): Indonesian villages randomized to a vitamin A program or control; ~20% of those assigned never received it (never-takers), and controls had no access (one-sided noncompliance). We rebuild the individual-level data from the published cell counts and compute the ITT — the effect of the program, comparing children as assigned.
suppressMessages(library(estimatr))
# Published Sommer-Zeger cell counts: assigned Z, received D, N, deaths
cells<-list(c(1,1,9675,12), c(1,0,2419,34), c(0,0,11588,74))
sz<-do.call(rbind, lapply(cells, function(r) data.frame(Z=r[1], D=r[2], death=c(rep(1,r[4]),rep(0,r[3]-r[4])))))
cat(sprintf("Sommer-Zeger vitamin A: %d children (%d vitamin-A, %d control); take-up when offered %.0f%%\n",
nrow(sz), sum(sz$Z==1), sum(sz$Z==0), 100*mean(sz$D[sz$Z==1])))
itt<-difference_in_means(death~Z, data=sz)
itt_y<-itt$coefficients[1]; itt_d<-mean(sz$D[sz$Z==1])-mean(sz$D[sz$Z==0])
cat(sprintf(" mortality: vitamin arm %.2f/1,000 control %.2f/1,000\n", 1000*mean(sz$death[sz$Z==1]), 1000*mean(sz$death[sz$Z==0])))
cat(sprintf(" ITT (program effect) = %+.2f per 1,000 (SE %.2f per 1,000); compliance gap ITT_D = %.2f\n",
1000*itt_y, 1000*itt$std.error[1], itt_d))
options(repr.plot.width=8, repr.plot.height=4.2)
bp<-barplot(c(control=1000*mean(sz$death[sz$Z==0]), `vitamin A`=1000*mean(sz$death[sz$Z==1])),
col=c("#2b6cb0","#2f855a"), ylab="mortality per 1,000", main=sprintf("ITT (as assigned) = %+.2f per 1,000", 1000*itt_y))
text(bp, c(1000*mean(sz$death[sz$Z==0]),1000*mean(sz$death[sz$Z==1]))-0.3, sprintf("%.2f",c(1000*mean(sz$death[sz$Z==0]),1000*mean(sz$death[sz$Z==1]))), col="white", font=2)
Warning message: "package 'estimatr' was built under R version 4.6.1"
Sommer-Zeger vitamin A: 23682 children (12094 vitamin-A, 11588 control); take-up when offered 80%
mortality: vitamin arm 3.80/1,000 control 6.39/1,000
ITT (program effect) = -2.58 per 1,000 (SE 0.93 per 1,000); compliance gap ITT_D = 0.80
2. As-treated is biased; CACE via iv_robust¶
Comparing children by what they took (as-treated) or compliers-to-all-controls (per-protocol) breaks randomization — refusers were sicker. The valid complier effect is CACE = ITT_Y / ITT_D, identical to 2SLS with assignment Z instrumenting receipt D. iv_robust(death ~ D | Z) computes it directly.
at <- mean(sz$death[sz$D==1]) - mean(sz$death[sz$D==0]) # as-treated
pp <- mean(sz$death[sz$Z==1 & sz$D==1]) - mean(sz$death[sz$Z==0]) # per-protocol
cace <- itt_y/itt_d # Bloom
iv <- iv_robust(death ~ D | Z, data=sz) # 2SLS
cat("Effect of vitamin A on mortality (per 1,000), four ways:\n")
cat(sprintf(" ITT (program) %+.2f valid, diluted\n", 1000*itt_y))
cat(sprintf(" as-treated (took vs didn't) %+.2f BIASED (refusers sicker -> overstates)\n", 1000*at))
cat(sprintf(" per-protocol %+.2f BIASED\n", 1000*pp))
cat(sprintf(" CACE = ITT_Y/ITT_D (Bloom) %+.2f compliers\n", 1000*cace))
cat(sprintf(" CACE via iv_robust (2SLS) %+.2f (SE %.2f) identical -- randomized instrument\n", 1000*iv$coefficients["D"], 1000*iv$std.error["D"]))
cat(sprintf("\n Among compliers, vitamin A cut mortality ~%.0f%% vs the control rate.\n", 100*(-cace)/mean(sz$death[sz$Z==0])))
options(repr.plot.width=8.5, repr.plot.height=4.4)
vals<-c(`ITT`=1000*itt_y, `as-treated`=1000*at, `per-protocol`=1000*pp, `CACE=IV`=1000*cace)
bp<-barplot(vals, col=c("#2b6cb0","#c53030","#dd6b20","#2f855a"), ylab="effect on mortality per 1,000",
main="ITT (diluted) vs CACE (compliers); as-treated exaggerates"); abline(h=0)
text(bp, vals-0.25, sprintf("%+.2f",vals), col="white", font=2)
Effect of vitamin A on mortality (per 1,000), four ways:
ITT (program) -2.58 valid, diluted
as-treated (took vs didn't) -6.47 BIASED (refusers sicker -> overstates)
per-protocol -5.15 BIASED
CACE = ITT_Y/ITT_D (Bloom) -3.23 compliers
CACE via iv_robust (2SLS) -3.23 (SE 1.16) identical -- randomized instrument
Among compliers, vitamin A cut mortality ~51% vs the control rate.
3. Cluster designs — Project STAR, cluster-robust standard errors¶
In Project STAR, class type is effectively assigned at the classroom/school level and reading scores are correlated within school (high ICC). lm_robust / difference_in_means with clusters = school give honest, cluster-robust (CR2) standard errors. The estimate is unchanged; the SE roughly doubles.
star<-read.csv("star_k.csv")
cat(sprintf("STAR kindergarten: %d students in %d schools (%d small, %d regular)\n",
nrow(star), length(unique(star$school)), sum(star$small), sum(1-star$small)))
mn <- lm_robust(read~small, data=star, se_type="HC2")
mc <- lm_robust(read~small, data=star, clusters=factor(school), se_type="CR2")
icc <- var(tapply(star$read, star$school, mean))/var(star$read)
cat(sprintf(" small-class effect on reading = %+.2f points\n", mn$coefficients["small"]))
cat(sprintf(" naive HC2 SE = %.3f 95%% CI [%.2f, %.2f]\n", mn$std.error["small"], mn$conf.low["small"], mn$conf.high["small"]))
cat(sprintf(" cluster-robust SE = %.3f 95%% CI [%.2f, %.2f] (%.2fx wider)\n", mc$std.error["small"], mc$conf.low["small"], mc$conf.high["small"], mc$std.error["small"]/mn$std.error["small"]))
cat(sprintf(" between-school variance share (ICC proxy) = %.3f\n", icc))
options(repr.plot.width=7.5, repr.plot.height=4.4)
est<-mn$coefficients["small"]
plot(c(0,1), c(est,est), ylim=c(0, est+2*mc$std.error["small"]), xlim=c(-.5,1.5), pch=19, col=c("#a0aec0","#2f855a"),
xaxt="n", xlab="", ylab="small-class effect (95% CI)", main=sprintf("Same estimate %+.1f, honest CI ~%.1fx wider", est, mc$std.error["small"]/mn$std.error["small"]))
axis(1, at=c(0,1), labels=c("naive HC2","cluster-robust"))
arrows(0, est-1.96*mn$std.error["small"], 0, est+1.96*mn$std.error["small"], angle=90, code=3, length=.08, col="#a0aec0", lwd=2)
arrows(1, est-1.96*mc$std.error["small"], 1, est+1.96*mc$std.error["small"], angle=90, code=3, length=.08, col="#2f855a", lwd=2)
abline(h=0, lty=3)
STAR kindergarten: 3743 students in 79 schools (1738 small, 2005 regular)
small-class effect on reading = +5.82 points
naive HC2 SE = 1.042 95% CI [3.78, 7.86]
cluster-robust SE = 1.856 95% CI [2.12, 9.52] (1.78x wider)
between-school variance share (ICC proxy) = 0.225
4. Why it matters — coverage under a cluster-randomized design (known truth)¶
Only a known-truth simulation shows the naive SE is wrong. We simulate a cluster-randomized experiment with a true effect of zero ($G=30$ clusters, $m=40$, ICC $=0.15$) and check 95% CI coverage. Naive intervals cover far below 95%; cluster-robust intervals are calibrated. The design effect $1+(m-1)\text{ICC}\approx6.8$ sets the SE inflation ($\sqrt{6.8}\approx2.6$).
set.seed(1)
one<-function(G=30,m=40,icc=0.15,tau=0){
su<-sqrt(icc); se<-sqrt(1-icc); tg<-sample(c(rep(1,G/2),rep(0,G-G/2)))
cl<-rep(1:G, each=m); a<-rep(rnorm(G,0,su), each=m); Tt<-rep(tg, each=m)
y<-tau*Tt + a + rnorm(G*m,0,se); d<-data.frame(cl,Tt,y)
nb<-lm_robust(y~Tt, data=d, se_type="HC2"); cb<-lm_robust(y~Tt, data=d, clusters=factor(cl), se_type="CR2")
c(est=as.numeric(nb$coefficients["Tt"]), cov_n=as.numeric(nb$conf.low["Tt"]<=0 & 0<=nb$conf.high["Tt"]),
cov_c=as.numeric(cb$conf.low["Tt"]<=0 & 0<=cb$conf.high["Tt"]), se_n=as.numeric(nb$std.error["Tt"]), se_c=as.numeric(cb$std.error["Tt"])) }
R<-t(sapply(1:400, function(i) one())); deff<-1+(40-1)*0.15
cat(sprintf("Cluster-randomized, true effect 0, 400 sims (G=30 x m=40, ICC=0.15):\n"))
cat(sprintf(" naive-SE 95%% coverage = %.0f%% <- should be 95: FALSE precision\n", 100*mean(R[,"cov_n"])))
cat(sprintf(" cluster-SE 95%% coverage = %.0f%% <- calibrated\n", 100*mean(R[,"cov_c"])))
cat(sprintf(" mean SE naive %.3f vs cluster %.3f (%.2fx); design effect %.1f, sqrt %.2f\n",
mean(R[,"se_n"]), mean(R[,"se_c"]), mean(R[,"se_c"])/mean(R[,"se_n"]), deff, sqrt(deff)))
options(repr.plot.width=7.5, repr.plot.height=4.2)
bp<-barplot(c(`naive SE`=100*mean(R[,"cov_n"]), `cluster SE`=100*mean(R[,"cov_c"])),
col=c("#c53030","#2f855a"), ylab="actual 95% CI coverage (%)", ylim=c(0,100), main="Naive SEs under-cover under clustering")
abline(h=95, lty=2); text(bp, c(100*mean(R[,"cov_n"]),100*mean(R[,"cov_c"]))-6, sprintf("%.0f%%",c(100*mean(R[,"cov_n"]),100*mean(R[,"cov_c"]))), col="white", font=2)
Cluster-randomized, true effect 0, 400 sims (G=30 x m=40, ICC=0.15):
naive-SE 95% coverage = 54% <- should be 95: FALSE precision
cluster-SE 95% coverage = 93% <- calibrated
mean SE naive 0.057 vs cluster 0.150 (2.61x); design effect 6.8, sqrt 2.62
5. Summary¶
estimatr reproduced both foundations on real experiments: difference_in_means gave the vitamin A ITT (−2.6 per 1,000, the diluted program effect), iv_robust gave the CACE = 2SLS complier effect (−3.2 per 1,000, ~50% mortality reduction) with randomized assignment as a textbook instrument, while as-treated (−6.5) was biased by who refuses; and lm_robust(clusters=) gave Project STAR's honest, ~1.8× wider cluster-robust SE, with a known-truth simulation confirming naive intervals cover only ~54%.
Guidance: report ITT as primary and CACE/IV as the complier complement (never as-treated); cluster SEs at the level of randomization. Cross-links: CACE is the IV/LATE estimator (subsection 3) with the cleanest possible instrument; cluster-robust inference recurs in panel data (subsection 5) and purged CV (Financial ML); the selection biasing as-treated is the confounding matching (subsection 2) tackles observationally.