Causal Inference IV — Regression Discontinuity (R companion)¶

rdrobust, rdplot, and rddensity — the reference implementation¶

Regression discontinuity is a rare case where the authors of the methodology maintain the canonical software in R, so this companion leads with it directly. The Calonico-Cattaneo-Titiunik suite is the field standard:

  • rdrobust — local-polynomial RD estimation with the MSE-optimal (CCT) bandwidth and bias-corrected, robust confidence intervals;
  • rdplot — the automatic, data-driven RD plot;
  • rddensity — the modern (Cattaneo-Jansson-Ma) manipulation test, a tuning-free successor to McCrary.

We reproduce the Python notebook's Lee (2008) incumbency-advantage results exactly, then close with the fuzzy-RD design that rdrobust estimates via its built-in 2SLS.

1. The RD plot — rdplot¶

The Lee (2008) data (from rddtools): running variable = Democratic margin of victory, outcome = Democratic vote share in the next election, cutoff = 0. rdplot produces the canonical figure — binned means of the outcome against the running variable with a fitted polynomial on each side — and the vertical gap at the cutoff is the incumbency advantage, visible by eye.

In [1]:
suppressMessages({library(rdrobust); library(rddensity)})
d<-read.csv("lee2008.csv"); x<-d$margin; y<-d$dem_voteshare_next
cat(sprintf("Lee (2008): n = %d U.S. House races; cutoff c = 0 (Democrat wins => incumbent)\n", length(x)))
options(repr.plot.width=8.5, repr.plot.height=5)
invisible(rdplot(y, x, c=0, p=4, x.label="Democratic margin of victory",
       y.label="Democratic vote share, next election",
       title="rdplot — incumbency advantage in U.S. House elections"))

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

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

Lee (2008): n = 6558 U.S. House races; cutoff c = 0 (Democrat wins => incumbent)

