Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
55,067 characters · 19 sections · 43 citation commands
rdhte: Conditional Average Treatment Effects in RD Designs
\inserttype[st0001]{article}
{ Sebastian Calonico\\UC Davis\\ Davis, CA\\[email removed] \and Matias D. Cattaneo\\Princeton University\\Princeton, NJ\\[email removed] \and Max H. Farrell\\UC Santa Barbara\\Santa Barbara, CA\\[email removed] \and Filippo Palomba\\Princeton University\\Princeton, NJ\\[email removed] \and Rocio Titiunik\\Princeton University\\Princeton, NJ\\[email removed] }
Studying causal heterogeneous treatment effects based on pretreatment covariates is an essential part of modern empirical work, as it helps uncover fairness concerns, differential impacts, and informs targeted policy interventions. While the regression discontinuity (RD) design has become a widely used tool for causal inference, existing methods primarily focus on estimating a single overall average treatment effect, leaving a gap in rigorous approaches for studying conditional average treatment effects. Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp addressed this gap by developing a unified, theoretically grounded methodology for heterogeneity analysis in RD designs, allowing researchers to systematically examine how treatment effects vary across subpopulations based on pretreatment characteristics. Cattaneo-Idrobo-Titiunik2019_book,Cattaneo-Idrobo-Titiunik2023_book provide a practical introduction to RD designs, estimation, and inference. For a review of recent literature on RD, see Cattaneo-Titiunik_2022_ARE.
This article discusses the general-purpose software package rdhte, which implements the main methodological results in Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp. The package is composed by three commands:
Section (ref) introduces the sharp RD setup, highlighting the role of pretreatment covariates, both heterogeneity analysis (which is our main focus) and for efficiency improvements. We detail how heterogeneous treatment effects are recovered using local polynomial regression methods with semi-linear interactions. The methods are adapted from Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp, and we refer to that paper for the relevant econometric theory. Cattaneo-Keele-Titiunik_2023_HandbookCh offer a broader discussion on the role of covariate-adjustments in RD designs.
Within the heterogeneity framework of Section (ref), there are two distinct cases depending on the type of covariates used: (i) dummy variables (or, more generally, factor variables) that identify orthogonal (mutually exclusive) subgroups; or (ii) generic covariates (discrete or continuous). The first case, treated in detail in Section (ref), arises when the covariates, and any included covariate interactions, correspond to non-overlapping subsets of the data, that is, when only binary orthogonal variables are used. This case brings the familiar subgroup analysis to RD designs, which naturally arises via indicator variables for different categories and from collapsing continuous variables into distinct bins. In this case, the treatment effect is different for each subgroup in an unrestricted way.
In the second case (Section (ref)), the covariates can be arbitrary, covering discrete, continuous, or mixed, and allowing for general interactions or transformations such as polynomials or other basis expansions thereof. For this generic setting, Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp provide formal identification, estimation, and robust bias-corrected inference methods for heterogeneity analysis based on standard local polynomial regression methods with semi-linear interactions. The treatment effect heterogeneity is defined as a varying coefficient linear function of the pretreatment covariates. As an alternative, see Reguly2021-wp and alcantara2025LearningConditionalAverage for heterogeneous treatment effect estimation leveraging machine learning methods, and hsu2019testing,Hsu-Shen_2021_JAE for specification testing methods.
The package rdhte complements the popular package rdrobust Calonico-Cattaneo-Farrell-Titiunik_2017_Stata, which focuses on estimation and robust bias-correction inference for the local average treatment effect at the cutoff in RD designs. Section (ref) compares the two RD packages, and explains precisely their differences. In certain specific cases, the package rdhte matches the analysis based on rdrobust.
Throughout Sections (ref), (ref), and (ref), we illustrate the features of {\tt rdhte} using the data from granzier2023coordination, who studied coordination behavior in French two-round elections. We add to their results by studying treatment effect heterogeneity. The running variable is the vote margin in the first round election, and the outcome of interest is a binary indicator for running in the second round or not. We explore heterogeneity by party ideology and party strength, as defined in more detail below.
Section (ref) concludes. Finally, the latest version of the package rdhte, replication files, and other related materials, are available at:
Extending RD analysis to heterogeneous treatment effects poses a nonparametric challenge, often requiring researchers to rely on semiparametric, yet parsimonious models. Following practice, Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp examines the common approach of using local least squares regression with linear interactions, clarifying the conditions under which this method yields meaningful causal interpretations. The authors established that when potential outcomes follow a local linear-in-parameters, functional coefficient model, heterogeneous effects are identifiable and interpretable, particularly for binary orthogonal covariates (subgroup effects). Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp further developed formal econometric methods for estimation and (robust bias-corrected) inference, including optimal bandwidth selection and standard error estimators robust to both heteroskedasticity and clustering. These results aim to add rigor and consistency to empirical practice, improving the applicability, and replicability, of RD heterogeneity analysis.
The observed random sample is $(Y_i, T_i, X_i, \mathbf{W}_i', \mathbf{Z}_i')'$, for $i=1,\dots,n$, where:
Without loss of generality, we set the cutoff $c=0$ to streamline the presentation. The canonical RD design corresponds to the case where neither $\mathbf{W}_i$ nor $\mathbf{Z}_i$ are included in the analysis.
The canonical sharp RD average treatment effect (at the cutoff $X_i=c=0$) is
It is common practice to employ the following least squares local polynomial RD estimator:
where $\mathbf{e}_{\ell}$ denotes the conformable unit vector with a $1$ in its $(\ell+1)$th element, and
with $\mathbf{r}_p(u) = (1,u,\dots,u^p)'$ denoting the $p$th order polynomial expansion, and $K_h(u) = K(u/h)/h$ for a kernel (i.e., weighting) function $K(\cdot)$ and bandwidth $h$.
For the local-linear case ($p=1$), the underlying implementation employing the Stata command reg is
properly localized to $c=0$ and weighted using the kernel $K(\cdot)$. Then, $\dot{\tau}$ corresponds to the coefficient estimate associated with t.
The classical RD treatment effect estimator $\dot{\tau}$ is consistent for $\tau$ under standard regularity conditions. Similarly, $\mathbf{e}_{1}'\dot{\boldsymbol{\beta}}$ is an estimator useful in the Kink RD design. Calonico-Cattaneo-Titiunik_2014_ECMA studied MSE-optimal point estimation and robust bias-corrected inference for this case, while Calonico-Cattaneo-Farrell_2020_ECTJ study optimal bandwidth selection for both point estimation and inference.
The covariates $\mathbf{Z}_i$ can be included in the regression estimation to improve efficiency in the estimation of RD treatment effect $\tau$. Calonico-Cattaneo-Farrell-Titiunik_2019_RESTAT studied MSE-optimal point estimation and robust bias-corrected inference for this case, and recommended the estimator
where
The covariate-adjusted, possibly more efficient RD estimator $\widetilde{\tau}$ is consistent for $\tau$ under regularity conditions. The estimated coefficients $\widetilde{\boldsymbol{\gamma}}$ do not have a causal interpretation; they are fitted to improve the precision of $\widetilde{\boldsymbol{\beta}}$. Standard Kink RD designs with covariate-adjustment for efficiency gains consider the estimator $\mathbf{e}_{1}'\widetilde{\boldsymbol{\beta}}$.
For the local-linear case ($p=1$), the underlying implementation employing the Stata command reg is
again properly localized to $c=0$ and weighted using the kernel $K(\cdot)$. Then, $\widetilde{\tau}$ corresponds to the coefficient associated with t.
The package rdrobust Calonico-Cattaneo-Farrell-Titiunik_2017_Stata provides software implementation for estimation and inference on the RD average treatment effects, both with and without adding covariates for efficiency gains.
To complement prior literature focusing on the RD average treatment effect $\tau$, Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp considers the (local to $X_i=c=0$) RD conditional average treatment effect (CATE) function
which employs the covariates $\mathbf{W}_i$ for RD treatment effect heterogeneity. The covariates $\mathbf{W}_i$ are distinct from those used for efficiency improvements, $\mathbf{Z}_i$, but importantly both are predetermined.
When $\mathbf{W}_i$ is continuous and/or high-dimensional, the RD CATE function $\kappa(\mathbf{w})$ can be difficult to estimate nonparametrically without further restrictions. Thus, it is common practice to employ the semilinear least squares local polynomial estimation procedure:
where $\mathbf{I}_d$ denotes the $(d\times d)$ identity matrix, $\mathbf{0}_{sd \times d}$ denotes the $(sd\times d)$ matrix of zeros, $\mathbf{w}$ takes values on the support of $\mathbf{W}_i$, and
with $\otimes$ denoting the Kronecker product. Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp study this estimation procedure and give easy-to-interpret sufficient conditions to ensure that $\widehat{\kappa}(\mathbf{w})$ is consistent for $\kappa(\mathbf{w})$. The most important of these is the assumption that the CATE function can be written as a (local) functional coefficient linear form: $\kappa(\mathbf{w}) = \theta(x) + \boldsymbol{\xi}(x)'\mathbf{w}$, see below. This is without loss of generality for binary orthogonal covariates (Section (ref)). Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp also establish MSE-optimal bandwidth selection and point estimation, and valid robust bias-corrected inference for uncertainty quantification.
The generic form of $\widehat{\kappa}(\mathbf{w})$, and its underlying least squares fit (ref), is notationally cumbersome but easy to understand. It arises from adding to the local polynomial estimation the interaction between the covariates $\mathbf{W}_i$ and the polynomial approximation $\mathbf{r}_s(X_i)$ (which may be a different polynomial than the main effect). For example, if $p=s=1$, so that $\mathbf{r}_s(X_i)=(1,X_i)'$, and $\mathbf{W}_i$ is a continuous variable, the {\tt Stata} implementation would be:
properly localized to $c=0$ and weighted using the kernel $K(\cdot)$. The estimate $\widehat{\kappa}(w_0)$ of Equation (ref), for a specific value $w_0$ is obtained as the coefficient on {\tt t} plus the product of $w_0$ and the coefficient on the interaction of {\tt t} and {\tt c.w}.
Covariate-adjustment based on $\mathbf{Z}_i$ for efficiency gains is also allowed: following Calonico-Cattaneo-Farrell-Titiunik_2019_RESTAT, in the regression (ref), $\mathbf{Z}_i$ and $\mathbf{Z}_i\otimes\mathbf{W}_i$ are included but without interaction with treatment assignment variable $T_i$. The resulting coefficients (on $\mathbf{Z}_i$ and $\mathbf{Z}_i\otimes\mathbf{W}_i$) do not have a causal interpretation, but the resulting estimator $\widehat{\kappa}(\mathbf{w})$ can exhibit efficiency gains. We omit further details for brevity, but the local-linear command is:
properly localized to $c=0$ and weighted using the kernel $K(\cdot)$. The CATE estimate is obtained exactly as before.
The causal interpretation of the probability limit of the estimator $\widehat{\kappa}(\mathbf{w})$ must be carefully considered and depends on the type of covariates as follows.
Regardless of the specific structure of $\mathbf{W}_i$, Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp develop MSE-optimal estimation and robust bias-correction inference for $\kappa(\mathbf{w})$, expanding prior results established for $\tau$. See Arai-Ichimura_2018_QE and Calonico-Cattaneo-Farrell_2020_ECTJ for more discussion on bandwidth selection for classical RD designs, Calonico-Cattaneo-Titiunik_2014_ECMA and Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2022_Bernoulli for foundational theoretical analysis of robust bias-correction inference. Hyytinen-Tukiainen-etal2018_QE and DeMagalhaes-etal_2025_PA offer comprehensive empirical validation of those estimation and inference methods.
We begin the illustration of heterogeneity analysis with {\tt rdhte} with the first case: binary orthogonal covariates (subgroup analysis). Recall the context of granzier2023coordination, where the outcome is $Y_i \in \{0,1\}$ indicates running in the second round election, and the running variable $X_i$ is the vote margin in the first round. We will show the features of {\tt rdhte} using different measures of party ideology and strength. Unless noted otherwise, all results will use standard errors clustered by district (stored as {\tt cluster_var}) following their original analysis. Because {\tt rdhte} is based on {\tt regress}, clustered standard errors are obtained using the standard {\tt vce} options: {\tt vce(cluster cluster_var}) yields HC1 clustered standard errors, while {\tt vce(hc2 cluster_var}) yields HC2, which is preferred. Currently, HC3 is not available, but it would provide a more robust option. We therefore recommend setting {\tt vce(hc2 cluster_var}), which we use throughout. In other software, HC3 is available and thus set as the default for both clustered and independent data.
The most basic form for such analysis is the case of a single binary variable, so that $d=1$ and $W_i \in \{0,1\}$. To illustrate, we create the dummy variable {\tt left}, taking the value one if the candidate belongs to the left or far left and zero otherwise. Running {\tt rdhte} yields a point estimate and robust bias corrected inference for each category.
{\fontsize{6}{7}\selectfont
}
We see that the RD treatment effect is larger for left-of-center candidates (0.089 versus 0.021) and is statistically significantly different from zero, which is not the case for the others. In this case, {\tt rdhte} automatically detects that {\tt w_left} is an indicator for subgroups because {\tt w_left} takes only the values zero and one. Any other coding (e.g., if it takes values 1 and 2) will require use of the {\tt i.}{\it varname} syntax. This output is omitted to save space but is available in the replication code.
The package, by default, automatically computes the optimal bandwidth for each group separately. The behavior can be overridden using the option {\tt bwjoint}. However, this is generally not advisable because the bias-variance trade-off that determines the optimal bandwidth need not be the same across subsets of the data, just as it need not be the same across different datasets. See the replication code for an example and the help file for {\tt rdbwhte} for more details on bandwidth selection.
A natural question following these results is whether the two subgroups are statistically significantly different. This is addressed using the featured post-estimation command {\tt rdhte_lincom}, which computes the point estimate of the difference (in this case $\widehat{\kappa}(1) - \widehat{\kappa}(0)$) and the robust bias-corrected confidence interval and p-value.
{\fontsize{8}{8}\selectfont
}
We see that, indeed, the effect for left-of-center candidates is significantly higher at the 5% level. In the next subsection, we also demonstrate post-estimation hypothesis testing using the Stata built-in command {\tt test}.
Moving beyond the case of a single binary variable, we now explore other instances where the set of covariates $\mathbf{W}_i$ are binary orthogonal variables and therefore identify subgroups of the data.
We obtain a more nuanced view of the heterogeneity by party ideology using the four-level factor variable {\tt w_ideology}. In this case, the four-level factor is “dummied out” to obtain an indicator for each group. This happens automatically using the {\tt i.} syntax.
{\fontsize{6}{7}\selectfont
}
For post-estimation, we use the built-in {\tt Stata} command {\tt test} to find that all the non-left categories are statistically indistinguishable from each other and from zero.
{\fontsize{8}{8}\selectfont
}
We bring in another pre-treatment variable for heterogeneity analysis: the strength at the national level, defined as the average of first-round vote shares of all candidates of the same orientation at the national level. The raw measure is stored as {\tt w_strength} and is considered in the next subsection. Here, we categorize strength into the four quartiles and obtain the effect for each. Compared to the case of {\tt w_ideology}, this variable has a clear ordering. Note here that the {\tt i.} syntax is required; otherwise {\tt rdhte} will yield a linear fit treating the variable as continuous (see below).
{\fontsize{5}{7}\selectfont
}
We find that the causal effect is increasing with strength, an interesting if not surprising finding. This analysis depends on binning, and can be viewed intuitively (though not formally) as a four-piece approximation to an unknown, nonparametric $\kappa(w)$, following Cattaneo-Crump-Farrell-Feng_2024_AER.
Finally, dummying out factor variables, and then interacting them, also yields binary orthogonal variables, one for each unique combination, provided the model is saturated. To illustrate, here we interact the dummy for left-of-center with an indicator for a “strong” candidate, defined as having above-median national strength.
{\fontsize{5}{7}\selectfont
}
As is typical for interaction effects, this adds nuance to the above findings, as we see that not only do left-of-center candidates have larger treatment effects, but this is even more pronounced when the candidates are stronger than the national median.
We now discuss the case of generic covariates: $\mathbf{W}_i$ is not a set of binary orthogonal variables, and so instead of obtaining heterogeneity by subgroups, we obtain a linear-in-parameters estimate of the function $\kappa(\mathbf{w})$. As discussed above, Calonico-Cattaneo-Farrell-Palomba-Titiunik_2025_wp use a functional coefficient model assumption to obtain a causal interpretation of the probability limit of $\widehat{\kappa}(\mathbf{w})$. In this case, instead of selecting an optimal bandwidth for each subset, a single bandwidth is selected.
To illustrate we use the measure of strength, {\tt w_strength}. Before, when grouped by quantiles, we saw that the treatment effect was increasing in strength. It may be that a linear fit is an appropriate and parsimonious way of capturing this relationship. Using {\tt w_strength}, we obtain the following. For later comparison with results from {\tt regress}, here we use the uniform kernel to obtain an unweighted local least squares fit. Note that the output of {\tt rdhte} is different relative to the subset analysis, reflecting that there is only one bandwidth and the coefficients here represent qualitatively different objects.
{\fontsize{8}{8}\selectfont
}
In Section (ref), each row of the output table showed the treatment effect estimate (and robust bias-corrected inference) for each subset of the data. Here, instead, we have the “intercept” and “slope” terms of the estimate $\widehat{\kappa}(\mathbf{w})$ of Equation (ref). That is, for a specific level of candidate strength $w_0$, we would obtain $\widehat{\kappa}(w_0) = -0.055 + 0.262 \times w_0$. The value -0.055 represents $\widehat{\kappa}(0)$, which may or may not have conceptual meaning, depending on the context and the definition of $\mathbf{W}_i$.
To build intuition, this interpretation is identical to interpreting coefficients in a linear regression, because here (by default) we are using {\tt rdhte} with a polynomial order of $p=s=1$. Recall that identical estimates can be obtained using {\tt regress}, properly localized and weighted, as shown in Equation (ref). Here, we have used the uniform kernel, so we must only localize using the same bandwidth in order to obtain the same point estimates. We extract the bandwidth from the ereturns of {\tt rdhte}, generate the treatment indicator, and then run the regression.
{\fontsize{8}{8}\selectfont
}
Here, the relevant coefficients are those on {\tt 1.T} and {\tt T\#c.w_strength}.
It is important to remember that this only replicates the point estimates, for inference, robust bias correction is required. The inference measures (standard errors, t-statistics, confidence intervals, and p-values) in this {\tt regress} output are not valid.
As a final illustration, consider the analog of the binary interaction at the close of Section (ref). There, we obtained an estimate for each of the four categories reflecting left or not and above/below median strength. Here, we consider the interaction of the binary indicator {\tt w_left} with the continuous measure of strength. By fully saturating the model, we obtain a separate intercept and slope for left-of-center and for center-and-right candidates.
{\fontsize{7}{8}\selectfont
}
This output reveals the same qualitative conclusion as the binary interaction, but expressed in a different way. It is worth noting that, because this model is fully saturated, the same results can be obtained studying {\tt w_strength} separately for each category of {\tt w_left}. This is shown in the replication code, but omitted here to save space.
In this section, we briefly illustrate how {\tt rdhte} does, and does not, match the output from the popular command {\tt rdrobust} Calonico-Cattaneo-Farrell-Titiunik_2017_Stata. For simplicity, all analyses in this section do not use clustered standard errors, and so for {\tt rdhte} the HC3 option is used. We show replication of average and of heterogeneous treatment effects. Both packages can use covariates for efficiency (using the {\tt covs_eff} option in {\tt rhdte}), but this is not a point of comparison here.
The discrepancies arise for three main reasons, all of which are driven by the fact that {\tt rdhte} (in all platforms) relies on built-in base commands for least squares regression.
First, consider the RD average treatment effect $\tau$ of Equation (ref).
Both {\tt rdrobust} and {\tt rdhte} can be used to obtain the corresponding estimator $\dot{\tau}$ of (ref) along with robust bias-corrected inference. For {\tt rdhte}, the average is estimated if no heterogeneity variables are specified. However, the default implementations differ, and so the output will not match, as shown here.
{\fontsize{8}{8}\selectfont
}
If all these settings are synchronized, the two commands report the same results (up to numerical differences in computation). We illustrate this point by setting $h=0.1$ (which automatically forces $h=b$ or $\rho=1$ in {\tt rdrobust}), and variance estimation to HC3.
{\fontsize{8}{8}\selectfont
}
The same equivalence can be obtained when conducting heterogeneity analysis by subgroup, because this amounts to estimating the average effect for a subset of the data. Here we illustrate that, again enforcing common settings, {\tt rdrobust} obtains the same results as {\tt rdhte} for left-of-center candidates (the other group is omitted to save space).
{\fontsize{6}{7}\selectfont
}
This article illustrated the main functionalities of the package {\tt rdhte} for heterogeneous RD treatment effects estimation and robust bias-corrected inference. We also discussed how this package complements the popular RD package {\tt rdrobust}. The methods implemented by the package {\tt rdhte} can also be used in the context of multi-cutoff and multi-score RD designs by discretizing the analysis along the multi-dimensional assignment rule; see Cattaneo-Titiunik-VazquezBare_2020_Stata for more discussion. Furthermore, {\tt rdhte} can also be used to implement other RD designs involving comparisons across subgroups, such as in difference-in-cutoffs designs Grembi-Nannicini-Troiano_2016_AEJ-Applied or dynamic designs Hsu-Shen_2024_QE. We do not discuss these connections further to conserve space.
Cattaneo and Titiunik gratefully acknowledge financial support from the National Science Foundation through grants SES-2019432 and SES-2241575.
Sebastian Calonico is an Assistant Professor at the University of California at Davis.
Matias D. Cattaneo is a Professor at Princeton University.
Filippo Palomba is a Ph.D. candidate at Princeton University.
Max H. Farrell is an Associate Professor at the University of California at Santa Barbara.
Rocio Titiunik is a Professor at Princeton University.