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.
88,886 characters · 17 sections · 38 citation commands
rd2d: Causal Inference in Boundary Discontinuity Designs
Keywords: treatment effects, regression discontinuity designs, nonparametric regression.
Regression Discontinuity (RD) designs are commonly used for treatment effect estimation and causal inference in quantitative sciences Cattaneo-Titiunik_2022_ARE. In their canonical form, each unit $i\in \{1,2,\cdots,n\}$ is assigned to control ($T_i=0$) or treatment ($T_i=1$) according to the discontinuous rule $T_i = \mathds{1}(X_i \geq c)$, where $X_i$ denotes a scalar score variable, $c$ denotes a scalar cutoff, and $\mathds{1}(\cdot)$ is the indicator function. The key idea underlying all RD designs is that units with a score near the cutoff determining treatment assignment are comparable in terms of all pretreatment observables and unobservable characteristics, the only difference being that some units are assigned to control ($X_i < c$) while other are assigned to treatment ($X_i \geq c$). Therefore, in the absence of score manipulation, units having a score near (but on different sides of) the cutoff can be used as counterfactual of each other to learn about causal treatment effects.
Boundary discontinuity designs generalize the canonical RD design to allow for a multi-dimensional score variable Papay-Willett-Murnane_2011_JoE,Reardon-Robinson_2012_JREE,Keele-Titiunik_2015_PA. The most common case is a bivariate score $\mathbf{X}_i = (X_{1i},X_{2i})^\top$ with an assignment boundary curve $\mathcal{B}$ on its support. For example, LondonoVelezRodriguezSanchez_2020_AEJ study the effect of a Colombian social program where $X_{1i}$ is a poverty index and $X_{2i}$ is an exam grade, for each student $i$, and where eligibility for receiving the treatment required either a minimum poverty index ($X_{1i} \geq c_1$) or a minimum exam grade ($X_{2i} \geq c_2$). Therefore, in their application, the boundary determining treatment is $\mathcal{B}=\{(X_{1i},X_{2i}): \{X_{1i} \geq c_1 \text{ and } X_{2i} = c_2\} \cup \{X_{1i}=c_1 \text{ and }X_{2i}\geq c_2\}\}$. This type of bivariate RD designs, where each score component has its own cutoff, is illustrated in Figure (ref). Another prototypical class of boundary discontinuity designs are the Geographic RD designs: for instance, Keele-Titiunik_2015_PA study the effect of political advertisements on voter turnout during a presidential campaign by leveraging sharp discontinuities in exposure to presidential ads induced across geographic media market boundaries. Figure (ref) illustrates a generic example of a Geographic RD design employing the boundary separating two US states. jardim2024local gives another recent empirical application of a boundary discontinuity design in the context of labor markets, and provides further references. See Cattaneo-Idrobo-Titiunik_2020_CUP,Cattaneo-Idrobo-Titiunik_2024_CUP for a two-part practical introductory monograph.
While classical RD designs based on a scalar score are well-understood in the literature, boundary discontinuity designs are surprisingly less studied. Methodological developments have lagged empirical practice for a while, leading to different approaches in practice, but without a foundational understanding of their relative merits. Cattaneo-Titiunik-Yu_2025_BoundaryRD address this gap in the literature by studying the properties of two leading approaches often used in empirical research leveraging boundary discontinuity designs:
Cattaneo-Titiunik-Yu_2025_BoundaryRD study the two methodologies, and establish novel identification, estimation and inference results, both pointwise for each point on $\mathcal{B}$ and uniformly over $\mathcal{B}$. Importantly, they demonstrate that the distance-based approach can exhibit a large bias near kinks or other irregularities in the assignment boundary $\mathcal{B}$, while the location-based approach remains valid even in those cases. Based on their findings, it is recommended to employ the bivariate location-based approach whenever possible, but their results provide foundational theoretical guidance for both empirical approaches.
This article introduces the R software package rd2d, which expands and implements the methodological results in Cattaneo-Titiunik-Yu_2025_BoundaryRD, thereby offering data-driven general-purpose methods for the analysis and interpretation of boundary discontinuity designs. The package includes the following four functions.
The methods print() and summary() are supported for objects returned by rd2d(), rdbw2d(), rd2d.dist(), and rdbw2d.dist(). We also demonstrate how to use the outputs to generate useful plots for empirical work, depicting treatment effect estimation, confidence intervals, and confidence bands, along the treatment assignment boundary $\mathcal{B}$.
In addition, the four functions in the package rd2d offer several practically relevant features, including (i) heteroskedasticity-robust and cluster-robust variance estimation, (ii) explicit regularization for the presence of mass points in the bivariate location score $\mathbf{X}_i$ or the univariate distance score $D_i(\mathbf{b}_j)$, and (iii) explicit regularization accounting for specific extreme shape features of the underlying unknown conditional expectation functions for bandwidth selection. In the case of rd2d(), the bivariate location-based approach, the data-driven point estimator is approximately MSE-optimal, while in the case of rd2d.dist(), the univariate distance-based approach, the rate-optimality of the data-driven point estimator depends on the underlying geometry of the assignment boundary $\mathcal{B}$. Putting aside the induced bias by the possibly non-smooth $\mathcal{B}$, both functions offer pointwise (for each $\mathbf{b}_j\in\mathcal{B}$) and uniform (over $\mathcal{B}$) robust bias-corrected inference Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2022_Bernoulli.
The main contribution of this article is to introduce and discuss the first general-purpose software implementation of (MSE-optimal) treatment effect estimation and (pointwise and uniform) uncertainty quantification methods for boundary discontinuity designs, given by the R package rd2d. To this end, the article develops second-generating DPI rules for bandwidth selection, along with principled regularization schemes for specific empirically relevant settings (e.g., mass points in $\mathbf{X}_i$) and practically relevant variance estimators (e.g., cluster-robust).
The rest of the paper proceeds as follows. Section (ref) reviews the main methodological contributions in Cattaneo-Titiunik-Yu_2025_BoundaryRD, and presents additional results related to bandwidth selection and regularized implementation. Section (ref) demonstrates the performance of the R package rd2d using a simulation study with a data generated process calibrated using the data in LondonoVelezRodriguezSanchez_2020_AEJ. Section (ref) concludes. Replication codes, background references, and other information related the software package rd2d can be found at: \url{https://rdpackages.github.io/rd2d/}.
We employ standard potential outcomes notation. Suppose that $(Y_i(0), Y_i(1), \mathbf{X}_i^\top)^\top$, $i = 1,2,\dots, n$, is a random sample, where $Y_i(0)$ and $Y_i(1)$ denote the scalar potential outcomes for unit $i$ under control and treatment assignment, respectively. Units are assigned to control group or treatment group according to their bivariate location score $\mathbf{X}_i= (X_{1i},X_{2i})^\top$ relative to a known one-dimensional boundary $\mathcal{B}$ splitting the support of $\mathbf{X}_i$ in two disjoint regions: $\mathcal{A}_0$ denotes the control region, and $\mathcal{A}_1$ denotes the treatment region. Thus, $\mathcal{B} = \mathtt{bd}(\mathcal{A}_0) \cap \mathtt{bd}(\mathcal{A}_1)$, where $\mathtt{bd}(\mathcal{A}_t)$ denotes the boundary of the set $\mathcal{A}_t$. The observed response variable is $Y_i = (1-T_i) \cdot Y_i(0) + T_i \cdot Y_i(1)$, where $T_i = \mathds{1}(\mathbf{X}_i \in \mathcal{A}_1)$. Without loss of generality, we assume that the boundary belongs to the treatment group, that is, $\mathtt{bd}(\mathcal{A}_1)\subset\mathcal{A}_1$ and $\mathcal{B}\cap\mathcal{A}_0=\emptyset$. Figure (ref) gives two graphical examples.
The causal parameter of interest is the average treatment effect curve along the boundary:
Identification follows directly by the usual continuity assumptions invoked in canonical RD designs Hahn-Todd-vanderKlaauw_2001_ECMA:
where $\mu_t(\mathbf{x})=\mathbb{E}[Y_i(t)|\mathbf{X}_i = \mathbf{x}]=\mathbb{E}[Y_i|\mathbf{X}_i = \mathbf{x}, T_i=t]$, $t\in\{0,1\}$, are assumed to be smooth functions. See Papay-Willett-Murnane_2011_JoE, Reardon-Robinson_2012_JREE, Keele-Titiunik_2015_PA, and references therein.
For implementation, the continuous assignment boundary $\mathcal{B}$ is first discretized into $J$ cutoff points $\mathbf{b}=(\mathbf{b}_1,\cdots,\mathbf{b}_J)^\top$ with $\mathbf{b}_j\in\mathcal{B}$ for all $j=1,\dots,J$. Then, the empirical analysis is conducted pointwise for each cutoff or uniformly over all cutoffs, employing either the bivariate location score $\mathbf{X}_i$ directly, or an induced univariate distance to each cutoff point. See Cattaneo-Idrobo-Titiunik_2024_CUP for an introductory discussion.
Following the methodological recommendations in Cattaneo-Titiunik-Yu_2025_BoundaryRD, most of the discussion focuses on location-based methods via bivariate local polynomial regression based on the data $(Y_i,\mathbf{X}_i^\top)^\top$, $i=1,\cdots,n$, which are implemented in the functions rd2d() and rdbw2d(). However, given their predominance in empirical work, Section (ref) also discusses distance-based methods via univariate local polynomial regression, which require a user-chosen scalar distance score to each cutoff point on the assignment boundary, and are implemented in the functions rd2d.dist() and rdbw2d.dist(). We omit assumptions and other technical details in the remaining of the paper, which can be found in the references given.
The location-based treatment effect curve estimator of $\tau(\mathbf{x})$ is
where, for $t \in \{0,1\}$,
with $\mathfrak{p}_p = (2+p)(1+p)/2-1$, $\mathbf{R}_p(\mathbf{u}) = (1,u_1,u_2,u_1^2,u_2^2,u_1u_2,\cdots,u_1^p,u_2^p)^\top$ denotes the $p$th order polynomial expansion of the bivariate vector $\mathbf{u}=(u_1,u_2)^\top$, $K_h(\mathbf{u})=K(u_1/h,u_2/h)/h^2$ for a bivariate kernel function $K(\cdot)$, and a bandwidth parameter $h$.
In practice, it is often important to first standardize each dimension of the bivariate location score $\mathbf{X}_i = (X_{1i},X_{2i})^\top$, and in some applications it may also be useful to allow for different bandwidths for control and treatment groups. Thus, the function rd2d() allows for four different bandwidths: $h_{01}$ is used for $X_{1i}$ and control units, $h_{02}$ is used for $X_{2i}$ and control units, $h_{11}$ is used for $X_{1i}$ and treatment units, and $h_{02}$ is used for $X_{2i}$ and treatment units. The data-driven bandwidth selection function rdbw2d() also allows for both standardization of each dimension of $\mathbf{X}_i$ (via the option stdvar=TRUE) and different bandwidth selection for control and treatment regions (via the options bwselect="msetwo" or bwselect="imsetwo"). Thus, the function rdbw2d() can report up to four distinct estimated bandwidths: $(\widehat{h}_{01},\widehat{h}_{02},\widehat{h}_{11},\widehat{h}_{11})$ corresponding to $(h_{01},h_{02},h_{11},h_{11})$. The discussion in this article focuses on a single common bandwidth $h$ for simplicity, but we explain how different MSE-optimal bandwidth are estimated as appropriate. In addition, while the notation does not explicitly reflect clustering, which is a common feature in geographic and other multidimensional RD designs, the package \texttt{rd2d} allows for cluster-robust inference as explained below.
Under minimal regularity conditions, the bivariate location-based treatment effect estimator $\widehat{\tau}(\mathbf{x})$ is pointwise and uniform consistent for the treatment effect curve $\tau(\mathbf{x})$ along the assignment boundary: $\widehat{\tau}(\mathbf{x}) \to_\mathbb{P} \tau(\mathbf{x})$ for each $\mathbf{x}\in\mathcal{B}$, and $\sup_{\mathbf{x} \in \mathcal{B}} |\widehat{\tau}(\mathbf{x}) - \tau(\mathbf{x})| \to_\mathbb{P} 0$, where $\to_\mathbb{P}$ denotes convergence in probability as $h\to0$ and $nh^2\to\infty$. Furthermore, under similar regularity conditions, precise (conditional) MSE expansions can be established along $\mathcal{B}$, which can then be used for principled bandwidth selection.
The pointwise (conditional) MSE expansion is
where $V_{\mathbf{x}} = V_{1,\mathbf{x}} + V_{0,\mathbf{x}}$ and $B_{\mathbf{x}} = B_{1,\mathbf{x}} - B_{0,\mathbf{x}}$ denote the fixed-$n$ conditional variance and the leading conditional bias of the treatment effect estimator, respectively, and $\approx_\mathbb{P}$ denotes equality in probability up to vanishing higher-order terms. More precisely, using standard multi-index notation and least squares algebra, the fixed-$n$ conditional variance and leading conditional bias for each group $t\in \{0,1\}$ are
respectively, where
and $\varepsilon_{i,t,\mathbf{x}} = Y_i - \mathds{1}(\mathbf{X}_i \in \mathcal{A}_t) \mu_t(\mathbf{x})$.
Therefore, for each $\mathbf{x}\in\mathcal{B}$, and noting that $\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}$, $\boldsymbol{\Sigma}_{t,\mathbf{x},\mathbf{x}}$, and $\widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k})$ converge in probability (as $nh^2\to\infty$) to well-defined limits independent on the bandwidth $h$, an MSE-optimal bandwidth choice is
provided that $B_{\mathbf{x}} \neq 0$.
Similarly, given a weighting function $w(\mathbf{x})$, an integrated MSE expansion is
using the notation already introduced. Therefore, along $\mathcal{B}$, an integrated MSE-optimal bandwidth choice is
provided that $\int_{\mathcal{B}} B_{\mathbf{x}}^2 w(\mathbf{x}) d \mathbf{x} \neq 0$.
The basic MSE-optimal and IMSE-optimal bandwidth choices, $h_{\operatorname{MSE},\mathbf{x}}$ and $h_{\operatorname{IMSE}}$, can be extended to accommodate different selections for each coordinate of the bivariate location score and/or for control and treatment groups separately. Different bandwidths for each coordinate in $\mathbf{X}_i = (X_{1i},X_{2i})^\top$ are obtained by first standardizing each component, then applying the basic bandwidth rules, and finally removing the standardization: letting $\tilde{\mathbf{X}}_i = (X_{1i}/\sigma_{X_1},X_{2i}/\sigma_{X_2})^\top$, where $\sigma_{X_l}^2 = \mathbb{V}[X_{li}]$ for $l\in\{1,2\}$, then for each coordinate $l\in\{1,2\}$ the bandwidth selectors are
where $h_{\operatorname{MSE},\mathbf{x}}$ and $h_{\operatorname{IMSE}}$ are computed using the standardized bivariate location score $\tilde{\mathbf{X}}_i$ (instead of using the original score $\mathbf{X}_i$). Different bandwidth selection for control and treatment groups are obtained by implementing $h_{\operatorname{MSE},\mathbf{x}}$ and $h_{\operatorname{IMSE}}$ for each group $t\in\{0,1\}$ separately:
assuming the denominators are not zero.
A combination of the two ideas gives four distinct MSE-optimal and IMSE-optimal bandwidth selection rules: $(\sigma_{X_1} \cdot h_{\operatorname{MSE},0,\mathbf{x}},\sigma_{X_2} \cdot h_{\operatorname{MSE},0,\mathbf{x}},\sigma_{X_1} \cdot h_{\operatorname{MSE},1,\mathbf{x}},\sigma_{X_2} \cdot h_{\operatorname{MSE},1,\mathbf{x}})$ and $(\sigma_{X_1} \cdot h_{\operatorname{IMSE},0},\sigma_{X_2} \cdot h_{\operatorname{IMSE},0},\sigma_{X_1} \cdot h_{\operatorname{IMSE},1},\sigma_{X_2} \cdot h_{\operatorname{IMSE},1})$, respectively. For implementation, the package rd2d employs the following options:
To implement the bandwidth selection procedures it is necessary to (i) estimate the unconditional variances $\sigma^2_{X_l}$ for standardization of each coordinate of $\mathbf{X}_i=(\mathbf{X}_{li}:l=1,2)^\top$, as needed; (ii) estimate residuals $\varepsilon_{i,t,\mathbf{x}}$ entering the variances $V_{t,\mathbf{x}}$ for each group $t\in\{0,1\}$; (iii) estimate higher-order derivatives $(\mu_t^{(\mathbf{k})}(\mathbf{x}) : |\mathbf{k}| = p + 1)$ entering the bias $B_{t,\mathbf{x}}$ quantities for each group $t\in\{0,1\}$; and (iv) select (preliminary) bandwidth entering the matrices $\widehat{\boldsymbol{\Gamma}}_{t,\mathbf{x}}$, $\boldsymbol{\Sigma}_{t,\mathbf{x},\mathbf{x}}$, and $\widehat{\boldsymbol{\vartheta}}_{t,\mathbf{x}}(\mathbf{k})$ for each group $t\in\{0,1\}$. (A data-driven regularization term is also added to the denominators to avoid near-zero bias, as explained below.) These unknown quantities are estimated as follows.
To assess uncertainty in the estimation of the causal effects $\tau(\mathbf{x})$ along the boundary $\mathcal{B}$, we consider the usual Wald-type test statistic
where $\widehat{\tau}(\mathbf{x})$ and $\widehat{V}_{\mathbf{x}}$ are constructed using the polynomial order $p$ and bandwidth $h$, after choosing the appropriate HC and CR variance estimator, as explained above in the context of bandwidth selection.
In practice, the bandwidth is chosen to be (I)MSE-optimal, and thus the sampling distribution of the statistic satisfies the following pointwise distributional approximation:
where $\overset{a}{\thicksim}$ denote an approximation in distribution as $nh^2\to\infty$, $\mathsf{N}(0,1)$ denotes the standard Gaussian distribution, and $\text{Bias}(\mathbf{x})$ denotes the standardized leading bias emerging whenever a “large” bandwidth is used (i.e., when the (I)MSE-optimal bandwidth is used, or any other bandwidth choice such that $\text{Bias}(\mathbf{x}) \not\to_\mathbb{P} 0$). The standard confidence interval estimator with $100(1-\alpha)\%$ nominal coverage is
where $\Phi_\alpha$ is the $\alpha$th quantile of the standard Gaussian distribution. However, for “large” bandwidths such as the (I)MSE-optimal choice, $\operatorname{CI}(\mathbf{x})$ will be invalid due to the bias $\text{Bias}(\mathbf{x})$, thereby delivering empirical coverage well below its nominal target. A solution to this problem is to employ ad-hoc undersmoothing, that is, to implement $\operatorname{CI}(\mathbf{x})$ with a “smaller” bandwidth $h$ relative to the (I)MSE-optimal choice. Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2022_Bernoulli showed that undersmoothing is sub-optimal (possibly invalid) under standard assumptions, while the robust bias-correction (RBC) methodology introduced by Calonico-Cattaneo-Titiunik_2014_ECMA enjoys validity and better (in some cases optimal) higher-order distributional properties. The core idea behind the RBC approach can be summarized as follows: (i) employ the (I)MSE-optimal bandwidth for constructing the point estimator $\widehat{\tau}(\mathbf{x})$, (ii) de-bias (bias correct) the numerator of the statistic $\operatorname{T}(\mathbf{x})$, and (iii) adjust the variance estimate to account for the variability introduced by the debiasing of $\widehat{\tau}(\mathbf{x})$.
The implementation of the RBC inference methodology is straightforward: given a chosen (I)MSE-optimal bandwidth for the $p$th order local polynomial point estimator $\widehat{\tau}(\mathbf{x})$, an adjsted test statistic is constructed using a $q$th order local polynomial point estimator $\widehat{\tau}_q(\mathbf{x})$ and its associated variance estimate $\widehat{V}_{\mathbf{x},q}$, with $q>p$. Thus, the approach employs the RBC statistic $T_{q}(\mathbf{x}) = (\widehat{\tau}_q(\mathbf{x}) - \tau(\mathbf{x}))/\sqrt{\widehat{V}_{\mathbf{x},q}/(nh^2)}$, instead of the original statistic $\operatorname{T}(\mathbf{x})$ above. The resulting RBC confidence interval estimator is
which is constructed with an (I)MSE-optimal bandwidth choice for the $p$th order local polynomial point estimator $\widehat{\tau}(\mathbf{x})$. (This amounts to a form of robust bias correction because $\widehat{\tau}_q(\mathbf{x}) = \widehat{\tau}(\mathbf{x}) - h^{p+1} \widehat{B}_\mathbf{x}$, where $\widehat{B}(\mathbf{x})$ is an “estimate” of $B_\mathbf{x}$.) The package rd2d employs $p=1$ and $q=2$ as defaults, which are standard choices for implementation.
Uniform inference and confidence bands along the boundary $\mathcal{B}$ also employ RBC methodology. Specifically, given an (I)MSE-optimal bandwidth for the $p$th order local polynomial point estimator $\widehat{\tau}(\mathbf{x})$, the associated RBC confidence band estimate is
where $\mathsf{q}_{\alpha}$ denotes a suitably chosen quantile to control false rejections uniformly over $\mathcal{B}$. In practice, the continuous assignment boundary is discretized to consider the $J$ cutoff in $\mathbf{b} = (\mathbf{b}_1,\cdots,\mathbf{b}_J)^\top$ jointly. Then, a feasible quantile choice is
where the $J$-dimensional standard Gaussian vector $\mathbf{Z}\thicksim\mathsf{N}(\mathbf{0}_J, \mathbf{I}_J)$ is independent of the data, and the $(J\times J)$ covariance matrix $\widehat{\mathbf{C}}_{\mathbf{b},q} = \big[\widehat{V}_{\mathbf{b}_j,\mathbf{b}_k} / (\widehat{V}_{\mathbf{b}_j,\mathbf{b}_j} \widehat{V}_{\mathbf{b}_k,\mathbf{b}_k})^{1/2}: 1\leq j,k \leq J \big]$ where $\widehat{V}_{\mathbf{b}_j,\mathbf{b}_k} = \widehat{V}_{1,\mathbf{b}_j,\mathbf{b}_k} + \widehat{V}_{0,\mathbf{b}_j,\mathbf{b}_k}$ with
and
with $\widehat{\varepsilon}_{i,t,\mathbf{x}} = Y_i - \mathbf{R}_q(\mathbf{X}_i - \mathbf{x})^\top\widehat{\boldsymbol{\beta}}_t(\mathbf{x})$, for $t\in\{0,1\}$. This formulation is robust to unknown conditional heteroskedasticity, while the clustered-robust analogue formula is omitted to save space; see Zeileis_2004_JSS and Zeileis-etal_2020_JSS.
The RBC method produces confidence intervals/bands that are not centered at the treatment effect point estimator because different polynomial orders are used for estimation and inference. As a result, the point estimates may lie outside the RBC confidence intervals/bands, particularly if the underlying treatment effect curve $\tau(\mathbf{x})$ exhibits high curvature at certain evaluation points. One possible solution is to increase the polynomial orders $p$ and $q$, or to use a bandwidth smaller than the (I)MSE-optimal one.
The package rd2d implements several regualization schemes to ensure robustness in applications.
For a weight function $w: \mathcal{B} \rightarrow [0,\infty)$, the aggregated average treatment effect (AATE) along the boundary $\mathcal{B}$ is
A choice of weight function recovers a scalar causal effect, aggregating potential heterogeneous treatment effects along the boundary. Estimation and inference methods for this class of causal parameters can be deduced from our results. For example, consider the generic plug-in “estimator” of $\tau_{\text{AATE},\mathcal{B}}$ given by
where the cutoff points $\mathbf{b}=(\mathbf{b}_1,\dots,\mathbf{b}_J)^\top$ along the assignment boundary are assumed to be evenly-spaced over $\mathcal{B}$, and $\widehat{\tau}(\mathbf{b}_j)$ is constructed using $p$th order local polynomial fits. It follows that the IMSE-optimal bandwidth choice is more natural. Feasible RBC confidence intervals are
where $\widehat{\tau}_{\text{AATE},\mathbf{b},q}$ is the point estimator $\widehat{\tau}_{\text{AATE},\mathbf{b}}$ constructed using a $q$th order polynomial basis, $\widehat{\mathbf{V}}_{\mathbf{b},q} = \big[\widehat{V}_{\mathbf{b}_j,\mathbf{b}_k} : 1\leq j,k \leq J \big]$ denotes its associated variance estimator, as explained previously, and $\mathbf{w} = (w(\mathbf{b}_1), \dots, w(\mathbf{b}_J))^\top/(\sum_{j = 1}^J w(\mathbf{b}_j))$ is the user-chosen vector of weights. In rd2d, AATE estimation and inference are implemented via the optional argument AATE = w in the summary() method for rd2d objects.
For each unit $i=1,\ldots,n$, their scalar distance-based score to the boundary point $\mathbf{x}=(x_1,x_2)^\top\in\mathcal{B}$ is $D_i(\mathbf{x}) = (2T_i-1)\mathcal{d}(\mathbf{X}_i, \mathbf{x})$, where $\mathcal{d}(\cdot,\cdot)$ denotes a distance function such as the Euclidean distance $\mathcal{d}(\mathbf{X}_i, \mathbf{x}) = \|\mathbf{X}_i - \mathbf{x}\| = \sqrt{(X_{1i} - x_1)^2 + (X_{2i} - x_2)^2}$. Therefore, for each $\mathbf{x}\in\mathcal{B}$, the setup reduces to a standard univariate RD design with distance score $D_i(\mathbf{x})\in\mathbb{R}$ and cutoff $c=0$, the observed data now being $(Y_1,D_1(\mathbf{x})),\cdots,(Y_n,D_n(\mathbf{x}))$ for each point on the assignment boundary $\mathbf{b}=(\mathbf{b}_1,\cdots,\mathbf{b}_J)^\top$.
The distance-based local polynomial treatment effect curve estimator of $\tau(\mathbf{x})$ is
where, for $t \in \{0,1\}$,
with $\mathbf{r}_p(u)=(1,u,u^2,\cdots,u^p)^\top$ the usual univariate polynomial basis, $k_h(u)=k(u/h)/h^2$ for univariate kernel function $k(\cdot)$ and bandwidth parameter $h$, and $\mathcal{I}_0 = (-\infty,0)$ and $\mathcal{I}_1 = [0,\infty)$. Cattaneo-Titiunik-Yu_2025_BoundaryRD studied the statistical properties of the distance-based approach in boundary discontinuity designs, and obtained the following main results (under regularity conditions).
Other implementation and regularization methods follow the same logic as for bivariate location-based estimation, taking into account the distance variable explicitly. In particular, $\text{\texttt{bwcheck}} = 50 + p + 1$ is used as default. See Cattaneo-Titiunik-Yu_2025_BoundaryRD for omitted technical and methodological details.
We illustrate the capabilities of the general-purpose R software package rd2d with a synthetic dataset of size $n = 20,000$ calibrated using the Ser Pilo Paga (SPP) dataset LondonoVelezRodriguezSanchez_2020_AEJ. We set $T_i = \mathds{1}(X_{1i} \geq 0, X_{2i} \geq 0)$, where $T_i = 1$ indicates unit $i$ is in the treatment group, and $T_i = 0$ indicates unit $i$ is in the control group. Covariates $\mathbf{X}_i = (X_{1i}, X_{2i})^\top$ are drawn from the product distribution $(100\mathsf{Beta}(3,4)-25,\mathsf{B}_2 = 100\mathsf{Beta}(3,4)-25)$ with independent components. Potential outcomes are generated by
where $\mathbf{X}_i$, $\varepsilon_{0,i}$ and $\varepsilon_{1,i}$ are mutually independent, and $\varepsilon_{t,i}\thicksim \mathsf{N}(0,\sigma_t^2)$, for $i = 1,2, \cdots, n$ and $t = 0,1$. We consider two DGPs as in Table (ref), where coefficients are estimated from the Ser Pilo Paga (SPP) dataset LondonoVelezRodriguezSanchez_2020_AEJ, and scaled by a factor of 2 to enhance signal strength.
Figure (ref) presents a scatterplot of a synthetic dataset constructed using the data generating process described above, which we use for numerical illustration of the main capabilities of the package rd2d. The gray assignment boundary separates the control (blue) and treatment (red) groups. The plot also includes forty grid points along the boundary, where the cutoff $\mathbf{b}_{21}$ is a kink point. Figure (ref) presents the population treatment effects corresponding to the two data generating processes. To demonstrate the variability in the outcome variable, the figure includes synthetic data points derived from 300 uniform draws along the boundary.
The function rd2d() provides point estimation, robust confidence intervals, and robust uniform confidence bands for (derivatives of) treatment effect function $\tau(\mathbf{x})$ based on bivariate location-based local polynomial regression. It takes as input an outcome vector y, a bivariate location score matrix X, a treatment indicator vector t, and a grid of evaluation points b along the boundary $\mathcal{B}$.
Optional arguments include bandwidth choices h, degrees of polynomial for point estimation (p) and inference (q), the partial derivative of treatment effect to be estimated deriv, and confidence level level. When optional arguments are not provided, the function defaults to estimate the value of treatment effect deriv = c(0,0) using the MSE-optimal bandwidth, with \texttt{p = 1} degree polynomial for point estimation and \texttt{q = 2} degree polynomial for robust bias-corrected confidence interval (and bands if requested), and \texttt{level = 95} percentage points confidence level. Additionally, \texttt{kernel_type} indicates whether a product kernel (\texttt{"prod"}) or a radial kernel (\texttt{"rad"}) is used for weighting. The default is \texttt{kernel_type = "prod"}.
Below is a demonstration of rd2d() applied on the synthetic dataset, with results for selected indices printed using the summary() method.
The first part of the output provides basic information on the options specified in the function. For example, the default estimand is the value of treatment effect, indicated by deriv = (0,0). The rest of the output gives estimation results, including (i) b1 and b2: First and second coordinate of the evaluation points; (ii) Coef.: Point estimation of (derivative) of treatment effect using $p=$ p polynomial order; and (iii) t-statistics, (iv) p-value, and (v) level% confidence intervals using $q=$ \texttt{q} polynomial order. When $q>p$, the resulting inference procedures correspond to robust bias correction Calonico-Cattaneo-Titiunik_2014_ECMA,Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2022_Bernoulli, which is the default and recommended method; $p=q$ corresponds to standard least squares methods. Point estimates, standard errors, and other information can be easily extracted for further statistical analysis. The output is stored in a standard matrix, and can be accessed with the following command,
where result.rd2d\$main.A0 contains results for the control group, result.rd2d\$main.A1 contains results for the treatment group, and result.rd2d\$main contains results for both.
The summary() method allows for three optional arguments. First, the option subset takes the indices of evaluation points to be presented, which should be a subset of c(1:nrow(eval)). The default is NULL, and thus all evaluation points are presented. Second, CBuniform is boolean variable for confidence bands construction, where FALSE indicates that pointwise confidence intervals are provided, and \texttt{TRUE} indicates uniform confidence bands are provided. The default is \texttt{CBuniform = FALSE}.
When the optional argument AATE $=\bigl(w_j\bigr)_{j=1}^J$, representing the weights for the evaluation points $\{\boldsymbol{b}_j\}_{j=1}^J$, is supplied to summary(), the output gains an extra row labeled “AATE.” This row reports the aggregated average treatment effect estimator \[ \widehat{\tau}_{\text{AATE}} =\frac{\sum_{j=1}^{J} w_j\,\widehat{\tau}\!\bigl(\boldsymbol{b}_j\bigr)} {\sum_{j=1}^{J} w_j}, \] along with its inference statistics.
Finally, the summary() method allows for presenting the underlying bandwidths used and associated effective sample sizes via the option output = "bw".
The function rdbw2d() is used for MSE (or IMSE) optimal bandwidth implementation, and is used internally in rd2d() when user does not specify bandwidth choices manually. The function takes the same input data as rd2d(), that is, an outcome vector y, a bivariate location score matrix X, a treatment indicator vector t, and a grid of evaluation points \texttt{b} along the boundary $\mathcal{B}$. In addition, the option \texttt{bwselect} encodes four options of bandwidth type: (i) \texttt{"mserd"} finds the MSE-optimal bandwidth for estimating (derivatives of) treatment effect, (ii) \texttt{"imserd"} finds the integreated MSE optimal bandwidth for estimating (derivatives of) treatment effect, (iii) \texttt{"msetwo"} finds the MSE optimal bandwidth for estimation (the derivatives of) conditional means of two potential outcome variables, (iv) \texttt{"imsetwo"} finds the integrated MSE optimal bandwidth for estimation (the derivatives of) conditional means of two potential outcome variables. The default is \texttt{bwselect = "mserd"}. An additional Boolean argument \texttt{stdvar} indicates whether the covariates are first standardized to unit standard deviation in each coordinate, in which case the optimal bandwidth is estimated and then converted back to the original scale. The default is \texttt{stdvar = TRUE}.
The first part of the summary output lists the options used for bandwidth selection. The second part of the summary output gives bandwidth selection results, including: (i) Boundary points, b1 for the first coordinate and b2 for the second coordinate; (ii) Bandwidths for control group, h01 for the first coordinate and h02 for the second coordinate; (iii) Bandwidths for treatment group, \texttt{h11} for the first coordinate and \texttt{h12} for the second coordinate.
The function rd2d.dist() provides point estimation and inference for boundary treatment effects using distance-based univariate local polynomial regression. It takes as input an outcome vector y, and a signed distance matrix D of distance to each boundary point, where each column of D corresponds to the signed distance from all observations to one evaluation point, with a positive sign indicating the unit is in the treatment group and a negative sign indicating the unit is in the control group.
Optional arguments include evaluation points b, bandwidth choices h, degrees of polynomial for point estimation (p) and inference (q), option for kink adjustment kink, and confidence level level, among other options. When not provided, the function defaults to estimate the value of treatment effect using the MSE-optimal bandwidth without kink adjustment (\texttt{kink = "off"}), using \texttt{p = 1} degree polynomial for point estimation and \texttt{q = 2} degree polynomial for robust bias-corrected confidence intervals and bands, providing \texttt{level = 95}% confidence interval and uniform confidence bands, and without displaying of evaluation points (\texttt{b = NULL}).
The first part of the summary output provides basic information on the options specified in the function. The rest of the summary output gives estimation and inference results, including: (i) b1 and b2 (when b is provided) report first and second coordinates of the evaluation points; (ii) Coef. reports the treatment effect estimate using a $p$th order polynomial regression; and (iii) the last three columns correspond to t-statistic, p-value and \texttt{level}% confidence intervals using $q$th order polynomial regression. The \texttt{summary()} method also has the option of displaying the uniform confidence bands (instead of the confidence intervals) in the last two columns as follows (numerical results omitted to conserve space).
When the optional argument AATE $=\bigl(w_j\bigr)_{j=1}^J$, representing the weights for the evaluation points $\{\boldsymbol{b}_j\}_{j=1}^J$, is supplied to summary(), the output gains an extra row labeled “AATE” for estimation and inference of aggregated average treatment effect.
In addition, The summary() method can also display the underlying bandwidths and effective samples sizes as follows (numerical results omitted to conserve space).
Point estimates, standard errors, and other information can be easily extracted for further statistical analysis. The output is stored in a standard matrix, and can be accessed with the following command.
result.dist\$main.A0 contains results for the control group, result.dist\$main.A1 contains results for the treatment group, and result.dist\$main contains results for both.
The function rdbw2d.dist() is used to implement MSE (or IMSE) rate-optimal ROT bandwidth selectors, and is used internally in rd2d.dist() when user does not provide bandwidths manually. As for rdbw2d(), four bandwidth types are allowed, that is, bwselect can be "mserd", "imserd", \texttt{"msetwo"} or \texttt{"imsetwo"}.
An additional argument kink, taking values "off" or "on", indicates whether a kink adjustment is made for estimation and inference. The default is kink = "off", but when kink = "on" is specified then the bandwidth is shrank to account for lack of smoothness of the assignment boundary $\mathcal{B}$.
The package rd2d provides an array of estimation and inference results that can be used for graphical presentation. Figure (ref) compares three point estimation methods: (i) bivariate method via rd2d; (ii) distance-based method via rd2d.dist, ignoring the kink (default argument kink = "off"); and (iii) distance-based method via rd2d.dist, adjusting for kink (kink = "on"). Figure (ref) plots point estimation using the bivariate method, along with its associated robust bias-corrected confidence intervals and confidence band. Figure (ref) presents a heatmap of treatment effects along the boundary, with high to low point estimation indicated by red to blue colors. Finally, Figure (ref) presents a heatmap of p-values along the boundary, with five colors assigned to five ranges of values. The codes for generating the graphical presentations are given in the replication \texttt{R} file.
It can be seen that distance-based estimation with kink = "off" overshoots compared to the bivariate estimation before the kink, and undershoots after the kink. This corresponds to the phenomena of getting a first-order bias using distance based method in the presence of a kink, despite using local polynomial regression of degree greater than or equal to $1$. See Cattaneo-Titiunik-Yu_2025_BoundaryRD for more methodological and theoretical discussions.
The discussion so far employed one realization of the data generating process to illustrate the main features of the package rd2d. In this final section, we conduct a Monte Carlo experiment to assess the performance of the package in repeated sampling. We consider $m = 1,000$ simulations of the two data generating processes defined at the beginning of this section, and report the simulation results in Table (ref) (DGP 1: linear model) and Table (ref) (DGP 2: quadratic model). Three methods are used and compared: the bivariate method rd2d, the distance-based method rd2d.dist ignoring the presence of the kink in the boundary (kink = "off"), and the distance-based method rd2d.dist adjusting for the kink (kink = "on"). Bandwidths are chosen automatically by the package, and their average across simulations is reported. We also report diagnostic measures including: bias, standard deviation of point estimator, root mean-squared error of point estimator, pointwise empirical coverage, pointwise interval length, uniform empirical coverage, and uniform interval length.
For both DPG 1 and DGP 2, all of the three methods give pointwise coverage around 95%, while the pointwise interval length for rd2d and rd2d.dist (kink = "off") are smaller compared to rd2d.dist (kink = "on"). This is likely due to bandwidth shrinkage for kink adjustment, which does results in a significantly smaller bias compared to the one ignoring the kink. Both rd2d and rd2d.dist (kink = "off") give around 95% uniform coverage and a relatively shorter interval length compared to rd2d.dist (kink = "on"), likely due to the same reason.
This paper introduced the R software package rd2d for causal inference in Boundary Discontinuity designs. The package provides pointwise and uniform (over the treatment assignment boundary) estimation and inference methods employing either a bivariate location score or a univariate distance score. In addition, the methods can be used for graphical presentation. From a methodological perspective, this paper introduced second generation bandwidth selection methods complementing the main results in Cattaneo-Titiunik-Yu_2025_BoundaryRD. Simulation evidence demonstrated a good performance of the package rd2d. Replication codes and related information are available at: \url{https://rdpackages.github.io/rd2d/}.
We thank Xinwei Ma and Filippo Palomba for insightful discussions. Cattaneo and Titiunik gratefully acknowledge financial support from the National Science Foundation (SES-2019432 and SES-2241575). Cattaneo gratefully acknowledge financial support from the Data-Driven Social Science initiative at Princeton University.