Causal Inference I — Randomized Experiments & Randomization Inference (R companion)¶
HistData, fisher.test, ri2, and estimatr on the same Fisher classics¶
This R companion reproduces the Python notebook's three foundations with the field's established packages, in the language where the randomization-inference and design-based-estimation traditions actually live:
fisher.test(base R) for the Lady Tasting Tea exact test;HistData::ZeaMays— Darwin's cross- vs self-fertilized maize, the very data in the package — analysed with the paired $t$-test and withri2(Alexander Coppock's randomization inference package), the modern implementation of Fisher's permutation logic on a declared experimental design;estimatr::difference_in_means, which computes Neyman's design-based estimate, standard error, and confidence interval for the matched-pair experiment directly.
The point of the R side, as throughout the portfolio, is to lead with the reference tools and confirm the from-scratch Python results against them. Data and conclusions match the Python notebook.
1. The Rubin Causal Model — randomization removes confounding¶
The same known-truth simulation as the Python notebook: a confounder $C$ drives both treatment and outcome, so the naive difference in means is biased; randomizing the same units removes it. Potential outcomes $Y(0),Y(1)$ are constructed with a known ATE of 4, and the randomization distribution of the difference-in-means estimator (re-randomizing with the potential outcomes held fixed) is centred on the true ATE — the design-based logic Neyman formalised.
set.seed(7); options(repr.plot.width=12, repr.plot.height=4.6)
N<-2000; ATE<-4
C<-rnorm(N); Y0<-50+5*C+rnorm(N,0,3); Y1<-Y0+ATE
p<-plogis(1.5*C); Wc<-rbinom(N,1,p); Yc<-Wc*Y1+(1-Wc)*Y0 # confounded
Wr<-rbinom(N,1,0.5); Yr<-Wr*Y1+(1-Wr)*Y0 # randomized
cat(sprintf("True ATE = %.1f\n", ATE))
cat(sprintf("(a) confounded diff-in-means = %.2f -> biased by %+.2f\n", mean(Yc[Wc==1])-mean(Yc[Wc==0]), (mean(Yc[Wc==1])-mean(Yc[Wc==0]))-ATE))
cat(sprintf("(b) randomized diff-in-means = %.2f -> unbiased\n", mean(Yr[Wr==1])-mean(Yr[Wr==0])))
draws<-replicate(4000,{w<-rbinom(N,1,0.5); y<-w*Y1+(1-w)*Y0; mean(y[w==1])-mean(y[w==0])})
par(mfrow=c(1,2))
hist(C[Wc==1],breaks=30,col=rgb(.77,.19,.19,.5),freq=FALSE,main="Confounded assignment imbalances C",xlab="confounder C")
hist(C[Wc==0],breaks=30,col=rgb(.17,.42,.69,.5),freq=FALSE,add=TRUE)
legend("topright",c("treated","control"),fill=c(rgb(.77,.19,.19,.5),rgb(.17,.42,.69,.5)),bty="n",cex=.8)
hist(draws,breaks=40,col="#2f855a",freq=FALSE,main="Randomization distribution of diff-in-means",xlab="estimate over re-randomizations")
abline(v=ATE,col="#c53030",lwd=2); abline(v=mean(draws),col="black",lty=2)
legend("topright",c(sprintf("true ATE=%.0f",ATE),sprintf("mean=%.2f",mean(draws))),col=c("#c53030","black"),lwd=2,lty=c(1,2),bty="n",cex=.8)
par(mfrow=c(1,1))
cat(sprintf("\nMean of estimator over 4000 re-randomizations = %.2f = ATE.\n", mean(draws)))
True ATE = 4.0
(a) confounded diff-in-means = 9.29 -> biased by +5.29
(b) randomized diff-in-means = 4.06 -> unbiased
Mean of estimator over 4000 re-randomizations = 4.00 = ATE.
2. The Lady Tasting Tea — fisher.test¶
Fisher's exact test is the base-R fisher.test. Muriel Bristol correctly identified all 4 milk-first cups out of 8 (told there were 4 of each). The 2×2 table of her guesses against the truth is perfectly diagonal; the one-sided exact $p$-value is $1/\binom{8}{4}=1/70$, computed from the hypergeometric distribution that fisher.test uses internally and that dhyper lets us plot.
options(repr.plot.width=7.5, repr.plot.height=4.2)
tea<-matrix(c(4,0,0,4),2,dimnames=list(guess=c("milk-1st","tea-1st"),truth=c("milk-1st","tea-1st")))
ft<-fisher.test(tea, alternative="greater")
cat("Lady's guesses vs truth:\n"); print(tea)
cat(sprintf("\nfisher.test one-sided p = %.4f = 1/70 = %.4f\n", ft$p.value, 1/70))
x<-0:4; pmf<-dhyper(x,4,4,4)
bp<-barplot(pmf,names.arg=x,col=ifelse(x==4,"#c53030","#a0aec0"),ylim=c(0,0.6),
main="Exact null distribution (hypergeometric) — Lady Tasting Tea",
xlab="cups correctly identified under guessing",ylab="probability under H0")
text(bp,pmf+0.03,sprintf("%.3f",pmf),cex=.8); text(bp[5],pmf[5]+0.09,"observed\np=1/70",col="#c53030",cex=.9)
Lady's guesses vs truth:
truth guess milk-1st tea-1st milk-1st 4 0 tea-1st 0 4
fisher.test one-sided p = 0.0143 = 1/70 = 0.0143
3. Darwin's maize — HistData::ZeaMays + ri2 randomization inference¶
Darwin's 15 matched pairs of Zea mays ship in HistData::ZeaMays (cross, self, and their diff, heights in inches). We test the sharp null of no fertilization effect two ways: the classical paired $t$-test, and ri2::conduct_ri, which declares the matched-pair assignment (declare_ra(blocks = pair) — one plant per pot randomly labelled "cross") and builds the exact randomization distribution by re-drawing that assignment under the sharp null. This is Fisher's permutation test, packaged as randomization inference on a declared design — and its $p$-value matches the from-scratch enumeration of all $2^{15}$ sign flips in the Python notebook (≈0.026 one-sided / ≈0.053 two-sided).
suppressMessages({library(HistData); library(ri2)})
options(repr.plot.width=8, repr.plot.height=4.6)
data(ZeaMays); cross<-ZeaMays$cross; self<-ZeaMays$self; d<-cross-self
cat(sprintf("Darwin's %d pairs: mean cross-fertilization advantage = %.3f in; %d/%d pairs favour cross.\n",
length(d), mean(d), sum(d>0), length(d)))
tt<-t.test(cross,self,paired=TRUE)
cat(sprintf("Paired t-test: t = %.4f, two-sided p = %.5f (Fisher's t=2.148)\n", tt$statistic, tt$p.value))
dat<-data.frame(height=c(cross,self), Z=c(rep(1,15),rep(0,15)), pair=factor(rep(1:15,2)))
decl<-with(dat, declare_ra(blocks=pair))
set.seed(1)
ri<-conduct_ri(height~Z, declaration=decl, assignment="Z", sharp_hypothesis=0, data=dat, sims=50000)
cat("\nri2 randomization inference (declared matched-pair design):\n"); print(summary(ri))
print(plot(ri) + ggplot2::ggtitle("ri2: exact randomization distribution — Darwin's maize"))
Warning message: "package 'HistData' was built under R version 4.6.1"
Warning message: "package 'ri2' was built under R version 4.6.1"
Darwin's 15 pairs: mean cross-fertilization advantage = 2.617 in; 13/15 pairs favour cross.
Paired t-test: t = 2.1480, two-sided p = 0.04970 (Fisher's t=2.148)
ri2 randomization inference (declared matched-pair design):
term estimate two_tailed_p_value 1 Z 2.616667 0.05267334
4. Neyman's estimate — estimatr::difference_in_means¶
estimatr::difference_in_means implements Neyman's design-based inference directly. Told the design is a matched pair (blocks = pair), it returns the block-adjusted average treatment effect — which for matched pairs is exactly the mean of the within-pair differences — together with the design-based standard error and 95% confidence interval. The estimate (~2.6 inches, CI excluding 0) matches the Python notebook: Fisher's test says an effect exists; Neyman's estimate says how big, and how sure.
suppressMessages(library(estimatr))
options(repr.plot.width=8, repr.plot.height=3.2)
dm<-difference_in_means(height~Z, blocks=pair, data=dat)
cat("Neyman design-based estimate (matched-pair):\n"); print(dm)
td<-tidy(dm); est<-td$estimate[1]; lo<-td$conf.low[1]; hi<-td$conf.high[1]
plot(est,1,xlim=c(-1,6),ylim=c(0.5,1.5),pch=19,cex=1.6,col="#2b6cb0",yaxt="n",ylab="",
xlab="average height advantage of cross-fertilization (inches)",
main="Neyman estimation: effect size with 95% CI")
arrows(lo,1,hi,1,angle=90,code=3,length=0.08,lwd=2,col="#2b6cb0")
abline(v=0,col="#c53030",lwd=2,lty=2)
legend("topleft",c(sprintf("ATE = %.2f in",est),sprintf("95%% CI [%.2f, %.2f]",lo,hi),"no effect"),
col=c("#2b6cb0","#2b6cb0","#c53030"),pch=c(19,NA,NA),lwd=c(NA,2,2),lty=c(NA,1,2),bty="n",cex=.85)
Warning message: "package 'estimatr' was built under R version 4.6.1"
Neyman design-based estimate (matched-pair):
Design: Matched-pair Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF Z 2.616667 1.218195 2.147987 0.04970294 0.003899165 5.229434 14
5. Summary¶
The R companion reproduces the foundations with the established tools: fisher.test gave the Lady Tasting Tea's exact $p=1/70$; HistData::ZeaMays supplied Darwin's maize, where the paired $t$-test ($t=2.148$) and ri2 randomization inference (a declared matched-pair design) agreed on a real cross-fertilization effect, matching the from-scratch $2^{15}$-enumeration of the Python notebook; and estimatr::difference_in_means delivered Neyman's design-based estimate and confidence interval (~2.6 inches, excluding zero).
ri2 is the package to remember for this whole arc — it makes randomization inference on an explicitly declared design the default, which is exactly the mindset (assignment mechanism first) that the observational subsections will have to assume rather than know. Next: Potential Outcomes & Matching, where MatchIt, WeightIt, and cobalt stand in for the randomization we no longer have.