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.
76,439 characters · 21 sections · 87 citation commands
Noise-Induced Randomization in Regression Discontinuity Designs
Regression discontinuity designs rely on known, discontinuous treatment assignment mechanisms to identify causal effects \citep*{hahn2001identification,imbens2008regression,thistlethwaite1960regression}: There is a running variable $Z_i \in \RR$ such that unit $i$ gets assigned treatment $W_i \in \cb{0, \, 1}$ whenever the running variable exceeds a cutoff $c \in \RR$, i.e., $W_i = 1(Z_i \geq c)$, and we estimate treatment effects by comparing units with $Z_i$ just above or below $c$. For example, in an educational setting where admission to a program hinges on a test score exceeding some cutoff, we could evaluate the effect of the program on marginal admits by comparing outcomes for students whose test scores fell right above and below the cutoff. Over the past decades, regression discontinuity designs have become one of the most widely used methods for causal inference, especially in the social sciences \citep*{currie2020technology}.
Explanations and qualitative justifications of identification in regression discontinuity designs often appeal to implicit, local randomization: Many factors outside the control of decision-makers determine the running variable $Z_i$ such that if some unit barely clears the eligibility cutoff for the intervention then the same unit could also plausibly have failed to clear the cutoff with a different realization of these chance factors lee2010regression. This is sometimes illustrated by reference to sampling error or other errors in measurement that cause units to have a measured running variable just above or just below the threshold. In our educational setting, there may be a group of marginal students who might barely pass or fail the test due to unpredictable variation in their test score, resulting in an effectively exogenous treatment assignment rule. Likewise, medical assays frequently involve a degree of random measurement error, whether because of sampling techniques or other sources of random variation bor2014regression.
Most formal and practical approaches to identification, estimation, and inference for treatment effects in regression discontinuity designs, however, do not use exogenous noise in the running variable to drive inference. Instead, following \citet*{hahn2001identification}, the dominant approach relies on a continuity argument. As in imbens2008regression, we assume potential outcomes $\cb{Y_i(0), \, Y_i(1)}$ such that $Y_i = Y_i(W_i)$. Then, we can identify a weighted treatment effect $\tau_c = \mathbb E\{Y_i(1) - Y_i(0) \mid Z_i = c\}$ via
provided that the conditional response functions \smash{$\mu_{(w)}(z) = \mathbb E\{Y(w) \mid Z = z\}$} are continuous. As we further explain in Section (ref), if we are willing to posit quantitative smoothness bounds on \smash{$\mu_{(w)}(\cdot)$}, then we can use this continuity-based argument to derive confidence intervals for $\tau_c$ with well understood asymptotics.
Despite its appeal and simple formulation, the continuity-based approach to regression discontinuity inference does not satisfy the criteria for credible causal inference as outlined by rubin2008objective. These criteria include outcome-free design and modeling randomness in the assignment mechanism rather than modeling the outcomes. rubin2008objective advocates for these principles to approximate the ideal of a design-based analysis of randomized controlled trials following neyman1923applications and rubin1974estimating. In contrast, the formal guarantees provided by continuity-based regression discontinuity analyses often take smoothness of \smash{$\mu_{(w)}(\cdot)$} as a primitive without explicitly modeling effective randomization in the assignment mechanism. While continuous measurement error in (or imprecise control of) the running variable implies continuity of \smash{$\mu_{(w)}(\cdot)$} lee2008randomized, this result is not used in estimation and inference.
Here we propose a new approach to regression discontinuity inference that goes back to the above qualitative argument used to justify regression discontinuity designs: Our approach directly exploits effectively random treatment assignment induced by noise in the running variable $Z_i$. Formally, we assume the existence of a latent variable $U_i$, and that the variation in the running variable $Z_i$ around $U_i$ is known and exogenous. For example, revisiting our educational setting, we can take $U_i$ to be a measure of the student's true ability; then the test score $Z_i$ is a noisy measurement of $U_i$ with well-documented psychometric properties. Likewise, in a medical setting, the running variable $Z_i$ may be a measurement of an underlying condition $U_i$ (e.g., CD4 counts); such diagnostic measurements often have well-studied test--retest reliability. In both cases, it is plausible that the measurements $Z_i$ are independent of relevant potential outcomes conditional on the underlying quantity $U_i$.
Our main result is that, under our assumptions, we can estimate weighted treatment effects that correspond to the effects of realistic changes to the existing treatment assignment rule. We then propose a practical approach to estimation and inference in regression discontinuity designs that builds on this result. Our approach is conceptually appealing, offers transparency on the key assumptions driving inference (noise-induced randomization and knowledge of noise mechanism), and allows for inference of policy-relevant estimands beyond (ref).
Our approach is motivated by settings where the researcher has limited understanding of the response variable $Y_i$ and the causal mechanism connecting $Y_i$ and the treatment $W_i$, but has substantive knowledge about the running variable $Z_i$. In this sense, our framework is akin to the model-X knockoff framework for controlled variable selection candes2018panning, which posits knowledge of the entire covariate distribution to facilitate inference of a poorly understood response variable conditionally on well-understood covariates.
We emphasize that, while this noise-induced randomization approach applies to many settings of interest, it does not apply to all regression discontinuity designs. Some running variables are not readily interpretable as having measurement error or other exogenous noise; or we may not have a-priori information on the distribution of this noise. For example, numerous studies have used geographic boundaries as discontinuities keele2014geographic, rischard2018bayesian, but it would be questionable to model the location of a household in space as having meaningful measurement error. Likewise, analyses of close elections, which are a central example of regression discontinuity designs in political science and economics caughey2011elections,lee2008randomized, may not allow for a natural noise model for $Z_i$ that would arise from, e.g., noisy counting of ballots, though perhaps there are other sources of exogenous noise gomez2007republicans,cooperman2017randomization. These considerations call attention to the limits of the proposed approach, but also highlight a difference in the foundational assumptions required for identification, estimation, and inference in regression discontinuity designs with a noisy running variable versus the assumptions required when the running variable is noiseless.
We consider the classical sharp regression discontinuity design with potential outcomes:
Our approach requires domain-specific knowledge about the distribution of the running variable.
Qualitatively, we interpret the latent variable $U_i$ in Assumption (ref) as a true measure of the property we want to use for treatment assignment, e.g., $U_i$ could capture ability in an educational setting or health in a medical one. The observed running variable $Z_i$ is then a noisy realization of $U_i$. The more noise there is in the running variable, the greater the effective randomization becomes, making our task easier. The assumption is flexible and can accommodate a wide range of noise models, including heteroskedastic noise and discrete running variables.
Assumption (ref) does not impose any restriction on $G$, which is a property of the studied population. While precise knowledge about the noise distribution $p(z \mid u)$ is a strong requirement, it enables us to obtain strong results (valid causal estimates justified via a form of effective randomization). In some applications, knowledge of the noise distribution may be available from test--retest data, prior modeling of item-level responses to tests, a physical model for the measurement device, biomedical knowledge, or direct control by the experimenter, e.g., in applications involving differential privacy dwork2013algorithmica. Any posited noise model should be carefully scrutinized since the credibility of the noise-induced randomization approach depends on the credibility of the noise model.
We also require for the additional noise to be exogenous. We formalize this requirement in terms of an unconfoundedness condition following rosenbaum1983central.
An implication of Assumption (ref) is that
where the \smash{$\alpha_{(w)}(u)$} are the response functions for the potential outcomes conditionally on the latent variable $u$. Following frangakis2002principal we can think of $u$ as indexing over unobserved principal strata; see also heckman2005structural.
A graphical illustration of our assumptions is presented in Fig. (ref). In view of Assumptions (ref) and (ref), the key argument for our strategy is captured by the following proposition.
We will apply this result by choosing functions $\gamma_+, \gamma_-$ and then averaging the response $Y_i(1)$ of treated units with weights $\gamma_+(Z_i)$ and the response $Y_i(0)$ of control units with weights $\gamma_-(Z_i)$. While there is no overlap between treated and control units in a sharp regression discontinuity design in terms of the running variable $Z_i$, Proposition (ref) establishes that by weighting treated units by $\gamma_+$ and control units by $\gamma_-$ we may achieve balance in the latent variable, if $\hkernel{\cdot}{\gamma_+} \approx \hkernel{\cdot}{\gamma_-}$.
As discussed above, the dominant approach to inference in regression discontinuity designs is via continuity-based arguments that build on (ref). Perhaps the most popular continuity-based approach is to use local linear regression to estimate the treatment effect (ref) at $Z_i = c$. This approach can be used for valid estimation and inference of $\tau_c$ provided the functions $\mu_{(w)}(z) = \mathbb E\{Y_i(w) \mid Z_i=z\}$ are smooth and the local linear regression bandwidth decays at an appropriate rate; the rate of convergence of $\htau_c$ and appropriate choice of bandwidth depend on the degree of smoothness assumed. Notable results in this line of work, covering topics such as robust confidence intervals and data-adaptive bandwidth choices, include armstrong2016simple, calonico2014robust and imbens2011optimal, as well as Bayesian approaches branson2019nonparametric,geneletti2015bayesian. More recently, extensions have been considered to the continuity-based approaches that improve over local linear regression by directly exploiting the assumed smoothness properties of $\mu_{(w)}(\cdot)$. Under the assumption that $\mu_{(w)}(\cdot)$ belongs to a convex class, e.g., \smash{$|\mu''_{(w)}(z)| \leq B$} for all $z \in \RR$, armstrong2018optimal and imbens2019optimized use numerical convex optimization to derive minimax linear estimators of \smash{$\tau_c$}.
One alternative approach to inference in regression discontinuity designs, which cattaneo2015randomization, li2015evaluating and mattei2016regression refer to as local randomization inference, starts by positing a non-trivial interval $\mathcal{I}$ with $c \in \mathcal{I}$, such that
They then focus on the subset of units with $Z_i \in \mathcal{I}$, and perform classical randomized study inference on this subset. Unlike continuity-based analysis, this approach is design-based in the sense of rubin2008objective. In practice, however, the assumption (ref) is often unrealistic and limits the applicability of methods relying on it sekhon2017interpreting. A testable implication of (ref) is that $\mu_{(w)}(z)$ should be constant over $\mathcal{I}$ for both $w = 0$ and $1$, but this structure rarely plays out in the data. One may try to fix this issue by de-trending outcomes and assuming (ref) on the residuals sales2020limitless; however, such an approach relies on correct specification of the trend removal, and is thus no longer justified by randomization. Furthermore, it is not clear how to choose the interval $\mathcal{I}$ used in (ref) via the types of methods typically used for regression discontinuity inference. There's no data-driven way of discovering an interval $\mathcal{I}$ over which (ref) holds that is itself justified by randomization; conversely, if the interval $\mathcal{I}$ is known a-priori, then the problem collapses to a basic randomized controlled trial where the regression discontinuity structure is not used for inference.
The idea that explicit structural modeling is valuable for causal inference has a long tradition in economics, going back to roy1951some and heckman1979sample, with recent developments by e.g., heckman2005structural, brinch2017beyond and mogstad2018using. At a high level, our work can be seen as connecting this tradition to the regression discontinuity design, and demonstrating how structural assumptions enable inference of policy-relevant causal estimands.
Knowledge of the presence of measurement error (or other noise) in running variables is often mentioned bor2014regression,bor2017treatment,harlow2020impact,lee2008randomized, yet this side-information is typically not directly used for inference. In a rare quantitative use of information about measurement error, fraga2016examining use margin of error statistics provided by the Census Bureau for the fraction or size of a voting-aged population that has limited English proficiency; they report some analyses using only units that are within a 90% margin of error of the cutoff.
Closer to our approach, rokkanen2015exam considers the regression discontinuity design under Assumptions (ref) and (ref). Instead of assuming prior knowledge of the noise distribution $p(\cdot \mid u$), rokkanen2015exam assumes that for each unit we observe at least two noisy measurements $Z_i', Z_i''$ of the underlying latent variable $U_i$ in addition to the running variable $Z_i$. While rokkanen2015exam provides conditions for the nonparametric identification of $\alpha_{(w)}(\cdot)$ in (ref) and consequently of treatment effects, the estimation and inference strategy posits strong parametric assumptions, namely joint normality of $(U_i, Z_i, Z_i', Z_i'')$ and linearity of $\alpha_{(w)}(u)$ as a function of $u$. In contrast, we assume knowledge of the noise distribution through, e.g., biomedical knowledge or test--retest data, however we impose no parametric restrictions on $G$ and $\alpha_{(w)}(u)$. Furthermore, we develop a practical and intuitive method for estimation and inference, that provides valid coverage even when treatment effects are only partially identified (e.g., when $p(\cdot \mid u)$ is finitely supported).
Our results are also connected to research on treatment effect estimation under biased or risk-based allocation robbins1989estimating, robbins1991estimating, finkelstein1996partB motivated from an empirical Bayes interpretation of the noise model in Assumption (ref). For example, robbins1989estimating study treatment effect estimation under what effectively amounts to our Assumptions (ref) and (ref) with Gaussian noise, \smash{$Z_i \mid U_i \sim \nn(U_i, \, \nu^2)$}, and control potential outcomes linked to $U_i$ via an additive shift, \smash{$\alpha_{(0)}(u) = \mathbb E\{Y_i(0) \mid U_i = u\} = u + s$} for \smash{$s \in \RR$}. These assumptions on $\alpha_{(0)}(u)$ are motivated by settings where measurements of the same quantity function as both the running variable and the outcome, such as the application in finkelstein1996partB where patients with high cholesterol are given a drug to lower cholesterol, and we are interested in the drug's effectiveness. However, such assumptions on control outcomes are not appropriate in the examples considered in this paper. Thus, while this line of work presents a notable yet largely overlooked chapter in the history of regression discontinuity designs cook2008waiting, it does not provide a methodological baseline for our approach.
Finally, we contrast our setup with a line of work that studies the regression discontinuity design when the running variable is unobserved, and instead a noisy measurement thereof is observed davezies2017regression, dong2021can; see the causal diagram in Supplementary Fig. (ref) for an illustration. Identification becomes subtle and estimation can be difficult because of the difficulties of nonparametrics with measurement error meister2009deconvolution. Instead, we use measurement error as our identifying assumption; the noise in our setup is beneficial for our estimation strategy rather than a barrier (and we observe the running variable).
Motivated by Proposition (ref), we consider ratio-form estimators,
where $\gamma_+, \gamma_-$ are pre-specified weighting functions such that $\gamma_+(z) = 0 \text{ for } z < c$, $\gamma_-(z) = 0 \text{ for } z\geq c$. The class (ref) is a broad and intuitive class of estimators that includes, for example, the difference-in-means of units that are close to the cutoff (with the choice $\gamma_+(z) = \ind\{z \in [c, c+h]\}$ and $\gamma_-(z) = \ind\{z \in [c-h, c)\}$ for $h>0$).
Our goal is to conduct inference for weighted treatment effects,
where $\tau(u)$ is the conditional average treatment effect (CATE) of the stratum with $U_i = u$,
and $w(\cdot)$ is a latent weighting (i.e., $w(\cdot)$ assigns weight to the latent $U$) chosen by the analyst.
In the following, we take $\gamma_+, \gamma_-$ as pre-specified by the researcher and seek to understand how to use the point estimate $\htau_{\gamma}$ in (ref) to form valid confidence intervals for $\tau_w$ in (ref) by accounting for potential bias. In Section (ref), we make a concrete recommendation for choosing $\gamma_+, \gamma_-$.
We first derive the asymptotic limit of $\htau_{\gamma}$ with fixed $\gamma_+(\cdot), \gamma_-(\cdot)$ given $n$ i.i.d. copies of $(U_i, Z_i, Y_i(0), Y_i(1))$ satisfying Assumptions (ref)-(ref).
In view of Theorem (ref) and the definition of $\tau_w$ in (ref), we derive an asymptotic decomposition of the bias in estimating $\tau_w$ through $\htau_{\gamma}$:
The bias decomposes into two terms. The first term (confounding bias) describes how well we are balancing units through their latent variable $u$ and will be small if $\hkernel{\cdot}{\gamma_+} \approx \hkernel{\cdot}{\gamma_-}$. The second term, which we call heterogeneity bias, is equal to zero when the conditional average treatment effect $\tau(u)$ is constant as a function of $u$, or when $\hkernel{u}{\gamma_+}=w(u)$ for all $u$.
We now provide examples of estimands that may be expressed as weighted treatment effects (ref).
In Supplement (ref), we also consider an estimand motivated by a policy intervention that involves reducing measurement error in the running variable.
In the previous section, we discussed the asymptotic limit of the ratio-form estimator in (ref) and the bias in estimating weighted treatment effects in regression discontinuity designs. To make use of such an estimator in practice, however, we also need to understand its sampling distribution and to control the bias. In this section, we describe our approach to inference.
We make the following additional assumption:
We start by studying the asymptotic distribution of the ratio-form estimator in (ref). We treat $\gamma_+, \gamma_-$ as deterministic but (in contrast to Theorem (ref)) allow them to vary with $n$, i.e., \smash{$\gamma_+ = \gamma_+^{(n)}$} and \smash{$\gamma_- = \gamma_-^{(n)}$}. Our first formal result is the following central limit theorem.
The assumption on $\gamma_+,\gamma_-$ is satisfied by the weighting functions proposed in Section (ref) and other choices. For example, the local difference-in-means estimator with $\gamma_+(z) = 1\{z \in [c, c+h_n]\}$, $\gamma_-(z) = \ind\{z \in [c-h_n, c)\}$ meets the assumption when $h_n^{-1} = O(n^{\beta})$ for $\beta \in (0,1/2)$, $\lambda$ is the Lebesgue measure and $\cb{p(\cdot \mid u)}_u$ are uniformly bounded and equicontinuous at $c$.
To construct confidence intervals for $\tau_w$ in (ref), we first estimate the asymptotic variance \smash{$V_\gamma$}.
Second, we account for the potential bias \smash{$|b_{\gamma}| = |\taugamma - \tau_w|$}. We do not assume negligible bias (i.e., undersmoothing) and accommodate settings wherein treatment effects are only partially identified, and bias does not decay to zero even asymptotically imbens2019optimized (e.g., when \smash{$Z_i \mid U_i$} has a binomial distribution). To do so, we derive an upper bound \smash{$\hB_{\gamma}$} for the bias $|b_{\gamma}|$. A challenge is that the expectations in Corollary (ref) involve integrals over the latent variable $U_i$ and the unknown functions $G$, $\tau(\cdot)$ and $\alpha_{(0)}(\cdot)$. Taking a clue from ignatiadis2019bias, we bound the worst-case bias over any data-generating distribution consistent with the observed data for the running variable $Z_i$. Define the marginal distribution function of $Z_i$ (marginalizing over $U_i \sim G$), $F_G(t)=\smallint 1(z \leq t) \smallint p(z\mid u)dG(u)d\lambda(z)$, and let $\mathcal{G}_n$ be the class of latent variable distributions with marginal distribution $F_G$ inside the Kolmogorov-Smirnov band massart1990tight centered at the empirical distribution $\widehat{F}_n(t) = \sum_{i=1}^n 1(Z_i \leq t)/n$,
We also consider a sensitivity model for treatment effect heterogeneity. For $M \in [0,\,1]$, we let
Above, $\mathcal{T}_0$ consists of all conditional average treatment effect (CATE) functions $\tau(\cdot)$ that are constant as a function of $u$. Under Assumption (ref), $\mathcal{T}_1 = \cb{\text{all CATE functions } \tau(\cdot)},\; \mathcal{T}_{1/2} \supset \cb{\text{all CATE functions } \tau(\cdot) \geq 0 },$ and so the choice $M=1$ avoids imposing any additional assumptions on heterogeneity, while $M=1/2$ is a conservative choice under the monotonicity restriction $\tau(\cdot) \geq 0$.
In Supplement (ref) we explain how to compute this bound on the bias. Finally, we build confidence intervals for $\tau$ that are robust to estimation bias up to \smash{$\hB_{\gamma,M}$} following imbens2004confidence, armstrong2018optimal, and imbens2019optimized.
Our method requires specifying the sensitivity model (ref). While one can adopt the unrestrictive model $\mathcal{T}_1$, we explore the robustness of our approach to misspecification of the sensitivity model: We suppose that $\tau(\cdot)$ is not constant as a function of $u$, yet we conduct inference using $\mathcal{T}_0$. In this case, our intervals attain the correct coverage for a convenience-weighted treatment effect.
The convenience-weighted treatment effect $\tau_{h,+}$ may be of interest if we are not directly interested in treatment heterogeneity crump2009dealing,li2016balancing,imbens2019optimized,kallus2020generalized. If we are interested in the null hypothesis of no treatment effects, $H_0: \tau(u) =0 \text{ for all }u,$ then we can form a valid test by forming confidence intervals for $\tau_w$ under the sensitivity model $\mathcal{T}_0$ and rejecting the null hypothesis when the resulting confidence interval does not include $0$.
In Supplement (ref), we study the asymptotic bias for our method when there is no effective randomization at all, but we proceed pretending Assumptions (ref) and (ref) hold for a given $p(\cdot \mid \cdot)$. Our result requires a strong functional form assumption according to which $\mu_{(w)}(z) = \mathbb E\{Y_i(w) \mid Z_i=z\}$ is a linear combination of $p(z \mid u)/f(z)$, where $f$ is the $d\lambda$-density of $Z$. The result implies consistency and valid inference (no matter the noise model we posit) when unbeknownst to us, $\mu_{(w)}(z)$ is constant as a function of $z$.
Given a choice of weighting functions $\gamma_{\pm}$ for (ref), Section (ref) provides a complete recipe for building valid confidence intervals justified by the noise-induced randomization framework. For example, one could take weighting functions implied by various regression discontinuity estimators. Existing weighting functions $\gamma_{\pm}$, however, were not designed for our framework, and so may not yield particularly short confidence intervals. Hence we turn to deriving weighting functions $\gamma_{\pm}$ with an eye towards making confidence intervals obtained via Corollary (ref) short.
Our (heuristic) strategy is to choose $\gamma_{\pm}$ by minimizing an approximate bound on the worst-case mean-squared error of the estimator in (ref). Let $w(\cdot)$ be the latent weighting of the estimand (ref) and suppose we posit the sensitivity model $\mathcal{T}_M$. Furthermore, let $\bar{F}(\cdot)$ be a guess or estimate of the marginal distribution $F_G(\cdot)$ of $Z_i$ under Assumption (ref) and let $\bar{w}(\cdot)$ be an estimate of the normalized latent weighting $w(\cdot)/\mathbb E_G\{w(U_i)\}$. We propose solving the following quadratic program (of which an appropriately discretized version can be solved using standard convex optimization software, e.g., MOSEK, mosek):
In choosing $\bar{F}(\cdot)$ and $\bar{w}(\cdot)$, we make use of the structure provided by Assumption (ref), and estimate $G$ as $\bar{G}$ via nonparametric maximum likelihood kiefer1956consistency and then we let $\bar{F}(\cdot) = F_{\bar{G}}(\cdot)$ and $\bar{w}(\cdot) = w(\cdot)/\mathbb E_{\bar{G}}\{w(U_i)\}$.
The first term in (ref) is a variance proxy for our estimator, motivated by the inequality $\operatorname{Var}(\gamma_{\diamond}(Z_i)Y_i) \leq \smallint \gamma_{\diamond}^2(z) \, dF(z)$ for $\diamond \in \cb{\pm}$. The second term, $(t_1 + t_2)^2$, approximately bounds the worst-case bias. The bias is decomposed through the triangle inequality into two terms resembling the bias decomposition of Corollary (ref); $t_1$ in (ref) seeks to bound the confounding bias and balances $\hkernel{\cdot}{\gamma_+}$ and $\hkernel{\cdot}{\gamma_-}$, while $t_2$ seeks to bound the heterogeneity bias and balances $h$ with the normalized $w(\cdot)$. Next, (ref) is a normalization constraint, and (ref) enforces that $\gamma_+$ ($\gamma_-$) assigns weight only to treated (control) units. Constraint (ref) ensures that no single observation has excessive influence (we omitted this constraint in our numerical implementation).
The following proposition shows that the weighting functions $\gamma_{\pm}$ derived from optimization problem (ref) satisfy the conditions of Theorem (ref) and thus enable valid inference.
In our implementation, we use all running variables $Z_i$ (but not the responses $Y_i$) to form estimates for $\bar{F}(\cdot)$ and $\bar{w}(\cdot)$; throughout our simulations we have not observed any undercoverage thereby. We summarize our approach in Algorithm (ref).
In this section, we apply our approach to a medical study. bor2017treatment study $11,306$ patients in South Africa (in 2011--2012) diagnosed with HIV aiming to understand whether immediate initiation of antiretroviral therapy (ART) helps retain patients in the medical system. The response $Y_i \in \cb{0,1}$ is an indicator of the $i$-th patient's retention, measured by the presence of a clinic visit, lab test, or ART initiation 6 to 18 months after the initial HIV diagnosis.
According to health guidelines used in South Africa at the time, an HIV-positive patient should receive immediate ART if their measured CD4 count was below $350\; \text{cells} /\mu L$ (a low CD4 count is indicative of poor immune function), lending itself to a natural regression discontinuity design for intention-to-treat effects. Figure (ref)(a) shows a histogram of the running variable $Z_i$, the log CD4 count (in cells/$\mu L$), with treatment cutoff $c = \log(350)$ denoted by a dashed line.
bor2017treatment emphasize that CD4 count measurements are noisy; causes of this noise include instrument imprecision and variability in the blood sample taken hughes1994within,wade2014multicenter. They then use the existence of such noise to qualitatively argue that treatment $W_i = \ind(Z_i < c)$ is effectively random close to the cutoff $c$, thus strengthening the credibility of the regression discontinuity analysis.
Here, we seek an approach to estimating the effect of ART on retention that is driven by the effective treatment randomization provided by the measurement error in $Z_i$ and (approximate) knowledge of the noise mechanism. To this end, we start by modeling this measurement error. venter provide pairs of repeated measurements \smash{$Z_i, \, Z_i'$} of the log CD4 count on 553 individuals (with measurements taken in the same laboratory). Figure (ref)(b) compares a histogram of the normalized differences \smash{$(Z_i-Z_i')/\surd{2}$} on the data of venter to a fitted Gaussian probability density function with noise $\hnu = 0.19$. We estimate the noise level $\hnu = 0.19$ using a robust method that ignores outliers by Winsorizing the smallest and largest 5% of the normalized differences \smash{$(Z_i-Z_i')/\surd{2}$} and rescales to maintain unbiasedness under Gaussian noise. The choice of Winsorization is motivated by the robustness of our method to underestimation of the noise level, as explained in Example (ref) and further demonstrated in the simulations of Section (ref) below.
The modeling choice \smash{$Z_i \mid U_i \sim \mathcal{N}(U_i, \,\hnu^2)$} illustrates our approach and serves as a starting point for analysis. We emphasize, however, that our method enables an epidemiologist to posit a more realistic model of the noise density $p(z \mid u)$ based on scientific understanding of CD4 counts and additional datasets with repeated measurements. Our causal identification strategy will be most credible when the scientist has considerable knowledge about the noise mechanism. Henceforth in applying our approach, we assume that measurement error in the log CD4 counts follows \smash{$Z_i \mid U_i \sim \mathcal{N}(U_i, \,\hnu^2)$}, where $U_i$ is the true underlying log CD4 count of patient $i$. Given this noise model, we apply our noise-induced randomization (NIR) approach, with sensitivity model $\mathcal{T}_0$ to test for the existence of any treatment effects (as explained after Corollary (ref)).
As a first comparison point, we consider treatment effect estimates obtained via the continuity-based approach proposed by calonico2014robust, which has recently become popular in applications. This approach involves first fitting the regression discontinuity parameter via local linear regression, and then estimating and correcting for its bias in a way that's asymptotically justified under higher-order smoothness assumptions calonico2014robust. We implement this approach via the R package rdrobust of calonico2015rdrobust with default tuning parameters.
As a second baseline, we consider the minimax linear inference approach developed by armstrong2018optimal,armstrong2016simple, imbens2019optimized and kolesar2018inference; we use the R package optrdd of imbens2019optimized. This approach posits a constant $B$ such that \smash{$|\mu''_{(w)}(z)| \leq B$} for all $w \in \cb{0, \, 1}$ and $z \in \RR$, and then provides intervals that are robust to the worst-case bias under the curvature bound. The main difficulty in using this approach is in choosing the curvature bound $B$. We use the heuristic considered in armstrong2016simple: We fit fourth-degree polynomials to $\mu_{(0)}(\cdot)$ and $\mu_{(1)}(\cdot)$, and take the largest estimated curvature obtained anywhere. Relative to rdrobust, the minimax linear inference approach makes explicit how smoothness is used for inference (i.e., if one believes in the proposed curvature bound $B$, one should also believe in the resulting intervals). In contrast, rdrobust relies more directly on asymptotics justified by higher-order smoothness; see calonico2018effect for further discussion.
We present the results in Table (ref). All displayed confidence intervals are significant at the 95% level. What differs is the assumptions we need to justify these confidence intervals. The baseline methods given here rely on quantifying the smoothness of the $\mu_{(w)}(\cdot)$ in a data-driven way; and the credibility of the resulting intervals hinges on how well we believe this task can be accomplished. In contrast, our NIR intervals are directly justified by the posited measurement error model for the running variable $Z_i$ and the induced effective treatment randomization.
Whether practitioners prefer the NIR intervals or the continuity-based alternative will likely depend on their intended use. Here, the continuity-based intervals are shorter than the NIR intervals, which is desirable in settings where precision is at a premium. (In the simulation study, we show examples where the NIR intervals are shorter.) On the other hand, the NIR intervals are explicitly justified using a form of effective treatment randomization, and thus may be seen as getting closer to best practices for credible causal inference as outlined by rubin2008objective. In some settings, practitioners may want to report both: One could see the NIR intervals as conservative intervals that may sustain stricter scrutiny in terms of identification (by scrutiny of the assumed noise model), and the continuity-based ones as sharper intervals if one is willing to rely on data-driven smoothness estimation.
Finally, for intuition, in Fig. (ref) we show the weighting functions $\gamma_\pm$ selected via quadratic programming and that were used by the NIR approach (Section (ref)), and the implied latent weighting $h(\cdot, \gamma_+)$, $h(\cdot, \gamma_-)$ as per (ref). Units with $Z_i$ close to the cutoff are strongly upweighted, and so we achieve approximate balance in terms of the latent $U_i$. The oscillations of the weighting functions $\gamma_{\pm}$ near the cutoff arise due to higher order bias corrections in nonparametric estimation and are common also for local linear regression estimates when represented as weighted averages (see, e.g, gelman2019high).
The two types of intervals discussed above may appear to rely on incomparable identification strategies. However, we can build a formal bridge connecting them. One can verify, that in the presence of Gaussian measurement error, the functions $\mu_{(w)}(z) = \mathbb E\{Y_i(w) \mid Z_i=z\}$ must be smooth. Under Assumptions (ref)--(ref), $\mu_{(w)}(z) = \textstyle \int \alpha_{(w)}(u)p(z\mid u)\,dG(u) \,\big/\,\int p(z\mid u)\,dG(u)$, so if $\alpha_{(w)}(\cdot)$ is bounded and $z \mapsto p(z \mid u)$ is continuous, then by the dominated convergence theorem we can show that $\mu_{(w)}(\cdot)$ is also continuous lee2008randomized. Furthermore, higher order differentiability of $p(\cdot \mid u)$ implies the same for $\mu_{(w)}(\cdot)$ dong2021can. Here, we will investigate this connection to gain further insights on the relationship between noise-induced-randomization and continuity-based methods shown above.
To this end, we define the worst-case possible curvature at $z$ among all data-generating distributions satisfying Assumptions (ref)--(ref) with conditional density $p(\cdot \mid \cdot)$ such that the marginal density of the running variable at $z$ is lower bounded by $\rho > 0$:
In (ref) we constrain ourselves to marginal densities such that $f_G(z) \geq \rho$ for $\rho>0$, because typically $\Curv(z,\,0, \, p) = \infty$. In Supplement (ref), we explain how the quantity (ref) may be computed numerically for any sufficiently regular $p$. One can then use the upper bounds on the second derivative of \smash{$\mu_{(w)}(\cdot)$} in (ref) in conjunction with, e.g., the estimators of imbens2019optimized and armstrong2016simple that provide uniform inference for the regression discontinuity parameter given a curvature bound on the response function.
To provide intuition for (ref), we provide analytic lower and upper bounds on (ref) in the case of Gaussian measurement error, i.e., with \smash{$Z_i \mid U_i \sim \nn(U_i, \, \nu^2)$} that quantify dependence on the noise level $\nu$ and the lower bound $\rho$ on the density.
We now return to the application of bor2017treatment. Recall that we assumed a measurement error with noise $\hnu = 0.19$. We estimate the density of the running variable at the cutoff as $\hat{f}(c)= 0.57$ using the nonparametric maximum likelihood estimator. Using optrdd with curvature parameter \smash{$\Curv\{c,\,\hat{f}(c), \, \mathcal{N}(\cdot,\, 0.19^2)\}$} (equal to $31.3$) yields intervals that are directly justified by our noise model, just like noise-induced randomization. The resulting 95% confidence interval is $\tau \in (0.071 \pm 0.130)$ which is much wider than any of the intervals reported in Table (ref).
The reason the optrdd intervals are wider than the continuity-based intervals in Table (ref) is that our noise model implies much less continuity than is discovered by the data-driven methods. For example, our noise model guarantees a curvature bound of $B = 31.3$, whereas the heuristic of armstrong2016simple gives a curvature bound of $B = 1.46$ (i.e., it finds the function to be 20x smoother than guaranteed by measurement error). This highlights the extent to which our proposal (and other methods justified by measurement error alone) can be seen as stricter than continuity-based alternatives in terms of the information used to estimate treatment effects.
We next consider the behavior of our method in a semi-synthetic regression discontinuity design built using data from the Early Childhood Longitudinal Study ECLS. This dataset has scaled mathematics test scores for $n = 18,174$ children from kindergarten to fifth grade. Furthermore, each test score is accompanied by a noise variance obtained via item response theory; see ECLS for further details.
Each sample $i = 1, \dotsc, n$ is built using the sequence of test scores from a single child. We set the running variable $Z_i$ to be the child's kindergarten spring semester score, and set treatment as $W_i = \ind(Z_i \geq c)$ for a cutoff $c = -0.2$. We set control potential outcomes $Y_i(0) \in \cb{0, \, 1}$ to indicate whether the child's score was above $a = 0.5$ in spring semester of their first grade, while $Y_i(1) \in \cb{0, \, 1}$ measures the same quantity in spring semester of their second grade; these are analogous to typically studied outcomes such as passing subsequent examinations. Thus, the treatment effect $Y_i(1) - Y_i(0)$ measures the child's improvement in passing the test (i.e., clearing the cutoff $a = 0.5$) between first and second grades.
As shown in Fig. (ref), there is considerable heterogeneity in the regression discontinuity parameter $\tau_{c'} = \mathbb E\{Y_i(1) - Y_i(0) \mid Z_i = c'\}$ as we vary $c'$ away from the cutoff: For children with either very good or very bad values of $Z_i$ the treatment effect is essentially 0 (since they will pass or, respectively, fail to pass the cutoff $a$ in both first and second grade with high probability), while for students with intermediate values of $Z_i$ there is a large treatment effect. We chose the parameters $a$ and $c$ in our semi-synthetic construction to accentuate this type of heterogeneity.
For this problem, the possibility and credibility of inference with noise-induced randomization builds on item response theory (IRT), a widely used model in educational testing. According to IRT, the $i$-th child's test score $Z_i$ is a noisy reflection of their true ability with (approximate) Gaussian measurement error with variance $\nu_i^2$ which is also determined by IRT. Below, we assume Gaussian errors in the running variable, i.e., \smash{$Z_i \mid U_i \sim \nn(U_i,\,\hnu^2)$}, and, following Example (ref), we set $\hnu = \min_i \{\nu_i\} = 0.2043$ to match the lowest noise estimate provided in the dataset. We run our method using sensitivity model $\mathcal{T}_{0.5}$. In this application, the monotonicity restriction $\tau(\cdot) \geq 0$ appears plausible, since the treatment effect measures the child's improvement between first and second grades, and as explained after (ref), $\mathcal{T}_{0.5}$ does not place further restrictions on treatment effect heterogeneity. We also construct confidence intervals centered at the same point estimates under $\mathcal{T}_{0.3}$; this sensitivity model is plausible based on the treatment heterogeneity in the ground truth individual treatment effects (Fig. (ref)).
Our main question is whether our procedure is able to estimate this heterogeneity, i.e., whether it can accurately recover variation in treatment effects away from the cutoff. To this end, we consider two statistical targets: First, we consider estimation of the regression discontinuity parameter (ref) at $c'$ away from the cutoff, and second, the policy-relevant parameter (ref) quantifying the effect of changing the cutoff from $c$ to \smash{$c'$}. Results for both targets are shown in Fig. (ref). Our method is able to recover heterogeneity. In both cases, the confidence intervals cover the ground truth. They are narrowest near the cutoff $c = -0.2$, and get wider as we move away from the cutoff.
To complement the picture given by our applications, we consider a simulation study to assess the performance of our method in terms of its accuracy and coverage. We first consider a data-generating distribution with null treatment effects $\tau(u)=0$ wherein $Z_i$ has discrete support, and has a binomial distribution conditionally on the latent $U_i$. For $i=1,\dotsc,n$, we generate,
where the number of trials $K$ and number of samples $n$ are simulation parameters and $c^*=0.6$.
We compare the following point estimates and $95\%$ confidence intervals for the (null) treatment effect:
We evaluate methods by computing the confidence interval coverage, the expected half-length of confidence intervals and the mean absolute error. These metrics are computed by averaging over 1,000 Monte Carlo replications.
The results of the simulation study are shown in Table (ref). All methods have approximately correct coverage, with optrdd and NIR always achieving the nominal $95\%$ level and rdrobust slightly undercovering. Although rdrobust and its distributional theory have been developed under the assumption of a continuous rather than discrete running variable, it performs reasonably well. For small $K$ and $n$, rdrobust sometimes returns an error, in which case we do not report its performance. NIR yields the shortest confidence intervals in most settings. The number of trials $K$ determines the effective noise level. Our method performs best when there is more effective noise in the running variable, that is, when $K$ is small ($K\leq 25$). This is in contrast to rdrobust, whose performance improves as $K$ increases and the running variable becomes less discrete, until at $K=200$ it leads to shorter confidence intervals than NIR. As expected, the confidence interval length decreases for all methods as the sample size $n$ increases.
This simulation experiment corroborates the claim that our method, NIR, can flexibly turn assumptions about exogenous noise in the running variable $Z_i$ into a practical procedure for inference in regression discontinuity designs. We achieve nominal coverage across simulation settings. Our results also point to the possibility that NIR may result in improved power in settings where running variables are discrete with known noise. This would not be unreasonable, as continuity-based approaches were not necessarily designed for this setting (although, as discussed in kolesar2018inference they can rigorously be used given appropriate interpretation).
We next explore the impact of misspecification of the noise model on the performance of NIR. We fix the sample size as $n=10,000$ and generate for $i=1,\dotsc,n$: $U_i \sim \nn(0, 1)$, $Z_i \mid U_i \sim p(\cdot \mid U_i)$, $W_i = \ind(Z_i \geq 0)$. We consider the following three location-scale models for \smash{$p(\cdot \mid U_i)$}: Gaussian, t with 6 degrees of freedom, and Laplace. In each case, the location is equal to $U_i$ and the scale is such that $\operatorname{Var}(Z_i \mid U_i) = 0.5^2$. The response is generated as in (ref) with $c^*=0$.
For each simulation setting, we compare the following methods: rdrobust, and noise-induced randomization (NIR) with noise model $\nn(U_i, \nu^2)$ for $\nu \in \cb{0.3, 0.5, 0.7, 0.9}$ and the sensitivity class $\mathcal{T}_0$. NIR is well-specified only in one case: when it is applied with noise level $\nu=0.5$ and data is generated according to the Gaussian location-scale model.
Our evaluation proceeds as in Section (ref) and the results are shown in Fig. (ref). Rdrobust performs well across all three scenarios and sets a benchmark (even if it has some undercoverage with Laplace noise). With this standard in mind, we discuss the robustness of our proposed approach when confronted with a misspecified noise model. We begin by examining the situation where the true noise model is Gaussian. Then, NIR has the correct (95$\%$) coverage for $\nu \in \cb{0.3, 0.5, 0.7}$. Our theoretical results provide justification for $\nu=0.5$ (well-specification), and $\nu=0.3$ (underestimated noise-level as described in Example (ref)). Coverage for $\nu=0.7$ in the simulation is not justified theoretically, but demonstrates some robustness of NIR to the specification of the noise level. On the other hand, for $\nu=0.9$, the coverage of NIR drops to roughly 65$\%$, showing that NIR is not robust to substantial overestimation of the noise level. The expected half-length of NIR confidence intervals is decreasing in $\nu$. The mean absolute error exhibits a trade-off behavior, being minimized at $\nu = 0.7$, decreasing before this point, and then increasing thereafter. This suggests that for point estimation (rather than inference), $\nu$ acts similarly to a standard bias-variance trade-off parameter. NIR is also moderately robust to misspecification of the shape of the noise distribution: NIR with $\nu=0.3$ attains 95$\%$ coverage with Laplace noise, and NIR with $\nu \in \cb{0.3,0.5}$ attains nominal coverage with t-noise. However, when both the noise level and the shape of the noise distribution are strongly misspecified, the coverage of NIR can be very low.
The results of this simulation study suggest that NIR is robust to moderate misspecification of the noise model. In applications, one should err toward underestimating the noise level if possible.
All numerical results in this paper are reproducible with the code in the following Github repository: \url{https://github.com/nignatiadis/noise-induced-randomization-paper}.\\ We provide an implementation of NIR as a package in the Julia programming language bezanson2017julia that depends, among others, on JuMP.jl DunningHuchetteLubin2017.