The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
rdhte: Empirical Illustration
Overview
rdhte estimates heterogeneous treatment
effects in sharp Regression Discontinuity (RD) designs along
one or more pretreatment covariates. The package exposes three main
commands:
rdhte() – point estimation and robust bias-corrected
inference of group-conditional RD effects.rdbwhte() – automatic per-cell bandwidth
selection.rdhte_lincom() – linear-combination tests over the
estimated CATEs.
This vignette illustrates each feature on a real-data extract from
Granzier, Pons, and Tricaud (2023, American Economic Journal:
Applied), bundled with the package as
rdhte_dataset.
The data and the design
Granzier, Pons, and Tricaud study coordination and bandwagon
effects in French two-round elections, where candidates must
clear a qualifying-vote threshold in the first round to advance to the
runoff. The institutional rule creates a sharp regression-discontinuity
design on every candidate’s first-round margin against that threshold:
candidates just above the cutoff advance, those just below do not, and
small differences in first-round support produce a discrete change in
runoff participation. The authors use the design to ask whether being
just-qualified causally affects subsequent voter and candidate behavior,
and document substantive heterogeneity across ideology and
candidate-strength dimensions.
The bundled extract has 39,534 candidate-race observations with the
following variables:
y – 0/1 indicator for advancing to the runoff (the
outcome).x – first-round margin against the qualifying threshold
(the running variable; cutoff at zero).cluster_var – district identifier, used for
cluster-robust inference.w_left – 0/1 indicator for left-of-center parties.w_ideology – unordered party-ideology bucket (4
categories).w_strength – continuous proxy for ex-ante candidate
strength (average prior national-level vote share).w_strong – 0/1 indicator for above-median
strength.w_strength_qrt – ordered quartile of
w_strength (4 levels).
These covariates support every covariate-incorporation pattern
rdhte exposes: binary cells, multi-level (unordered and
ordered) factor cells, factor-by-factor interactions, continuous slopes,
and binary x continuous interactions.
library(rdhte)
library(rdrobust)
library(broom)
data(rdhte_dataset)
attach(rdhte_dataset)
Single binary variable
When covs.hte is binary (or a factor with two levels),
rdhte reports one CATE per cell.
summary(rd_left <- rdhte(y = y, x = x, covs.hte = factor(w_left),
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Subgroups., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#>
#> ==================================================================================================================
#> Point Robust Inference
#> factor(w_left) Estimate z Pr(>|z|) [ 95% C.I. ] Nh- Nh+ h- h+
#> ------------------------------------------------------------------------------------------------------------------
#> 0 0.021 1.523 0.128 [-0.003 , 0.027] 4352 4422 0.077 0.077
#> 1 0.089 6.821 0.000 [0.062 , 0.112] 5036 4840 0.117 0.117
#> ==================================================================================================================
How large is the advantage for left-of-center candidates? Test the
difference with rdhte_lincom:
rdhte_lincom(model = rd_left,
linfct = c("`factor(w_left)1` - `factor(w_left)0` = 0"))
#> $individual
#> hypothesis estimate z_stat p_value conf.low
#> 1 `factor(w_left)1` - `factor(w_left)0` 0.068 5.056 0 0.046
#> conf.high
#> 1 0.105
#>
#> $joint
#> statistic X1L p_value
#> 1 25.563 1 0
Forcing a common bandwidth across cells
(bw.joint = TRUE) can distort the cell-specific inference –
here it makes the right-of-center effect (incorrectly) statistically
significant:
summary(rdhte(y = y, x = x, covs.hte = w_left, cluster = cluster_var,
bw.joint = TRUE))
#> Sharp RD Heterogeneous Treatment Effects: Subgroups., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#>
#> ================================================================================================================
#> Point Robust Inference
#> w_left Estimate z Pr(>|z|) [ 95% C.I. ] Nh- Nh+ h- h+
#> ----------------------------------------------------------------------------------------------------------------
#> 0 0.023 2.272 0.023 [0.002 , 0.029] 5208 5364 0.097 0.097
#> 1 0.088 6.437 0.000 [0.062 , 0.116] 4325 4183 0.097 0.097
#> ================================================================================================================
A 0/1 numeric variable is auto-promoted to a factor; coefficient
names then carry the variable name as a prefix:
summary(rd_left2 <- rdhte(y = y, x = x, covs.hte = w_left,
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Subgroups., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#>
#> ================================================================================================================
#> Point Robust Inference
#> w_left Estimate z Pr(>|z|) [ 95% C.I. ] Nh- Nh+ h- h+
#> ----------------------------------------------------------------------------------------------------------------
#> 0 0.021 1.523 0.128 [-0.003 , 0.027] 4352 4422 0.077 0.077
#> 1 0.089 6.821 0.000 [0.062 , 0.112] 5036 4840 0.117 0.117
#> ================================================================================================================
rdhte_lincom(model = rd_left2,
linfct = c("`w_left1` - `w_left0` = 0"))
#> $individual
#> hypothesis estimate z_stat p_value conf.low conf.high
#> 1 w_left1 - w_left0 0.068 5.056 0 0.046 0.105
#>
#> $joint
#> statistic X1L p_value
#> 1 25.563 1 0
Single categorical variable – unordered
A multi-level factor produces one CATE per level. The reference
category is the factor’s first level.
summary(rd_ideology <- rdhte(y = y, x = x,
covs.hte = factor(w_ideology),
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Subgroups., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#>
#> ======================================================================================================================
#> Point Robust Inference
#> factor(w_ideology) Estimate z Pr(>|z|) [ 95% C.I. ] Nh- Nh+ h- h+
#> ----------------------------------------------------------------------------------------------------------------------
#> 1 0.089 6.819 0.000 [0.062 , 0.112] 5036 4840 0.117 0.117
#> 2 -0.024 -0.853 0.393 [-0.195 , 0.077] 181 156 0.056 0.056
#> 3 0.026 1.897 0.058 [-0.001 , 0.032] 3816 4202 0.086 0.086
#> 4 0.007 0.884 0.377 [-0.012 , 0.032] 713 384 0.090 0.090
#> ======================================================================================================================
Joint test that the three non-reference ideology cells are all
indistinguishable from zero:
rdhte_lincom(model = rd_ideology,
linfct = c("`factor(w_ideology)2` = 0",
"`factor(w_ideology)3` = 0",
"`factor(w_ideology)4` = 0"))
#> $individual
#> hypothesis estimate z_stat p_value conf.low conf.high
#> 1 factor(w_ideology)2 -0.024 -0.853 0.393 -0.195 0.077
#> 2 factor(w_ideology)3 0.026 1.897 0.058 -0.001 0.032
#> 3 factor(w_ideology)4 0.007 0.884 0.377 -0.012 0.032
#>
#> $joint
#> statistic X3L p_value
#> 1 5.116 3 0.163
Single categorical variable – ordered
For an ordered categorical variable, wrapping in
factor() keeps the per-level CATE interpretation:
summary(rdhte(y = y, x = x, covs.hte = factor(w_strength_qrt),
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Subgroups., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#>
#> ==========================================================================================================================
#> Point Robust Inference
#> factor(w_strength_qrt) Estimate z Pr(>|z|) [ 95% C.I. ] Nh- Nh+ h- h+
#> --------------------------------------------------------------------------------------------------------------------------
#> 0 0.016 1.803 0.071 [-0.001 , 0.032] 2583 2283 0.094 0.094
#> 1 0.039 2.219 0.026 [0.003 , 0.055] 2227 2300 0.085 0.085
#> 2 0.054 3.490 0.000 [0.025 , 0.090] 2040 2329 0.105 0.105
#> 3 0.094 6.652 0.000 [0.069 , 0.127] 3436 3486 0.143 0.143
#> ==========================================================================================================================
The per-quartile CATE increases monotonically: candidates with higher
average strength have larger treatment effects.
A bare numeric variable (no factor()) is treated as
continuous – useful when the ordering is meaningful but you want the
linear-in-W parametrization:
summary(rdhte(y = y, x = x, covs.hte = w_strength_qrt,
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Continuous., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Eff. Number of Obs. 9533 9547
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.097 0.097
#>
#> ==========================================================================
#> Point Robust Inference
#> Estimate z Pr(>|z|) [ 95% C.I. ]
#> --------------------------------------------------------------------------
#> T 0.014 1.389 0.165 [-0.005 , 0.029]
#> T:w_strength_qrt 0.026 3.997 0.000 [0.013 , 0.037]
#> ==========================================================================
Two binary variables – interaction
Interactions of two factors produce one CATE per joint cell. The
ordering of cells follows R’s interaction() convention.
summary(rdhte(y = y, x = x,
covs.hte = factor(w_left):factor(w_strong),
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Subgroups., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#>
#> ===================================================================================================================================
#> Point Robust Inference
#> factor(w_left):factor(w_strong) Estimate z Pr(>|z|) [ 95% C.I. ] Nh- Nh+ h- h+
#> -----------------------------------------------------------------------------------------------------------------------------------
#> 0:1 0.003 0.004 0.997 [-0.016 , 0.016] 2348 2257 0.070 0.070
#> 0:2 0.044 3.014 0.003 [0.014 , 0.064] 2316 2643 0.104 0.104
#> 1:1 0.061 3.838 0.000 [0.029 , 0.090] 2066 1960 0.097 0.097
#> 1:2 0.114 6.462 0.000 [0.082 , 0.153] 3107 3055 0.143 0.143
#> ===================================================================================================================================
The same model can be expressed by pre-building the joint cell
variable:
interactions <- 1*(w_left==0)*(w_strong==1) +
2*(w_left==0)*(w_strong==2) +
3*(w_left==1)*(w_strong==1) +
4*(w_left==1)*(w_strong==2)
summary(rdhte(y = y, x = x, covs.hte = factor(interactions),
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Subgroups., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#>
#> ========================================================================================================================
#> Point Robust Inference
#> factor(interactions) Estimate z Pr(>|z|) [ 95% C.I. ] Nh- Nh+ h- h+
#> ------------------------------------------------------------------------------------------------------------------------
#> 1 0.003 0.004 0.997 [-0.016 , 0.016] 2348 2257 0.070 0.070
#> 2 0.044 3.014 0.003 [0.014 , 0.064] 2316 2643 0.104 0.104
#> 3 0.061 3.838 0.000 [0.029 , 0.090] 2066 1960 0.097 0.097
#> 4 0.114 6.462 0.000 [0.082 , 0.153] 3107 3055 0.143 0.143
#> ========================================================================================================================
Single continuous variable
When covs.hte is continuous, rdhte switches
to a linear-in-W parametrization of the CATE function:
\[
\tau(w) = \beta_T + \beta_{T:W}\,w.
\]
summary(rd_continuous <- rdhte(y = y, x = x, covs.hte = w_strength,
kernel = "uni",
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Continuous., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Uniform
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Eff. Number of Obs. 9514 9529
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.096 0.096
#>
#> ========================================================================
#> Point Robust Inference
#> Estimate z Pr(>|z|) [ 95% C.I. ]
#> ------------------------------------------------------------------------
#> T -0.055 -2.507 0.012 [-0.129 , -0.016]
#> T:w_strength 0.262 3.927 0.000 [0.151 , 0.451]
#> ========================================================================
To aid interpretation, the slope coefficient is precisely a
partial-linear regression slope. With the uniform
kernel and a fixed bandwidth, the rdhte point
estimates match plain lm() on the same in-bandwidth
subset:
trt <- (x > 0)
new.data <- data.frame(y, x, w_strength, trt)
using.lm <- coef(lm(y ~ trt * x * w_strength,
data = new.data[abs(new.data$x) < rd_continuous$h[1], ]))
rd_continuous$Estimate[1] - using.lm["trtTRUE"]
#> T
#> 0
rd_continuous$Estimate[2] - using.lm["trtTRUE:w_strength"]
#> T:w_strength
#> 0
(Inference, however, requires the robust bias correction in
rdhte and cannot be obtained from the lm() fit
alone.)
Interaction effect: binary x continuous
A fully interacted specification gives a separate intercept and slope
for each level of the categorical covariate.
summary(rd_interaction <- rdhte(y = y, x = x,
covs.hte = "w_left*w_strength",
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Continuous., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Eff. Number of Obs. 9533 9547
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.097 0.097
#>
#> ===============================================================================
#> Point Robust Inference
#> Estimate z Pr(>|z|) [ 95% C.I. ]
#> -------------------------------------------------------------------------------
#> T -0.033 -1.364 0.173 [-0.094 , 0.017]
#> T:w_left -0.059 -0.275 0.783 [-0.306 , 0.230]
#> T:w_strength 0.141 1.779 0.075 [-0.014 , 0.295]
#> T:w_left:w_strength 0.274 0.742 0.458 [-0.395 , 0.876]
#> ===============================================================================
Wrapping the binary in factor() is also allowed but
shifts the baseline category and therefore the reported intercepts:
summary(rdhte(y = y, x = x, covs.hte = "factor(w_left)*w_strength",
cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Continuous., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Eff. Number of Obs. 9533 9547
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.097 0.097
#>
#> =======================================================================================
#> Point Robust Inference
#> Estimate z Pr(>|z|) [ 95% C.I. ]
#> ---------------------------------------------------------------------------------------
#> T -0.092 -0.571 0.568 [-0.338 , 0.186]
#> T:factor(w_left)0 0.059 0.275 0.783 [-0.230 , 0.306]
#> T:w_strength 0.141 1.779 0.075 [-0.014 , 0.295]
#> T:factor(w_left)1:w_strength 0.274 0.742 0.458 [-0.395 , 0.876]
#> =======================================================================================
Each cell-specific coefficient may look insignificant individually
while the joint test still rejects:
rdhte_lincom(model = rd_interaction,
linfct = c("`T` = 0",
"`T:w_left` = 0",
"`T:w_strength` = 0",
"`T:w_left:w_strength` = 0"))
#> $individual
#> hypothesis estimate z_stat p_value conf.low conf.high
#> 1 T -0.033 -1.364 0.173 -0.094 0.017
#> 2 T:w_left -0.059 -0.275 0.783 -0.306 0.230
#> 3 T:w_strength 0.141 1.779 0.075 -0.014 0.295
#> 4 T:w_left:w_strength 0.274 0.742 0.458 -0.395 0.876
#>
#> $joint
#> statistic X4L p_value
#> 1 47.392 4 0
With the bandwidth fixed, the fully interacted model matches results
from category-specific estimation. Use rdhte_lincom to read
off the per-cell intercept and slope:
summary(rd_interaction <- rdhte(y = y, x = x,
covs.hte = "w_left*w_strength",
h = 0.1, cluster = cluster_var))
#> Sharp RD Heterogeneous Treatment Effects: Continuous., Cluster-Adjusted
#>
#> Number of Obs. 39534
#> BW type Manual
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 19723 19811
#> Eff. Number of Obs. 9829 9843
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.100 0.100
#>
#> ===============================================================================
#> Point Robust Inference
#> Estimate z Pr(>|z|) [ 95% C.I. ]
#> -------------------------------------------------------------------------------
#> T -0.031 -1.434 0.152 [-0.095 , 0.015]
#> T:w_left -0.061 -0.266 0.791 [-0.303 , 0.231]
#> T:w_strength 0.138 1.875 0.061 [-0.007 , 0.297]
#> T:w_left:w_strength 0.281 0.727 0.467 [-0.397 , 0.866]
#> ===============================================================================
summary(rdhte(y = y[w_left == 0], x = x[w_left == 0],
covs.hte = w_strength[w_left == 0], h = 0.1,
cluster = cluster_var[w_left == 0]))
#> Sharp RD Heterogeneous Treatment Effects: Continuous., Cluster-Adjusted
#>
#> Number of Obs. 21938
#> BW type Manual
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 10287 11651
#> Eff. Number of Obs. 5371 5530
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.100 0.100
#>
#> ===================================================================================
#> Point Robust Inference
#> Estimate z Pr(>|z|) [ 95% C.I. ]
#> -----------------------------------------------------------------------------------
#> T -0.031 -1.434 0.152 [-0.095 , 0.015]
#> T:w_strength[w_left == 0] 0.138 1.875 0.061 [-0.007 , 0.297]
#> ===================================================================================
summary(rdhte(y = y[w_left == 1], x = x[w_left == 1],
covs.hte = w_strength[w_left == 1], h = 0.1,
cluster = cluster_var[w_left == 1]))
#> Sharp RD Heterogeneous Treatment Effects: Continuous., Cluster-Adjusted
#>
#> Number of Obs. 17596
#> BW type Manual
#> Kernel Triangular
#> VCE method CR1
#>
#> Number of Obs. 9436 8160
#> Eff. Number of Obs. 4458 4313
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.100 0.100
#>
#> ===================================================================================
#> Point Robust Inference
#> Estimate z Pr(>|z|) [ 95% C.I. ]
#> -----------------------------------------------------------------------------------
#> T -0.093 -0.571 0.568 [-0.337 , 0.185]
#> T:w_strength[w_left == 1] 0.419 1.213 0.225 [-0.233 , 0.992]
#> ===================================================================================
rdhte_lincom(model = rd_interaction,
linfct = c("`T` + `T:w_left` = 0",
"`T:w_strength` + `T:w_left:w_strength` = 0"))
#> $individual
#> hypothesis estimate z_stat p_value conf.low
#> 1 T + `T:w_left` -0.093 -0.571 0.568 -0.337
#> 2 `T:w_strength` + `T:w_left:w_strength` 0.419 1.213 0.225 -0.233
#> conf.high
#> 1 0.185
#> 2 0.992
#>
#> $joint
#> statistic X2L p_value
#> 1 42.13 2 0
Standalone bandwidth selection (rdbwhte)
rdbwhte runs the same data-driven bandwidth selectors as
rdhte but returns the bandwidths without estimating the
CATEs. Useful when you want to inspect or fix h and then
explore alternative variance estimators or compare specifications at a
common bandwidth.
bw <- rdbwhte(y = y, x = x, covs.hte = factor(w_ideology),
cluster = cluster_var)
bw$h
#> [,1] [,2]
#> [1,] 0.11696532 0.11696532
#> [2,] 0.05645238 0.05645238
#> [3,] 0.08564972 0.08564972
#> [4,] 0.09037549 0.09037549
bw.joint = TRUE forces a single shared bandwidth across
cells:
bw_joint <- rdbwhte(y = y, x = x, covs.hte = factor(w_ideology),
cluster = cluster_var, bw.joint = TRUE)
bw_joint$h
#> [,1] [,2]
#> [1,] 0.09650206 0.09650206
#> [2,] 0.09650206 0.09650206
#> [3,] 0.09650206 0.09650206
#> [4,] 0.09650206 0.09650206
Efficiency-improving covariates (covs.eff)
covs.eff adds covariates to the local-polynomial
regression to shrink standard errors without changing
the identification of the CATE. They enter additively (and as
covs.eff x W interactions in the heterogeneity-aware paths)
but never with the treatment indicator or the running-variable
polynomial. Useful when you have strong pretreatment predictors of the
outcome.
m_no_eff <- rdhte(y = y, x = x, covs.hte = factor(w_ideology),
cluster = cluster_var)
m_eff <- rdhte(y = y, x = x, covs.hte = factor(w_ideology),
covs.eff = w_strength, cluster = cluster_var)
data.frame(cell = m_no_eff$W.lev,
se_no_eff = round(m_no_eff$se.rb, 4),
se_eff = round(m_eff$se.rb, 4),
pct_change = round(100 * (m_eff$se.rb - m_no_eff$se.rb)
/ m_no_eff$se.rb, 1))
#> cell se_no_eff se_eff pct_change
#> factor(w_ideology)1 1 0.0128 0.0129 0.6
#> factor(w_ideology)2 2 0.0691 0.0689 -0.3
#> factor(w_ideology)3 3 0.0083 0.0084 0.6
#> factor(w_ideology)4 4 0.0113 0.0117 3.0
Plotting
plot() is a post-estimation method for categorical
covs.hte: one point per cell at the conventional point
estimate, with the robust bias-corrected CI as an error bar.
sort = TRUE reorders cells along the x-axis by their
point estimate – often the more useful default for ranking-style
narratives:
plot(rd_ideology, sort = TRUE)
The returned object is a ggplot, so additional layers
compose naturally. For example, a custom title with horizontal
orientation:
p <- plot(rd_ideology, sort = TRUE,
title = "Heterogeneity by ideology bucket",
ylab = "Sharp RD ITT")
if (requireNamespace("ggplot2", quietly = TRUE)) {
p + ggplot2::theme_minimal(base_size = 11) +
ggplot2::coord_flip()
}
Building publication-ready tables
rdhte exposes the per-cell results through
tidy() and a one-row summary through glance()
(both broom-compatible). Combined with knitr::kable() /
xtable::xtable() the same numbers can be turned into
Markdown, HTML, or LaTeX tables for papers.
(A) per-cell table from a single call
tidy() returns one row per cell with the conventional
point estimate, robust BC standard error, z / p-value, BC confidence
interval, and per-side h and Nh. Renaming or
formatting columns before kable() gives a clean publication
panel:
tab_A <- tidy(rd_ideology)
tab_A
#> term estimate std.error statistic p.value conf.low conf.high
#> 1 1 0.088652653 0.01279375 6.8186555 9.189647e-12 0.0621608620 0.11231142
#> 2 2 -0.024026038 0.06914491 -0.8533574 3.934611e-01 -0.1945268631 0.07651622
#> 3 3 0.025632358 0.00831174 1.8974042 5.777461e-02 -0.0005199808 0.03206144
#> 4 4 0.007008009 0.01133009 0.8840432 3.766729e-01 -0.0121902781 0.03222286
#> estimate.bc h.left h.right n.eff.left n.eff.right
#> 1 0.08723614 0.11696532 0.11696532 5036 4840
#> 2 -0.05900532 0.05645238 0.05645238 181 156
#> 3 0.01577073 0.08564972 0.08564972 3816 4202
#> 4 0.01001629 0.09037549 0.09037549 713 384
knitr::kable(
tab_A[, c("term", "estimate", "std.error", "conf.low", "conf.high",
"n.eff.left", "n.eff.right", "h.left", "h.right")],
digits = c(NA, 3, 3, 3, 3, 0, 0, 3, 3),
col.names = c("Cell", "Estimate", "SE", "CI lo", "CI hi",
"Nh-", "Nh+", "h-", "h+"),
caption = "Per-cell RD effects by ideology bucket"
)
Per-cell RD effects by ideology bucket
| 1 |
0.089 |
0.013 |
0.062 |
0.112 |
5036 |
4840 |
0.117 |
0.117 |
| 2 |
-0.024 |
0.069 |
-0.195 |
0.077 |
181 |
156 |
0.056 |
0.056 |
| 3 |
0.026 |
0.008 |
-0.001 |
0.032 |
3816 |
4202 |
0.086 |
0.086 |
| 4 |
0.007 |
0.011 |
-0.012 |
0.032 |
713 |
384 |
0.090 |
0.090 |
glance() complements tidy() with a one-row
description of the fit (N, polynomial orders, kernel, VCE, bandwidth
selector):
glance(rd_ideology)
#> n n.left n.right cutoff p q kernel vce vce_select bwselect level
#> 1 39534 19723 19811 0 1 2 Triangular CR1 cr1 mserd 95
#> n.terms covs.continuous
#> 1 4 FALSE
(B) cell x specification comparison
Fix the covs.hte argument and vary vce
across columns to see how the choice of variance estimator moves the
per-cell standard errors.
specs <- list(HC0 = list(vce = "hc0"),
HC2 = list(vce = "hc2"),
HC3 = list(vce = "hc3"),
CR1 = list(vce = "cr1", cluster = cluster_var))
fit_one <- function(args) {
args <- c(list(y = y, x = x, covs.hte = factor(w_ideology)), args)
do.call(rdhte, args)
}
cells <- sort(unique(w_ideology))
mat_pt <- sapply(specs, function(s) tidy(fit_one(s))$estimate)
mat_se <- sapply(specs, function(s) tidy(fit_one(s))$std.error)
rownames(mat_pt) <- rownames(mat_se) <- as.character(cells)
mat_pt
#> HC0 HC2 HC3 CR1
#> 1 0.088623908 0.088625020 0.088625979 0.088652653
#> 2 -0.024247190 -0.023878639 -0.023577737 -0.024026038
#> 3 0.025784637 0.025791260 0.025797429 0.025632358
#> 4 0.007022683 0.007006824 0.006991369 0.007008009
mat_se
#> HC0 HC2 HC3 CR1
#> 1 0.011736372 0.01174108 0.011746082 0.01279375
#> 2 0.069234017 0.07084168 0.072492654 0.06914491
#> 3 0.008274419 0.00827868 0.008283159 0.00831174
#> 4 0.011338034 0.01136898 0.011400526 0.01133009
A common publication layout interleaves the point estimate and its SE
in parentheses below it:
n_cells <- nrow(mat_pt)
out <- data.frame(matrix(NA_character_, nrow = 2 * n_cells,
ncol = ncol(mat_pt)))
for (i in seq_len(n_cells)) {
out[2 * i - 1, ] <- sprintf("%.3f", mat_pt[i, ])
out[2 * i, ] <- sprintf("(%.3f)", mat_se[i, ])
}
out <- cbind(Cell = rep(rownames(mat_pt), each = 2), out)
out$Cell[seq(2, nrow(out), by = 2)] <- ""
colnames(out)[-1] <- colnames(mat_pt)
knitr::kable(out, caption = "Per-cell estimates by variance option (SE in parentheses)")
Per-cell estimates by variance option (SE in
parentheses)
| 1 |
0.089 |
0.089 |
0.089 |
0.089 |
|
(0.012) |
(0.012) |
(0.012) |
(0.013) |
| 2 |
-0.024 |
-0.024 |
-0.024 |
-0.024 |
|
(0.069) |
(0.071) |
(0.072) |
(0.069) |
| 3 |
0.026 |
0.026 |
0.026 |
0.026 |
|
(0.008) |
(0.008) |
(0.008) |
(0.008) |
| 4 |
0.007 |
0.007 |
0.007 |
0.007 |
|
(0.011) |
(0.011) |
(0.011) |
(0.011) |
(C) LaTeX export
xtable or kableExtra can render the same
data frame as LaTeX. The snippet below produces a paper-ready table;
uncomment the writeLines() call to write to a file:
tex <- xtable::xtable(
tab_A[, c("term", "estimate", "std.error", "conf.low", "conf.high")],
digits = c(0, 0, 3, 3, 3, 3),
caption = "Per-cell RD effects by ideology bucket",
label = "tab:rdhte-ideology"
)
print(tex, include.rownames = FALSE, booktabs = TRUE,
caption.placement = "top")
#> % latex table generated in R 4.6.0 by xtable 1.8-8 package
#> % Tue May 26 09:50:57 2026
#> \begin{table}[ht]
#> \centering
#> \caption{Per-cell RD effects by ideology bucket}
#> \label{tab:rdhte-ideology}
#> \begin{tabular}{lrrrr}
#> \toprule
#> term & estimate & std.error & conf.low & conf.high \\
#> \midrule
#> 1 & 0.089 & 0.013 & 0.062 & 0.112 \\
#> 2 & -0.024 & 0.069 & -0.195 & 0.077 \\
#> 3 & 0.026 & 0.008 & -0.001 & 0.032 \\
#> 4 & 0.007 & 0.011 & -0.012 & 0.032 \\
#> \bottomrule
#> \end{tabular}
#> \end{table}
# writeLines(capture.output(print(tex, include.rownames = FALSE,
# booktabs = TRUE,
# caption.placement = "top")),
# con = "rdhte_table_A.tex")
Replicating rdrobust
Average effects from rdhte() and rdrobust()
will not match under default settings, because the two packages use
different defaults for rho and vce:
summary(rdhte(y = y, x = x))
#> Sharp RD Average Treatment Effect.
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method HC3
#>
#> Number of Obs. 19723 19811
#> Eff. Number of Obs. 9536 9550
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.097 0.097
#>
#> ========================================================================
#> Point Robust Inference
#> Estimate z Pr(>|z|) [ 95% C.I. ]
#> ------------------------------------------------------------------------
#> T 0.051 7.071 0.000 [0.035 , 0.062]
#> ========================================================================
summary(rdrobust(y = y, x = x))
#> Call: rdrobust
#>
#> Sharp RD estimates using local polynomial regression.
#>
#> Number of Obs. 39534
#> BW type mserd
#> Kernel Triangular
#> VCE method NN
#>
#> Left Right
#> Number of Obs. 19723 19811
#> Eff. Number of Obs. 9546 9560
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.097 0.097
#> BW bias (b) 0.159 0.159
#> rho (h/b) 0.608 0.608
#> Unique Obs. 19660 19746
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 0.051 8.894 0.000 [0.039 , 0.062]
#> =====================================================================
To replicate exactly with rdrobust, set
rho = 1, choose a matching vce, and enforce
the same bandwidth:
summary(rdhte(y = y, x = x, h = 0.1, vce = "hc3"))
#> Sharp RD Average Treatment Effect.
#>
#> Number of Obs. 39534
#> BW type Manual
#> Kernel Triangular
#> VCE method HC3
#>
#> Number of Obs. 19723 19811
#> Eff. Number of Obs. 9829 9843
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.100 0.100
#>
#> ========================================================================
#> Point Robust Inference
#> Estimate z Pr(>|z|) [ 95% C.I. ]
#> ------------------------------------------------------------------------
#> T 0.051 7.194 0.000 [0.035 , 0.062]
#> ========================================================================
summary(rdrobust(y = y, x = x, h = 0.1, rho = 1, vce = "hc3"))
#> Call: rdrobust
#>
#> Sharp RD estimates using local polynomial regression.
#>
#> Number of Obs. 39534
#> BW type Manual
#> Kernel Triangular
#> VCE method HC3
#>
#> Left Right
#> Number of Obs. 19723 19811
#> Eff. Number of Obs. 9829 9843
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.100 0.100
#> BW bias (b) 0.100 0.100
#> rho (h/b) 1.000 1.000
#> Unique Obs. 19723 19811
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 0.051 7.194 0.000 [0.035 , 0.062]
#> =====================================================================
The same trick reproduces the per-cell rdrobust
estimates from a single rdhte call with per-cell
bandwidths:
summary(rdhte(y = y, x = x, covs.hte = w_left,
h = c(0.078, 0.116), vce = "hc3"))
#> Sharp RD Heterogeneous Treatment Effects: Subgroups.
#>
#> Number of Obs. 39534
#> BW type Manual
#> Kernel Triangular
#> VCE method HC3
#>
#> Number of Obs. 19723 19811
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#>
#> ================================================================================================================
#> Point Robust Inference
#> w_left Estimate z Pr(>|z|) [ 95% C.I. ] Nh- Nh+ h- h+
#> ----------------------------------------------------------------------------------------------------------------
#> 0 0.021 1.522 0.128 [-0.003 , 0.027] 4388 4453 0.078 0.078
#> 1 0.089 7.435 0.000 [0.064 , 0.110] 5010 4799 0.116 0.116
#> ================================================================================================================
summary(rdrobust(y = y[w_left == 1], x = x[w_left == 1],
h = 0.116, rho = 1, vce = "hc3"))
#> Call: rdrobust
#>
#> Sharp RD estimates using local polynomial regression.
#>
#> Number of Obs. 17596
#> BW type Manual
#> Kernel Triangular
#> VCE method HC3
#>
#> Left Right
#> Number of Obs. 9436 8160
#> Eff. Number of Obs. 5010 4799
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 0.116 0.116
#> BW bias (b) 0.116 0.116
#> rho (h/b) 1.000 1.000
#> Unique Obs. 9436 8160
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 0.089 7.435 0.000 [0.064 , 0.110]
#> =====================================================================
See also
rdhte(), rdbwhte(),
rdhte_lincom(), plot.rdhte().- The companion Stata package ships an
rdhte_plot command
that mirrors the R plot() method.
Reference
Calonico, S., Cattaneo, M.D., Farrell, M.H., Palomba, F., and
Titiunik, R. (2025). Treatment Effect Heterogeneity in Regression
Discontinuity Designs. arXiv:2503.13696.
Granzier, R., Pons, V., and Tricaud, C. (2023). Coordination and
Bandwagon Effects: How Past Rankings Shape the Behavior of Voters and
Candidates. American Economic Journal: Applied
Economics 15(4): 177-217.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.