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.
77,463 characters · 18 sections · 116 citation commands
Using Forests in Multivariate Regression Discontinuity Designs
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 { \\ Department of Economics, Cornell University\\ and \\ Yuan Qi\\ [email removed]}\\ Department of Economics, University of Washington} } \fi
\if10 {
} \fi
\spacingset{1.45}
Regression discontinuity (RD) designs are widely used to estimate causal effects in observational settings, where individuals receive a binary treatment if an observed running variable is above some known threshold. Under the potential outcomes framework, HTV01 show that the treatment is as good as random in the subpopulation arbitrarily close to the threshold, if the conditional expectations of the two potential outcomes are continuous at the cutoff.\footnote{See, e.g., IL08 and CT22, for reviews on RD designs.} Estimation and inference of treatment effects in this context are, therefore, a local problem in the neighborhood of the threshold. Following HTV01 and P03, local polynomial regressions remain the most popular estimation strategy in the empirical RD literature. At the heart of these local methods is bandwidth selection for the weighting kernel. In designs with one running variable, IK12 propose an optimal bandwidth from minimizing the asymptotic mean squared error (MSE) at the cutoff. This procedure, however, leads to an asymptotically non-vanishing bias. CCT14 put forth a bias-corrected estimator that first subtracts an estimated bias term and then adjusts the standard error accordingly.
In many situations, the treatment status is determined by more than one criterion. Yet, univariate local linear regressions do not readily generalize to the multivariate case: kernel-based methods are less useful when the dimension far exceeds two HTF16, and standard bias correction quickly becomes algebraically involved and computationally difficult. Existing methods for multiple scores include multiple bandwidths selected from cross-validation PWM11, CL18, and Z12 extends the MSE-optimal bandwidth of IK12 to the bivariate case. But, as pointed out by IW19, these methods yield oversmoothing bandwidths and asymptotically biased estimates;\footnote{CL18 also consider rule-of-thumb bandwidths that scale with $n^{-1/6}$, which are again oversmoothing in the sense that $nh^5\not\to0$.} to avoid selecting a multivariate bandwidth, following the literature on minimax linear estimation armstrong2018optimal, kolesar2018inference, they propose an estimator that is minimax-optimal over a class of conditional expectations with second derivatives bounded by a given constant $\mathcal{B}$. This method, however, requires estimating $\mathcal{B}$ and, as we show in Section (ref), its performance depends crucially on accurate estimation of this nuisance parameter.
In this paper, we conduct a systematic Monte Carlo simulation to evaluate existing approaches to multivariate RD designs, including the minimax-optimal estimator of IW19 and the commonly used local linear regression approach with a transformed univariate score: to sidestep the complications of running local linear regressions in higher dimensions, empirical researchers often reduce the problem to a univariate one. Our contribution is twofold. First, we point out a potential issue of the common practice of implementing univariate local linear regressions when the running variable is in fact multivariate, give a formal characterization of when such problem arises, and show in simulations that it can produce biased estimates even for very large sample sizes. Second, we demonstrate two additional valid, easy-to-implement estimators that can be leveraged in multivariate RD designs, and together with a systematic Monte Carlo evaluation of existing estimators and an empirical application, provide researchers with a menu of nonparametric methods that they can tailor to different empirical settings.
Despite the prevalence of multi-score RD designs, empirical researchers often adopt a univariate approach. For example, unemployment insurance programs usually set multiple eligibility criteria, such as age and work history in determining the maximum duration of benefits, but SVB12 subset the data to consider only individuals with eligible employment histories within their respective age groups, and use age as the effective running variable NW17. A similar subsetting approach is taken by JL04 and M08 to evaluate education policies where treatment eligibility is determined by multiple test scores. Another common approach is to transform the multivariate score into a univariate one, especially when the treatment boundary takes arbitrary shapes and finding a suitable subset where a univariate score applies is infeasible. For example, to exploit geographic discontinuities, B99 considers the shortest distance to the boundary of a geographic region; KT15 transform bivariate coordinates by taking their rescaled Euclidean (chordal) distance to a given boundary point and then run univariate local linear regressions, and BCI14 take a similar Euclidean transformation. Other examples include BM22 that use the average of two scores, and BBRW09 and CM14 that use the minimum of multiple scores.
These approaches, however, have several drawbacks. Considering only subsets of the data leaves unexplored a rich set of treatment effects that could be identified under a valid multivariate RD design; see Section (ref). Heterogeneous treatment effects of this kind can be important to answer policy questions, as exemplified by our empirical application in Section (ref). Besides, subsetting does not fully exploit signals from multiple dimensions and could suffer a loss of precision IW19. Perhaps more precarious is the case where observations with the same distance to the cutoff are collapsed to the same point on the real line: transforming a multivariate score may incur unexpected behaviors of the resulting density of the effective univariate score, violating key assumptions needed for both identification and estimation. In particular, the density of the transformed univariate score may be zero at the new cutoff. In Section (ref), we formally discuss the issue of zero density and provide a characterization of bivariate densities subject to this problem. We also show in simulations that such problem can make local linear regressions biased and undercover the true treatment effect, even when the sample size gets large. Nevertheless, these aforementioned issues that could arise when reducing a multivariate RD design to a univariate one have not received much attention in empirical RD research and are difficult to detect in practice without knowing the true data generating process (DGP).
What should an empirical researcher do in the presence of multiple scores, then? We consider, as alternatives to these existing approaches, two estimators that have gained popularity in estimating heterogeneous treatment effect under unconfoundedness, but have not been fully exploited in the empirical RD literature---honest regression forests of WA18 and local linear forests of FTAW20. Like local linear regressions, these ensemble-tree methods construct a neighborhood (i.e., leaf) for a given cutoff, weighting observations by how frequently they fall into the same leaf as the cutoff across all the trees in the forest. Regression forests return the weighted average of outcomes within that neighborhood as the estimate, whereas local linear forests use these forest weights to run a weighted linear regression to capture smooth signals. Thus, regression forests and local linear forests are, respectively, analogues of local constant and local linear regressions, with forest weights replacing kernel weights. Since the identified conditional treatment effect at a given cutoff is simply the difference between two conditional means, we mimic the standard local linear regression approach by building two forests using treated and control observations, respectively, and estimate the conditional treatment effect by subtracting the control forest estimate from the treated counterpart, both evaluated at the same cutoff. These forest-based estimators are easy to implement, allow valid inference, and can accommodate any fixed number of scores that does not grow with the sample size.
To evaluate these nonparametric estimators, we design a systematic Monte Carlo study with DGPs built both from functional forms that we specify, and, following the proposal of AIMM21, from Wasserstein Generative Adversarial Networks (WGANs) with a high degree of complexity that can closely mimic the observed data. We find that no single estimator dominates across all simulations: (i) local linear regressions perform well in univariate settings, but can undercover when multivariate scores are transformed into a univariate index---a common practice---possibly due to the “zero-density" issue detailed in Section (ref); (ii) the minimax-optimal estimator has good bias-variance properties when the bound on the second derivative of the conditional expectation function (CEF) is accurately estimated, but its performance is sensitive to the estimation of this nuisance parameter and its current implementation only accepts up to two scores; (iii) forest-based estimators are not designed for estimation at boundary points and can suffer from bias in finite sample---though local linear forests improve upon regression forests when the underlying signal exhibits strong linear trends---but their flexibility in modeling multivariate scores opens the door to a wide range of empirical applications in multivariate RD designs.
The rest of the paper proceeds as follows. In Section (ref), we begin with the identification of conditional treatment effects in multivariate RD designs. In Section (ref), we review current estimation and inference practices using local linear regressions and the issue thereof, as well as the minimax-optimal estimator, followed by an introduction to honest regression forests and local linear forests. We present simulation results in Section (ref) and an empirical application in Section (ref). Section (ref) concludes. Appendix (ref) collects the proof of Proposition \hyperref[prop:1]{1}, and Appendix (ref) provides additional information on the Monte Carlo study.
Suppose we have a sample of independent and identically distributed observations $\mathcal{S}_n :=\{(Y_i, X_i)\}_{i=1}^n$, where $Y_i \in \mathbb{R}$ is the outcome of interest and $X_i \in \mathcal{X} \subseteq \mathbb{R}^{d_X}$ is the $d_X$-dimensional running variable that determines the binary treatment status, $D_i\in\{0,1\}$, of the $i$-th unit according to some exogenous treatment assignment rule, $A:\mathcal{X} \to \{0,1\}$. If $d_X=1$, then $D_i=1$ if and only if $A(X_i)=\mathds{1}\{X_i\geq x^c\}=1$, where $x^c$ is the cutoff as in the standard sharp RD design with a univariate score. In higher dimensions, the treatment assignment rule is more complex. For example, KT15 study the effect of television advertising on election turnout, where the treatment status is determined by a geographic media market boundary. Following Z12, let $$ \mathscr{B} := \{x\in\mathcal{X}: \forall \epsilon > 0, \ \exists x_0, x_1 \in \mathrm{Ball}(x, \epsilon) \text{ such that } A(x_0)=0\ \text{and}\ A(x_1)=1\} $$ denote the treatment boundary, i.e., a point $x\in\mathcal{X}$ lies on the treatment boundary if every $\epsilon$-ball around it contains both treated and control scores. In univariate designs, $\mathscr{B}$ is simply a singleton containing the one-dimensional cutoff. Let $\bigl(Y_i(1), Y_i(0)\bigr)$ denote the pair of potential outcomes with and without the binary treatment. We are interested in the following conditional treatment effect at any $x^c\in \mathscr{B}$:
For each unit $i$, we only ever observe one of the potential outcomes. However, if we assume (i) the CEFs $\mathbb{E}\left[Y_i(1)|X_i=x^c\right]$ and $\mathbb{E}\left[Y_i(0)|X_i=x^c\right]$ are continuous at all $x^c \in \mathscr{B}$, and (ii) $X$ has a density, $f_X(\cdot)$, that is bounded away from 0 in some neighborhood around $x^c$ for all $x^c \in \mathscr{B}$, then we can identify $\tau(x^c)$ using neighboring observations around $x^c$:
where $\tau(x^c)$ is identified as the difference between the conditional expectation of the observed outcome for the treated and control populations, respectively, both evaluated at $x^c$.\footnote{A detailed proof can be found in Z12.} Thus, under a valid multivariate RD design, we can identify the conditional treatment effect as a function of scores along the treatment boundary $\mathscr{B}$. Such heterogeneity may be of interest to the researcher. Alternatively, Z12 and KT15 discuss a summary measure of the average treatment effect along the boundary, i.e., $\mathbb{E}[Y_i(1)-Y_i(0)|X_i \in \mathscr{B}]$, which can be estimated given a sample of boundary points by averaging the estimated conditional treatment effects at those sample boundary points. However, as pointed out by NY23, the Lebesgue measure of the boundary $\mathscr{B}$ is typically zero, and hence they consider the Hausdorff measure instead and propose a two-stage least squares (2SLS) estimator for the average treatment effect along the treatment boundary under the Hausdorff measure. In what follows, we focus on estimating conditional treatment effects at a given treatment boundary point $x^c\in\mathscr{B}$ and leave how exactly to aggregate these heterogeneous treatment effects to the empirical researcher.
In this section, we review the current estimation and inference practices using local linear regressions and the minimax-optimal estimator of IW19 in multivariate RD designs. We highlight the problem of the common practice of transforming multivariate scores to a univariate one and then running local linear regressions. We conclude this section with an introduction to the honest regression forests of WA18 and the local linear linear forests of FTAW20, and discuss how they can be leveraged to estimate heterogeneous treatment effects in RD designs with multiple scores.
Given a boundary point $x^c \in \mathscr{B}$, estimation of $\tau(x^c)$ boils down to estimating two CEFs,
Consider univariate designs first. The standard estimation approach is to run a local linear regression on each side of the cutoff using some weighting kernel $K_h(\cdot)$ with bandwidth $h$,
where in (ref), the subscript “llr” refers to local linear regression, $\bm{e}_1:=[\,1\quad 0\,]^\intercal$ is the standard basis in $\mathbb{R}^2$ that picks out the first element of the minimizer (i.e., the fitted intercept term $\beta_0$, as the slope $\beta_1$ goes away when we evaluate the fitted regression at $x^c$); $X_+$ is the $n_1\times 2$ design matrix, where $n_1$ is the number of treated observations. The first column of $X_+$ contains all ones and the second column stacks all $\{X_i-x^c: X_i\geq x^c\}$. $W(X_+,x^c;h)$ is the $n_1\times n_1$ diagonal weighting matrix with the $i$-th diagonal element being the kernel weight $K_h(X_i-x^c)$ for $X_i\geq x^c$. $Y_+$ is the $n_1\times 1$ outcome vector that stacks the outcomes for all treated observations. Notations in (ref) are defined analogously for the control observations. The treatment effect at $x^c$ is estimated as
Following IK12, the choice of the bandwidth $h$ is usually obtained by minimizing the asymptotic expansion of
with respect to $h$. However, CCT14 note that the MSE-optimal bandwidth obtained from optimizing (ref), denoted as $h^*$, is oversmoothing in the sense that $n(h^{{*}})^{5}$ does not converge to 0, leading to a bias term that must be corrected in order to conduct valid inference. They propose subtracting from $\widehat{\tau}_{\texttt{llr}}(x^c;h^*)$ an estimated bias term and then adjusting the standard error to account for the uncertainty in the bias estimation. This procedure, however, does not readily generalize to multivariate RD designs. Although Z12 proposes a similar MSE-optimal bandwidth selector in the bivariate case, valid inference still requires bias correction due to oversmoothing, which is technically challenging beyond one dimension and to our knowledge there is no such proposal in the existing literature. In empirical research involving multiple scores, the multivariate problem is usually reduced to one dimension so that univariate methods can apply.
As mentioned earlier, common practices of reducing a multivariate RD design to a univariate one include subsetting the data and reducing the multivariate score to one dimension, and each has its own pitfalls. In particular, transforming a multivariate score in some cases incurs unexpected behaviors of the resulting density of the effective univariate running variable. In this section, we detail a situation that may arise under such transformation: the transformed univariate score may have zero density at the new cutoff, violating a key assumption for both the identification and estimation of treatment effects in RD designs.
We illustrate this problem for a transformation used by KT15, who transform a bivariate score to a univariate one using rescaled Euclidean (chordal) distance to a given treatment boundary point and then estimate the treatment effect at the new cutoff $x^c=0$ using local linear regression with robust bias correction as in CCT14. The following proposition gives a general characterization of bivariate densities that have a zero density at the effective univariate cutoff $x^c=0$, once the bivariate score is collapsed to a univariate variable by taking the Euclidean distance between the bivariate score and a given boundary point.
Proposition 1. Let $(X_1,X_2)$ be a continuous bivariate random variable with joint density $f_{(X_1,X_2)}(x_1,x_2)$ over the support $\Omega(X_1,X_2)$. Given any point $(c_1,c_2) \in \Omega(X_1,X_2)$, consider the bivariate transformation $(X_1,X_2)\mapsto(E,V)$ where $E:=\sqrt{(X_1-c_1)^2+(X_2-c_2)^2}$ and $V:=X_1-c_1$ with joint support
where $l(e)$ and $u(e)$ are the lower and upper bound of $v$ that can depend on $e$.\footnote{To simplify notation, we assume without loss of generality that the lower and upper bounds of $(e,v)$ are attained so that the joint support of $(E,V)$ is a closed subset of $\mathbb{R}^2$; otherwise the closed brackets need to be changed to open brackets accordingly.} Then the Euclidean distance $E$ to the given point $(c_1, c_2)$ has zero density at the new cutoff $e=0$, $f_E(0)=0$, if and only if
}
Intuitively, such transformation involves a loss of information: many distinct multivariate realizations can map to the same distance on the real line. In particular, for the collapsed variable to be zero, all coordinates of the original variable must simultaneously have zero distance, which is a rare event. We give the proof of Proposition \hyperref[prop:1]{1} and two examples under which (ref) holds in Appendix (ref). This characterization depends on the population joint density of the scores and is specific to the transformation chosen by the researcher.
{
}
To illustrate, the first panel of Figure (ref) is the true density of a transformed bivariate score that is uniformly distributed over $[-1,1]^2$, and the point to which we calculate the Euclidean distance is $(0,0)$; see Example \hyperref[appendix:eg]{A.1} for the derivation of this density. As in KT15, we sign the Euclidean distance so that treated (control) observations receive a positive (negative) transformed score, and the new cutoff after transformation becomes $0$, at which the density of the transformed score is exactly $0$. In practice, the researcher has access to only finite samples, and problems of this kind are difficult to detect without knowing the true DGP. In the last two panels of Figure (ref), we plot the raw histograms of the transformed bivariate score with respect to a given treatment boundary point using the exact same transformation employed by KT15; each histogram uses data simulated from a different DGP built using the real data from KT15, and thus the simulated bivariate scores do not follow any known distributions.\footnote{We provide details on how the simulation data is generated in Section (ref).} In all three panels, the empirical density of the transformed bivariate score, be it true or empirical, is zero at the effective univariate cutoff $x^c=0$. When such zero-density problem arises, not only is the conditional treatment effect at the transformed cutoff not identified because the conditional expectations are not defined, the problem is also reflected in estimation, especially for very large sample sizes, as we show in the simulations in Section (ref).
To avoid the problem of choosing a multivariate bandwidth, IW19 propose an estimator based on numerical optimization that is minimax-optimal among all estimators that are linear in the outcome variable over all problems with finite $\sigma_i^2:=Var[Y_i|X_i]$ and a second derivative bound $\mathcal{B}<\infty$ on the CEF of potential outcomes. Specifically, their estimator for $\tau(x^c)$ takes the form
where
for $\mu_{A_i}(X_i)=\mu^+(X_i)$ if $A_i:=A(X_i)=1$ and $\mu_{A_i}(X_i)=\mu^-(X_i)$ otherwise. The norm $\lVert\cdot\rVert$ denotes the Euclidean norm when $d_X=1$ and the operator norm when $d_X>1$. In words, the estimator ((ref)) is a weighted average of the outcomes, where the weights $\widehat{\gamma}(x^c;\mathcal{B})$ are constructed by minimizing the resulting variance plus the worst-case bias over a class of CEFs with second derivatives bounded by $\mathcal{B}$. Inference follows from obtaining confidence intervals $\mathcal{I}_{\alpha}$ that achieve uniform coverage over this class of CEFs:
Implementing the minimax-optimal estimator ((ref)) requires estimates of $\sigma_i^2$ and $\mathcal{B}$. IW19 recommend estimating $\sigma_i^2$ by averaging the residuals from an ordinary least-squares regression of $Y_i$ on the interaction of $X_i$ and $A(X_i)$. Estimating $\mathcal{B}$, however, is problem-specific and less straightforward. IW19 recommend estimating the CEFs globally using second-order polynomials and then multiplying the maximal estimated curvature by 2 or 4. In one of our simulations in Section (ref) using data from L08, we find that using second-order polynomials yields estimates of $\mathcal{B}$ that are too small, even after multiplying by $4$; only after multiplying the maximal estimated curvature by $30$ do we get an estimate that is close to the true (least upper) bound $\mathcal{B}$.
We now introduce two additional theoretically valid estimators that can be leveraged in multivariate RD designs but have been largely overlooked in the literature---the honest regression forests of WA18 and the local linear forests of FTAW20. Both are variants of random forests originally proposed by B01 and have received popularity in estimating heterogeneous treatment effects under unconfoundedness.
A regression forest evaluated at some point $x \in \mathcal{X}$, denoted as $RF(x)$, is an average of $B$ trees, $\{T_b(x)\}_{b=1}^B$, each fitted using a bootstrapped sample or a random subsample of $\mathcal{S}_n:=\{(Y_i, X_i)\}_{i=1}^n$. Following ATW19, FTAW20, and WA18, we focus on forests built using random subsamples, for which central limit theorems are established. Given a subsample of size $s$ from $\mathcal{S}_n$, $\{(Y_{i_k}, X_{i_k})\}_{k=1}^s$, a regression tree recursively divides the sample into mutually exclusive and collectively exhaustive convex subsets called leaves. For a point $x \in \mathcal{X}$, let $L(x)$ denote the leaf that contains $x$. A tree $T$ estimates the CEF, $\mathbb{E}[Y_i|X_i=x]$, by averaging the outcomes of observations that fall into the same leaf as $x$:
Similar to selecting an MSE-optimal bandwidth for local linear regressions, we look for partitions of the feature space that are MSE-optimal. Unlike local linear regressions, trees do not target the MSE at a particular treatment boundary point, nor do they have an analytical form for the asymptotic expansion of the overall MSE,
and searching for all possible partitions is computationally infeasible. Instead, standard splitting rules rely on greedy algorithms that immediately improve the in-sample fit of leaf-wide averages after each split. As such, estimates from a single tree can be highly unstable without regularization: in the extreme case, one can achieve perfect in-sample fit by partitioning the sample data so that each leaf contains only one observation.
Regression forests stabilize predictions from individual trees by averaging many of them: for a collection of $b=1,...,B$ trees, we fit each tree $T_{b}(\cdot)$ using a random subsample of size $s$, and obtain the regression forest estimate by averaging the predictions across all $B$ trees, $\frac{1}{B}\sum_{b=1}^B T_{b}(x)$. To achieve a sizable reduction in variance, we may want to decorrelate the trees as much as possible by, for example, randomly selecting a dimension to split at each node. Following WA18, we assume throughout that $B$ is large enough so that the Monte Carlo variability of which subsamples are drawn is negligible and adopt their definition of regression forests given the training sample $\mathcal{S}_n$ and subsample size $s$:
which averages across trees fitted using all possible draws of a size-$s$ subsample, marginalized over an auxiliary random term $\xi$ that serves as a decorrelation tool encoding, for example, the random split information. To establish asymptotic normality of regression forests, WA18 leverage a specific subsampling scheme satisfying a property called honesty to reduce bias, which requires that the subsample used to make the splits is independent of the subsample used to estimate the within-leaf conditional means. We refer the reader to WA18 for a discussion on honesty and other regularity conditions needed to establish asymptotic normality of honest regression forests.
We now define the estimator based on honest regression forests for the conditional treatment effect at some point $x^c \in \mathscr{B}$ along the treatment boundary, $\tau(x^c)$:
Just like estimation using local linear regressions, $\widehat{\tau}_{\texttt{rf}}(x^c)$ takes the difference between a regression forest fitted using observations that fall into the treated region and a forest fitted using observations that fall into the control region, both evaluated at the same treatment boundary point $x^c\in\mathscr{B}$.
\noindentRemark 1. Although it seems that the local linear estimator ((ref)) depends on one more parameter, $h$, than the regression forest estimator ((ref)), we note that this is not the case. Each tree $T$ in the regression forest implicitly defines a neighborhood $L(x^c)$ for the treatment boundary point $x^c$. These neighborhoods are then aggregated, producing a weight for each observation that is proportional to the fraction of time it falls into the same leaf as $x^c$ across all trees in the regression forest. These forest weights are analogous to the weights generated by a kernel with bandwidth $h$. Thus ((ref)) implicitly depends on some “pseudo-bandwidth” generated by tree partitions across the forest. Analogously, the $\gamma$ weights of the minimax-optimal estimator ((ref)) play the role of kernels and bandwidths. We can equivalently express regression forests as solutions to weighted least-squares constant regressions:
where the elements in ((ref)) and ((ref)) are analogous to the terms in ((ref)) and ((ref)), except $X_+$ and $X_-$ are now $n_1\times 1$ and $n_0\times 1$ vectors of ones, denoted respectively by $\mathbf{1}_{n_1}$ and $\mathbf{1}_{n_0}$, and we replace the kernel weights with forest weights,
where $L^+_b(x^c)$ denotes the leaf containing $x^c$ in the $b$-th tree fitted using the treated observations; $W_{\texttt{rf}}^-$ is defined similarly for the control observations. More generally, the least-square representation of forests in ((ref)) and ((ref)) effectively views $\widehat{\mu}_{\texttt{rf}}^+(x^c)$ and $\widehat{\mu}_{\texttt{rf}}^-(x^c)$ as solutions to $\mathbb{E}\left[(Y_i-\theta)^2|X_i=x^c, D_i=1\right]=0$ and $\mathbb{E}\left[(Y_i-\theta)^2|X_i=x^c, D_i=0\right]=0$, and solves for $\theta$ using their empirical analogs with forest weights reflective of the importance of each observation $i$ in terms of minimizing the objective function. ATW19 use forests as weighting kernels to solve parameters identified by local moments in a more general setting.
Though a powerful nonparametric tool and popular in applications, regression forests resemble local constant regressions with convex weights, and hence cannot extrapolate and are known to have limited ability to capture smoothness or linear trends in the CEF, as noted by FTAW20 and the references therein, especially at points near the boundary of the support where we observe data asymmetrically from only one side of a given boundary point. To improve the performance of forests in the presence of strong linear signals, FTAW20 build on the honest regression forests of WA18 and propose a linear regression adjustment by running a ridge local linear regression with forest weights, under which we have the following estimator for $\tau(x^c)$:
where, similar to the expressions in (ref) and (ref) for local linear regressions, $\bm{e}_1$ is the standard basis in $\mathbb{R}^{(d_X+1)}$ with the first element being $1$ and $0$ elsewhere to pick out the fitted intercept $\beta_0$; $W_{\texttt{rf}}^+(X_i,x^c)$ and $W_{\texttt{rf}}^-(X_i,x^c)$ are the forest weights defined in (ref),\footnote{To account for the fact that local linear forests will use linear regression adjustments to estimate conditional means (as opposed to using simple weighted means, as in the case of regression forests), FTAW20 also modify the standard splitting rules described in Section (ref) to “ridge residual splitting;" see Section 2.1 of FTAW20.} effectively replacing the kernel weights in (ref) and (ref); $\lambda\|\beta_1\|_2^2$ is the ridge penalty term that prevents overfitting to local trends, and as explained in FTAW20, is needed both to improve finite-sample performance and for proving convergence rates to establish a central limit theorem for local linear forests. Just like local linear regressions and regression forests, local linear forests (ref) and (ref) can also be expressed as solutions to weighted least-square problems; see Equation (5) of FTAW20.
The assumptions needed for the asymptotic normality of local linear forests largely overlap with those required for honest regression forests to be asymptotically normal. However, with an additional assumption on the smoothness of the CEF $\mathbb{E}[Y_i|X_i=x]$---that it is differentiable with a Lipschitz continuous derivative, as opposed to the WA18 condition that only requires it to be Lipschitz continuous---FTAW20 show an improved convergence rate of local linear forests relative to the rate of regression forests, effectively exploiting the additional smoothness assumption with a linear adjustment. We refer the reader to WA18 and FTAW20 for details on the technical assumptions and discussions thereof. Pointwise confidence intervals for both honest regression forests and honest local linear forests are constructed using variance estimates from the bootstrap of little bags algorithm of SL09, with consistency results established by the random forest delta method of ATW19. In our Monte Carlo study, we find that local linear forests are in general less biased than regression forests and have better coverage rates whenever the simulation DGP exhibits strong linear trends.
\noindentRemark 2. One caveat of using local linear forests in multivariate RD designs is that the point at which we evaluate local linear forests, $x^c\in\mathscr{B}$, is necessarily a boundary point of the support of $X_i$ for both treated and control populations. On the other hand, the central limit theorem of FTAW20 only applies to points in the interior of the support, as the proof starts with a usual Taylor expansion of $\mathbb{E}[Y_i|X_i=x]$ around $x^c$.\footnote{Asymptotic normality of local linear regressions also relies on a similar Taylor expansion, but it is the CEFs of potential outcomes that get expanded (as opposed to the CEFs of observed outcomes of either the treated or the control population), and so in this case $x^c\in\mathscr{B}$ is an interior point of the support of $X$.} We sidestep this problem by adding an arbitrarily small “buffer" to $x^c\in\mathscr{B}$: we pick another point $\Tilde{x}^c$ in the interior of the treated score region such that $\|\Tilde{x}^c-x^c\|=\epsilon$ for an infinitesimal $\epsilon>0$ and evaluate the treated local linear forest $\widehat{\mu}^+_{\texttt{llf}}$ at $\Tilde{x}^c$ to approximate $\widehat{\mu}^+_{\texttt{llf}}(x^c)$; we do the same for the control local linear forest $\widehat{\mu}^-_{\texttt{llf}}$. Under the assumptions in FTAW20, both $\mu^+(x^c)$ and $\mu^-(x^c)$ are Lipschitz in $x^c$ and therefore very close to $\mu^+(\Tilde{x}^c)$ and $\mu^-(\Tilde{x}^c)$, respectively, for an infinitesimal $\epsilon$. Ideally, one should adjust the confidence interval accordingly to account for such approximation. We leave this for future research, and in our simulations and empirical application we choose $\epsilon=10^{-30}$ so that such correction is essentially negligible up to machine precision.
In this section, we conduct a systematic Monte Carlo study to evaluate the two existing estimation and inference practices in multivariate RD designs, as well as the two forest-based estimators that have not received much attention in the empirical RD literature before. All replication files can be accessed at \href{https://github.com/yqi3/GRF-RD}{github.com/yqi3/RDForest}.
We first consider a univariate design simulated using data from L08, which is a popular choice in the RD literature to conduct simulations. L08 studies the effect of being the current incumbent party in a district on the votes obtained in the district’s next election. The running variable is univariate and measures the difference in vote share between the Democratic candidate and their strongest opponent in the current election, and the outcome variable is the Democratic vote share in the next election. This univariate design normalizes the assignment cutoff to $0$, and units with $X_i\geq 0$ are treated (i.e., winning the current election).
We then consider a bivariate design simulated using data from KT15, who study the effect of television advertising on election turnout, where the running variable is the two-dimensional geographic coordinate of a voter's location and the outcome of interest is election turnout and takes binary values. Treatment status is determined by a geographic media market boundary. Their data also contains several other outcome variables for placebo analysis. In addition to election turnout, we include age and housing price per square foot in our simulation to cover both discrete and continuous outcome variables. KT15 consider three equally spaced points along the treatment boundary, depicted as black squares in Figure (ref), to be the points at which they estimate the treatment effect. For simplicity, we pick the middle point, circled in green in Figure (ref), to be the boundary point of interest in our simulation.
We consider designs up to two scores in our Monte Carlo study because the current R implementation of the minimax-optimal estimator $\widehat{\tau}_{\texttt{mm}}(x^c;\mathcal{B})$ in (ref) provided by IW19 only accepts up to two scores. In an empirical application in Section (ref), we illustrate how to use the forest-based estimators in an RD design with more than two scores to explore heterogeneous treatment effects along the treatment boundary.
We consider data generating processes built both from functional forms that we specify, and, following the proposal of AIMM21, from Wasserstein Generative Adversarial Networks (WGANs) with a high degree of complexity that can closely mimic the observed data. We first fit the observed data with a polynomial: IK12, CCT14, and CCFT19 all use the L08 data to generate simulation data by fitting two $5$-th order global polynomials separately for the treated and control populations, which we also use in our simulation. For the KT15 data, we fit a $3$-rd order interacted polynomial for continuous outcomes (i.e., age and housing price) and a logistic regression using a $3$-rd order interacted polynomial inside the logit link function for the binary outcome turnout, separately for the treated and control populations. Details on these polynomial DGPs can be found in Appendix (ref).
As argued by AIMM21, Monte Carlo studies are more credible when the data used to evaluate the performance of econometric methods does not come from researcher-chosen functional forms. Following their proposal, we also consider simulation data generated from WGAN, which trains two neural networks---one called the generator and the other called the critic---in an adversarial setting to minimize the Wasserstein distance between a real dataset and the artificial dataset generated from the neural nets. We refer the reader to AIMM21 for an introduction to WGAN.
Using WGAN to generate data to evaluate estimators for treatment effects, however, is task-specific, since we never directly observe the treatment effect in the data, but rather we “identify" it under a set of identification assumptions. AIMM21 provide an example of how to generate simulation data to evaluate estimators of treatment effects under unconfoundedness, where they train two conditional WGANs to learn $X_i|D_i$ and $Y_i|(X_i, D_i)$ from the real sample. They then draw an artificial sample of treatment status and covariates from the WGAN-generated distribution of $X_i|D_i$, and for each $(X_i, D_i)$ drawn, an outcome is drawn from the WGAN-generated distribution of $Y_i|(X_i, D_i)$. Under unconfoundedness, $Y_i|(X_i, D_i=d)$ identifies $Y_i(d)|X_i$ and hence for each $Y_i$ drawn, its counterfactual outcome is generated by mutating the corresponding treatment status $D_i$ to $(1-D_i)$ but keeping the same $X_i$. In other words, the counterfactual outcome is drawn from the WGAN-generated distribution of $Y_i|(X_i, 1-D_i)$. The true average treatment effect is simulated by averaging the difference between the pair of generated outcome and its counterfactual over a large WGAN-generated sample.
Our setting, with a different design and identification assumption for the treatment effect of interest, requires a different way of using WGAN to generate suitable simulation data. Specifically, the treatment effects are only identified for $X_i\in\mathscr{B}$, so we need to simulate the treatment effect at each treatment boundary point of interest, instead of simply taking averages as in the unconfoundedness setting. To do so, we first train two conditional WGANs to learn $X_i|D_i$ and $Y_i|(X_i, D_i)$ from the real sample, as in AIMM21. But for a given treatment boundary point $x^c\in\mathscr{B}$, we draw a separate artificial sample to simulate the true treatment effect at $x^c$: after training the WGAN to learn the distribution of $Y_i|(X_i, D_i)$, we generate 100,000 outcomes from $Y_i|(X_i=x^c, D_i=1)$ and another 100,000 outcomes from $Y_i|(X_i=x^c, D_i=0)$, and take the average of their differences as the true treatment effect at $x^c$.\footnote{We note that WGANs do not directly yield closed-form estimates of the fitted distribution, but rather generate a new, artificial sample of arbitrary size that closely mimics the real sample. This is why we need to simulate the true treatment effect at a given treatment boundary point $x^c\in\mathscr{B}$, unlike in the polynomial case where we can directly calculate the true treatment effects from the coefficients of the fitted polynomials evaluated at $x^c$.} Figures (ref) and (ref) visualize the actual data and the simulated data from the polynomial specification and WGAN. In Appendix (ref), we include additional descriptions of the WGAN-generated DGPs in Figures (ref)--(ref), and we report in Table (ref) the Wasserstein distance between the real and WGAN data. We use the default options available from the wgan package of AIMM21 to train the WGANs; see Appendix (ref) for details. Following Section 5 of AIMM21, we conduct a thorough robustness check to assess the variability of the Monte Carlo results with respect to sampling variation and WGAN tuning parameters, with robustness results reported in Appendix (ref).
We report the simulation results in Table (ref). All simulations are implemented in R, where we use the rdrobust package of CCT14 for local linear regressions, the optrdd package of IW19 for the minimax-optimal estimator, and the grf package of ATW19 for the two forest estimators---all with default options without tuning, except we specify pilot values for a few hyperparameters of the two forest estimators, for reasons explained in Appendix (ref).
Local linear regressions perform well in univariate settings, although they slightly undercover in the polynomial DGP (Panel A, Columns 1-5 of Table (ref)).\footnote{This undercoverage is consistent with the simulation result ($90.6\%$ coverage) of CCFT19, where they use the same polynomial DGP but with $n=1000$ for a total of $5000$ replications; see Table SA-1, Row 2 and Column 3 in their supplemental appendix. We have done our simulation using their STATA code and obtained similar results to those reported in Table (ref).} However, they are generally noisier in bivariate settings, where we follow KT15 to first collapse the bivariate score using rescaled Euclidean (chordal) distance and then run rdrobust, especially in small samples. Although the variance of local linear regressions shrinks as the sample size gets large, their bias does not shrink fast enough relative to the diminishing variance in those bivariate WGAN DGPs (Panels F-H, Columns 1-5 of Table (ref)), and therefore the coverage rate drops as the sample size increases. This problem arises possibly due to the zero-density issue of the collapsed multivariate score at the new univariate cutoff, as we point out in Section (ref) and empirically illustrate in the last two histograms of Figure (ref), cautioning against such univariate transformations.
The minimax-optimal estimator of IW19 has a relatively low bias but tends to overcover, which is expected since it is designed for minimizing the worst-case MSE and inference aims at achieving uniform coverage over a class of CEFs explained in Section (ref). Nevertheless, its performance is contingent on how precisely the second derivative bound $\mathcal{B}$ is estimated: for the two artificial datasets built using the L08 data (Panels A-B of Table (ref)), the minimax-optimal estimator undercovers severely, especially for the polynomial DGP in Panel A. There, we follow the recommendation of IW19 to estimate the CEFs globally using a second-order polynomial and multiply the maximal estimated curvature by $4$. But this turns out to yield an estimate of $\mathcal{B}$ that is too small; Table (ref) in Appendix (ref) shows how the simulation results change with the constant multiplied by $\mathcal{B}$, and it is only after multiplying the maximal estimated curvature by $30$ do we get an estimate that is close to the true (least upper) bound $\mathcal{B}$.
Of course, one could argue that estimating $\mathcal{B}$ by a global second-order polynomial is a bad choice to start with. For the bivariate design, we follow their working paper IW19w and estimate $\mathcal{B}$ using the $95$-th percentile (over all sample points) of the operator norm curvature via a cross-validated ridge regression with interacted $7$-th order natural splines, which yields good coverage (Panels C-H, Columns 6-10 of Table (ref)). However, without knowledge about the true DGP, choosing a “good” model is not straightforward, and estimating second derivatives is more difficult when the dimension of the running variable increases. In addition, the current implementation of optrdd accepts up to two scores only, limiting its application more broadly.
Honest regression forests and local linear forests have comparable performance to local linear regressions in the univariate design (Panels A-B, Columns 11-20 of Table (ref)), though their variances are slightly higher. On the other hand, regression forests are generally more biased than local linear forests, especially in small samples, which is a known boundary problem for regression forests. In addition, the variance of regression forests tends to be smaller than that of local linear forests, exacerbating the undercoverage problem. As we discussed in Section (ref), local linear forests are better at capturing smooth signals via a linear correction, and its advantage over regression forests is especially pronounced when the underlying DGP exhibits strong linear trends near the boundary point at which we estimate the treatment effect. In Appendix (ref), we examine the robustness of the results in Table (ref) to the number of training epochs, which is a tuning parameter of WGAN, and increasing the number of training epochs leads to stronger linear trends near the boundary point (see Figure (ref)); in these cases, the performance of regression forests deteriorates, whereas local linear regressions have robust performance. With that said, these two forest-based methods are theoretically valid alternatives to existing methods for multivariate RD designs and can flexibly model any given number of scores, and empirical researchers may leverage and tailor them accordingly to fit different empirical settings.
We revisit the empirical analysis in NY23, who study the effect of hospital relief funding during COVID-19.
The pandemic inflicted substantial financial pressure on U.S. hospitals, compounding revenue losses with increasing operational expenses. To alleviate this, the U.S. government enacted the Coronavirus Aid, Relief, and Economic Security (CARES) Act in March 2020 that provides emergency financial support to healthcare providers. A \$10 billion portion of this funding was allocated to “safety net" hospitals, as determined by a set of eligibility criteria that target hospitals serving vulnerable populations. Hospitals were eligible for safety net funding if they met all of the following three conditions:\footnote{Source: \href{https://www.hrsa.gov/sites/default/files/hrsa/coronavirus/provider-relief-fund-faq-complete.pdf}{Provider Relief Fund Frequently Asked Questions} from the Health Resources and Services Administration, U.S. Department of Health and Human Services.}
To determine safety net funding eligibility, we use data from the Healthcare Cost Report Information System (HCRIS) for the 2019 financial year TL-303. The 2019 measurements are the closest available approximation of the values hospitals likely reported when they applied for funding in 2020.\footnote{Although using the 2020 data may seem more temporally aligned with the application timeline, it is unclear whether these records were actually used by the hospitals at the time of application. Since the coverage of the 2019 data was low at the time of NY23's study, they use the 2018 HCRIS data instead. We thank them for their suggestion to use the updated data and help with data preparation.} For the outcome variable, following NY23, we use the number of confirmed or suspected COVID patients hospitalized in each hospital during the week spanning August 2 to August 8, 2020, drawn from the most recent (June 2024) version of the COVID Reported Patient Impact and Hospital Capacity by Facility dataset HHS2024. We note that, due to changes in the reporting format, our weekly outcome is based on a slightly different weakly reporting window from that of NY23, which starts from July 31 to August 6, 2020 using an earlier version of the outcome dataset.
We follow the steps in Appendix H of NY23 for data cleaning and missing value imputation.\footnote{In addition, we remove one hospital that appears twice in the outcome data, as the two occurrences have considerably different values for the outcome variable.} The resulting treated and control sample sizes are $n_1 = 656$ and $n_0 = 3$,199, respectively. Figure (ref) visualizes the sample data as a scatter plot, where red (black) colored dots correspond to treated (control) hospitals; the treatment boundary is consisted of three partially bounded planes colored in pink. Under the assumptions in Section (ref), we can identify the conditional treatment effect of funding eligibility at any point on the pink surface.
Safety net funding eligibility criteria create a three-dimensional RD setup, where our framework applies. Acknowledging that directly estimating the conditional treatment effects nonparametrically at each treatment boundary point $X_i$ “is hard to use in practice, however, when $X_i$ has many elements (p. 12)," NY23 propose a 2SLS estimator that is consistent for a weighted average of conditional treatment effects, nesting both propensity score matching and regression discontinuity in the case where the overlap condition fails in the first stage. In their empirical application, they use funding eligibility as an instrument for the amount of funding received and estimate the average effect of received funding amount on multiple outcomes, including the number of COVID patients hospitalized. Their findings suggest little to no effect. We find a similar null effect, but our analysis views funding eligibility as a binary treatment and directly estimates its causal effect. Moreover, we estimate conditional treatment effects as a function of treatment boundary point, allowing investigation of treatment effect heterogeneity along the treatment boundary. We find the estimated treatment effects, albeit with confidence intervals (CIs) covering zero, exhibit sizable heterogeneity along the treatment boundary.
We apply the honest regression forest and local linear forest estimators in (ref) and (ref) to the CARES Act data and, for ease of illustration, we plot how the treatment effect varies along each of the three running variables, while holding the other two fixed. This corresponds to estimating the treatment effect along the three “edges" of the treatment boundary, shown in blue in Figure (ref). Due to the potential zero-density issue of local linear regressions with a collapsed univariate score, and the inability of the minimax-optimal estimator to accommodate more than two scores, we omit analyses using these two estimators.
For simplicity and better generalizability, we do not tune the hyperparameters of the forest estimators, but instead use the default values and set a few to specific pilot values as detailed in Section (ref). We use 50,000 trees for more stable results, as the grf package is known to be irreproducible across platforms.\footnote{See the grf \href{https://grf-labs.github.io/grf/REFERENCE.html\#forests-predict-different-values-depending-on-the-platform-even-though-the-seed-is-the-same}{package manual}.} Figure (ref) presents the estimation results with 95% CIs. Comparing across the top and bottom panels, we see little qualitative difference between the two estimators, and the CIs for the point estimates generally cover zero, supporting NY23's conclusion that the relief funding was not well targeted toward hospitals most in need. There is, however, noticeable treatment effect heterogeneity as DPP and profit margin vary, suggesting that more vulnerable hospitals---those with higher DPP and lower profit margin---tend to utilize the relief funding more efficiently by providing care to more patients.
This paper studies RD designs with multiple scores. We first highlight the critical but often overlooked zero-density issue with running local linear regressions using a collapsed univariate score, and then explore the use of nonparametric estimators based on forests that directly model multivariate scores. Our simulation results suggest that the zero-density issue, while hard to detect from an empirical dataset, is likely present and would impair the performance of local linear regressions. Regression forests and local linear forests, on the other hand, generally require larger sample sizes to achieve good performance at boundary points, though local linear forests are better at capturing linear trends. Although this paper does not offer a simple, universal solution to multivariate RD problems, it adds two forest-based estimators to empirical researchers' toolkit. Overall, we recommend using the two forest-based methods as the primary approach when the sample size is large and the number of scores exceeds two because of their flexibility and theoretical validity. In cases where the sample size is small, local linear regressions with a collapsed univariate score can serve as a starting point, but researchers should proceed with caution and perhaps compare the results against other nonparametric alternatives introduced in this paper as a way of robustness checks.
\singlespacing