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.
94,163 characters · 28 sections · 70 citation commands
Flexible Covariate Adjustments in Regression Discontinuity DesignsFirst version: July 16, 2021. This version: . We thank Sebastian Calonico, Michal Koles\'ar, Thomas Lemieux, Jonathan Roth, Vira Semenova, Stefan Wager, Daniel Wilhelm, Andrei Zeleneev, and numerous conference and seminar participants for helpful comments and suggestions. We thank Tobias Grob\"olting and Merve \"Ogretmek for excellent research assistance. The authors gratefully acknowledge financial support by the European Research Council (ERC) through grant SH1-77202. The second author also gratefully acknowledges support from the European Research Council ERC through grant SH-1852332. Author contact information: Claudia Noack, Department of Economics, University of Bonn. E-Mail: [email removed]. Website: https://claudianoack.github.io. Tomasz Olma, Department of Statistics, Ludwig Maximilian University of Munich. E-Mail: [email removed]. Website: https://tomaszolma.github.io. Christoph Rothe, Department of Economics, University of Mannheim. E-Mail: [email removed]. Website: http://www.christophrothe.net.
\pagestyle{plain} \onehalfspacing
\onehalfspacing
Regression discontinuity (RD) designs are widely used for estimating causal effects from observational data in economics and other social sciences. The design exploits the fact that in many contexts a unit's treatment status is determined by whether its realization of a running variable exceeds some known cutoff value. For example, students might qualify for a scholarship if their GPA is above some threshold. Under continuity conditions on the distribution of potential outcomes, the average treatment effect at the cutoff is identified by the jump in the conditional expectation of the outcome given the running variable at the cutoff. Estimation and inference methods based on local linear regression are widely used and their properties are by now well understood hahn2001identification, imbens2012optimal, calonico2014robust, armstrong2020simple.
While an RD analysis generally does not require data beyond the outcome and the running variable, additional covariate information can be used to reduce the variance of empirical estimates. A common approach is to include the covariates linearly and without separate localization in a local linear RD regression calonico2019regression. This conventional linear adjustment estimator is consistent without functional form assumptions on the underlying conditional expectations if the covariates are unaffected by the treatment in some appropriate sense. It does not exploit the available covariate information efficiently though, and is also not well-suited for settings with many covariates.
In this paper, we propose a novel class of flexible covariate-adjusted RD estimators. Our approach involves running a standard local linear RD regression after subtracting an (estimated) function of the covariates from the original outcome variable. We characterize the function that leads to the RD estimator with the smallest asymptotic variance and show how this function can be estimated with modern machine learning techniques. We also show that the conventional RD framework gives rise to a special robustness property which implies our final RD estimator is very insensitive to estimation uncertainty about the estimated adjustment function. This implies, for example, that existing methods for bandwidth choice and inference can directly be used with our adjusted outcome variable. Our approach is thus easily implemented with existing software packages.
To motivate our proposed procedure, let $Y_i$ and $Z_i$ denote the outcome and covariates, respectively, of observational unit $i$. Then the conventional linear adjustment RD estimator is asymptotically equivalent to a local linear RD estimator with the modified outcome variable $Y_i -Z_i^\top \gamma_0$, where $\gamma_0$ is a vector of projection coefficients. We consider generalizations of such estimators which replace the linearly adjusted outcome with a flexibly adjusted outcome of the form $Y_i - \eta(Z_i)$, for some generic function $\eta$. Such estimators are easily seen to be consistent for any fixed $\eta$ if the distribution of the covariates varies smoothly around the cutoff in some appropriate sense, which is in line with the notion of the covariates being unaffected by the treatment in a causal sense.\footnote{We note that the adjustment function must be the same for observations on either side of the cutoff, as using different adjustments on either side of the cutoff would generally yield inconsistent RD estimates calonico2019regression.} We show that the asymptotic variance in this class of estimators is minimized over $\eta$ by the average of the two conditional expectations of the outcome variable given the running variable and the covariates just above and below the cutoff. This optimal adjustment function $\eta_0$ is generally nonlinear and unknown in practice but can be estimated from the data.
Our proposed estimators hence take the form of a local linear RD regression with generated outcome $Y_i -\widehat\eta(Z_i)$, where $\widehat\eta$ is some estimate of $\eta_0$ obtained in a preliminary stage. We implement such estimators with cross-fitting chernozhukov2018double, which is an efficient form of sample splitting that removes the bias induced by overfitting in the first stage and can accommodate a wide range of estimators of the optimal adjustment function. In particular, one can use modern machine learning methods like lasso regression, random forests, deep neural networks, or ensemble combinations thereof, to estimate the optimal adjustment function. In low-dimensional settings, researchers can also use classical nonparametric approaches like local polynomials or series regression, or estimators based on parametric specifications.
Our theory does not require that $\eta_0$ is consistently estimated for valid inference on the RD parameter. We only require that in large samples the first-stage estimates concentrate in a mean-square sense around some deterministic function $\bar\eta$, which could in principle be different from $\eta_0$. The rate of this convergence can be arbitrarily slow. Our setup can allow for this kind of potential misspecification because our proposed RD estimators are highly insensitive to estimation errors in the preliminary stage. Specifically, they are constructed as sample analogues of a moment function that contains $\eta_0$ as a nuisance function, but does not vary with it: as discussed above, our parameter of interest is equal to the jump in the conditional expectation of $Y_i - \eta(Z_i)$ given the running variable at the cutoff for any fixed function $\eta$. This insensitivity property is related to Neyman orthogonality, which features prominently in many modern two-stage estimation methods chernozhukov2018double, but here it is a global rather than a local property and has therefore substantially stronger implications.\footnote{A moment function is Neyman orthogonal if its first functional derivative with respect to the nuisance function is zero. In contrast, the (conditional) moment function on which our estimates are based is fully invariant with respect to the nuisance function. chernozhukov2018double give examples of setups with (unconditional) moment functions in which such a global insensitivity property occurs. These include optimal instrument problems, certain partial linear models, and treatment effect estimation under unconfoundedness with known propensity score. This property also occurs generally if one of the two nuisance functions in a doubly robust moment condition robins01rotnitsky is known.}
Our theoretical analysis shows that, under the conditions outlined above, our proposed RD estimator is first-order asymptotically equivalent to a local linear “no covariates” RD estimator with $Y_i -\bar\eta(Z_i)$ as the dependent variable. This result is then used to study its asymptotic bias and variance, and to derive an asymptotic normality result. The asymptotic variance of our estimator depends on the function $\bar\eta$ and achieves its minimum value if $\bar\eta = \eta_0$ (that is if $\eta_0$ is consistently estimated in the first stage), but the variance can be estimated consistently irrespective of whether or not that is the case. As our result does not require a particular rate of convergence for the first step estimate of $\eta_0$, our RD estimator can be seen as shielded from the “curse of dimensionality” to some degree, and can hence be expected to perform well in settings with many covariates. We show that these results also extend to fuzzy RD designs in a straightforward manner.
Practical issues like bandwidth choice and construction of confidence intervals can be addressed in a straightforward manner. Specifically, we show that standard methods remain valid if they are applied to a data set in which the outcome $Y_i$ is replaced with the generated outcome $Y_i -\widehat\eta(Z_i)$ and one ignores that $\widehat\eta$ has been estimated. Our approach can therefore easily be integrated into existing software packages.
Our theoretical results are qualitatively similar to those that have been obtained for efficient influence function (EIF) estimators of the population average treatment effect in simple randomized experiments with known and constant propensity scores wager2016high. Such parallels arise because EIF estimators are also based on a moment function that is globally invariant with respect to a nuisance function. In fact, we argue that our RD estimator is in many ways a direct analogue of the EIF estimator, and that the variance it achieves under the optimal adjustment function is similar in structure to the semiparametric efficiency bound in simple randomized experiments.
Our proposed flexible covariate adjustments can lead to substantial efficiency gains in practice. To illustrate this, we collect data from empirical papers recently published in leading journals that use RD estimation with covariates. In total, we reanalyze 56 specifications from 16 papers, and study how different types of covariate adjustments affect confidence interval lengths. Including covariates in the RD regression does not meaningfully reduce the length of the confidence intervals in about half of the specifications we consider irrespective of the specific method used for adjustment, but our proposed flexible adjustments also achieve a reduction of more than 30% in one setting. To put this into perspective, obtaining this reduction would require roughly increasing the sample size by a factor of 2.4 if the covariates were not used. We also observe that linear adjustments alone are often unable to exhaust all the available covariate information: the largest reduction in the confidence interval length from using our flexible adjustment relative to linear adjustments exceeds 20%.
We also conduct simulations based on the data set from one of the papers from our empirical literature survey. In order to cover all types of settings from our empirical literature survey, we consider simulation setups of large sample sizes and a moderate number of covariates as well as small sample sizes and a varying number of covariates. Our proposed RD estimators perform very well in all these settings, in the sense that their standard errors are close to their standard deviations and the associated confidence intervals have simulated coverage rate close to the nominal one.
Our paper contributes to an extensive literature on estimation and inference in RD designs; see, e.g., imbens2008regression and lee2010regression for surveys, and cattaneo2019practical for a textbook treatment. Different ad-hoc methods for incorporating covariates into an RD analysis have long been used in empirical economics lee2010regression. Following calonico2019regression, it has become common practice to include covariates without localization into the usual local linear regression estimator. Our approach nests this estimator as a special case, but is generally more efficient. Other closely related papers are kreiss2021regression and arai2024regression, who extend the approach in calonico2019regression to settings with high-dimensional covariates under sparsity conditions using the lasso. In contrast, our approach allows for flexible use of other machine learning methods in the spirit of double-debiased machine learning of chernozhukov2018double. Moreover, even if one commits to lasso-based adjustments, there are two ways in which our approach can improve upon the methods of kreiss2021regression and arai2024regression. First, we propose a different variant of (post-) lasso adjustments that can be more stable in finite samples (see “global adjustments” in Section (ref)). Second, cross-fitting yields a more precise standard error even if the number of selected covariates is not small relative to the effective sample size. froelich2019including propose to incorporate covariates into an RD analysis in a fully nonparametric fashion, but their approach is generally affected by the curse of dimensionality, and is thus unlikely to perform well in practice.
Our results are also related in a more general sense to the vast literature on two-step estimation problems with infinite-dimensional nuisance parameters andrews1994asymptotics,newey1994variance, especially the recent strand that exploits Neyman orthogonal (or debiased) moment functions and cross-fitting belloni2017program, chernozhukov2018double. The latter literature focuses mostly on regular (root-$n$ estimable) parameters, while our RD treatment effect is a non-regular (nonparametric) quantity. Some general results on non-regular estimation based on orthogonal moments are derived in chernozhukov2019double, and specific results for estimating conditional average treatment effects in models with unconfoundedness are given, for example, in kennedy2017nonparametric, kennedy2020optimal and fan2020estimation. Our results are qualitatively different because, as explained above, our estimator is based on a moment function that satisfies a property that is stronger than Neyman orthogonality.
Finally, our work is linked to the literature on inference in randomized experiments with covariates Freedman2008regressionadjustmentsseveraltreatments, freedman2008regression, lin2013agnostic, wager2016high, lei2021regression, chiang2023regression, chang2024exact.
The remainder of this paper is organized as follows. In Section (ref), we introduce the setup and review existing procedures. In Section (ref), we describe our proposed covariate-adjusted RD estimator, and we present our main theoretical results in Section (ref). Further extensions are discussed in Section (ref). We present our empirical results in Section (ref) and simulation results in Section (ref). Section (ref) concludes. The proofs of our main results are given in Appendix (ref). Appendix (ref) formally studies the proposed inference procedures and Appendix (ref) gives details on our literature survey. The Online Supplement contains additional empirical and simulation results.
We begin by considering sharp RD designs. The data $\{W_i\}_{i\in [n]}= \{(Y_i,X_i,Z_i)\}_{i\in [n]}$, where $[n]= \{1, \dots, n\}$, are an i.i.d.\ sample of size $n$ from the distribution of $W=(Y,X,Z)$. Here, $Y_i\in\mathbb{R}$ is the outcome variable, $X_i\in\mathbb{R}$ is the running variable, and $Z_i \in \mathbb{R}^d$ is a (possibly high-dimensional) vector of covariates.\footnote{Throughout the paper, we assume that the distribution of the running variable $X_i$ is fixed, but we allow the conditional distribution of $(Y_i,Z_i)$ given $X_i$ to change with the sample size in our asymptotic analysis. In particular, we allow the dimension of $Z_i$ to grow with $n$ in order to accommodate high-dimensional settings, but we generally leave such dependence on $n$ implicit in our notation. } Units receive the treatment if and only if the running variable exceeds a known threshold, which we normalize to zero without loss of generality. We denote the treatment indicator by $T_i$, so that $T_i=\mathbf{1}\{X_i \geq 0\}$. The parameter of interest is the height of the jump in the conditional expectation function of the observed outcome variable given the running variable at zero:
where we use the notation that $f(0^+) =\lim_{x\downarrow 0}f(x)$ and $f(0^-) =\lim_{x\uparrow 0}f(x)$ are the right and left limit, respectively, of a generic function $f(x)$ at zero. In a potential outcomes framework, the parameter $\tau$ coincides with the average treatment effect of units at the cutoff under certain continuity conditions hahn2001identification.
Without the use of covariates, the parameter $\tau$ is typically estimated by running separate local linear regressions fan1996local on each side of the cutoff. That is, the baseline no covariates RD estimator takes the form
where $S_i =(T_i, X_i, T_i X_i,1)^\top$, $K_h(v)=K(v/h)/h$ with $K(\cdot)$ a kernel function and $h>0$ a bandwidth, and $e_1 = (1,0,0,0)^\top$ is the first unit vector. This estimator is a linear smoother that can also be written as a weighted sum of the realizations of the outcome variable, $$ \widehat\tau_{\scriptscriptstyle base}(h) = \sum_{i=1}^n w_i(h) Y_i,$$ where the $w_i(h)$ are local linear regression weights that depend on the data through the realizations of the running variable only; see Appendix (ref) for an explicit expression.
Under standard conditions hahn2001identification, which include a continuously distributed running variable and that the bandwidth $h$ tends to zero at an appropriate rate, the estimator $\widehat\tau_{\scriptscriptstyle base}(h)$ is approximately normally distributed in large samples under conventional pointwise asymptotics, with bias of order $h^2$ and variance of order $(nh)^{-1}$:
Here “$\stackrel{a}{\sim}$” indicates a finite-sample distributional approximation justified by an asymptotic normality result, and the bias and variance terms are given, respectively, by
The terms $\bar{\nu}$ and $\bar{\kappa}$ are kernel constants, defined as $\bar{\nu}= (\bar{\nu}_2^2 - \bar{\nu}_1 \bar{\nu}_{3})/(\bar{\nu}_2\bar{\nu}_0-\bar{\nu}_1^2)$ for $\bar{\nu}_{j}= \int_{0}^{\infty}v^j K(v)dv$ and $\bar{\kappa}= \int_0^\infty(K(v)(\bar{\nu}_1 v - \bar{\nu}_2))^2dv/ (\bar{\nu}_2\bar{\nu}_0-\bar{\nu}_1^2)^2$, and $f_X$ denotes the density of $X_i$. Practical methods for inference based on approximations like (ref) are discussed, for instance, by calonico2014robust and armstrong2020simple.
If covariates are available, they can be used to improve the accuracy of empirical RD estimates.\footnote{Throughout the paper, we focus on settings in which covariates are included to improve estimation efficiency and not to restore identification of the RD parameter by making the design more plausible. } The arguably most popular strategy calonico2019regression is to include them linearly and without kernel localization in the local linear regression (ref):
By simple least squares algebra, this “linear adjustment” estimator can be written as a no covariates estimator with covariate-adjusted outcome $Y_i-Z_i^{\top}\widehat\gamma_h$, where $\widehat\gamma_h$ is the minimizer with respect to $\gamma$ in (ref):
The linear adjustment estimator is consistent for the RD parameter without functional form assumptions on the underlying conditional expectations if the covariates are predetermined, in the sense that they are not causally affected by the treatment, and thus their conditional expectation given the running variable varies smoothly around the cutoff. Moreover, if $\mathbb{E}[Z_i|X_i=x]$ is twice continuously differentiable around the cutoff, then
under pointwise asymptotics and regularity conditions analogous to those for the no covariates estimator. Here the bias term $B_{\scriptscriptstyle base}$ is the same as that of the no covariates estimator and the new variance term is
where $\gamma_0$, a non-random vector of projection coefficients, is the probability limit of $\widehat\gamma_h$. The first-order asymptotic properties of $\widehat\tau_{\scriptscriptstyle lin}(h)$ are thus the same as that of its infeasible counterpart $\widetilde\tau_{\scriptscriptstyle lin}(h)=\sum_{i=1}^n w_i(h)(Y_i-Z_i^{\top}\gamma_0)$ that uses the population projection coefficients $\gamma_0$ instead of their estimates $\widehat\gamma_h$ to create the adjusted outcome variable. As $V_{\scriptscriptstyle lin}\leq V_{\scriptscriptstyle base}$ under standard conditions kreiss2021regression, including a fixed number of covariates generally increases the precision of the estimator in large samples. To construct standard errors and confidence intervals, one can then use methods developed for the no covariates case, replacing the original outcome $Y_i$ with the adjusted outcome $Y_i-Z_i^{\top}\widehat\gamma_h$ in the respective formulas calonico2019regression,armstrong2018optimal. For instance, a nearest-neighbor standard error $\widehat{\textnormal{se}}_{\scriptscriptstyle lin}(h)$ of $\widehat\tau_{\scriptscriptstyle lin}(h)$ is
Here $\widehat\sigma_{i, \scriptscriptstyle lin}^2$ is an estimate of $\sigma_{i, \scriptscriptstyle lin}^2=\mathbb{V}(Y_i-Z_i^{\top}\gamma_0|X_i)$, $R\geq 1$ is a (small) fixed integer, and $\mathcal{R}_i$ is the set that contains the indices of the $R$ nearest neighbors of unit $i$ in terms of their realization of the running variable among units on the same side of the cutoff.
While linear adjustments are easy to implement, they might not exploit the available covariate information efficiently. Inference might also not be reliable with linear adjustments if the number of covariates is large relative to the effective sample size.\footnote{If there are many covariates relative to the number of observations that receive positive kernel weights in (ref), the standard error in (ref) is generally downward biased. This bias occurs because the local empirical variances $\widehat\sigma^2_{i,lin}$ are typically smaller than their population counterparts $\sigma^2_{i,lin}$ in such cases due to overfitting. If the number of covariates exceeds the number of observations with positive kernel weights, the estimator in equation (ref) is of course not even well-defined in the first place.} In this paper, we propose a new method to address these issues. It allows for general nonlinear covariate adjustments and can accommodate regularization methods in the estimation of the adjustment terms.
To motivate our flexible covariate adjustments, recall that the linear adjustment estimator is asymptotically equivalent to a no covariates RD estimator of the form in (ref) that uses the covariate-adjusted outcome $Y_i-Z_i^{\top}\gamma_0$ instead of the original outcome $Y_i$. We generalize this by considering a class of estimators with covariate-adjusted outcomes based on potentially nonlinear adjustment functions $\eta$:
We note that the adjustment function must be the same for observations on either side of the cutoff, as using different adjustments on either side of the cutoff would generally yield inconsistent RD estimates calonico2019regression.
If the covariates are predetermined, their conditional distribution given the running variable should vary smoothly around the cutoff. We formalize this notion by assuming that for every adjustment function $\eta$, the function $\mathbb{E}[\eta(Z_i)|X_i=x]$ is twice continuously differentiable around the cutoff.\footnote{Our analysis only rules out adjustment functions that do not satisfy certain technical regularity conditions, such as functions for which the respective conditional expectation does not exist in the first place. Assuming smoothness of $\mathbb{E}[\eta(Z_i)|X_i=x]$ for (essentially) all $\eta$ is stronger than only assuming smoothness of $\mathbb{E}[Z_i|X_i=x]$, as in calonico2019regression. Our stronger assumption, however, is still very much in line with the notion of covariates being predetermined.} This assumption implies that
The estimator $\widehat\tau(h;\eta)$ can thus be seen as a sample analog estimator based on the moment condition (ref), which identifies $\tau$ and is globally invariant with respect to the adjustment function $\eta$. Because of this global invariance, we expect that
for every $\eta$ under standard regularity conditions. Under these pointwise asymptotics, the bias term in (ref) is again that of the baseline no covariates estimator as it does not depend on the adjustment function due to the assumed smoothness of $\mathbb{E}[\eta(Z_i)|X_i=x]$. On the other hand, the variance term in (ref) does depend on $\eta$, and is given by
To maximize the precision of the estimator $\widehat\tau(h;\eta)$ for any particular bandwidth $h$, we want to choose $\eta$ such that $V(\eta)$ is as small as possible. Our analysis below shows that using the equally-weighted average of the left and right limits of the “long” conditional expectation function $\mathbb{E}[Y_i|X_i=x,Z_i=z]$ at the cutoff achieves this goal. That is, we show that $$V(\eta)\geq V(\eta_0) \textnormal{ for all } \eta,$$ where
As the optimal adjustment function $\eta_0$ is generally unknown in practice, we propose to estimate the RD parameter $\tau$ by a feasible version of $\widehat \tau(h;\eta_0)$.
Our proposed estimator requires a first-stage estimate of the optimal adjustment function, which does not have to be of a particular type: practitioners can use classical nonparametric or modern machine learning methods to reduce the risk of model misspecification, or choose suitable parametric methods (conventional linear adjustments can be seen as a special case of the latter type). Our proposed estimator also uses cross-fitting, which is an efficient form of sample splitting that prevents overfitting of the estimated adjustment function and avoids unrealistic empirical process conditions in the theoretical analysis chernozhukov2018double. Specifically, our estimator is computed in two steps:
Our theoretical analysis below shows that under weak conditions the estimator $\widehat{\tau}(h; \widehat\eta)$ is asymptotically equivalent to the infeasible estimator $\widehat{\tau}(h; \bar\eta) = \sum_{i=1}^n w_i(h) M_i(\bar\eta)$, where $\bar\eta$ is a deterministic approximation of $\widehat\eta$ whose error vanishes in large samples in some appropriate sense. Importantly, our approach does not require the first-stage estimator of $\eta_0$ to be consistent, in the sense that we allow for the possibility that $\bar\eta\neq\eta_0$. The first-stage estimator also does not have to converge with a particularly fast rate. In view of (ref), it then holds that
As mentioned above, the variance $V(\bar\eta)$ is minimized if $\bar\eta =\eta_0$. However, the distributional approximation is also valid if $\bar\eta \neq\eta_0$ because the moment condition (ref) holds for all adjustment functions, and not just the optimal one. In that sense, our procedure is robust to misspecification or over-regularized estimation of the optimal adjustment function. Moreover, we argue that $ V(\bar\eta)$ is typically smaller than $V_{\scriptscriptstyle base}$ or $V_{\scriptscriptstyle lin}$ even if $\bar\eta \neq\eta_0$.
We also show that other common steps in an empirical RD analysis can easily be implemented by applying existing methods that are devised for settings without covariates to the generated data set $\{(X_i,M_i(\widehat\eta_{s(i)}))\}_{i\in[n]}$. For example, we can construct an estimator of the bandwidth that minimizes the asymptotic mean squared error of $\widehat{\tau}(h; \widehat\eta)$ by using the procedures proposed by calonico2014robust or imbens2012optimal. Similarly, we can generalize the standard error (ref) and construct a valid nearest-neighbor standard error $\widehat{\textnormal{se}}(h; \widehat\eta)$ as
and construct “robust bias correction” and “bias-aware” confidence intervals as in calonico2014robust and armstrong2020simple, respectively. To reduce the sensitivity of empirical findings to the particular data split in the cross-fitting step, we can proceed as in chernozhukov2018double by repeating the respective procedure several times and reporting a summary measure of the results, such as the median. We recommend proceeding like this especially when working with smaller sample sizes.
We now discuss some implementation details for the first-stage estimator of the optimal adjustment function $\eta_0$. Our theoretical analysis allows for a variety of different methods to be used in this context. If one wishes to maintain the simplicity of the conventional linear adjustment, one can obtain a “cross-fitted” version of $\widehat\tau_{\scriptscriptstyle lin}(h)$ by setting $\widehat\eta_s(z) = z^\top\widehat\gamma_{s,h}$, for $s \in \{1,\ldots,S\}$, where $\widehat\gamma_{s,h}$ is the minimizer w.r.t.\ $\gamma$ in the minimization problem
We refer to this procedure as the cross-fitted “localized” linear adjustment. The adjustment coefficients, however, do not need to be estimated using the kernel weights from the second-stage regression. As an important variant of cross-fitted linear adjustments, we consider $\widehat\eta_s(z) = z^\top\widehat\gamma_{s,\infty}$, where $\widehat\gamma_{s,\infty}$ is obtained via a version of (ref) without kernel weights. Since all the observations outside of fold $s$ are used to obtain $\widehat\gamma_{s,\infty}$, we refer to this procedure as the cross-fitted “global” linear adjustment. In finite samples, the global version can outperform the localized one in terms of the resulting standard deviation of the RD estimator due to the increased stability of the first-stage estimates. This approach can be naturally extended to other parametric models where the components involving $S_i$ and $Z_i$ are additively separable, and it can be combined with lasso regularization. The cross-fitted post-lasso adjustments are then obtained via (ref) with the set of covariates restricted to those “selected” by the lasso.
More generally, our approach allows for any parametric, classical nonparametric as well as generic modern machine learning methods. To allow for such generality, we consider adjustment functions of the form $$ \widehat\eta_s (z) = \frac{1}{2}(\widehat\mu_s^+(z)+\widehat\mu_s^-(z)),\quad s \in \{1,\ldots,S\}, $$ where $\widehat\mu_s^+(z)$ and $\widehat\mu_s^-(z)$ are separate estimators of $\mu_0^+(z)=\mathbb{E}[Y_i|X_i=0^+,Z_i=z]$ and $\mu_0^-(z)=\mathbb{E}[Y_i|X_i=0^-,Z_i=z]$, respectively, using the data outside of fold $s$. With appropriate subject knowledge, one can then, for example, specify parametric models for $\mathbb{E}[Y_i|X_i=x,Z_i=z]$. If the number of covariates is small, the functions $\mu_0^+$ and $\mu_0^-$ can be also estimated using classical nonparametric methods under smoothness conditions, with local polynomial regression being particularly suitable due to its good boundary properties. If the number of covariates is large, however, we recommend using modern machine learning methods, such as lasso or post-lasso regression, random forests, deep neural networks, boosting, or ensemble combinations thereof.
One issue to consider is that the default implementations of generic machine learning estimators of $\mathbb{E}[Y_i|X_i=x,Z_i=z]$ will not automatically produce an estimate with a jump at the cutoff. As having this feature is potentially important in our context, we consider two simple variations of generic machine learning estimators. First, let
be a generic machine learning estimator of $\mathbb{E}[Y_i|T_i=t, X_i=x, Z_i=z]$, computed by minimizing some empirical loss function $L(f) =\sum_{i \in I_s^c} l(Y_i,f(T_i, X_i, Z_i))$ over a set of candidate functions $\mathcal{F}$. We can then estimate $\mu^+(z)$ by $\widehat \mathbb{E}_s[Y_i|T_i=1, X_i=0, Z_i=z]$ and $\mu^-(z)$ by $\widehat \mathbb{E}_s[Y_i|T_i=0, X_i=0, Z_i=z]$. Here including the seemingly superfluous treatment indicator $T_i=\mathbf{1}\{X_i\geq 0\}$ as a predictor allows the machine learner to create the “jump” in the estimated function at the cutoff value. We refer to this type of implementation as “global”, as it uses all available observations.
To define the second type of implementation of machine learning we consider in this paper, let
be a generic machine learning estimator of $\mathbb{E}[Y_i|T_i=t, Z_i=z]$, where $b>0$ is some positive bandwidth and $K$ is again a kernel function. We can then estimate $\mu^+(z)$ as $\widehat \mathbb{E}_s[Y_i|T_i=1, Z_i=z]$ and $\mu^-(z)$ as $\widehat \mathbb{E}_s[Y_i|T_i=0, Z_i=z]$. We refer to this type of implementation as “localized”, as it effectively only uses data points whose realization of the running variable is close to the cutoff. The idea is to produce an estimate with small empirical loss in the relevant area around the cutoff rather than one with small “overall” loss. The downside of this approach is the reduced effective sample size and the need to choose the tuning parameter $b$.\footnote{The choice of $b$ involves a bias-variance trade-off similar to the one encountered in classical nonparametric kernel regression problems. We are not aware of generic theoretical results for such estimators in settings with $b\to 0$ as $n\to\infty$. Specific results are given by su2019non for the lasso, and by colangelo2022double for series estimators and deep neural networks. In our applications below, we simply use $b=h$. To make this simultaneous choice feasible, we use an iterative procedure. We first choose a reasonable preliminary first-stage bandwidth, like the one that would be optimal for RD estimation without covariates, and generate preliminary versions of the adjustment terms as described above. Next, we use the preliminary covariate-adjusted outcomes to pick an optimized second-stage bandwidth. Finally, we rerun both stages with this last bandwidth to obtain our empirical results.}
The specific flexible covariate adjustment that we propose and implement in our empirical analysis and simulations is an ensemble of the following methods:
All four methods are implemented in localized and global versions discussed above. We use the cross-fitted linear and post-lasso adjustments specified in the discussion following (ref), and the boosted trees and random forest adjustments are based on the formulations in (ref) and (ref). Our proposed flexible covariate adjustment is a convex combination of these eight adjustment functions and the trivial no-adjustment function that minimizes the mean squared error for predicting the outcome close to the cutoff. Specifically, we employ the super learning approach of laan2007super, where the optimal weights are chosen via cross-validation.\footnote{We implemented our procedure in the {\tt R} programming language. The boosted trees adjustments are obtained using the package {\tt xgboost} with trees of depth 2 and shrinkage rate 0.1. The number of boosting iterations is chosen via cross-validation with a maximum of 1000 iterations, separately for the localized and global versions. The random forest with 1000 trees is implemented using the package {\tt ranger} with the minimal node size set to the maximum of 10 and 0.1% of the sample size. All other parameters are set to the default values in the respective packages. For post-lasso estimation, we use the function {\tt rlasso} from the package {\tt hdm} with a data-driven penalty parameter. We use the package {\tt SuperLearner} to choose the optimal weights via cross-validation.}
We study the theoretical properties of our proposed estimator under a number of conditions that are either standard in the RD literature, or concern the general properties of the first-stage estimator $\widehat\eta$. To describe them, we denote the support of $Z_i$ by $\mathcal{Z}$, and the support of $X_i$ by $\mathcal{X}$. We write $\mathcal{X}_h=\mathcal{X} \cap [-h,h]$, and $\mathcal{Z}_h$ denotes the support of $Z_i$ given $X_i\in \mathcal{X}_h$. We also define the following class of admissible adjustment functions:
The class $\mathcal E$ implicitly depends on the underlying conditional distribution of the covariates given the running variable. If this conditional distribution changes smoothly around the cutoff, the class $\mathcal E$ contains essentially all functions of the covariates, subject only to technical integrability conditions.\footnote{For example, if the conditional distribution of $Z_i$ given $X_i$ admits a density $f_{Z|X}(z|x)$ that is twice continuously differentiable in $x$ and $ |\partial_x^jf_{Z|X}(z|x)|\leq g_j(z)$ for all $x$ in a neighborhood of the cutoff, some integrable functions $g_j$, and $j \in \{0,1,2\}$, then $\mathcal E$ contains all bounded Borel functions. The class $\mathcal E$ also contains all polynomials if the corresponding conditional moments of $Z_i$ exist and are twice continuously differentiable.}
Assumption (ref) states that with high probability the first-stage estimator belongs to some realization set $\mathcal{T}_n \subset \mathcal E$. As discussed above, this requirement seems weak as we generally expect the class $\mathcal E$ to be very large. The assumption also states that the sets $\mathcal{T}_n$ contract around a deterministic sequence of functions in a particular $L_2$-type sense. Note that taking the supremum in Assumption (ref) over $\mathcal X_h$ instead of $\mathcal X$ suffices as the properties of the first-stage estimator are only relevant for observations with non-zero kernel weights in the second-stage local linear regression. The assumption does not impose any restrictions on the speed at which $\widehat\eta$ concentrates around $\bar\eta$. It also allows the function $\bar\eta$ to be different from the target function $\eta_0$, so that $\widehat\eta$ could be inconsistent for $\eta_0$.
Mean-square error consistency as in Assumption (ref) follows under classical conditions for the parametric and nonparametric procedures for settings in which the number of covariates is fixed. For the “localized” versions of the machine learning methods described in Section (ref), existing results imply that for fixed $b>0$ and $K$ the uniform kernel,
with $\bar\eta(z) = ( \mathbb{E}[Y_i|X_i\in(- b,0), Z_i=z] + \mathbb{E}[Y_i|X_i\in(0,b), Z_i=z])/2$ and some $ r_{n}=o(1)$, under general conditions. For example, if $\bar\eta(z)$ is contained in a Hölder class of order $s$, then (ref) can hold with $r_n^2 = n^{-2s/(2s+d)}$ for estimators that exploit smoothness. If $\bar\eta(z)$ is $s$-sparse, then (ref) can hold with $r_n^2 = s\log(d)/n$ for estimators that exploit sparsity. Assumption (ref) then follows from (ref) if the conditional distribution of the covariates does not change “too quickly” when moving away from the cutoff. For example, if the covariates are continuously distributed conditional on the running variable, having that
for some constant $C$ and all $n$ sufficiently large, suffices. Similar conditions can be given for discrete conditional covariate distributions, or intermediate cases. If $\mathbb{E}[Y_i|X_i=x, Z_i=z]$ is sufficiently smooth in $x$ on both sides of the cutoff, we can also expect that $\bar\eta$ is “close” to $\eta_0$ for “small” values of $b$. Formal rate results with $b\to 0$ are given by su2019non for the Lasso, and by colangelo2022double for series and deep neural networks.
Assumption (ref) also concerns the first-stage estimator, and requires the first and second derivatives of $\mathbb{E}\left[ \eta(Z_i)-\bar\eta(Z_i) |X_i=x\right]$ to be close to zero in large samples for all $\eta\in\mathcal{T}_n$. We generally expect this condition to hold with $v_{1,n}=v_{2,n}=r_n$ with $r_n$ as in Assumption (ref).\footnote{For example, this can easily be seen to be the case if $\widehat\eta$ converges to $\bar\eta$ uniformly over $\mathcal Z$ with rate $r_n$ and the smoothness conditions for $f_{Z|X}(z|x)$ given in footnote (ref) hold. Similarly, under regularity conditions on $\mathbb{E}[Z_i|X_i=x]$, these three rates coincide if $\mathcal{T}_n$ contains only linear functions. Without any additional restrictions on first-stage estimators or $\mathcal T_n$, except that it contains only bounded functions, Assumption (ref) also follows from Assumption (ref), again with $v_{1,n}=v_{2,n}=r_n$, under restrictions concerning solely the conditional density $f_{Z|X}(z|x)$. Specifically, it suffices that $\mathbb{E}\big[\big(\partial_x^jf_{Z|X}(Z_i|x)/f_{Z|X}(Z_i|x)\big)^2|X_i=x\big]$ is bounded for $j \in \{1,2\}$ uniformly in $x$ and the conditions from footnote (ref) hold.}
Assumption (ref) is a standard condition from the RD literature. Continuity of the running variable's density $f_X$ around the cutoff is strictly speaking not required for an RD analysis. However, a discontinuity in $f_X$ is typically considered to be an indication of a design failure that prevents $\tau$ from being interpreted as a causal parameter mccrary2008manipulation, gerard2020bounds. For this reason, we focus on the case of a continuous running variable density in this paper.
The conditions on the kernel and the bandwidth that are imposed in Assumption (ref) are standard in the RD literature.
Assumption (ref) collects standard conditions for an RD analysis with $M_i(\bar\eta)$ as the outcome variable. Part (i) imposes smoothness conditions on $\mathbb{E}[M_i(\bar\eta)|X_i=x]$, and parts (ref) and (ref) impose restrictions on conditional moments of the outcome variable. Throughout, we use constants $C$ and $L$ independent of the sample size to ensure asymptotic normality of the infeasible estimator $\widehat\tau(h;\bar\eta)$ even in settings where the distribution of the data, and thus $\bar\eta$, might change with $n$.
We give four main results in this subsection. The first shows that our proposed estimator $\widehat{\tau}(h; \widehat\eta)$ is asymptotically equivalent to an infeasible analog $\widehat{\tau}(h;\bar\eta)$ that replaces the estimator $\widehat\eta$ with the deterministic sequence $\bar\eta$; the second shows the asymptotic normality of the estimator; the third characterizes how the asymptotic variance changes with the adjustment function and shows that $\eta_0$ is indeed the optimal adjustment; and the fourth shows the impact of flexible covariate adjustments on the optimal bandwidth and the corresponding mean squared error.
Theorem (ref) is easiest to interpret in what is arguably the standard case that $v_{1,n}=v_{2,n}=r_n$, in which it holds that $$\widehat{\tau}(h; \widehat\eta) = \widehat{\tau}(h;\bar\eta) +O_P(r_n(h^2+(nh)^{-1/2})) = \widehat{\tau}(h;\bar\eta) + O_P(r_n|\widehat{\tau}(h;\bar\eta)-\tau|).$$ The accuracy of the approximation that $\widehat{\tau}(h; \widehat\eta) \approx \widehat{\tau}(h;\bar\eta)$ thus increases with the rate at which $\widehat\eta$ concentrates around $\bar\eta$, but first-order asymptotic equivalence holds even if the first-stage estimator converges arbitrarily slowly. This insensitivity of $\widehat{\tau}(h; \widehat\eta)$ to sampling variation in $\widehat\eta$ occurs because $\widehat{\tau}(h; \widehat\eta)$ is based on the moment condition $$\tau = \mathbb{E}[M_i(\eta)|X_i=0^+] - \mathbb{E}[M_i(\eta)|X_i=0^-], $$ which is insensitive to variation in $\eta$. Moment conditions with a local form of insensitivity with respect to a nuisance function, often called Neyman orthogonality, are used extensively in the recent literature on two-stage estimators that use machine learning in the first stage belloni2017program, chernozhukov2018double. The global insensitivity that arises in our RD setup is stronger, and allows us to work with weaker conditions on the first-stage estimates than those used in papers that work with Neyman orthogonality. Similarly, globally insensitive moment function exists, for example, in certain types of randomized experiments, and our proposed estimator is in many ways analogous to efficient estimators in such setups; see Section (ref) for further discussion.
Theorem (ref) follows from Theorem (ref) under the additional regularity conditions of Assumption (ref). It shows that our estimator is asymptotically normal, gives explicit expressions for its asymptotic bias and variance, and justifies the distributional approximation given in Section (ref).
Theorem (ref) introduces a function class $\mathcal V$ that, similarly to the class $\mathcal E$ above, enforces some technical integrability conditions. The theorem shows that $V(\eta^{(a)}) < V(\eta^{(b)})$ for generic adjustment functions $\eta^{(a)}$ and $\eta^{(b)}$ if and only if $\mathbb{V}[\eta_0(Z_i)-\eta^{(a)}(Z_i)|X_i=0] < \mathbb{V}[\eta_0(Z_i)-\eta^{(b)}(Z_i)|X_i=0]$. That is, the “closer” (in a particular $L_2$-sense) the adjustment function is to the optimal one, the smaller the asymptotic variance. In consequence, the lowest possible value of $V(\bar\eta)$ is achieved for $\bar\eta=\eta_0$. Even if $\bar\eta \neq \eta_0$, our flexible covariate adjustment RD estimators typically still have smaller asymptotic variances than the no covariates and linear adjustment RD estimators. Specifically, $V(\bar\eta) < V_{base}$ if and only if $\mathbb{V}[\eta_0(Z_i)-\bar\eta(Z_i)|X_i=0] < \mathbb{V}[\eta_0(Z_i)|X_i=0]$, i.e.\ whenever $\bar\eta(Z_i)$ captures some of the variance of $\eta_0(Z_i)$ among units near the cutoff; and $V(\bar\eta) < V_{lin}$ if and only if $\mathbb{V}[\eta_0(Z_i)-\bar\eta(Z_i)|X_i=0] < \mathbb{V}[\eta_0(Z_i)-Z_i^\top\gamma_0|X_i=0]$, i.e.\ whenever $\bar\eta$ is “closer” to $\eta_0$ in our particular $L_2$-sense than the population linear adjustment function.
Theorem (ref) implies that flexible covariate adjustments can reduce the approximate mean squared error of our estimator not only directly through a smaller asymptotic variance term but also indirectly through a change in the optimal bandwidth and a corresponding reduction in bias. That is, if $V(\eta^{(a)}) < V(\eta^{(b)})$ for generic adjustment functions $\eta^{(a)}$ and $\eta^{(b)}$, then the optimal bandwidth $h_{\scriptscriptstyle{AMSE}}(\eta^{(a)})$ is smaller than $h_{\scriptscriptstyle{AMSE}}(\eta^{(b)})$, and the corresponding estimator $\widehat{\tau}(h_{\scriptscriptstyle{AMSE}}(\eta^{(a)}); \eta^{(a)})$ has both smaller asymptotic bias and smaller asymptotic variance than $\widehat{\tau}(h_{\scriptscriptstyle{AMSE}}(\eta^{(b)}); \eta^{(b)})$.
We formally show in Appendix (ref) that standard methods for bandwidth choice and confidence interval construction based on the no covariates RD estimator maintain their general asymptotic properties when they are applied to the generated data set $\{(X_i,M_i(\widehat\eta_{s(i)}))\}_{i\in[n]}$ without any adjustment for the sampling uncertainty about the estimated adjustment function. Specifically, we derive three groups of results, all under conditions that are rather weak and analogous to those commonly imposed in setups without covariates.
First, we show that the nearest neighbor standard error (ref) is consistent, in the sense that $$nh\,\widehat{\text{se}}^2( h; \widehat \eta) /V(\bar\eta)\overset{ p}{\to}1.$$ Second, we show that commonly used methods for confidence interval construction achieve correct asymptotic coverage. For example, assuming a bound on $|\partial^2_x \mathbb{E}[M_i(\bar\eta)|X_i=x]|$, the absolute value of the second derivative of the conditional expectation of the adjusted outcome given the running variable, one can construct a “bias-aware” confidence interval as in armstrong2020simple as $$ CI^{ba}_{1-\alpha} = \left[\widehat{\tau}(h; \widehat \eta) \pm z_\alpha(\bar b(h)/\widehat{\text{se}}( h; \widehat \eta)) \,\widehat{\text{se}}( h; \widehat \eta)\right]. $$ Here $z_\alpha(r)$ is the $(1-\alpha)$-quantile of $|N(r,1)|$, the folded normal distribution with mean $r$ and variance one, and $\bar b(h)$ is an explicit bound on the finite-sample bias of $\widehat\tau(h,\bar\eta)$ given in the appendix. Alternatively, one can also construct a “robust bias correction” confidence interval as in calonico2014robust by subtracting a local quadratic estimate of the first-order bias of $\widehat{\tau}(h; \widehat \eta)$ from the estimator, and adjusting the standard error appropriately. This yields a confidence interval of the form \[ CI^{rbc}_{1-\alpha} = \big[\widehat{\tau}^{rbc}(h; \widehat \eta) \pm z_\alpha \widehat{\text{se}}^{rbc}( h; \widehat \eta)\big], \] where $z_\alpha = z_\alpha(0)$ and the other terms are formally defined in the appendix. Third, we show that the MSE-optimal bandwidth selector $\widehat h_{n}$ of calonico2014robust, which is similar to that of imbens2012optimal, consistently estimates the AMSE-optimal bandwidth $h_{\scriptscriptstyle{AMSE}}(\bar \eta)$ defined in Theorem (ref), in the sense that $$ \widehat h_{n} / h_{\scriptscriptstyle{AMSE}}(\bar \eta) \overset{p}{\to} 1.$$ RD estimation and inference with flexible covariate adjustments are thus easy to implement with existing software packages.
The results in Section (ref) are qualitatively similar to ones obtained for efficient influence function (EIF) estimators of the population average treatment effect (PATE) in randomized experiments with known and constant propensity scores wager2016high,chernozhukov2018double. To see this, consider a randomized experiment with unconfounded treatment assignment and known constant propensity score $p$. Using our notation in an analogous fashion, the EIF of the PATE in such a setup is typically written in the form $$\psi_i(m_0^0, m_0^1) = m_0^1(Z_i) - m_0^0(Z_i) + \frac{T_i(Y_i-m_0^1(Z_i))}{p} - \frac{(1-T_i)(Y_i-m_0^0(Z_i))}{1-p}, $$ where $m_0^t(z) = \mathbb{E}[Y_i|Z_i=z,T_i=t]$ for $t \in \{0,1\}$ hahn1998role. The minimum variance any regular estimator of the PATE can achieve is thus $V_{\scriptscriptstyle \textnormal{PATE}}=\mathbb{V}(\psi_i(m_0^0, m_0^1))$. By randomization, it also holds that $\tau_{\scriptscriptstyle \textnormal{PATE}} = \mathbb{E}[\psi_i(m^0, m^1 )]$ for all (suitably integrable) functions $m^0 $ and $m^1$. The PATE is thus identified by a moment function that satisfies a global invariance property. A sample analog estimator of $\tau_{\scriptscriptstyle \textnormal{PATE}}$ based on this moment function has asymptotic variance $V_{\scriptscriptstyle \textnormal{PATE}}$ if $\widehat m^t$ is a consistent estimator of $m_0^t$ for $t\in\{0,1\}$, but remains consistent and asymptotically normal with asymptotic variance $\mathbb{V}(\psi_i(\bar m^0, \bar m^1))$ if $\widehat m^t$ is consistent for some other function $\bar{m}^t$, $t\in\{0,1\}$. The convergence of $\widehat m^t$ to $\bar{m}^t$ can be arbitrarily slow for these results wager2016high,chernozhukov2018double.
The qualitative parallels between these findings and ours in Section (ref) arise because our covariate-adjusted RD estimator is in many ways a direct analog of such EIF estimators. To show this, write $ m(z)= (1- p){m}^1(z) + p\,{m}^0(z)$ for any two functions $m^0$ and $m^1$, so that $ m_0(z)= (1- p){m}_0^1(z) + p\,{m}_0^0(z)$. The PATE's influence function can then be expressed as $$\psi_i({m}_0^0, {m}_0^1) = \frac{T_i(Y_i- m_0(Z_i))}{p} - \frac{(1-T_i)(Y_i-m_0(Z_i))}{1-p}, $$ and it holds that $$\mathbb{E}[\psi_i(m^0, m^1 )]=\mathbb{E} [Y_i- m(Z_i)|T_i=1] - \mathbb{E}[Y_i- m(Z_i)|T_i=0],$$ which is the difference in average covariate-adjusted outcomes between treated and untreated units. This last equation is fully analogous to our equation (ref), with $p=1/2$, and conditioning on $T_i=1$ and $T_i=0$ replaced by conditioning on $X_i$ in infinitesimal right and left neighborhoods of the cutoff (the value $p=1/2$ is appropriate here because continuity of the running variable's density implies that an equal share of units close to the cutoff can be found on either side). An EIF estimator of $\tau_{\scriptscriptstyle \textnormal{PATE}}$ is thus analogous to our estimator $\widehat{\tau}(h; \widehat\eta)$, as they are both sample analogs of a moment function with the same basic properties.\footnote{We note that neither in our setting nor with EIF estimation of the PATE would replacing the known propensity score with some empirical estimate result in any efficiency gains. The finding from hahn1998role that using an estimated propensity score can be more efficient than using the true one refers to inverse probability weighting (IPW) type estimators and does not apply to the EIF type estimators we consider here.}
In fuzzy RD designs, units are assigned to treatment if their realization of the running variable falls above the threshold value, but might not comply with their assignment. The conditional treatment probability given the running variable hence changes discontinuously at the cutoff, but in contrast to sharp RD designs it does not jump from zero to one. The parameter of interest in fuzzy RD designs is $$ \theta = \frac{\tau_Y}{\tau_T} \equiv \frac{\mathbb{E} [Y_i|X_i=0^+] - \mathbb{E} [Y_i|X_i=0^-]}{\mathbb{E} [T_i|X_i=0^+]- \mathbb{E} [T_i|X_i=0^-]}, $$ which is the ratio of two sharp RD estimands (throughout this subsection, the notation is analogous to that used before, with the subscripts $Y$ and $T$ referencing the respective outcome variable). Under standard conditions hahn2001identification, dong2014alternative, one can interpret $\theta$ as the average causal effect of the treatment among units at the cutoff whose treatment decision is affected by whether their value of the running variable is above or below the cutoff.
Similarly to sharp RD designs, predetermined covariates can be used in fuzzy RD designs to improve efficiency. Building on our proposed method, we consider estimating $\theta$ by the ratio of two generic flexible covariate-adjusted sharp RD estimators: \[ \widehat{\theta}(h; \widehat \eta_Y,\widehat \eta_T) = \frac{ \widehat{\tau}_{Y} (h; \widehat \eta_Y )}{\widehat{\tau}_{T} (h; \widehat \eta_T ) } =\frac{ \sum_{i=1}^n w_{i}(h) (Y_i- \widehat \eta_{Y, s(i)}(Z_i)) }{ \sum_{i=1}^n w_{i}(h) (T_i- \widehat \eta_{T, s(i)}(Z_i)) }. \]
The first part of the proposition shows that our flexible covariate-adjusted fuzzy RD estimator is asymptotically normal, with asymptotic variance that depends on the population counterparts $\bar\eta_Y$ and $\bar\eta_T$ of the two estimated adjustment functions. This result can then be used to construct a confidence interval for $\theta$ based on the t-statistic. Alternatively, confidence sets for $\theta$ can be constructed via an Anderson-Rubin-type approach, which circumvents certain problems of ratio estimators noack2021bias.
The second part of the proposition shows that the asymptotic variance of our estimator is minimized if the estimated adjustment functions concentrate around $\bar\eta_Y=\eta_{Y,0}$ and $\bar\eta_T=\eta_{T,0}$, respectively. That is, the optimal adjustment functions for fuzzy RD designs can be obtained by separately considering two covariate-adjusted sharp RD problems with outcomes $Y_i$ and $T_i$, respectively. This holds because for fixed adjustment functions $\eta_Y$ and $\eta_T$ we have that $\widehat{\theta}(h; \eta_Y, \eta_T) -\theta$ is first-order asymptotically equivalent to a sharp RD estimator with the infeasible outcome $U_i(\eta_Y,\eta_T)=\left(Y_i- \theta T_i - (\eta_Y(Z_i) - \theta \eta_T(Z_i))\right)/\tau_T.$ By our Theorem (ref), the asymptotic variance of $\widehat{\theta}(h; \eta_Y, \eta_T)$ is minimized if $(\eta_Y(Z_i) - \theta \eta_T(Z_i))/\tau_T$ equals the optimal adjustment function for the outcome $\left(Y_i- \theta T_i\right)/\tau_T$. By linearity of conditional expectations, this holds if $\eta_Y=\eta_{Y,0}$ and $\eta_T=\eta_{T,0}$.
We note that instead of the type of cross-fitting described in Section (ref), which is analogous to the DML2 method in chernozhukov2018double, one could also consider an analog of their DML1 method, which creates an overall estimate by averaging separate estimates from each data fold. In our context, this would yield an estimator of the form $$ \widehat{\tau}_{alt}( h; \widehat \eta) = \frac{1}{S}\sum_{s \in[S]} \sum_{i \in I_s} w_{i, s}(h) M_i(\widehat \eta_s),$$ where $w_{i, s}(h)$ is the local linear regression weight of unit $i$ using only data from the $s$-th fold; see Appendix (ref). Under the conditions of Theorem (ref), we see from its proof that
The estimators $\widehat{\tau}(h; \widehat \eta)$ and $\widehat{\tau}_{alt}(h; \widehat \eta)$ thus have the same first-order asymptotic distribution. However, comparing the rate in (ref) to that in Theorem (ref) shows that the alternative implementation removes a term of order $O_P(v_{1,n}h(nh)^{-1/2})$. We still prefer our proposed implementation of cross-fitting despite this improvement in second-order asymptotic properties because it allows existing routines for bandwidth selection and confidence interval construction to be applied directly to the generated data set $\{(X_i,M_i(\widehat\eta_{s(i)}))\}_{i\in[n]}$, as discussed in Section (ref).\footnote{velez2024asymptotic also argues in the context of a setting with regular parameters DML2 should be preferred DML1 due to better bias and mean-squared error properties under particular asymptotic regimes.}
To illustrate the scope for efficiency improvements that flexible covariate adjustments can achieve in practically relevant settings, we applied our method to a number of recent RD studies. Specifically, we collected data from all articles that appeared between 2018 and 2023 in the main AEA journals for applied microeconomic research, fit into our general framework, use covariates, and have directly available public replication data. We found 16 such papers with a total of 56 main specifications. For each of these specifications, we computed the length of bias-aware 95% confidence intervals for the respective RD parameter based on local linear estimators that use our flexible covariate adjustments, linear adjustments, and no covariate adjustments, respectively.\footnote{We implement the linear adjustment using cross-fitting to ensure a fair benchmark for our flexible adjustment. In Section (ref), we illustrate in simulations that the standard error based on the conventional linear adjustment may be downward-biased in settings where the number of covariates is large relative to the effective sample size.} For illustration purposes, we use the same “no covariates” smoothness bound for all estimators within each specification.\footnote{Specifically, we use the rule of thumb for the smoothness bound of imbens2017optimized, which equals twice the maximal second derivative of a second-order global polynomial fitted on each side of the cutoff. The smoothness bound is calibrated based on the original outcomes. The results in the main text are not very sensitive to specific choices of the smoothness bounds.} This analysis captures the effect of covariate adjustments on the confidence interval length that is due to a reduction in variance, abstracting from other issues. In Section (ref) of the Online Supplement we also provide results where the smoothness bound is calibrated based on the adjusted outcomes as well as results based on the robust bias correction. Appendix (ref) contains further details on the implementation and the data collection process, and Table (ref) in the Online Supplement provides the complete list of papers and specifications used.
Figure (ref) shows the distribution of the ratio of confidence interval lengths for flexible adjustments relative to no adjustments in its left panel, and for flexible adjustments relative to linear adjustments in its right panel. We first note that the confidence intervals with flexible adjustments are never noticeably wider than the ones of the no covariates or linear adjustment RD estimators. From the left panel, we can see that in nearly half of our specifications the flexible covariate adjustments yield confidence intervals that are not noticeably shorter than the ones obtained without covariate adjustments. Given the flexibility of our methods, this suggests that the covariates are not informative about the outcome in these specifications, and hence, there is no scope for efficiency gains. In many specifications, however, flexible covariate adjustments lead to substantially shorter confidence intervals, with the biggest reduction being exceeding 30%. To put this into perspective, note that one would have to increase the sample size used by the no covariates RD estimator by a factor of about 2.4 to achieve a similar reduction in the length of the confidence interval. From the right panel of Figure (ref), we can see that linear adjustments are in general unable to exhaust all the available covariate information. Indeed, the confidence intervals based on flexible adjustments can be substantially shorter, with the biggest reduction reaching 20% in our empirical exercise.
In this section, we investigate the finite-sample properties of our proposed flexible covariate adjustment RD estimators under realistic conditions in two simulation studies. The first study's purpose is to show that our theoretical results provide accurate approximations to our estimator's actual finite properties, whereas the second study's purpose is to document how the properties of our estimator and that of existing methods are affected if the number of covariates becomes large relative to the effective sample size.
Our simulations are based on real data from londono2020upstream, who study the impact of merit-based college financial aid for low-income students in a sharp RD design. Their data contain the outcome variable, a dummy for immediate enrollment in any post-secondary education, the running variable, a test score,\footnote{londono2020upstream consider two different test scores as running variables. We focus on the SABER 11 test score in this section as it is available for a larger number of data points.} and 21 covariates, namely age, family size, indicators for gender, ethnicity, employment status, parent's education, household residential stratum, high school schedule, and high school type. Our simulations involve repeatedly drawing random samples from a version of the data that is restricted to the $n=259,419$ observations with test scores below the original treatment threshold (so that none of the students remaining in the data set are actually assigned to treatment), and then estimating the effect of a placebo treatment “received” by students with test scores above the median test score value. We use either the original outcome (enrollment in any post-secondary education) or age (one of the original covariates) as the dependent variable. These two dependent variables correspond to settings in which covariate adjustments achieve almost no and quite substantial efficiency gains, respectively; see the RD estimates from the entire restricted data in Table (ref) in the Online Supplement for details. In the main text, we conduct inference using the bias-aware approach. All simulations are based on 10,000 Monte Carlo draws.
\afterpage{
}
In this simulation study, we evaluate the finite-sample performance of our methods in a typical RD setting with a moderate number of covariates and a relatively large number of observations. Specifically, we consider estimation with the original baseline covariates and samples of size 5000, which is around the median of the sample sizes of the empirical applications of our literature survey. We apply our flexible covariate adjustment discussed in Section (ref). Additionally, we consider the deterministic approximations of all the feasible adjustment methods, which were obtained by running the respective method on the full restricted dataset. By comparing the feasible adjustments and their respective deterministic approximation, we can assess the quality of the approximation in our equivalence result of Theorem (ref). For comparison, we also report the results without covariate adjustments and with conventional linear adjustments. For each adjustment method, we select the bandwidth and construct a confidence interval using the bias-aware approach with the smoothness bound calibrated using the adjusted outcomes via the rule of thumb of imbens2017optimized in each Monte Carlow draw.\footnote{The number of effective observations used in the second stage is on average around 3447 for the original outcome and around 3622 for age as the dependent variable.} The results are based on 5-fold cross-fitting with one random data split.
Table (ref) reports the main results of this simulation study. Our methods work very well for both dependent variables and all types of adjustments in that the mean simulated standard errors are close to the simulated standard deviations and the confidence intervals have simulated coverage rates close to the nominal one. The confidence intervals are slightly conservative, which is typical in bias-aware inference. The changes in the mean bias for different types of adjustments are negligible relative to the standard deviation, which is consistent with our conjecture that covariate adjustments should typically have no first-order effect on the leading bias constant.
In Panel A, the covariates have essentially no impact on the dependent variable, and so none of the methods leads to noticeable reductions in the standard deviation. In Panel B, where the covariates have some explanatory power for the dependent variable, the cross-fitted RD estimator with localized linear adjustment yields a confidence interval that is on average 12% shorter than the no covariates confidence interval. The flexible adjustment improves this performance even further. As can be expected in a setting with a small number of covariates relative to the sample size, the conventional and cross-fitted localized linear covariate adjustments yield similar results.
In Appendix (ref) of the Online Supplement, we present additional simulation results for all individual adjustment methods described in Section (ref). We further investigate the asymptotic equivalence result presented in Theorem (ref) within this simulation design and we illustrate that the estimation uncertainty of the adjustment functions is indeed asymptotically negligible relative to the overall estimation uncertainty of the RD estimators. Additionally, we present estimation and inference results based on robust bias corrections. The qualitative conclusions remain very similar to those presented above.
This simulation setting is motivated by the fact that in some empirical applications of our literature survey, researchers estimated specifications where the ratio of the effective sample size to the number of covariates is relatively small. Indeed, in Table (ref), the minimal value of this ratio is 2.7 and it falls below 20 in about one fifth of the specifications.
To mimic settings where there are many covariates relative to the effective sample size, we sample 500 observations without replacement within a distance of 25 from the placebo cutoff and use all of them in the RD regressions, i.e. we use a fixed bandwidth $h=25$. We create additional covariates by generating all second-order interaction terms of the original covariates, and we consider different settings by including the first 2, 10, 50, 100, and 150 of these covariates.\footnote{The first 21 of the technical covariates correspond to the original covariates, followed by the interaction terms. Since the covariates have essentially no explanatory power for the original outcome, the exact order of inclusion does not affect the results in this section.} In this setting, the ratio of the effective sample size to the number of covariates lies in the range between 250 and 3.33, which corresponds to the settings in our literature analysis with small values of this ratio. For each subsample, we estimate the RD parameter using the no covariates, the conventional linear adjustment, and our cross-fitted RD estimator with localized linear adjustments and localized random forest adjustments.\footnote{We chose the random forest to represent the machine learning adjustments here, but the qualitative results are similar when employing other methods. In this section, we focus on the individual adjustment methods, rather than on the flexible ensemble, to offer more direct insights into the mechanics of the linear and regularized adjustments.} In this simulation, we calibrated the smoothness constant for each estimator and number of covariates using the full sample and we kept them fixed across simulation draws. The results are based on the bias-aware inference approach and $B=11$ data splits for each Monte Carlo draw.\footnote{In Simulation II, we calibrated the smoothness bound for each method and number of covariates only once, using the rule of thumb of imbens2017optimized and oracle adjusted outcomes based on the full restricted sample described in Section (ref).}
Figure (ref) shows the bias and the standard deviation of the four estimation methods, normalized by the standard deviation of the no covariates RD estimator, for a varying number of covariates. We note that the simulated bias is insensitive to including many covariates, and we therefore focus on the standard deviation. In this setting, the covariates seem to have essentially no explanatory power for the dependent variable, and so adjustments based on them cannot lead to a reduction in the asymptotic variance of the RD estimator; see estimation results in Table (ref). As predicted by our theory, when the number of covariates remains moderate, all estimators perform very similarly, meaning that all the adjustments concentrate around the optimal function of no adjustment. However, as the number of covariates increases, the standard deviations of both the conventional and cross-fitted localized linear adjustment estimators become substantially larger than that of the no covariates RD estimator. The reason for that is that the linear regression with a large number of covariates is very variable, such that the estimated adjustments are no longer close to a constant\footnote{Such finite-sample behavior renders our asymptotic theory as well as the results of calonico2019regression inapplicable in this setting.} and the high-dimensional linear adjustments effectively add non-negligible noise to the outcome variable in this setting.
In contrast, the RD estimator with random forest adjustments, due to built-in regularization, does not become more variable as the number of covariates increases, meaning that the estimated adjustment function remains close to the optimal function of no adjustment. In general, it is therefore advisable to always rely on regularized adjustments in high-dimensional settings.
We now turn to the standard error and coverage of the confidence intervals for the respective methods. The left panel of Figure (ref) shows that the standard error of the conventional linear adjustment estimator exhibits a downward bias that increases substantially with the number of covariates, with its ratio to the estimator's standard deviation reaching less than 70% for 150 covariates. This effect is due to overfitting: with many covariates, the regression residuals that enter the standard error formula become “too close to zero”, and standard errors therefore become “too small”. With cross-fitting, this issue occurs neither for the linear adjustment nor for the random forest adjustment.
The right panel of Figure (ref) shows that, due to increasingly biased standard errors, the coverage of the conventional linear adjustment bias-aware confidence intervals with nominal level $95\%$ falls below 85% for 150 covariates. With cross-fitting, bias-aware confidence intervals have close to nominal coverage for both adjustment methods and all numbers of covariates under consideration. These simulation results demonstrate that cross-fitting yields consistent standard errors and valid inference even in high-dimensional settings where conventional methods may fail.
We have proposed a novel class of estimators that can make use of covariate information more efficiently than the conventional linear adjustment estimators that are currently used widely in practice. In particular, our approach allows the use of modern machine learning tools to adjust for covariates, and is at the same time largely unaffected by the “curse of dimensionality”. Our estimator is also easy to implement in practice, and can be combined in a straightforward manner with existing methods for bandwidth choice and the construction of confidence intervals. In our reanalysis of the literature, we show that our proposed estimator yields shorter confidence intervals in almost all empirical applications, and in some cases, these reductions can be substantial. We therefore expect our proposed estimator to be very attractive for a wide range of future economic applications.