Causal Inference II(b) — Sensitivity Analysis (R companion)¶
rbounds for Rosenbaum bounds and EValue for the E-value¶
R has the reference packages for both frameworks:
rbounds::psens— Rosenbaum bounds (Wilcoxon signed-rank) for matched-pair data, tabulating the worst-case p-value across values of $\Gamma$;EValue::evalues.RR— the VanderWeele-Ding E-value for a risk ratio and its confidence interval.
Using the exact matched-pair differences exported from the Python notebook, R reproduces every number: the strong effect tolerates $\Gamma\approx7.8$, the weak one breaks at $\Gamma\approx1$, the LaLonde matching estimate is fragile at $\Gamma\approx1.25$, and a smoking-scale risk ratio has an E-value of ~7.
1. Rosenbaum bounds — rbounds::psens¶
psens(x, y, Gamma, GammaInc) runs the Wilcoxon signed-rank sensitivity analysis on matched pairs; passing the vector of within-pair differences against a zero vector tabulates the worst-case (upper-bound) p-value as $\Gamma$ increases. We read off $\Gamma^\star$ — the point where the upper-bound p-value crosses 0.05 — for a strong and a weak simulated effect (same noise). The strong effect stays significant up to $\Gamma\approx7.8$; the weak one loses significance immediately.
suppressMessages({library(rbounds); library(EValue)})
d<-read.csv("sens_sim_diffs.csv")
gstar<-function(x){b<-psens(x, rep(0,length(x)), Gamma=10, GammaInc=0.05)$bounds
hi<-b[b[,3]>0.05,1]; if(length(hi)) hi[1] else Inf}
gs<-gstar(d$strong); gw<-gstar(d$weak)
cat(sprintf("STRONG effect: Rosenbaum Gamma* = %.2f (robust)\n", gs))
cat(sprintf("WEAK effect : Rosenbaum Gamma* = %.2f (fragile)\n", gw))
bs<-psens(d$strong, rep(0,nrow(d)), Gamma=8, GammaInc=0.1)$bounds
bw<-psens(d$weak, rep(0,nrow(d)), Gamma=8, GammaInc=0.1)$bounds
options(repr.plot.width=8.5, repr.plot.height=4.6)
plot(bs[,1], bs[,3], type="l", col="#2f855a", lwd=2.5, ylim=c(0,0.3), xlab="Gamma (hidden-bias sensitivity)",
ylab="worst-case p-value", main="Rosenbaum bounds: Gamma* = where significance is lost")
lines(bw[,1], bw[,3], col="#dd6b20", lwd=2.5); abline(h=0.05, col="#c53030", lty=2)
abline(v=gs, col="#2f855a", lty=3); abline(v=gw, col="#dd6b20", lty=3)
legend("topleft", c(sprintf("strong (Gamma*=%.1f)",gs), sprintf("weak (Gamma*=%.2f)",gw), "p=0.05"),
col=c("#2f855a","#dd6b20","#c53030"), lwd=2, lty=c(1,1,2), bty="n")
Warning message: "package 'EValue' was built under R version 4.6.1"
STRONG effect: Rosenbaum Gamma* = 7.75 (robust)
WEAK effect : Rosenbaum Gamma* = 1.00 (fragile)
2. The LaLonde matching estimate — fragile¶
The 185 matched-pair earnings differences from the Matching notebook (which recovered the experimental benchmark of ~\$1,792). psens shows the estimate is robust to hidden bias only up to $\Gamma\approx1.25$: an unmeasured confounder making one member of a pair 1.25× more likely to enrol would overturn it. Recovering the truth did not make it robust — the quantified version of "unconfoundedness is untestable".
lal<-read.csv("sens_lalonde_diffs.csv")$lalonde
cat(sprintf("LaLonde matched pairs: n=%d, mean difference $%.0f, Rosenbaum Gamma* = %.2f (FRAGILE)\n", length(lal), mean(lal), gstar(lal)))
b<-psens(lal, rep(0,length(lal)), Gamma=2.5, GammaInc=0.05)$bounds
options(repr.plot.width=8, repr.plot.height=4.4)
plot(b[,1], b[,3], type="l", col="#2b6cb0", lwd=2.5, xlab="Gamma", ylab="worst-case p-value",
main="LaLonde matching estimate: robust only to Gamma ~ 1.25")
abline(h=0.05, col="#c53030", lty=2); abline(v=gstar(lal), col="#2f855a", lty=3)
legend("bottomright", c(sprintf("Gamma* = %.2f",gstar(lal)),"p=0.05"), col=c("#2f855a","#c53030"), lty=c(3,2), lwd=2, bty="n")
LaLonde matched pairs: n=185, mean difference $1792, Rosenbaum Gamma* = 1.25 (FRAGILE)
3. The E-value — EValue::evalues.RR¶
evalues.RR(est, lo, hi) returns the E-value for a risk ratio and for its confidence-interval limit. A smoking-scale $RR=3.9$ gives an E-value of ~7.3 (point) and ~5.5 (CI): a confounder would need ~7-fold associations with both treatment and outcome to explain the effect, and ~5.5-fold to erase its significance — implausibly strong, so the finding is robust. A weak $RR=1.3$ gives ~1.9 and ~1.3 — a modest confounder could explain it away.
strong_ev<-evalues.RR(est=3.9, lo=3.0, hi=5.1)
weak_ev <-evalues.RR(est=1.3, lo=1.05, hi=1.61)
cat("STRONG (smoking-scale) RR=3.9:\n"); print(round(strong_ev,2))
cat("\nWEAK RR=1.3:\n"); print(round(weak_ev,2))
ev<-function(rr){rr<-ifelse(rr>=1,rr,1/rr); rr+sqrt(rr*(rr-1))}
rrs<-seq(1.05,6,length.out=100)
options(repr.plot.width=8, repr.plot.height=4.4)
plot(rrs, sapply(rrs,ev), type="l", col="#6b46c1", lwd=2.5, xlab="observed risk ratio", ylab="E-value",
main="E-value: bigger effects need stronger confounders to explain away")
points(c(3.9,1.3), c(ev(3.9),ev(1.3)), pch=19, cex=1.4, col=c("#2f855a","#dd6b20"))
text(3.9, ev(3.9), sprintf(" smoking-scale (E=%.1f)",ev(3.9)), pos=2, cex=.85)
text(1.3, ev(1.3), sprintf(" weak (E=%.1f)",ev(1.3)), pos=4, cex=.85)
STRONG (smoking-scale) RR=3.9:
point lower upper RR 3.90 3.00 5.1 E-values 7.26 5.45 NA
WEAK RR=1.3:
point lower upper RR 1.30 1.05 1.61 E-values 1.92 1.28 NA
4. Summary¶
rbounds and EValue reproduced the Python notebook's sensitivity analyses exactly: Rosenbaum $\Gamma^\star\approx7.8$ (strong, robust) vs $\approx1$ (weak, fragile), the LaLonde matching estimate fragile at $\Gamma\approx1.25$, and E-values of ~7 (smoking-scale, robust) vs ~1.9 (weak). These are the reference implementations to reach for: psens/binarysens for matched designs, evalues.RR/evalues.OLS for ratio and continuous estimates.
The discipline this notebook adds to the whole arc: an observational effect estimate is incomplete without a sensitivity analysis. Report $\Gamma^\star$ (matched) or an E-value (ratio), for both the point estimate and the interval, and judge it against the confounders that could plausibly exist. It is the honest companion to matching, IPW, doubly-robust, DML, and causal-survival estimation — every method that assumes unconfoundedness — and the empirical face of the DAG's unmeasured-confounder arrow.