Warning message in rdplot(y, x, c = 0, p = 4, x.label = "Democratic margin of victory", :
"Mass points detected in the running variable."
No description has been provided for this image

2. The RD estimate — rdrobust with the CCT bandwidth¶

rdrobust fits local-linear regressions with a triangular kernel, selects the MSE-optimal bandwidth automatically, and reports three inference rows: conventional, bias-corrected, and robust (the recommended one, valid at the optimal bandwidth). The point estimate — the jump in next-election vote share at the winning threshold — matches the from-scratch Python result (≈ 0.06), and the robust confidence interval excludes zero: a real, sizeable incumbency advantage.

In [2]:
r<-rdrobust(y, x, c=0)
cat(sprintf("CCT optimal bandwidth  = %.3f\n", r$bws[1,1]))
cat(sprintf("RD estimate (conventional) = %.4f\n", r$coef[1]))
cat(sprintf("Robust 95%% CI = [%.4f, %.4f],  robust p = %.4g\n\n", r$ci[3,1], r$ci[3,2], r$pv[3]))
print(summary(r))

Warning message in rdrobust(y, x, c = 0):
"Mass points detected in the running variable."

CCT optimal bandwidth  = 0.136

RD estimate (conventional) = 0.0637

Robust 95% CI = [0.0348, 0.0839],  robust p = 2.155e-06


Call: rdrobust

Sharp RD estimates using local polynomial regression.

Number of Obs.                 6558
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

                               Left        Right
Number of Obs.                 2740         3818
Eff. Number of Obs.             789          817
Order est. (p)                    1            1
Order bias (q)                    2            2
BW est. (h)                   0.136        0.136
BW bias (b)                   0.240        0.240
rho (h/b)                     0.565        0.565
Unique Obs.                    2108         2581

=====================================================================
                   Point    Robust Inference
                Estimate         z     P>|z|      [ 95% C.I. ]       
---------------------------------------------------------------------
     RD Effect     0.064     4.738     0.000     [0.035 , 0.084]     
=====================================================================

3. Manipulation test — rddensity¶

RD requires that units cannot precisely sort across the cutoff. rddensity implements the Cattaneo-Jansson-Ma density test: a discontinuity in the density of the running variable at the cutoff would signal manipulation. For elections it should pass — a razor-thin win cannot be engineered — and it does (p well above 0.05). rdplotdensity shows the estimated density with its confidence bands meeting smoothly at zero.

In [3]:
options(repr.plot.width=8, repr.plot.height=4.6)
dens<-rddensity(X=x, c=0)
cat(sprintf("Cattaneo-Jansson-Ma manipulation test:  T = %.3f,  p = %.3f\n", dens$test$t_jk, dens$test$p_jk))
cat("p > 0.05 -> fail to reject continuous density -> no evidence of manipulation -> design valid.\n")
invisible(rdplotdensity(dens, X=x, plotRange=c(-0.5,0.5),
            xlabel="Democratic margin of victory", ylabel="density",
            title="rddensity — running-variable density is smooth at the cutoff"))

Cattaneo-Jansson-Ma manipulation test:  T = 1.435,  p = 0.151

p > 0.05 -> fail to reject continuous density -> no evidence of manipulation -> design valid.
No description has been provided for this image

4. Fuzzy RD — rdrobust with a treatment indicator¶

When crossing the cutoff changes only the probability of treatment, rdrobust accepts a treatment variable via the fuzzy = argument and returns the fuzzy estimate — the outcome jump divided by the treatment jump, i.e. 2SLS with "above the cutoff" as the instrument (the connection to the Instrumental Variables notebook). We reproduce the Python simulation with a known effect of 3: the sharp jump understates it, and the fuzzy estimate recovers it.

In [4]:
set.seed(2); m<-8000; xf<-runif(m,-1,1); TRUE_EFF<-3
pr<-ifelse(xf>=0,0.85,0.25); Df<-rbinom(m,1,pr); Yf<-2+TRUE_EFF*Df+1.0*xf+rnorm(m,0,0.7)
sharp<-rdrobust(Yf, xf, c=0)                        # ITT: jump in outcome only
fuzzy<-rdrobust(Yf, xf, c=0, fuzzy=Df)             # fuzzy: outcome jump / treatment jump
cat(sprintf("True effect                      = %d\n", TRUE_EFF))
cat(sprintf("Sharp (ITT) RD jump in outcome   = %.3f   (biased toward 0)\n", sharp$coef[1]))
cat(sprintf("Fuzzy RD (rdrobust, 2SLS at cutoff) = %.3f   (recovers the truth)\n", fuzzy$coef[1]))
options(repr.plot.width=7.5, repr.plot.height=4)
barplot(c(`sharp/ITT`=sharp$coef[1], `fuzzy RD`=fuzzy$coef[1]), col=c("#a0aec0","#2f855a"),
        ylab="estimated effect", main="Fuzzy RD recovers the effect; the sharp jump does not")
abline(h=TRUE_EFF, col="#c53030", lwd=2, lty=2); legend("topleft", sprintf("true effect %d",TRUE_EFF), col="#c53030", lwd=2, lty=2, bty="n")

True effect                      = 3

Sharp (ITT) RD jump in outcome   = 1.709   (biased toward 0)

Fuzzy RD (rdrobust, 2SLS at cutoff) = 2.999   (recovers the truth)
No description has been provided for this image

5. Summary¶

The rdrobust suite reproduced the Lee (2008) analysis end to end: rdplot drew the discontinuity, rdrobust estimated the incumbency advantage at ≈ 0.06 (a ~6-point jump in next-election vote share) with the CCT bandwidth and robust bias-corrected inference — matching the from-scratch Python estimate — rddensity confirmed no manipulation, and the fuzzy = argument recovered a known effect through the RD-as-IV design.

This is the arc's cleanest illustration that the reference implementation and the from-scratch estimator agree, because RD's estimator is transparent: a local line on each side of a cutoff. The design's strength is the mildness of its assumption (continuity, not unconfoundedness or an excluded instrument); its limitation is that the estimate is local to the threshold — the incumbency advantage for near-tied races need not equal it for landslides. Next: Panel Data & Fixed Effects, where plm and fixest difference away time-invariant confounders across repeated observations.