Existing functionality — bandwidth selection, covariates, kernels,
fuzzy and kink designs — is unchanged. See ?rdrobust and
?rdbwselect for full documentation.
1. The data = argument
Previously, callers either had to attach() the data
frame or pass each column with df$ prefixing. In 4.0.0,
rdrobust(), rdbwselect(), and
rdplot() accept a data = argument so
y, x, covs, cluster,
fuzzy, weights, and subset can be
given as bare names, resolved against the supplied data frame:
r <- rdrobust(y = vote, x = margin, data = df)
round(r$coef, 3)
#> Coeff
#> Conventional 7.414
#> Bias-Corrected 7.507
#> Robust 7.507
It also works for expressions built from columns
(e.g. cbind() for multiple covariates, comparisons for
subset):
r_cov <- rdrobust(y = vote, x = margin,
covs = cbind(class, termshouse, termssenate),
data = df)
round(r_cov$coef, 3)
#> Coeff
#> Conventional 6.850
#> Bias-Corrected 6.988
#> Robust 6.988
# subset via an expression evaluated inside `data`
r_sub <- rdrobust(y = vote, x = margin,
subset = year >= 1950, data = df)
round(r_sub$coef, 3)
#> Coeff
#> Conventional 9.223
#> Bias-Corrected 9.295
#> Robust 9.295
The same pattern works for rdbwselect() and
rdplot(). For the rest of the vignette we pass variables
through data =.
2. Heteroskedasticity-robust variance (HC family)
The full HC0–HC3 family is available via
vce =. HC0 is the plug-in sandwich;
HC1 adds a finite-sample scalar; HC2 and
HC3 reweight residuals by leverage, with HC3
being the most conservative.
hc_fits <- lapply(c("nn", "hc0", "hc1", "hc2", "hc3"), function(v) {
rdrobust(y = vote, x = margin, vce = v, data = df)
})
names(hc_fits) <- c("NN", "HC0", "HC1", "HC2", "HC3")
hc_tbl <- do.call(rbind, lapply(names(hc_fits), function(nm) {
r <- hc_fits[[nm]]
data.frame(vce = nm,
coef = round(r$coef[2], 3), # bias-corrected
se_rb = round(r$se[3], 3), # robust SE
pval = round(r$pv[3], 3))
}))
print(hc_tbl, row.names = FALSE)
#> vce coef se_rb pval
#> NN 7.507 1.741 0
#> HC0 7.506 1.739 0
#> HC1 7.504 1.744 0
#> HC2 7.502 1.747 0
#> HC3 7.497 1.755 0
As expected, the point estimate is identical across HC variants (the
variance estimator does not change it); only the standard errors and
p-values move.
3. Cluster-Robust Variance Estimators (CR family)
When observations are grouped into clusters, standard errors should
account for within-cluster correlation. Version 4.0.0 adds three
cluster-robust variance estimators:
"cr1" |
CR1 — cluster sandwich with degrees-of-freedom correction; the
default when cluster is supplied |
"cr2" |
CR2 — Bell & McCaffrey (2002); accounts for within-cluster
leverage |
"cr3" |
CR3 — Pustejovsky & Tipton (2018); leave-one-cluster-out
jackknife; most conservative |
cat("Number of clusters:", length(unique(df$state)), "\n")
#> Number of clusters: 50
r_cr1 <- rdrobust(y = vote, x = margin, cluster = state, vce = "cr1", data = df)
r_cr2 <- rdrobust(y = vote, x = margin, cluster = state, vce = "cr2", data = df)
r_cr3 <- rdrobust(y = vote, x = margin, cluster = state, vce = "cr3", data = df)
cl_tbl <- data.frame(
vce = c("CR1", "CR2", "CR3"),
coef = c(r_cr1$coef[2], r_cr2$coef[2], r_cr3$coef[2]),
se_rb = c(r_cr1$se[3], r_cr2$se[3], r_cr3$se[3]),
pval = c(r_cr1$pv[3], r_cr2$pv[3], r_cr3$pv[3])
)
cl_tbl[, 2:4] <- round(cl_tbl[, 2:4], 3)
print(cl_tbl, row.names = FALSE)
#> vce coef se_rb pval
#> CR1 7.505 1.796 0
#> CR2 7.501 1.800 0
#> CR3 7.494 1.826 0
The point estimate (bias-corrected) is the same across estimators —
only the standard errors differ. CR3 produces the widest intervals and
is the safest default when the number of clusters is moderate.
Automatic validation
When cluster is supplied, only
cr1/cr2/cr3 are valid; the HC
variants and nn are silently remapped (or warned and
remapped, for HC) to the corresponding cluster-robust estimator:
vce = "nn" + cluster → cr1 (silent
default)vce = "hc0" + cluster → cr1 (with
warning)vce = "hc1" + cluster → cr1 (with
warning)vce = "hc2" + cluster → cr2 (with
warning)vce = "hc3" + cluster → cr3 (with
warning)
Conversely, asking for cr1/cr2/cr3 without supplying a
cluster falls back to the corresponding HC variant with a
warning.
Full output for CR3
summary(r_cr3)
#> Call: rdrobust
#>
#> Sharp RD estimates using local polynomial regression. Std. errors are clustered (100 clusters).
#>
#> Number of Obs. 1297
#> BW type mserd
#> Kernel Triangular
#> VCE method CR3
#>
#> Left Right
#> Number of Obs. 595 702
#> Eff. Number of Obs. 369 328
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 18.259 18.259
#> BW bias (b) 27.844 27.844
#> rho (h/b) 0.656 0.656
#> Clusters (g) 50 50
#> Unique Obs. 595 665
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 7.386 4.105 0.000 [3.916 , 11.072]
#> =====================================================================
4. plot.rdrobust: Visual Summary of RD Results
Version 4.0.0 adds a plot S3 method for
rdrobust objects. It produces a binscatter with polynomial
fits and confidence bands on a common scale.
Signature:
plot(x, y, x_run, nbins = 20, ci = TRUE, show_effect = FALSE,
title = NULL, x.label = "Running Variable", y.label = "Outcome", ...)
Basic plot
plot(r_cr1, df$vote, df$margin,
title = "Senate: incumbency advantage",
x.label = "Vote margin (t)",
y.label = "Vote share (t+2)")
Binscatter with local polynomial fits and 95% CI
bands.
Adding an effect panel
show_effect = TRUE appends a second panel with the
conventional point estimate, the robust bias-corrected confidence
interval, and significance stars.
plot(r_cr1, df$vote, df$margin,
show_effect = TRUE,
title = "Senate: incumbency advantage",
x.label = "Vote margin (t)",
y.label = "Vote share (t+2)")
Main plot with RD effect panel.
Clustered plot
The method works identically for any rdrobust object,
regardless of how it was estimated.
plot(r_cr3, df$vote, df$margin,
show_effect = TRUE,
title = "Senate: CR3 clustered SE",
x.label = "Vote margin (t)",
y.label = "Vote share (t+2)")
CR3 clustered variance, with effect panel.
5. Custom RD plots from rdplot()
rdplot() now accepts a hide = TRUE flag and
returns the bin-level data used to draw the plot, so you can build any
visualization on top of ggplot2, lattice, or
base graphics without re-implementing the binning logic.
The returned list contains:
vars_bins — one row per bin: bin midpoint
(rdplot_mean_bin), within-bin mean of y
(rdplot_mean_y), within-bin standard error
(rdplot_se_y), and pointwise CI endpoints
(rdplot_ci_l / rdplot_ci_r) when
ci = is set.vars_poly — a fine evaluation grid
(rdplot_x) with the fitted polynomial values on each side
of the cutoff (rdplot_y).coef — the per-side polynomial coefficient matrix
(Left / Right columns) that produced
vars_poly.
out <- rdplot(y = vote, x = margin, c = 0,
hide = TRUE, ci = 95, data = df)
str(out$vars_bins, give.attr = FALSE)
#> 'data.frame': 49 obs. of 9 variables:
#> $ rdplot_mean_bin: num -96.7 -90 -83.3 -76.7 -70 ...
#> $ rdplot_mean_x : num -98.4 -87.7 -83.2 -79.3 -67.9 ...
#> $ rdplot_mean_y : num 25.4 57 42.3 43.5 43.9 ...
#> $ rdplot_min_bin : num -100 -93.3 -86.7 -80 -73.3 ...
#> $ rdplot_max_bin : num -93.3 -86.7 -80 -73.3 -66.7 ...
#> $ rdplot_se_y : num 10.12 0 6.77 0 3.86 ...
#> $ rdplot_N : int 4 1 5 1 2 4 9 11 30 30 ...
#> $ rdplot_ci_l : num -6.76 57.01 23.54 43.5 -5.07 ...
#> $ rdplot_ci_r : num 57.6 57 61.1 43.5 92.9 ...
str(out$vars_poly, give.attr = FALSE)
#> 'data.frame': 1000 obs. of 2 variables:
#> $ rdplot_x: num -100 -99.8 -99.6 -99.4 -99.2 ...
#> $ rdplot_y: num 29.8 30 30.2 30.5 30.7 ...
out$coef
#> Left Right
#> [1,] 43.9372949735 5.334437e+01
#> [2,] -0.3118100526 5.642182e-02
#> [3,] -0.0371985718 8.489156e-03
#> [4,] -0.0007184850 -5.409779e-05
#> [5,] -0.0000039186 -4.578999e-09
A minimal custom plot using ggplot2:
library(ggplot2)
poly_l <- subset(out$vars_poly, rdplot_x < 0)
poly_r <- subset(out$vars_poly, rdplot_x >= 0)
ggplot() +
geom_vline(xintercept = 0, linewidth = 0.4) +
geom_errorbar(data = out$vars_bins,
aes(x = rdplot_mean_bin,
ymin = rdplot_ci_l, ymax = rdplot_ci_r),
color = "grey60", width = 0) +
geom_point(data = out$vars_bins,
aes(x = rdplot_mean_bin, y = rdplot_mean_y),
color = "grey40", size = 1.5) +
geom_line(data = poly_l, aes(x = rdplot_x, y = rdplot_y)) +
geom_line(data = poly_r, aes(x = rdplot_x, y = rdplot_y)) +
theme_bw() +
labs(x = "Vote margin (t)", y = "Vote share (t+2)",
title = "Senate: regression function fit")
Custom RD plot built from rdplot() output.
You can also reconstruct the fitted polynomial on a custom x-range
directly from coef, which is convenient when the default
evaluation grid doesn’t span the range you want:
eval_fit <- function(xs, coef, c)
sapply(xs, function(xi) sum(coef * (xi - c) ^ (0:(length(coef) - 1))))
xs <- seq(-50, 50, length.out = 5)
data.frame(
x = xs,
y_left = ifelse(xs < 0, eval_fit(xs, out$coef[, "Left"], 0), NA),
y_right = ifelse(xs >= 0, eval_fit(xs, out$coef[, "Right"], 0), NA)
)
#> x y_left y_right
#> 1 -50 31.85074 NA
#> 2 -25 38.17906 NA
#> 3 0 NA 53.34437
#> 4 25 NA 59.21357
#> 5 50 NA 70.59751
The full set of worked examples (default plot, error-bar variant,
shaded CI ribbons, coefficient-based reconstruction) is in
rdplot_illustration.R at the project root, with Python and
Stata twins (rdplot_illustration.py,
rdplot_illustration.do).
6. broom Integration
rdrobust 4.0.0 registers tidy and
glance methods so that results integrate seamlessly with
the broom package and tidyverse workflows.
tidy(): Coefficient-level output
tidy() returns one row per estimator (Conventional,
Bias-Corrected, Robust), with standard broom column names:
library(broom)
tidy(r_cr3)
#> term estimate std.error statistic p.value conf.low conf.high
#> 1 Conventional 7.386333 1.559265 4.737061 2.168401e-06 4.330230 10.44244
#> 2 Bias-Corrected 7.494018 1.559265 4.806122 1.538859e-06 4.437915 10.55012
#> 3 Robust 7.494018 1.825625 4.104905 4.044813e-05 3.915858 11.07218
For most applications, the Robust row (row 3) is the
recommended inference — it combines the bias-corrected point estimate
with robust standard errors.
glance(): Model-level summaries
glance() returns a single-row data frame with key
estimation details: sample sizes, bandwidths, VCE method, and whether
covariates or clusters were used.
glance(r_cr3)
#> nobs.left nobs.right nobs.effective.left nobs.effective.right h.left h.right
#> 1 595 702 369 328 18.2592 18.2592
#> b.left b.right cutoff order.regression order.bias kernel vce bwselect
#> 1 27.84356 27.84356 0 1 2 Triangular CR3 mserd
#> covariates clusters.left clusters.right
#> 1 FALSE 50 50
Collecting results across specifications
The main payoff of broom integration is easy comparison across
specifications. Here we compare all three cluster-robust estimators:
specs <- list(CR1 = r_cr1, CR2 = r_cr2, CR3 = r_cr3)
# Robust BC estimates side-by-side
res <- do.call(rbind, lapply(names(specs), function(nm) {
cbind(spec = nm, tidy(specs[[nm]])[3, ]) # row 3 = robust
}))
res[, c("spec", "estimate", "std.error", "conf.low", "conf.high")]
#> spec estimate std.error conf.low conf.high
#> 3 CR1 7.504956 1.795945 3.984970 11.02494
#> 31 CR2 7.501257 1.800322 3.972691 11.02982
#> 32 CR3 7.494018 1.825625 3.915858 11.07218
We can also compare model-level details:
gl <- do.call(rbind, lapply(names(specs), function(nm) {
cbind(spec = nm, glance(specs[[nm]]))
}))
gl[, c("spec", "nobs.effective.left", "nobs.effective.right",
"h.left", "vce", "clusters.left", "clusters.right")]
#> spec nobs.effective.left nobs.effective.right h.left vce clusters.left
#> 1 CR1 366 325 18.08429 CR1 50
#> 2 CR2 367 326 18.14594 CR2 50
#> 3 CR3 369 328 18.25920 CR3 50
#> clusters.right
#> 1 50
#> 2 50
#> 3 50
Broader comparison: HC3 vs clustered
A common workflow is comparing heteroskedasticity-robust and
cluster-robust results to assess sensitivity to within-cluster
correlation:
r_hc3 <- rdrobust(y = vote, x = margin, vce = "hc3", data = df)
all_specs <- list("HC3" = r_hc3,
"CR1" = r_cr1,
"CR2" = r_cr2,
"CR3" = r_cr3)
comparison <- do.call(rbind, lapply(names(all_specs), function(nm) {
g <- glance(all_specs[[nm]])
t <- tidy(all_specs[[nm]])[3, ] # robust row
data.frame(spec = nm,
estimate = round(t$estimate, 3),
std.error = round(t$std.error, 3),
ci = paste0("[", round(t$conf.low, 3), ", ", round(t$conf.high, 3), "]"),
h = round(g$h.left, 3),
clusters = g$clusters.left)
}))
print(comparison, row.names = FALSE)
#> spec estimate std.error ci h clusters
#> HC3 7.497 1.755 [4.058, 10.936] 17.766 NA
#> CR1 7.505 1.796 [3.985, 11.025] 18.084 50
#> CR2 7.501 1.800 [3.973, 11.03] 18.146 50
#> CR3 7.494 1.826 [3.916, 11.072] 18.259 50
7. Reporting tables: conventional point + RBC CI
A common workflow is to present a small panel of specifications
(sharp, covariate-adjusted, cluster-robust, etc.) in a single table. The
canonical Calonico-Cattaneo-Titiunik presentation pairs the
conventional point estimate with the robust
bias-corrected (RBC) confidence interval — the point estimate
is from row 1 of coef, the CI is row 3 of ci.
The same pattern appears in the Stata illustration
(rdrobust_illustration_new_features.do), built there with
esttab; the R analog uses knitr::kable().
We fit four specifications side-by-side: sharp, sharp+covs, cluster,
cluster+covs.
covs_form <- ~ class + termshouse + termssenate
m1 <- rdrobust(y = vote, x = margin, data = df)
m2 <- rdrobust(y = vote, x = margin, covs = covs_form, data = df)
m3 <- rdrobust(y = vote, x = margin, vce = "cr1", cluster = state, data = df)
m4 <- rdrobust(y = vote, x = margin, vce = "cr1", cluster = state,
covs = covs_form, data = df)
A small helper extracts the conventional point estimate, the RBC
confidence interval, and a few model-level stats from a fitted
object:
rdr_row <- function(fit, label) {
data.frame(
spec = label,
estimate = fit$coef[1, 1], # tau_conv
rbc_ci = sprintf("[%.3f, %.3f]", fit$ci[3, 1], fit$ci[3, 2]),
rbc_pval = fit$pv[3, 1],
h = fit$bws[1, 1], # h_left
N_eff = fit$N_h[1] + fit$N_h[2],
clusters = if (!is.null(fit$n_clust)) sum(fit$n_clust) else NA_integer_
)
}
tab <- rbind(
rdr_row(m1, "sharp"),
rdr_row(m2, "sharp + covs"),
rdr_row(m3, "cluster"),
rdr_row(m4, "cluster + covs")
)
Render with knitr::kable() for the report:
knitr::kable(
tab,
digits = 3,
col.names = c("Specification", "Estimate", "Robust 95% CI",
"RBC p-value", "h", "Eff. N", "Clusters"),
caption = "RD estimates: conventional point estimate with RBC confidence interval."
)
RD estimates: conventional point estimate with RBC confidence
interval.
| sharp |
7.414 |
[4.094, 10.919] |
0 |
17.754 |
683 |
NA |
| sharp + covs |
6.850 |
[3.728, 10.249] |
0 |
18.033 |
598 |
NA |
| cluster |
7.395 |
[3.985, 11.025] |
0 |
18.084 |
691 |
100 |
| cluster + covs |
6.832 |
[3.726, 10.217] |
0 |
18.431 |
612 |
98 |
The Estimate column reports tau_conv (the standard
local-polynomial point estimate); the Robust 95% CI is
the bias-corrected, robust-SE confidence interval — the inference
recommended by CCT 2014 and the same one printed by
summary(fit) in the Robust row. The CI is centered on
tau_bc, not on tau_conv, which is why a CI may
not be symmetric about the reported point estimate. This is a feature of
robust bias-corrected inference, not a typo.
Specs-as-columns variant (paper style)
For a typeset table closer to the Stata esttab layout
(specifications as columns, statistics as rows), transpose:
paper <- t(tab[, c("estimate", "rbc_ci", "h", "N_eff")])
colnames(paper) <- tab$spec
rownames(paper) <- c("RD effect (conv.)", "Robust 95% CI", "h", "Eff. N")
knitr::kable(
paper,
caption = "Senate incumbency: RD estimates across specifications."
)
Senate incumbency: RD estimates across
specifications.
| RD effect (conv.) |
7.414131 |
6.849875 |
7.394838 |
6.832365 |
| Robust 95% CI |
[4.094, 10.919] |
[3.728, 10.249] |
[3.985, 11.025] |
[3.726, 10.217] |
| h |
17.75440 |
18.03336 |
18.08429 |
18.43127 |
| Eff. N |
683 |
598 |
691 |
612 |
LaTeX export
To write a publication-ready LaTeX table to a file (the R analog of
the Stata esttab ... using rdrobust_table.tex):
tex <- knitr::kable(
paper,
format = "latex",
booktabs = TRUE,
caption = "RD estimates for Senate incumbency.",
label = "tab:rdrobust"
)
writeLines(tex, "rdrobust_table.tex")
For richer formatting (multi-row headers, footnotes, grouping), the
kableExtra and gt packages build on top of
kable() / gt::gt() and accept the same input
data frame.
