Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
55,587 characters · 16 sections · 55 citation commands
Semiparametric Inference for Regression-Discontinuity Designs
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Causal inference; Regression discontinuity designs; Partially linear model; $P$-spline; Two-Stage Estimator.
\spacingset{1.8}
\pagenumbering{arabic} \setcounter{page}{1}
The regression discontinuity design (RDD) was initially introduced by TC:60 to examine the influence of merit awards on students' future academic outcomes. Since then, it has evolved into one of the most widely used strategies for estimating treatment effects across disciplines such as economics, political science, biology, and medicine. In an RDD, units receive scores, and treatment allocation is determined by whether the score exceeds a predetermined cutoff value: units scoring above the cutoff are designated to the treatment condition, while those scoring below are assigned to the control condition. This treatment assignment rule creates a discontinuity in the probability of receiving treatment, allowing researchers to estimate the treatment effect by comparing units barely above and barely below the cutoff.
The local linear or quadratic approximation Fan:92,FG:96 is a commonly employed approach in RDDs Hahn:01. Specifically, researchers typically estimate the treatment effect locally around the cutoff, with the estimates depending on the selected bandwidth. Studies have investigated optimal bandwidth selection for local linear and quadratic regression specifications IK:12, CCT:14, CCT:15. Another approach explored by researchers is the use of a global polynomial regression approach LL:10. In practice, researchers often employ high-order polynomials (up to the fifth or sixth degree), with the polynomial degree selected using statistical information criteria or cross-validation. However, GI:19 argued against the use of high-order global polynomial approximations, instead advocating for inference based on local low-order polynomials, such as local linear or quadratic models. For comprehensive reviews of the RDD literature, refer to IL:08 and CIT:19.
Hahn:01 demonstrated that identifying the average treatment effect at the cutoff in RDDs relies on the unconfoundedness assumption. They also showed that, in the absence of this assumption, only the local average treatment effect at the cutoff can be identified. In this article, we delve into a semiparametric framework for RDDs, proposing a global approach that identifies the average treatment effect when the unconfoundedness assumption is violated. Our primary focus is on identifying the average treatment effect within a partially linear model framework. This entails capturing the treatment effect through the coefficient of a specified polynomial function of the assignment probability (or propensity score) in the linear component of the model.
This result, stemming from the unique identification structure, underpins our semiparametric inference for RDDs, offering a more efficient estimation of the treatment effect at the cutoff compared to the earlier global polynomial regression approach LL:10. By employing $P$-spline to approximate the nonparameteric component, we reformulate the semiparametric RDD framework into a regression based on a linear mixed-effects model. We then derive the treatment effect estimate within this linear mixed-effects model framework. Theoretically, we establish the asymptotic normality of the estimate under regular conditions, highlighting the advantages of our approach, and propose an inference procedure. Finally, leveraging the asymptotic normality, we optimize the form of the polynomial function of the propensity score, which enhances the accuracy of our inference.
Empirically, our experiments on simulated datasets support the assertion that our global approach provides a viable and superior alternative to local estimation methods. Furthermore, we apply our global approach to evaluate the effect of antihypertensive treatment on reducing the risk of cardiovascular disease calonico2024regression, and to assess the impact of incumbent advantage in U.S. House and Senate elections Lee:08,CFT:15. Our analysis of these real-world datasets demonstrates that our global estimation approach yields more reliable inference.
The rest of the article proceeds as follows. Section (ref) outlines the problem setting, and Section (ref) introduces a semiparametric framework for identifying the average treatment effect at the cutoff. Section (ref) presents estimation and inference procedures, demonstrating that both the bias and variance are approximately negligible compared to those of the local method. Section (ref) shows the simulation results comparing our approach with the local method, while Section (ref) evaluates its performance on real datasets. Concluding remarks are provided in Section (ref). Proofs are gathered in the appendix.
Throughout this paper, we consider the following problem setup. We have a random sample $(X_i, W_i, Y_i)$ of individuals that are independently distributed, with the individuals in the sample indexed by $i=1,\cdots,n$. Using the potential outcome approach, let $(Y_i(0),Y_i(1))$ denote the pair of potential outcomes for individual $i$, and $W_i\in\{0,1\}$ represent the treatment received. Each individual $i$ is assigned only one treatment $W_i$, so we observe only the realized outcome $Y_i =Y_i(W_i)$.
Let $X_i$ be the running or forcing variable, and $c$ be the cutoff for $X_i$ at which the probability of treatment changes. Let $W_i^*$ be a dummy that indicates whether $X_i$ is above or below the threshold $c$, i.e., $W_i^*=\boldsymbol{1}(X_i \geq c)$. Define the potential treatment status $W_i(w^*)$ as what an individual's treatment status would be if $W_i^*=w^*$. We have $W_i=W_i(1)W_i^*+W_i(0)(1-W_i^*)$. If $W_i^*=W_i$, we call it the sharp RDD, where the probability of treatment assignment jumps from 0 to 1 when the forcing variable $X_i$ crosses the cutoff $x=c$. Otherwise, we call it the fuzzy RDD which allows a smaller jump in the probability of treatment assignment at the cutoff and requires $$\lim\limits_{x\rightarrow c^-}p(x)\neq\lim\limits_{x\rightarrow c^+}p(x),$$ where $p(x)=\mathbb{P}(W_i=1 | X_i=x)$. Note that the propensity score $p(\cdot)$ is generally unknown and is estimated using the dataset.
Our primary concern throughout this paper is the identification and inference under the fuzzy RDD. The sharp RDD, as a special case, can be applied with minor adjustments, which are detailed as needed.
In fuzzy RDD settings, two causal quantities are of interest at the cutoff: average treatment effect (ATE) and local average treatment effect (LATE), defined as follows:
Identifying the average treatment effect $\text{ATE}$ typically requires the unconfoundedness assumption. When this assumption does not hold, the LATE can be identified instead if the monotonicity assumption holds. However, the LATE only captures the treatment effect for compliers, defined as individuals having $W_i(0)<W_i(1)$. As a result, the LATE provides a limited perspective on the overall treatment effect, as it is restricted to the subset of compliers.
In this section, we propose a semiparametric framework that identifies the ATE without requiring the unconfoundedness assumption.
Define the conditional expectations of the potential outcomes $Y_i(1)$ and $Y_i(0)$, and the treatment assignment $W_i$, given $X_i$, as follows:
where $\mu_1(x)$ and $\mu_0(x)$ represent the expected outcomes under treatment and control, respectively, and $p(x)$ denotes the propensity score. Define the following deviations:
where $\epsilon_i(0)$ and $\epsilon_i(1)$ represent deviations of the potential outcomes from their conditional expectations, and $\epsilon_i$ denotes the deviation of the treatment assignment from the propensity score. By definition, the expectations of $ \epsilon_i(0), \epsilon_i(1), \text{and} \ \epsilon_i $ are all zero and uncorrelated with $X_i$. It is worth noting that we do not assume $X_i$ to be independent of $(\epsilon_i(0), \epsilon_i(1), \epsilon_i )$. The average treatment effect at the cutoff is defined as $$\tau_{c}=\mu_1(c)-\mu_0(c).$$
We impose the following assumptions:
Assumption (ref) imposes the standard requirement in RDDs that $\mu_1(x)$ and $\mu_0(x)$ are continuous over the entire domain of the running variable. Assumption (ref) reflects a fundamental condition in fuzzy RDDs, requiring the propensity score $p(x)$ to be discontinuous at the cutoff $x=c$ while remaining continuous elsewhere. Additionally, we assume that the propensity score is right-continuous at the cutoff. This is a minimal requirement, as in many applications, the propensity score between $c$ and $c+e$ is nearly identical when $e$ is very small. It is worth noting that Assumption (ref) holds for sharp RDDs by setting $p(x)=\boldsymbol{1}(x\geq c)$.
Assumption (ref), a novel and critical condition proposed in this paper, posits that the covariance between $\epsilon_i(1) - \epsilon_i(0)$ and $\epsilon_i$, conditional on $X_i=x$, exhibits continuity properties. Specifically, let $C(x)=\mathbb{E}[(\epsilon_i(1)-\epsilon_i(0)) \epsilon_i \mid X_i=x]$. Assumption (ref) requires that $C(x)$ is continuous when $x\neq c$ while allowing for a potential discontinuity at $x=c$ due to variations in confounding variables around the cutoff. This condition enhances the flexibility and applicability of our approach by accommodating settings where confounding factors differ near the cutoff. Importantly, it does not require a discontinuity at the cutoff but permits its existence if warranted by the data.
Let us take an example to illustrate this assumption. Consider a confounding variable $Z_i$ correlated with both $W_i$ and $(Y_i(0), Y_i(1))$. Suppose $\epsilon_i=b_wZ_i+\tilde{\epsilon}_i$ and $\epsilon_i(1)-\epsilon_i(0) =b_yZ_i+\breve{\epsilon}_i$, where $b_w$ and $b_y$ are the coefficients, $\breve{\epsilon}_i \perp \!\!\! \perp \tilde{\epsilon}_i \mid X_i$, and $\mathbb{E}[Z_i \mid X_i]=0$. In this case, $C(x)=b_wb_y\text{Var}[Z_i \mid X_i=x]$. Assumption (ref) is satisfied when $\text{Var}[Z_i \mid X_i=x]$ is continuous when $x \neq c$ and right-continuous at $x=c$. This example illustrates that Assumption (ref) is quite mild, as it only requires the continuity of specific moments of the confounding variable conditional on $X_i=x$ for $x\neq c$, while permitting its discontinuity at the cutoff $x=c$.
We now present a semiparametric framework for identifying the average treatment effect at the cutoff in RDDs.
Theorem (ref) demonstrates that the average treatment effect at the cutoff $x=c$ can be expressed as the coefficient of the propensity score $p(X_i)$ in the linear component of the partial linear model (ref).
However, the identification of $\tau_c$ in the partially linear model (ref) is infeasible without imposing any restriction on the nonparametric function $f(x)$ Robinson:88. It is important to note that the identification condition in Robinson:88 does not hold in this model. To ensure the identification, we impose a structural restriction on $f(x)$. Specifically, we approximate $f(x)$ using a basis representation. This restriction facilitates the identification of $\tau_c$.
Corollary (ref) further demonstrates that the average treatment effect at the cutoff $x=c$ in RDDs is identified by the coefficient of the propensity score $p(X_i)$ in the linear component. The restriction corresponds to approximating $f(\cdot)$ with splines. In the estimation step, we employ $P$-spline approach, which achieves a high level of accuracy, with the bias introduced by the restriction being negligible Ruppert:02,Ruppert:03,CKO:09. We further discuss this in Section (ref).
We further generalize in the following corollary that the coefficient of a specified polynomial function of the propensity score in the linear component captures the average treatment effect at the cutoff in fuzzy RDDs.
Corollary (ref) generalizes $p(X_i)$ to $g(p(X_i))$, where $g(\cdot)$ is a polynomial function. This generalization can enhance estimation by selecting an optimal $g(\cdot)$, as discussed in Section (ref). The restriction $\mathbf{a}^\top \mathbf{G} \mathbf{a} =1$ serves as a normalization to facilitate the determination of $\mathbf{a}$. The restriction $a_1>0$ ensures consistency with the case that $g(\cdot)$ degenerates to the identical function when $m=1$.
In this section, we first propose a two-stage estimator for the treatment effect at the cutoff within the semiparametric identification framework. We then establish its asymptotic normality for conducting inference and determine an optimal form of $g(\cdot)$.
In the first-stage estimation, we estimate $p(X_i) = \mathbb{E}[W_i \mid X_i]$, which is flexible and allows for a jump at the cutoff. To achieve this, we apply a nonparametric logistic regression HTF:11 with two segments, allowing for a jump at the cutoff. Specifically, we approximate $p(X_i)$ by a linear combination of a set of basis functions within the following model:
where $\text{logit}(x)= \log\frac{x}{1-x}$. Here, for $j=0,1$, the vector of basis functions evaluated at $X_i$ is $$\mathbf{S}_{j}(X_i) = \left(1, S_{j,1}(X_i), \cdots, S_{j,K_j}(X_i) \right)^\top, $$ where $S_{j,k}(X_i)$ denotes the $k$-th basis function for $k=1, \cdots, K_j$, and $K_j$ represents the number of basis functions. The vectors $\boldsymbol{\alpha}_{0}$ and $\boldsymbol{\alpha}_{1}$ are of dimensions $(K_0+1)$ and $(K_1+1)$, respectively.
We estimate $\boldsymbol{\alpha}_{0}$ and $\boldsymbol{\alpha}_{1}$ by fitting the logistic regression model, obtaining the propensity score estimate $\hat{p}(X_i)$ of $p(X_i)$ for $i=1,\cdots,n$. Once the propensity scores are estimated, the corresponding values of $g(\hat{p}(X_i))$ for $i=1,\cdots,n$ are subsequently computed.
We employ the P-spline approach to estimate $f(x)$. Specifically, $f(x)$ is approximated using splines with a set of radial basis functions in the following form:
where $q$ denotes the polynomial and $\kappa_1 < \kappa_2 < \cdots < \kappa_K$ represent a set of knots along the running variable. The radial basis functions are the default splines in the R package SemiPar SemiPar:05, as described in FKW:01 and Ruppert:03. Specifically, the terms $1, x, \cdots, x^q$ form the polynomial, capturing global trends, while $|x - \kappa_1|^{2q+1}, \cdots, |x - \kappa_K|^{2q+1}$ are radial basis functions, representing local deviations at the specified knots $\kappa_1, \cdots, \kappa_K$. As the number of knots increases, the approximation is of a high level of accuracy. However, it may lead to severe overfitting. To mitigate this issue, we employ the $P$-spline approach to estimate $f(x)$ by imposing a penalty on the parameters $\gamma_k$, $k= 1, \cdots, K$, assuming that they are identically and independently distributed according to a normal distribution $\mathcal{N}(0, \sigma_\gamma^2)$.
The choice of $K$ is a critical consideration in nonparametric estimation. Ruppert:02 provides an algorithm for selecting the optimal number of truncated polynomial basis functions by minimizing the generalized cross-validation. A simple alternative is to select $K$ as $\max \{n/4, 20\}$, ensuring that there are four or five points in each sub-interval consecutive knots Ruppert:02,Ruppert:03.
Denoting $\mathbf{z} = ( |x- \kappa_1|^{2q+1}, \cdots, |x- \kappa_K|^{2q+1} )^\top$ and $\boldsymbol{\gamma} = (\gamma_1, \cdots, \gamma_K)^\top$, Equation (ref) is reformulated as the following form:
where $\boldsymbol{\gamma} \sim \mathcal{N}(\boldsymbol{0}_K, \sigma_\gamma^2 \mathbf{I}_K)$. Let $\boldsymbol{\beta}=( \beta, \beta_0, \cdots, \beta_q)^\top$, $\mathbf{g} = ( g(\hat{p}(X_i)), \cdots, g(\hat{p}(X_n)))^\top$, $\mathbf{x}=( \boldsymbol{1}(x \geq c), 1, x, \cdots, x^q)^{\top}$, $\mathbf{y} = (y_1, \cdots, y_n)^\top$, $\mathbf{X} = (\mathbf{x}_1, \cdots, \mathbf{x}_n)^\top$, $\mathbf{Z} = (\mathbf{z}_1, \cdots, \mathbf{z}_n)^\top$ and $\boldsymbol{\varepsilon} = ({\varepsilon}_1, \cdots, {\varepsilon}_n)^\top$. From Equation (ref), we rewrite the model (ref) in Corollary (ref) as the following linear mixed-effects model:
where $(\boldsymbol{\gamma}^\top, \boldsymbol{\varepsilon}^\top )^\top \sim \mathcal{N}(\boldsymbol{0}_{K+n}, \text{diag}(\sigma_\gamma^2 \mathbf{1}_K, \sigma_1^2, \cdots, \sigma_n^2 ))$. The model (ref) disregards the estimation error of $\hat{p}(X_i)$, which is negligible due to the consistency of nonparametric logistic regression GS:93.
From the model, the covariance matrix of $\mathbf{y}$ is given by $\mathbf{V}_0 = \sigma_\gamma^2 \mathbf{Z}\mathbf{Z}^\top + \sigma^2 \text{diag}(\lambda_1^2, \cdots, \lambda_n^2)$, where $\lambda_i^2 = \sigma_i^2/\sigma^2$ for $i=1, \cdots, n$. However, to simplify the estimation step, we disregard the heteroscedasticity and instead use the simplified covariance term $\mathbf{V} = \sigma_\gamma^2 \mathbf{Z}\mathbf{Z}^\top + \sigma^2 \mathbf{I}_n$. Thus, we estimate $\boldsymbol{\theta} = (\tau_c, \boldsymbol{\beta}^\top)^\top$ by applying generalized least squares, yielding the following estimate:
where $\mathbf{U} = (\mathbf{g}, \mathbf{X})$. It is clear that the estimate $\tilde{\boldsymbol{\theta}}$ remains consistent and unbiased. We obtain the estimates of $\sigma_\gamma^2, \sigma^2$ by applying the maximum likelihood or restricted maximum likelihood method corbeil1976restricted. Then, we calculate the estimated covariance matrix as $\hat{\mathbf{V}} = \hat{\sigma}_\gamma^2 \mathbf{Z}\mathbf{Z}^\top + \hat{\sigma}^2 \mathbf{I}_n $. Substituting these estimates to Equation (ref), the estimate of $\boldsymbol{\theta}$ is obtained as
Finally, the two-stage estimator of $\tau_c$ is given by
where $\mathbf{e}_1$ is a $(q+3)$-dimensional vector whose first element is $1$ and others are $0$.
We present the asymptotic normality of the two-stage estimator $\hat{\tau}_c$ within the model (ref) for inference. Before introducing the property, we first impose the following assumptions.
Part (a) in Assumption (ref) requires that $U$ has full column rank, which is a mild condition as $q$ is typically small, often taking values in $\{1,2,3\}$. Parts (b) and (c) impose constraints on the number of knots $K$, requiring that $K$ approaches infinity as $n\rightarrow \infty$, and that $K$ is of a smaller or the same order as $n$. Part (d) allows the rank of the design matrix of $\mathbf{Z}$ to be smaller but of the same order as $n$, and of the same order as $K$.
Theorem (ref) is proven in the case where the unknown variance components $\sigma_{\gamma}^2$ and $\sigma^2$ are estimated by maximum likelihood. If the restricted maximum likelihood method is employed instead, the above asymptotic normality still holds, as demonstrated by Jiang:98. The condition that $\varepsilon_1,\cdots,\varepsilon_n$ are Gaussian is necessary to leverage the asymptotic results from Miller:77. The condition that $\lambda_1,\cdots,\lambda_n$ are known simplifies the analysis and facilitates the application of the results in Miller:77.
We employ heteroskedasticity-consistent (HC) standard errors to account for heteroskedasticity, as suggested by white1980heteroskedasticity and huang2022accounting. Specifically, the variance $\mathbf{V}_0$ is estimated as $\hat{\mathbf{V}}_0 = \text{diag}( \hat{v}^2_1,\cdots, \hat{v}^2_n)$, where $\hat{v}_i = \frac{e_i}{1-h_i}$, with $e_i$ denoting the marginal residual and $h_i$ the leverage value for unit $i$.
By substituting $\hat{\mathbf{V}}$ and $\hat{\mathbf{V}}_0$ into the definitions of $\mathbf{R}$ and $\mathbf{S}$ as stated in Corollary (ref), we obtain the estimates $\hat{\mathbf{R}}$ and $\hat{\mathbf{S}}$, respectively. Consequently, from Corollary (ref), an $1-\alpha$ confidence interval for $\tau_{c}$ is given by
where $\hat{V}_{\tau} = {\mathbf{g}^{\top} \hat{\mathbf{R}}\mathbf{g} }/{ (\mathbf{g}^{\top}\hat{\mathbf{S}}\mathbf{g})^2 }$ and $z_{\alpha/2}$ is the upper $\alpha/2$-quantile from a standard Gaussian distribution.
Intuitively, our global approach utilizes all data points to estimate the treatment effect, whereas the local approach replies on data near the cutoff Hahn:01. As a result, our global approach is expected to be more efficient than the local approach. We theoretically compare our global approach with the local approach and demonstrates that our approach exhibits lower bias and variance, achieving a faster convergence rate.
Theorem (ref) has demonstrated that the asymptotic variance of $\hat{\tau}_c$ is $O(n^{-1/2})$, while the asymptotic bias term is negligible, i.e., $o(n^{-1/2})$ as discussed in Section (ref). For local linear regression, FG:96 shows that the asymptotic bias is of order $h^2$, where $h$ is the bandwidth, and the asymptotic variance is of order $1/(nh)$. Hahn:01 discusses that these properties specifically apply to the local estimator in the context of RDDs. Typically, the optimal bandwidth is chosen as $h=O(n^{-1/5})$. Therefore, we have the following proposition.
Proposition (ref) demonstrates that our global approach is more effective than the local approach for RDDs.
Until now, we have investigated the estimation and inference procedure of ${\tau}_c$ given a specified form of $g(\cdot)$. We now turn to the method for determining the optimal form of $g(\cdot)$ that minimizes the variance $V_{\tau}$.
Recall that $g(t) = a_1t+\cdots+a_m t^m$. Let $\mathbf{p}_k = (\hat{p}^k(X_1), \cdots, \hat{p}^k(X_n) )^\top$ for $k =1, \cdots, m$, and denote $\mathbf{P} = (\mathbf{p}_1, \cdots, \mathbf{p}_m)$. We then express $\mathbf{g}=(g(\hat{p}(X_1)),\cdots,g(\hat{p}(X_n)))^{\top}$ as $\mathbf{g} = \mathbf{P}\mathbf{a}$. The form of $g(\cdot)$ is determined by finding a coefficient vector $\mathbf{a}$ for the polynomial function $g(\cdot)$ that minimizes the variance $V_{\tau}$.
Define $\mathbf{Q}_{\text{R}} = m \mathbf{P}^\top \mathbf{R} \mathbf{P} / \text{tr}(\mathbf{P}^\top \mathbf{R} \mathbf{P})$ and $\mathbf{Q}_{\text{S}} = m \mathbf{P}^\top \mathbf{S} \mathbf{P} / \text{tr}(\mathbf{P}^\top \mathbf{S} \mathbf{P})$. We specify $\mathbf{G}$ in Corollary (ref) as $\mathbf{Q}_{\text{S}}$, which implies that the optimal coefficient vector is determined under the conditions $\mathbf{a}^\top \mathbf{Q}_{\text{S}} \mathbf{a} = 1$ and $a_1 \geq 0$. Therefore, based on Equation (ref), the variance minimization problem is equivalent to solve the following optimization problem:
Note that $\mathbf{Q}_{\text{S}}$ is positive semi-definite and might be singular. To address this issue, we first perform an orthogonal decomposition of it:
where $\mathbf{\Lambda}$ is a diagonal matrix containing all positive eigenvalues, and $\mathbf{S}_1$ is the matrix consisting of the corresponding eigenvectors. In practice, eigenvalues below $10^{-5}$ are treated as zero. Let $\mathbf{a} = \mathbf{S}_1 \mathbf{b}$ and then we just need to solve the following problem:
Since $\mathbf{S}^\top_1 \mathbf{Q}_{\text{S}} \mathbf{S}_1$ is positive definite, Problem (ref) is a classical optimization problem for minimizing a generalized Rayleigh quotient, with the optimal solution given by $\mathbf{b}_{opt}$. The optimal coefficient vector $\mathbf{a}_{opt} = \mathbf{S}_1 \mathbf{b}_{opt}$ is uniquely determined. Consequently, the optimal function $\mathbf{g}$ is $\mathbf{g}_{opt} = \mathbf{P} \mathbf{a}_{opt}$.
We conclude this section by discussing the choice of $m$. The form of $g(\cdot)$ that we use can be viewed as an approximation using a set of polynomial functions derived from $p(X_i)$. When $m$ is large, the search space for $g(\cdot)$ expands, which, however, also increases the variance. This is analogous to nonparametric fitting with a polynomial basis, where selecting a high-degree polynomial can lead to excessive variance. Thus, a practice approach is to select $m$ from a moderate range, such as $\{3,4,5,6,7\}$. In our calculations, we set $m=5$ and also evaluate performance across various $m$ values, observing that the superiority of our approach is not sensitive to this choice. To optimally balance bias and variance, cross-validation can be employed.
In this section, we conduct Monte Carlo experiments to assess the performance of our method. We consider two scenarios: one where the unconfoundedness assumption holds and the other where it does not. Note that in this section we focus on the performance for fuzzy RDDs. However, we also include a simulation study for sharp RDDs and report the results in Table (ref) of the Appendix (ref), where our method demonstrates superior performance.
We set the cutoff at 0 and generate the running variable $X_i$ from a uniform distribution, i.e. $X_i \sim \mathcal{U}(-1, 1)$. We consider three regression functions for the potential outcomes, labeled M1, M2, and M3, respectively:
Each of the three models above represents a distinct pattern of change. In M1, the average treatment effect increases monotonically throughout the running variable. In M2, the treatment effect remains 0 when the running variable is negative and increases for $x \geq 0$. M3 corresponds to the function in Lee:08.
Based on the setup outlined above, we now consider the two scenarios to generate the observed treatments $W_i$ and outcomes $Y_i$.
Scenario 1: UA holds. We first consider the scenario where the unconfoundedness assumption (UA) holds. The treatment variable is given as follows:
where $\text{expit}(x) = 1/(1+e^{-x})$. The observed outcome for unit $i$ is generated as follows:
where $\eta_i \sim \mathcal{N}(0, \sigma^2_{\eta})$. We set $\sigma^2_{\eta}$ to ensure that $R^2$ in Equation (ref) is maintained as 0.75.
Scenario 2: UA is violated. We introduce an extra exogenous random error $\varepsilon_i$, which is defined as $\varepsilon_i \sim \boldsymbol{1}(X_i < 0) \mathcal{N}(0, \sigma^2_\varepsilon) + \boldsymbol{1}(X_i \geq 0) \mathcal{N}(0, 2\sigma^2_\varepsilon)$, and add it to the propensity score model, as follows:
We set the variance $\sigma^2_\varepsilon = \text{Var}\left(0.5X_i+0.2X_i^2 + 2\boldsymbol{1}(X_i \geq 0) -1 \right)/3$. We then generate $W_i$ according to Equation (ref) as well. $Y_i(0)$ and $Y_i(1)$ are generated as follows:
where $\varepsilon_i(0) = c_0 \varepsilon_i$ and $\varepsilon_i(1) = c_1 \varepsilon_i$. Here, $c_0$ and $c_1$ are constants, determined by setting $R^2$ to 0.75 in Equations (ref) and (ref), respectively. The observed outcome is then generated by $Y_i = Y_i(0) + (Y_i(1) - Y_i(0))W_i$. Clearly, in this case, the unconfoundedness assumption is violated, while our Assumption 3 holds.
In the first-stage regression, we use natural splines for $\mathbf{S}_0(\cdot)$ and $\mathbf{S}_1(\cdot)$ in Equation (ref). Following the recommendation of harrell2001regression, the number of knots is set to either 3 or 5, with the final selection based on the largest $R^2$ zhang2017coefficient.
In the second-stage regression, we use the R package SemiPar SemiPar:05 to implement our global estimation. The degree $q$ is set to 1 in Equation (ref), which results in basis functions corresponding to cubic thin plate splines. The number and placement of knots are set by default in the SemiPar package. The parameter $m$ is chosen to be $5$, two degrees higher than the order of $f$ and natural splines, consistent with the configuration used in Lee:08. We evaluate the performance of our method for different values of $m$ and report the results in Table (ref) of the Appendix (ref). The results demonstrate that our method maintains its superiority across varying values of $m$. In the implementation of local estimators, the bandwidth is selected using two methods: those proposed by IK:12 and CCT:14. A triangular kernel is employed in both cases.
Simulations are conducted for two sample sizes, 500 and 1000, with 10,000 Monte Carlo replications for each case. The results for the case where the unconfoundedness assumption holds are presented in Table (ref), while the results for the case where the unconfoundedness assumption does not hold are shown in Table (ref). We compare our method, denoted as “PL" for the partially linear model, with two local estimators from IK:12 and CCT:14, referred to as “IK" and “CCT", respectively.
From both tables, we make the following observations. First, our method consistently achieves significantly smaller RMSE compared with the local methods across all cases in both scenarios. In terms of bias, our method shows much smaller values in M3 and comparable values in M1 and M2. Second, the empirical coverage of our method is close to nominal coverage, with generally smaller confidence interval lengths. Third, when comparing the two sample sizes, $n=500$ and $n=1000$, the advantage of our method over the local methods is more pronounced when the sample size is smaller. Finally, comparing the two scenarios, our method shows a larger advantage when the unconfoundedness assumption is violated than when it holds. In the scenario where the unconfoundedness assumption does not hold, both the bias and RMSE values of our method are consistently much smaller than those of the local methods across all cases.
\end{threeparttable} \end{table} \fi
In this section, we apply our global approach to study two problems: evaluating the effect of antihypertensive treatment on reducing the risk of cardiovascular disease calonico2024regression, and assessing the effect of incumbent advantage in U.S. House and Senate elections Lee:08,CFT:15.
calonico2024regression provided an insightful article in which they used the teaching version of the Framingham Heart Study dataset to explore whether antihypertensive medication could reduce the risk of cardiovascular disease, treating systolic blood pressure as the running variable. The Framingham Heart Study is a long-term prospective study on the etiology of cardiovascular disease among a population of free-living individuals in Framingham, Massachusetts tsao2015cohort. The teaching version of the dataset is available upon request from the Biologic Specimen and Data Repository Information Coordinating Center.
We follow the steps outlined in calonico2024regression and continue to use the teaching version of the dataset to explore the treatment effect of antihypertensive medication on the occurrence of cardiovascular disease. However, we exclude all observations from examination cycles 1 and 2 to ensure the independence of the samples. Since the diagnostic criteria for hypertension are a systolic blood pressure greater than 140 mm Hg or a diastolic blood pressure greater than 90 mm Hg, we select these as the running variables, respectively. The treatment variable indicates whether the antihypertensive medication was used during the exam.
In our study, Cardiovascular Disease is defined as Angina Pectoris, Myocardial infarction (hospitalized and silent or unrecognized), coronary insufficiency (unstable Angina), or Fatal Coronary Heart Disease. Thus, the observed outcome is a binary variable which equals 1 if the condition is present during follow-up and 0 otherwise. The covariates include sex, age, the number of cigarettes smoked per day, BMI, heart rate (ventricular rate) in beats/min, and the presence of diabetes and prevalent coronary heart disease on exam. Diastolic blood pressure will serve as a covariate if systolic blood pressure is the running variable and vice versa. After excluding all missing values, the final study dataset includes 2,799 units. It is worth noting, however, that any findings should not be used to interpret actual results, as the data have been redacted for teaching purposes.
Table (ref) presents the results, where “Systolic Blood Pressure" indicates that systolic blood pressure is used as the running variable, and “Diastolic Blood Pressure" conveys a similar meaning. Similar to the results from the local methods calonico2024regression, the results from our approach are not statistically significant. However, compared with the local methods, our approach yields treatment effect estimates with much lower standard errors, suggesting that it may be more reliable.
Moreover, we restrict the scope of Cardiovascular Disease to Hospitalized Myocardial Infarction, and re-evaluate the effect of antihypertensive treatment on reducing the risk of cardiovascular disease. The results are reported in Table (ref). From the table, the results from our method are similar to those in Table (ref), with no significant change and the same sign. However, there is a substantial difference in magnitude and the direction reverses for the two local estimates. Given the assumption that defining Cardiovascular Disease as Hospitalized Myocardial Infarction largely overlaps with its definition based on multiple measures, this study demonstrates that our global approach is more reliable and consistent.
We demonstrate the performance of our method in sharp RDDs by analyzing the effects of incumbent advantage in U.S. House and Senate elections by analyzing two datasets. The first dateset, consisting of 6558 samples, is sourced from U.S. House elections Lee:08, while the other dataset, comprising 1390 samples, examines party-level advantage in U.S. Senate elections spanning from 1914 to 2010 CFT:15. In these datasets, the forcing variable $x_i$ represents the margin of victory of the Democratic party in a given election, with the outcome variable being the Democratic vote share in the subsequent election.
Our identification result can directly be applied sharp RDDs, and the estimation method is also applicable to sharp RDDs. Unlike in fuzzy RDDs, there is no need to perform a first-stage regression or account for $g(\cdot)$. We utilize cubic thin plate splines for global estimation and compare our approach with the local estimators. The results are reported in Table (ref). For both datasets, the PL method performs comparably to the local estimator, suggesting that our method could serve as a viable alternative option for sharp RDDs in practice.
We lack knowledge of the true treatment effects in the real datasets. To further assess the performance of our method relative to local estimators, we conduct a simulated experiment by setting artificial cutoffs at $c=\pm 0.1$, retaining the original control and treated units, and generating corresponding virtual datasets. Using our global estimator and the local estimators on these virtual datasets, we infer the treatment effects at $c=\pm 0.1$ and summarize the results in Table (ref).
Under the conditions of the sharp RDD, the true estimands at the cutoffs at $c=\pm 0.1$ are $0$. From Table (ref), our method fails to reject the null hypothesis in all situations, while IK significantly rejects it when $c=0.1$ in the U.S. House election dataset. Additionally, CCT only fails to reject it when $c=-0.1$ in the U.S. Senate election dataset. Therefore, our method offers more reliable inferences in this context. This investigation underscores the potential advantage of our method in providing more reliable inferences compared to the local estimators for RDDs.
The estimation of the average treatment effect at the cutoff is generally considered infeasible in RDDs without the unconfoundedness assumption; instead, only the local treatment effect is identified. In this paper, we have propose a semiparametric inference framework that identifies the average treatment without relying on the unconfoundedness assumption.
Moreover, this framework offers an efficient global estimator for the average treatment effect without replying on the unconfoundedness assumption. Compared to local estimators commonly used for RDDs, we have demonstrated the superiority of our global estimator through theoretical analysis and simulation experiments. We have also applied our proposed method to three datasets, yielding reliable inference. We believe that our proposed global estimator is a viable alternative to local estimators for RDDs.
There are several potential extensions to consider. This paper focuses on the use of $P$-splines to represent the nonparametric part, but other nonparametric methods could also be efficient. Further investigation in this area is warranted. Another potential avenue is the regression kink design Card:15, where CCT:14 and GI:19 recommend using local quadratic approach, building on the argument presented by Hahn:01. Similar to our study on RDDs, our semiparametric inference framework could offer a robust alternative to local estimators in the regression kink design. Additionally, extending our framework to RDDs with covariates C:16 and RDDs with a discrete running variable KR:18 represents worthwhile directions for future research.