A Nonparametric Bayesian Methodology for Regression Discontinuity Designs

TL;DR

This paper introduces a Bayesian Gaussian process regression approach for regression discontinuity designs, offering a flexible and uncertainty-aware alternative to local linear regression, with proven consistency and promising empirical performance.

Contribution

It develops a Gaussian process-based Bayesian methodology for regression discontinuity, improving flexibility and uncertainty quantification over traditional local linear regression.

Findings

Method is consistent in estimating treatment effects.

Simulation shows better coverage and mean squared error.

Applied to NBA draft data with meaningful insights.

Abstract

One of the most popular methodologies for estimating the average treatment effect at the threshold in a regression discontinuity design is local linear regression (LLR), which places larger weight on units closer to the threshold. We propose a Gaussian process regression methodology that acts as a Bayesian analog to LLR for regression discontinuity designs. Our methodology provides a flexible fit for treatment and control responses by placing a general prior on the mean response functions. Furthermore, unlike LLR, our methodology can incorporate uncertainty in how units are weighted when estimating the treatment effect. We prove our method is consistent in estimating the average treatment effect at the threshold. Furthermore, we find via simulation that our method exhibits promising coverage, interval length, and mean squared error properties compared to standard LLR and state-of-the-art LLR methodologies. Finally, we explore the performance of our method on a real-world example by studying the impact of being a first-round draft pick on the performance and playing time of basketball players in the National Basketball Association.

Tables

Table 1. Table 1 : Treatment effect estimation for LLR, robust LLR, and GPR on NBA data

	LLR		Robust LLR		GPR
Outcome	Estimate	95% CI	Estimate	95% CI	Estimate	95% CI
Box Plus-Minus	-3.06	[-5.20, -0.92]	-3.37	[-5.86, -0.88]	-2.42	[-5.02, 0.01]
Win Shares	-1.58	[-4.09, 0.93]	-1.73	[-4.76, 1.29]	-1.08	[-3.04, 0.88]
Minutes Played	-446.29	[-1042.81, 150.23]	-406.63	[-1123.99, 310.72]	-640.75	[-1259.65, 5.51]
Games Played	-29.71	[-40.20, -19.22]	-28.16	[-40.81, -15.51]	-32.00	[-45.69, -18.14]

Table 2. Table 2 : Simulation for n = 500 𝑛 500 n=500 , shown in Figure 2

	Method
	LLR				Robust LLR				GPR
Dataset	EC	$\overline{\text{IL}}$	Bias	RMSE	EC	$\overline{\text{IL}}$	Bias	RMSE	EC	$\overline{\text{IL}}$	Bias	RMSE
Lee	89.2%	0.21	0.02	0.06	91.7%	0.25	0.01	0.07	82.3%	0.19	0.05	0.07
Quad	92.8%	0.21	0.00	0.06	93.9%	0.25	0.00	0.07	97.0%	0.18	0.00	0.04
Cubic	92.5%	0.22	-0.01	0.06	93.0%	0.25	0.00	0.07	96.3%	0.20	-0.01	0.05
Cate 1	92.4%	0.24	-0.01	0.07	92.4%	0.28	0.00	0.08	94.4%	0.25	-0.01	0.07
Cate 2	92.4%	0.24	-0.01	0.07	92.4%	0.28	0.00	0.08	94.4%	0.25	-0.01	0.07
Ludwig	85.1%	0.32	0.05	0.10	93.1%	0.35	0.01	0.09	75.0%	0.25	0.08	0.10
Curvature	91.2%	0.22	-0.01	0.06	92.7%	0.25	0.00	0.07	96.3%	0.21	-0.01	0.05

Table 3. Table 3 : Simulation Results for GPR, GPR (Cut), GPR (Zero Mean), and GPR (MLE)

	Method
	GPR				GPR (Cut)				GPR (Zero Mean)				GPR (MLE)
Data	EC	$\overline{\text{IL}}$	Bias	RMSE	EC	$\overline{\text{IL}}$	Bias	RMSE	EC	$\overline{\text{IL}}$	Bias	RMSE	EC	$\overline{\text{IL}}$	Bias	RMSE
Lee	82.3%	0.19	0.05	0.07	94.3%	0.22	0.03	0.06	81.6%	0.18	0.05	0.07	76.4%	0.15	0.04	0.06
Quad	97.0%	0.18	0.00	0.04	97.2%	0.23	-0.00	0.05	97.0%	0.19	0.01	0.04	95.9%	0.18	0.00	0.05
Cubic	96.3%	0.20	-0.01	0.05	96.0%	0.35	-0.01	0.08	96.7%	0.20	0.00	0.05	91.1%	0.19	-0.03	0.06
Cate 1	94.4%	0.25	-0.01	0.07	96.9%	0.29	-0.01	0.07	95.2%	0.26	0.00	0.07	93.8%	0.25	0.00	0.07
Cate 2	94.4%	0.25	-0.01	0.07	97.2%	0.29	-0.01	0.07	95.4%	0.26	0.00	0.07	94.2%	0.25	0.00	0.07
Ludwig	75.0%	0.25	0.08	0.10	88.9%	0.34	0.06	0.10	64.5%	0.24	0.10	0.12	52.0%	0.24	0.11	0.13
Curvature	96.3%	0.21	-0.01	0.05	96.1%	0.26	-0.01	0.06	96.9%	0.20	0.00	0.05	89.6%	0.20	-0.03	0.06

Table 4. Table 4 : Simulation for n = 500 𝑛 500 n=500 , using the IK bandwidth

	Method
	LLR				Robust LLR				GPR
Dataset	EC	$\overline{\text{IL}}$	Bias	RMSE	EC	$\overline{\text{IL}}$	Bias	RMSE	EC	$\overline{\text{IL}}$	Bias	RMSE
Lee	82.7%	0.15	0.04	0.05	91.2%	0.27	0.01	0.08	82.3%	0.19	0.05	0.07
Quad	94.6%	0.15	0.00	0.04	93.2%	0.24	0.00	0.07	97.0%	0.18	0.00	0.04
Cubic	93.7%	0.19	-0.01	0.05	93.8%	0.24	0.01	0.06	96.3%	0.20	-0.01	0.05
Cate 1	91.2%	0.21	-0.01	0.06	92.9%	0.26	0.01	0.07	94.4%	0.25	-0.01	0.07
Cate 2	92.1%	0.22	-0.01	0.06	92.8%	0.26	0.01	0.07	94.4%	0.25	-0.01	0.07
Ludwig	31.5%	0.23	0.15	0.16	89.8%	0.27	0.04	0.08	75.0%	0.25	0.08	0.10
Curvature	84.4%	0.19	-0.03	0.06	94.8%	0.23	0.00	0.06	96.3%	0.21	-0.01	0.05

Table 5. Table 5 : Simulation for n = 500 𝑛 500 n=500 , using the CER bandwidth

	Method
	LLR				Robust LLR				GPR
Dataset	EC	$\overline{\text{IL}}$	Bias	RMSE	EC	$\overline{\text{IL}}$	Bias	RMSE	EC	$\overline{\text{IL}}$	Bias	RMSE
Lee	91.1%	0.25	0.01	0.07	91.8%	0.27	0.01	0.08	82.3%	0.19	0.05	0.07
Quad	92.3%	0.24	-0.00	0.07	92.7%	0.27	-0.00	0.08	97.0%	0.18	0.00	0.04
Cubic	92.2%	0.25	-0.00	0.07	92.1%	0.27	0.00	0.08	96.3%	0.20	-0.01	0.05
Cate 1	91.9%	0.28	-0.00	0.08	92.3%	0.30	0.00	0.08	94.4%	0.25	-0.01	0.07
Cate 2	92.0%	0.28	-0.00	0.08	92.3%	0.30	0.00	0.08	94.4%	0.25	-0.01	0.07
Ludwig	91.5%	0.39	0.02	0.10	92.9%	0.41	0.00	0.10	75.0%	0.25	0.08	0.10
Curvature	91.4%	0.26	-0.00	0.07	93.3%	0.28	0.00	0.08	96.3%	0.21	-0.01	0.05

Equations

$$y_i=W_i{Y}_i(1)+(1-W_i)Y_i(0)$$

$$0<P(W_i=1|X_i=x)<1$$

$$W_i=\{\begin{array}{ll}1& \text{mbox}if{X}_i\ge b\\ 0& \text{mbox}otherwise.\end{array}$$

$${\mu }_T(x)\equiv \mathbb{E}[Y_i(1)|X_i=x],{\mu }_C(x)\equiv \mathbb{E}[Y_i(0)|X_i=x],$$

$$\tau =\underset{x\downarrow b}{\lim }\mathbb{E}[y_i|X_i=x]-\underset{x\uparrow b}{\lim }\mathbb{E}[y_i|X_i=x]={\mu }_T(b)-{\mu }_C(b).$$

$${\hat{\mu }}_T(b)=X_b(X_T^T{W}_T{X}_T{)}^{-1}{X}_T^T{W}_T{Y}_T$$

$$(W_T{)}_{ii}\equiv K\left(\frac{x_i-b}{h}\right),i=1,\dots ,n_T$$

$$\mathbf{y}=f(\mathbf{x})+ϵ$$

$$f(\mathbf{x})\sim \text{GP}(m(\mathbf{x}),K(\mathbf{x},{\mathbf{x}}^{\prime }))$$

$$K(x_i,x_j)\equiv {\sigma }_{GP}^2\exp \left(-\frac{1}{2{\ell }^2}(x_i-x_j{)}^2\right).$$

$$f(x^{\ast })|\mathbf{x},\mathbf{y}$$

$${\mu }^{\ast }$$

$${\Sigma }^{\ast }$$

$$Y_i(1)={\mu }_T(x_i)+{ϵ}_{i1},\text{and}{Y}_i(0)={\mu }_C(x_i)+{ϵ}_{i0},\text{where}$$

$$\displaystyle\mu_{T}(x_{i})\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\mu_{C}(x_{i}),\hskip 7.22743pt\epsilon_{i1}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sigma_{y1}^{2}),\hskip 7.22743pt\text{and }\epsilon_{i0}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sigma_{y0}^{2})$$

$${\mu }_T(x)$$

$${\mu }_C(x)$$

$${\mu }_T(b)|\mathbf{x},\mathbf{y}$$

$${\mu }_C(b)|\mathbf{x},\mathbf{y}$$

$${\mu }_{b|T}$$

$${\mu }_{b|C}$$

$${\Sigma }_{b|T}$$

$${\Sigma }_{b|C}$$

$$\tau \equiv {\mu }_T(b)-{\mu }_C(b)|\mathbf{x},\mathbf{y}\sim N\left({\mu }_{b|T}-{\mu }_{b|C},2{\sigma }_{GP}^2-{\Sigma }_{b|T}-{\Sigma }_{b|C}\right).$$

$${\mu }_{b|T}$$

$${\mu }_{b|C}$$

$${\Sigma }_{b|T}$$

$${\Sigma }_{b|C}$$

$$m_T(\mathbf{x})$$

$$m_C(\mathbf{x})$$

$$\prod \left(\tau :h(\tau ,{\tau }^{\ast })\ge M{ϵ}_n|x_1,\dots ,x_n\right)\stackrel{P_{{\tau }^{\ast }}}{\to }0$$

$$m_T(\mathbf{x})$$

$$m_C(\mathbf{x})$$

$$y_i$$

$$C_1{a}^p\exp (-D_{1}a{\log }^{q}a)\le \kappa (a)\le C_2{a}^p\exp (-D_{2}a{\log }^{q}a)$$

$$\prod \left({\mu }_T:h({\mu }_T,{\mu }_T^{\ast })\ge M_T{ϵ}_{n_T}|x_1,\dots ,x_{n_T}\right)\stackrel{P_{{\mu }_T^{\ast }}}{\to }0,\text{and}$$

Full text

Keywords: Regression discontinuity designs, Gaussian process regression, Bayesian nonparametrics, coverage, posterior consistency

1 Introduction

Recently there has been a renewed interest in regression discontinuity designs (RDDs), which originated with Thistlethwaite & Campbell (1960). In an RDD, the treatment assignment is discontinuous at a certain covariate value, or “threshold,” such that only units whose covariate is above the threshold will receive treatment. There are many examples of RDDs in the applied econometrics literature: the United States providing additional funding to only the 300 poorest counties for the Head Start education program (Ludwig & Miller, 2007); schools mandating students to attend summer school if their exam scores are below a threshold (Matsudaira, 2008); colleges offering financial aid to students whose academic ability is above a cutoff (Van der Klaauw, 2002); and Medicare increasing insurance coverage after age 65 (Card, Dobkin, & Maestas, 2004). The main goal of an RDD is to estimate a treatment effect while addressing likely confounding by the covariate that determines treatment assignment.

One of the most popular methodologies for estimating the average treatment effect at the threshold in an RDD is local linear regression (LLR), which places larger weight on units closer to the threshold. Implementation of LLR is straightforward and there is a wide literature on its theoretical properties. However, recent works have found that LLR can exhibit poor inferential properties—such as confidence intervals that tend to undercover—which has motivated a strand of literature started by Calonico et al. (2014) that modifies LLR to improve coverage and other inferential properties.

Adding to this literature, we propose a nonparametric regression approach that acts as a Bayesian analog to LLR for sharp regression discontinuity designs. Our approach utilizes Gaussian process regression (GPR) to provide a flexible fit for treatment and control responses by placing a general prior on the mean response functions. While GPR has been widely used in the machine learning and statistics literature, it has not previously been proposed for estimating treatment effects in RDDs. Thus, our main contribution is outlining how to use Gaussian processes to make causal inferences in RDDs and assess how such a methodology compares to current LLR methodologies.

In the remainder of this section, we review RDDs and LLR methodologies for estimating the average treatment effect at the threshold. In Section 2, we outline GPR for sharp RDDs and note various analogies to LLR, which builds intuition for implementing our method. In Section 3, we establish that our method is consistent in estimating the average treatment effect at the boundary. In Section 4, we show via simulation that our method exhibits promising coverage, interval length, and mean squared error properties compared to standard LLR and state-of-the-art LLR methodologies. In Section 5, we use GPR on data from the National Basketball Association (NBA) to estimate the effect of being a first-round versus a second-round pick on basketball player performance and playing time, and we find that GPR detects treatment effects that are more in line with previous results in the sports literature than do LLR methodologies. In Section 6, we conclude by discussing extensions to our methodology to tackle problems beyond sharp RDDs.

1.1 Overview of Regression Discontinuity Designs

We follow the notation of Imbens & Lemieux (2008) and discuss RDDs within the potential outcomes framework: For each unit $i$, there are two potential outcomes, $Y_{i}(1)$ and $Y_{i}(0)$, corresponding to treatment and control, respectively, and a covariate $X_{i}$. Only one of these two potential outcomes is observed, but $X_{i}$ is always observed. Let $W_{i}$ denote the assignment for unit $i$, where $W_{i}=1$ if $i$ is assigned to treatment and $W_{i}=0$ if unit $i$ is assigned to control. The observed outcome for unit $i$ is then

[TABLE]

Often, one wants to estimate the average treatment effect $\mathbb{E}[Y_{i}(1)-Y_{i}(0)]$, but usually this treatment effect is confounded by $X$ (and possibly other unobserved covariates), such that examining the difference in mean response between treatment and control is not appropriate. In these cases, methods such as stratification, matching, and regression are often employed to address covariate confounding when estimating the average treatment effect. However, such methods are only appropriate when there is sufficient overlap in the treatment and control covariate distributions, i.e.,

[TABLE]

For example, this overlap assumption is essential for propensity score methodologies, where the relationship between $W_{i}$ and $X_{i}$ is estimated and then accounted for during the analysis.

In an RDD, the relationship between the treatment assignment $W_{i}$ and the covariates is known. More specifically, it is known that the function $P(W=1|X=x)$ is discontinuous at some threshold or boundary $X=b$. In this paper we focus on a special case, sharp RDDs, where the treatment assignment for a unit $i$ is

[TABLE]

The covariate $X$ that determines treatment assignment in an RDD is often called the “running variable” in order to not confuse it with other background covariates that do not necessarily determine treatment assignment. The more general case, where $P(W=1|X=x)$ is discontinuous at $X=b$ but is not necessarily 0 or 1, is known as a fuzzy RDD.

Rubin (1977) states that when treatment assignment depends on one covariate, estimating $\mathbb{E}[Y_{i}(1)|X_{i}]$ and $\mathbb{E}[Y_{i}(0)|X_{i}]$ is essential for estimating the average treatment effect. Furthermore, treatment effect estimates are particularly sensitive to model specification of $\mathbb{E}[Y_{i}(1)|X_{i}]$ and $\mathbb{E}[Y_{i}(0)|X_{i}]$ when there is not substantial covariate overlap, as in a sharp RDD. Aware of this sensitivity, RDD analyses typically do not attempt to estimate the average treatment effect. Instead, they focus on the estimand that requires the least amount of extrapolation to overcome this lack of covariate overlap: the average treatment effect at the boundary $b$. Defining the conditional expectations for treatment and control as

[TABLE]

the treatment effect at the boundary $b$ is (Imbens and Lemieux 2008)

[TABLE]

The notation $\mu_{T}(x)$ and $\mu_{C}(x)$ emphasizes that the goal of an RDD requires estimating two unknown mean functions. Hahn et al. (2001) showed that sufficient conditions for $\tau$ to be identifiable in a sharp RDD are that $\mu_{C}(x)$ and $\mu_{T}(x)-\mu_{C}(x)$ are continuous at $b$. They further state that “we can use any nonparametric estimator to estimate” $\mu_{T}(x)$ and $\mu_{C}(x)$, and recommended local linear regression (hereafter called LLR), which is currently the most popular methodology for estimating the average treatment effect at the boundary in an RDD.

1.2 Review of Local Linear Regression

The goal of an RDD is to estimate $\tau$ defined in (5), i.e., to estimate $\mu_{T}(b)$ and $\mu_{C}(b)$. LLR estimates $\mu_{T}(b)$ as

[TABLE]

where $X_{b}\equiv\begin{pmatrix}1&b\end{pmatrix}$, $X_{T}$ is the $n_{T}\times 2$ design matrix corresponding to the intercept and running variable $X$ for treated units, $Y_{T}$ is the $n_{T}$-dimensional column vector of treated units’ responses, and $W_{T}$ is a $n_{T}\times n_{T}$ diagonal weight matrix whose entries are

[TABLE]

for some kernel $K(\cdot)$ and bandwidth $h$. The estimator $\hat{\mu}_{C}(b)$ is analogously defined for the control. Then, the estimator for the treatment effect is $\hat{\tau}=\hat{\mu}_{T}(b)-\hat{\mu}_{C}(b)$. Here, $\hat{\mu}_{T}(b)$ and $\hat{\mu}_{C}(b)$ are weighted least squares estimators, where the weights depend on units’ closeness to the boundary $b$. To perform local polynomial regression, one appends higher orders of $X$ to the design matrices $X_{T}$ and $X_{C}$.

The RDD literature has focused on LLR largely because of its boundary bias properties. For example, Hahn et al. (2001) recommend LLR over alternatives like kernel regression because Fan (1992) showed that LLR exhibits better bias properties for boundary points than kernel regression. For more details on the bias comparison between kernel regression and LLR, see Imbens & Lemieux 2008 (p. 624-625) as well as Fan & Gijbels (1992) and Porter (2003).

Furthermore, LLR’s implementation is straightforward once a kernel $K(\cdot)$ and bandwidth $h$ are chosen. The most common choice of $K(\cdot)$ is the rectangular or triangular kernel; Imbens & Lemieux (2008) argue that more complicated kernels rarely make a difference in estimation. Much more attention has been given to the bandwidth choice $h$, largely because the bias is characterized by $h$. In the 2000s, choosing an appropriate $h$ for LLR in RDDs was an open problem: For example, Ludwig & Miller (2007) stated that “there is currently no widely-agreed-upon method for selection of optimal bandwidths…so our strategy is to present results for a broad range of candidate bandwidths.” One widely-used bandwidth selection method is that of Imbens & Kalyanaraman (2012), who derived a data-driven, MSE-optimal bandwidth for LLR estimators. This provided practitioners with clear guidelines for implementing LLR for RDDs, which made its use very popular.

The bandwidth is arguably the most important choice to be made in the LLR methodology for RDDs, because the treatment effect is often sensitive to the bandwidth choice. This motivates sensitivity checks such as that in Ludwig & Miller (2007), where the treatment effect is estimated several times with different bandwidths to ensure that estimates do not vary too greatly. Some have noted that confidence intervals from LLR have a tendency to undercover when a single bandwidth is chosen for inference when the treatment effect is sensitive to the bandwidth choice (Armstrong & Kolesár, 2017; Gelman & Imbens, 2018). A recent extension to the LLR methodology—that of Calonico et al. (2014), hereafter called “robust LLR”—was one of the first methods to address the undercoverage issue of LLR by incorporating a bias correction and inflated confidence intervals corresponding to the uncertainty in estimating the bias correction. Because of its promising inference properties, Calonico et al. (2014) has arguably become the state-of-the-art for conducting inference for the average treatment effect at the boundary in an RDD.

1.3 Other Methods Besides LLR

Other methodologies besides LLR have been proposed for estimating the average treatment effect at the boundary in a sharp RDD. For example, many practitioners have used high-order global polynomials to estimate $\mu_{T}(x)$ and $\mu_{C}(x)$: Matsudaira (2008) argued for a global third-order polynomial regression, and also considered fourth- and fifth-order polynomials as a sensitivity check; similarly, Van der Klaauw (2002) used a global third-order polynomial and noted that LLR could have been an alternative; finally, Card et al. (2004) argued for using a global third-order polynomial regression instead of LLR because the running variable, age, was discrete. However, in recent years many have argued against the use of high-order polynomials in RDDs because of their tendency to yield point estimates and confidence intervals that are highly sensitive to the order of the polynomial (Calonico et al., 2015; Gelman & Imbens, 2018).

Others have focused on local randomization methodologies, where units within a window around the boundary are viewed as-if randomized to treatment and control. For example, Cattaneo et al. (2015) recommends a series of covariate balance tests to decide the window around the boundary such that the as-if randomized assumption is most plausible. Li et al. (2015) extended these ideas to develop a notion of a local overlap assumption and used a Bayesian hierarchical modeling approach for deciding the window around the boundary where this assumption is most plausible. Cattaneo et al. (2017) compared the local randomization approach to local polynomial estimators, and they extended the local randomization approach to incorporate adjustments via parametric models as well.

Finally, others have developed Bayesian methodologies for RDDs. Li et al. (2015) propose a principal stratification approach that provides alternative identification assumptions based on a formal definition of local randomization. Geneletti et al. (2015) propose a Bayesian methodology that incorporates prior information in the treatment effect. Chib & Greenberg (2014) use Bayesian splines to estimate treatment effects in RDDs. Chib & Jacobi (2015) propose a Bayesian methodology specific to fuzzy RDDs.

1.4 Our Proposal: Gaussian Process Regression for Sharp RDDs

We propose a methodology that utilizes Gaussian process regression (GPR), which is one of the most popular nonparametric methodologies in the machine learning and Bayesian modeling literature for estimating unknown functions (Rasmussen & Williams, 2006). The notion of using GPR for RDDs is very much in line with the claim in Hahn et al. (2001) that any nonparametric estimator can be used to estimate the treatment and control response in an RDD. However, to our knowledge, GPR has not been previously proposed for RDDs.

Similar to LLR, our GPR methodology provides a flexible fit to the mean functions $\mu_{T}(x)$ and $\mu_{C}(x)$. Furthermore, our methodology can incorporate both prior knowledge and uncertainty in various parameters in the RDD problem—such as how units are weighted near the boundary—which is not necessarily as straightforward with current LLR methodologies. Finally, our GPR methodology can be used in conjunction with a local randomization perspective. Our methodology adds to the strand of literature started by Calonico et al. (2014) that addresses the undercoverage of standard LLR, as well as the strand of literature on Bayesian methodologies for RDDs.

2 GPR Models to Estimate the Average Treatment Effect at the Boundary

First we review notation for GPR and how GPR is used to estimate a single unknown function. We then discuss GPR models that estimate the two unknown mean functions in sharp RDDs.

2.1 Notation for Estimating One Unknown Function

Define a dataset $\{x_{i},y_{i}\}_{i=1}^{n}$ of responses $\mathbf{y}=(y_{1},\dots,y_{n})$ that varies around some unknown function of the covariate $\mathbf{x}=(x_{1},\dots,x_{n})$:

[TABLE]

where $\epsilon\sim\mathcal{N}_{n}(0,\sigma^{2}_{y}\mathbf{I}_{n})$ and $f(\mathbf{x})\equiv(f(x_{1}),\dots,f(x_{n}))$. If the goal is to well-estimate $\mathbb{E}[f(x^{*})]$ for a particular covariate value $x^{*}$, one option is to specify a functional form for $f(\mathbf{x})$ and then predict $\mathbb{E}[f(x^{*})]$ from this specified model, such as local linear regression, as discussed in Section 1. Instead of specifying a functional form for $f(\mathbf{x})$, we consider nonparametrically inferring $f(\mathbf{x})$ by placing a prior on $f(\mathbf{x})$. A Gaussian process is one such prior:

[TABLE]

for some unknown mean function $m(\mathbf{x})$ and covariance function $K(\mathbf{x},\mathbf{x}^{\prime})$. The notation $f(\mathbf{x})\sim GP(m(\mathbf{x}),K(\mathbf{x},\mathbf{x}^{\prime}))$ denotes a Gaussian process prior on the unknown function $f(\mathbf{x})$, which states that, for any $(x_{1},\dots,x_{n})$, the joint distribution $(f(x_{1}),\dots,f(x_{n}))$ is an $n$-dimensional multivariate normal distribution with mean vector $(m(x_{1}),\dots,m(x_{n}))$ and covariance matrix $K(\mathbf{x},\mathbf{x}^{\prime})$ whose $(i,j)$ entries are $K(x_{i},x_{j}^{\prime})$.

There are many choices one could make for the mean and covariance functions. A common choice for the mean function is $m(\mathbf{x})=\mathbf{0}$; a common choice for the covariance function is the squared exponential covariance function, whose entries are

[TABLE]

Placing a Gaussian process prior with a squared exponential covariance function on $f(\mathbf{x})$ assumes that $f(\mathbf{x})$ is infinitely differentiable, which is similar to other assumptions in the RDD literature (e.g., Assumption 3.3 of Imbens & Kalyanaraman 2012 and Assumption 1 of Calonico et al. 2014). The covariance parameters $\sigma^{2}_{GP}$ and $\ell$ are called the variance and lengthscale, respectively. The variance determines the amplitude of $f(\mathbf{x})$, i.e., how much the function varies from its mean. The lengthscale determines the smoothness of the function: Small lengthscales correspond to $f(\mathbf{x})$ changing rapidly. Most importantly, the covariance function assumes that the response at a particular covariate value $f(x^{*})$ will be similar to the response at covariate values close to $x^{*}$.

Given the Gaussian process prior with mean function $m(\mathbf{x})$ and covariance function $K(\mathbf{x},\mathbf{x}^{\prime})$, as well as their parameters, the posterior for $f(x^{*})$ at any particular covariate value $x^{*}$ can be obtained via standard conditional multivariate normal theory (for an exposition, see Rasmussen & Williams 2006, Pages 16-17):

[TABLE]

[TABLE]

The above posterior can thus be used to obtain point estimates and credible intervals for the value of an unknown function at a particular covariate value $x^{*}$. In practice, the covariance parameters are estimated from the data, such as through maximum likelihood or cross-validation Rasmussen & Williams 2006 (Chapter 5). In Section 2.2 we assume that the covariance parameters are fixed, and in Section 2.3 we extend to a full Bayesian approach that places priors on $\ell,\sigma^{2}_{GP}$, and $\sigma^{2}_{y}$.

2.2 GPR Models for Sharp RDDs

The notion of using GPR to estimate the average treatment effect at the boundary in an RDD suggests a class of models that has not previously been considered in the RDD literature. We focus on two GPR models, which correspond to different assumptions placed on the unknown response functions $\mu_{T}(x)$ and $\mu_{C}(x)$. For each model we show the resulting posterior for the average treatment effect at the boundary and compare it to its analogous LLR model. In Section 2.3 we discuss how—unlike LLR methodologies—the uncertainty in how units are weighted can be incorporated into these GPR models.

For both GPR models, we assume that the treatment response $Y_{i}(1)$ and the control response $Y_{i}(0)$ have the following relationship with the running variable $x_{i}$

[TABLE]

Thus, intuitively, the procedure outlined in Section 2.1 can simply be performed twice—once for $\mu_{T}(x)$ and once for $\mu_{C}(x)$. However, there are assumptions on $\mu_{T}(x)$ and $\mu_{C}(x)$ that, if true, can simplify our GPR models and make inference more precise.

In particular, assumptions can be placed on the covariance structure of $\mu_{T}(x)$ and $\mu_{C}(x)$. The two models we present correspond to two different sets of assumptions—the first assumes that the covariance structure of $\mu_{T}(x)$ and $\mu_{C}(x)$ are the same, while the second allows them to be different. In both models, we assume that $\mu_{T}(x)$ and $\mu_{C}(x)$ are stationary processes, i.e., the covariance parameters of $\mu_{T}(x)$ and $\mu_{C}(x)$ do not vary with the running variable $X$. We discuss cases when this stationarity assumption is inappropriate in Sections 4 and 6.

Same Covariance Assumption: $\text{Cov}(\mu_{T}(x))=\text{Cov}(\mu_{C}(x))$, and $\mu_{T}(x)$ and $\mu_{C}(x)$ are stationary processes.

If the Same Covariance Assumption holds, a natural LLR procedure is to fit local linear regressions on both sides of the boundary with the same bandwidth but different intercepts and slopes. This is largely the standard practice in the RDD literature (Imbens & Lemieux, 2008). Analogously, we place Gaussian process priors on $\mu_{T}(x)$ and $\mu_{C}(x)$ for given mean functions $m_{T}(\mathbf{x})$ and $m_{C}(\mathbf{x})$ and covariance function $K(\mathbf{x},\mathbf{x}^{\prime})$:

[TABLE]

Then, estimates $\hat{\mu}_{T}(b)$ and $\hat{\mu}_{C}(b)$ are obtained, which results in a treatment effect estimate $\hat{\tau}=\hat{\mu}_{T}(b)-\hat{\mu}_{C}(b)$. We now outline how such estimates $\hat{\mu}_{T}(b)$ and $\hat{\mu}_{C}(b)$ are obtained. Using standard conditional multivariate normal theory as in Section 2.1, we have the following posteriors for $\mu_{T}(b)$ and $\mu_{C}(b)$:

[TABLE]

[TABLE]

Here, $\mu_{b|T}$ and $\mu_{b|C}$ denote the posterior mean for $\mu_{T}(b)$ and $\mu_{C}(b)$, respectively, which are in the definition of the treatment effect $\tau$ defined in (5). Note that $\mu_{b|T}$ and $\mu_{b|C}$ are weighted averages of the observed response, where the weights $K(b,\mathbf{x})[K(\mathbf{x},\mathbf{x})+\sigma^{2}_{y}\mathbf{I}]^{-1}$ depend on the covariance parameters $\ell$, $\sigma^{2}_{GP}$, and $\sigma^{2}_{y}$, as well as $\mathbf{x}$. For more discussion on the behavior of the weights $K(b,\mathbf{x})[K(\mathbf{x},\mathbf{x})+\sigma^{2}_{y}\mathbf{I}]^{-1}$, see Rasmussen & Williams (2006, Section 2.6).

The posterior for the treatment effect under the Same Covariance Assumption is then

[TABLE]

where we have also used the independence of $\mu_{T}(x)$ and $\mu_{C}(x)$ stated in (13). If the Same Covariance Assumption does not hold, one can still assume that the mean treatment and control response processes are stationary, but allow both the mean and covariance to vary on either side of the boundary.

Stationary Assumption: $\mu_{T}(x)$ and $\mu_{C}(x)$ are stationary processes.

The posterior in this case would be identical to (17), except the means $\mu_{b|T}$ and $\mu_{b|C}$ and covariances $\Sigma_{b|T}$ and $\Sigma_{b|C}$ are instead defined as

[TABLE]

i.e., the shared covariance $K(\cdot,\cdot)$ is replaced with $K_{T}(\cdot,\cdot)$ for units receiving treatment and $K_{C}(\cdot,\cdot)$ for units receiving control. The analogous LLR methodology would be to allow different intercepts, slopes, and bandwidths on either side of the boundary. However, using different bandwidths on either side of the boundary is rarely done in practice. For example, Imbens & Lemieux (2008) argue that if the curvature of $\mu_{T}(x)$ and $\mu_{C}(x)$ is the same, then the large-sample optimal bandwidths should be the same; and, furthermore, there is additional variance in estimating two optimal bandwidths rather than one, due to the smaller sample used to estimate each bandwidth. Thus, a benefit of the Same Covariance Assumption is that it allows researchers to use the entire data to estimate one set of covariance parameters, instead of estimating two separate sets of covariance parameters for treatment and control. However, when fitting our GPR model, we do not recommend sharing information between $\mu_{T}(x)$ and $\mu_{C}(x)$ beyond estimating their covariance structure—this follows the general practice in the RDD literature to fit separate regression functions (that may nonetheless share the same bandwidth) for the treatment and control groups.

The above posteriors for these two GPR models assume fixed mean and covariance parameters. In practice, maximum-likelihood or cross-validation can be used for estimating these parameters. In Section 2.3, we extend the above GPR models to a full-Bayesian approach that incorporates uncertainty in the mean and covariance parameters.

2.3 Accounting for Mean and Covariance Function Uncertainty

The GPR models in Section 2.2 assume that $\mu_{T}(b)$ and $\mu_{C}(b)$ will be similar to $\mu_{T}(x)$ and $\mu_{C}(x)$, respectively, for $x$ near $b$. The extent of this similarity is determined by the covariance function $K(\mathbf{x},\mathbf{x}^{\prime})$ and its parameters. In particular, recall that in Section 2.2 we showed that the posterior mean of the average treatment effect for GPR is characterized by a difference of two weighted averages, where the weights are of the form $K(b,\mathbf{x})[K(\mathbf{x},\mathbf{x})+\sigma^{2}_{y}\mathbf{I}]^{-1}$. Thus, incorporating uncertainty in the covariance parameters in turn incorporates uncertainty in how units are weighted when estimating the average treatment effect.

Denote the mean function parameters by $\boldsymbol{\theta}_{m}$ and the covariance function parameters by $\boldsymbol{\theta}_{K}$. For example, consider the mean functions

[TABLE]

where $\mathbf{h}(x)=(1,x,\dots,x^{p-1})$, and $\boldsymbol{\beta}_{T}$ and $\boldsymbol{\beta}_{C}$ are the corresponding $p$-dimensional column vectors. In this case, $\boldsymbol{\theta}_{m}=(\boldsymbol{\beta}_{T},\boldsymbol{\beta}_{C})$. For the squared exponential covariance function defined in (10), $\boldsymbol{\theta}_{K}=(\ell,\sigma^{2}_{GP},\sigma^{2}_{y})$.

In order to incorporate uncertainty in $\boldsymbol{\theta_{m}}$ and $\boldsymbol{\theta_{K}}$, one can first obtain draws $1,\dots,D$ from the joint posterior of $(\boldsymbol{\theta}_{m},\boldsymbol{\theta}_{K})$, rather than obtaining point-estimates $\hat{\boldsymbol{\theta}}_{m}$ and $\hat{\boldsymbol{\theta}}_{K}$. Then, for each draw $(\boldsymbol{\theta}_{m},\boldsymbol{\theta}_{K})_{1},\dots,(\boldsymbol{\theta}_{m},\boldsymbol{\theta}_{K})_{D}$, one draws from the posterior for $\tau$, defined in (17).

Section 2.2 already defines the likelihood for the GPR models, so all that remains is to specify priors for $(\boldsymbol{\theta}_{m},\boldsymbol{\theta}_{K})$ in order to obtain draws from the joint posterior of $(\boldsymbol{\theta}_{m},\boldsymbol{\theta}_{K})$. Priors for $(\boldsymbol{\theta}_{m},\boldsymbol{\theta}_{K})$ will depend on the choice of mean and covariance functions. For example, for the mean functions defined in (19), we recommend $\mathcal{N}(\mathbf{0},B)$ as a prior for $\boldsymbol{\beta}_{T}$ and $\boldsymbol{\beta}_{C}$, where $B$ is a $p\times p$ diagonal matrix with reasonably large entries (Rasmussen & Williams 2006, Pages 28-29). For the squared exponential covariance function given by (10), we recommend half-Cauchy priors for the covariance parameters $\frac{1}{\ell^{2}}$ and $\sigma^{2}_{GP}$ and noise $\sigma^{2}_{y}$, following advice from Gelman (2006) about stable priors for variance parameters.

Now we prove that our GPR methodology is consistent in estimating the average treatment effect at the boundary. First we establish consistency when the covariance parameters are fixed, and then we consider the case where priors are placed on the covariance parameters.

3 Posterior Consistency of our GPR Models

Gaussian processes are known to exhibit posterior consistency under minimal assumptions. Ghosal & Roy (2006) proved posterior consistency of binary GPR for fixed covariance parameters, and Choi & Schervish (2007) proved posterior consistency of GPR when the response is continuous. More generally, van der Vaart & van Zanten (2008) studied the contraction rate for Gaussian process priors for density estimation and regression problems, and van der Vaart & van Zanten (2009) extended these results to when a prior is placed on the lengthscale of a Gaussian process.

Here we evaluate our Bayesian methodology from a frequentist point-of-view, which assumes a fixed treatment effect at the boundary. The GPR models in Section 2.2 estimate the treatment effect as the difference between two Gaussian process regressions; thus, our posterior of the treatment effect is consistent if the separate GPRs on either side of the discontinuity are consistent. First we establish posterior consistency assuming the covariance parameters are fixed, as in Section 2.2, and then we extend these results to when a prior is placed on the covariance parameters, as in Section 2.3.

We prove posterior consistency assuming the Stationary Assumption in Section 2.2, but the results also hold for the Same Covariance Assumption. Furthermore, we assume the mean functions $m_{T}(x)=0$, $m_{C}(x)=0$ and the squared exponential covariance function defined in (10). Discussion about the extent to which these results extend to other choices for the mean and covariance functions is in the Appendix. The other assumptions necessary for Theorems 1 and 2 below follow van der Vaart & van Zanten (2009) and are also given in the Appendix.

Theorem 1: Assume that the Stationary Assumption holds, the covariance functions $K_{T}(x,x)$ and $K_{C}(x,x)$ are fixed, and Assumptions A1, A2, and A3 given in the Appendix hold. Denote the true average treatment effect at the boundary as $\tau^{*}=\mu_{T}^{*}(b)-\mu_{C}^{*}(b)$, where $\mu_{T}^{*}(x)$ and $\mu_{C}^{*}(x)$ are the true mean treatment and control response functions in the model (13). Let $\prod(\tau|x_{1},\dots,x_{n})$ denote the posterior distribution of the average treatment effect at the boundary, defined in (17). Then, this posterior is consistent, in the sense that

[TABLE]

for sufficiently large $M$, where $h$ is the Hellinger distance, and $\epsilon_{n}$ is the rate at which the posterior of $\tau$ contracts to the true $\tau^{*}$.

The proof of Theorem 1, as well as a discussion about the nature of the contraction rate, is given in the Appendix. Theorem 2 establishes posterior consistency when a prior is placed on the lengthscale parameter $\ell$, instead of being held fixed (see Section 2.3). A discussion about posterior consistency when an additional prior is placed on $\sigma^{2}_{GP}$ is in the Appendix.

Theorem 2: Assume that the Stationary Assumption holds, the $\sigma^{2}_{GP}$ parameters in $K_{T}(x,x)$ and $K_{C}(x,x)$ are fixed, and Assumptions A1, A2, A3, and A4 given in the Appendix hold. Then, Theorem 1 holds.

The proof of Theorem 2 is given in the Appendix. A corollary follows from the proofs of Theorem 1 and 2.

Corollary 1: Theorems 1 and 2 hold if the Same Covariance Assumption holds instead of the Stationary Assumption.

4 Simulation Results

Here we investigate how our Gaussian Process model compares to standard LLR and the robust LLR method introduced in Calonico et al. (2014). We choose these two methods because the former is the standard in both applied work and the RDD literature at large, and the latter is a recent method that attempts to solve the undercoverage issue of standard LLR. We focus on the GPR model assuming the Same Covariance Assumption in Section 2.2 and use the mean functions

[TABLE]

and the squared exponential covariance function given by (10). These assumptions are analogous to the LLR procedure of fitting separate local linear regressions in treatment and control with differing slopes but the same bandwidth. Specification of the mean function in the Gaussian process prior is typically not consequential for estimation; however, as discussed in Rasmussen & Williams (2006, Section 2.7), there can be some benefits to specifying a mean function, as we do here. In particular, the above specification allows GPR predictions to pull towards a global linear trend instead of a global mean (which would be the case if we used a zero mean function—a common choice in the literature—for the Gaussian process prior). This can be useful within the context of extrapolation towards the boundary, as in an RDD. In the Appendix in Table 3, we present simulation results for our GPR methodology using a zero mean function instead of the above linear mean function for the Gaussian process prior. The results for that case are largely the same as the results presented here, which suggests that our results are insensitive to specification of the mean function in the Gaussian process prior.

As discussed in Section 2.3, we took a full-Bayesian approach to our GPR methodology and placed independent $\mathcal{N}(0,100^{2})$ priors on the mean function parameters in (21) and independent $\text{half-Cauchy}(0,5)$ priors on the covariance parameters. These choices for the prior distributions are in line with common recommendations in the Bayesian data analysis literature: The choice of Normal priors on the mean function parameters follows recommendations from Rasmussen & Williams (2006), and the choice of half-Cauchy priors on the covariance parameters follows recommendations from Gelman (2006), Polson & Scott (2012), and Gelman et al. (2013, Chapter 5). We used the R package rstan (Carpenter et al., 2016) to sample from the posterior of these parameters. In the Appendix in Table 3, we discuss simulation results for GPR when we instead plug in the MLE for the covariance parameters; however, we found that the full-Bayesian approach is preferable in terms of inferential properties, which suggests that it is beneficial to incorporate uncertainty in the covariance parameters for our GPR method.

We conduct a simulation study based on simulations from Imbens & Kalyanaraman (2012) and Calonico et al. (2014). In all simulations, we generated 1,000 datasets of 500 observations $\{(x_{i},y_{i},\epsilon_{i}):i=1,\dots,500\}$, where $x_{i}\sim 2\text{Beta}(2,4)-1$, $\epsilon_{i}\sim N(0,0.1295^{2})$, and $y_{i}=\mu_{j}(x_{i})+\epsilon_{i}$ for different mean functions $\mu_{j}(x_{i})$.

We consider seven different mean functions (see Figure 1), which we call Lee, Quad, Cate 1, Cate 2, Ludwig, Curvature, and Cubic. Lee, Quad, Cate 1, and Cate 2 were used in Imbens & Kalyanaraman (2012), and Lee, Ludwig, and Curvature were used in Calonico et al. (2014); details about these datasets can be found in Imbens & Kalyanaraman 2012 (Page 18) and Calonico et al. 2014 (Page 20). We also introduce the Cubic mean function as a comparison to the Quad mean function, because in the Quad mean function the linear trends on either side of $b$ are the opposite sign, whereas those for the Cubic mean function are the same sign.

The boundary for each dataset is $b=0$. The treatment effect is $\tau=0.04$ for Lee and Curvature, $\tau=0$ for Quad and Cubic, $\tau=0.1$ for Cate 1 and 2, and $\tau=-3.35$ for Ludwig. Also displayed in Figure 1 are a set of sample points $\{(x_{i},y_{i}),i=1,\dots,500\}$ for each mean function, which shows what one dataset looks like for each mean function. One can see that—although $\epsilon_{i}\sim N(0,0.1295^{2})$ for all mean functions—the relative noise varies across mean functions.

For standard LLR and robust LLR we used the rdrobust R package (Calonico et al., 2017). For both methods, we used an MSE-optimal bandwidth that is the default in rdrobust. Simulation results using the bandwidth introduced in Imbens & Kalyanaraman (2012)—also known as the IK bandwidth, which has been widely used in practice—are provided in the Appendix in Table 4, and results using the coverage error rate optimal bandwidth—an alternative bandwidth choice within the rdrobust R package that is also discussed in Calonico et al. (2018a)—are provided in the Appendix in Table 5. The results using those bandwidths are largely the same as the results presented here. Simulation results for other bandwidth choices appear in Calonico et al. (2014).

Imbens & Kalyanaraman (2012) ran a simulation study that focused on the Lee, Quad, Cate 1, and Cate 2 mean functions, and they compared different bandwidth selectors for LLR in terms of bias and root mean squared error (RMSE). Similarly, Calonico et al. (2014) ran a simulation study that focused on the Lee, Ludwig, and Curvature mean functions, and they compared different bandwidth selectors for LLR and their methodology in terms of coverage and mean interval length (IL). We synthesize these simulation studies and report in Figure 2 how LLR, robust LLR, and two versions of our GPR methodology perform on the seven mean functions in Figure 1 in terms of coverage, IL, absolute bias, and RMSE. The numbers plotted in Figure 2 are in Tables 2 and 3 in the Appendix. Point estimates and 95% confidence intervals for LLR and robust LLR were obtained from rdrobust. Point estimates and 95% credible intervals for our methodology corresponded to the mean and 2.5% and 97.5% quantiles, respectively, of the posterior of the average treatment effect, shown in (17).

Robust LLR is meant to improve the coverage of LLR, and indeed it does for all datasets. The better coverage is in part due to wider confidence intervals (see the systematically higher mean interval length at the top right of Figure 2) and in part due to better bias properties (see the bottom left of Figure 2). Robust LLR also tends to exhibit worse RMSE than LLR (see the bottom right of Figure 2).

Our primary method (“GPR”) tends to exhibit narrower intervals than both LLR and robust LLR. GPR also exhibits better coverage than both methods, except for the Lee and Ludwig datasets. Furthermore, our method tends to exhibit lower RMSE than both LLR and robust LLR. However, our method always exhibits more bias than robust LLR, which explicitly uses a bias correction.

For the Lee and Ludwig datasets, our method did worse than robust LLR in terms of coverage and bias, and this may be because it is inappropriate to assume that $\mu_{T}(x)$ and $\mu_{C}(x)$ are stationary processes—i.e., that the covariance parameters do not vary across the running variable $X$—in these cases. By assuming stationarity, our GPR model uses data both close to and far from the boundary to estimate the single variance $\sigma^{2}_{GP}$ and lengthscale $\ell$. This assumption is related to the stability of the second derivative of $\mu_{T}(x)$ and $\mu_{C}(x)$, because the covariance parameters of a Gaussian process dictate their derivative processes; see Wang (2012) for a further discussion of this relationship. Figure 3 displays the absolute second derivative of the seven mean functions (in blue). The Lee and Ludwig datasets are characterized by the absolute second derivative rapidly increasing near the boundary. Our GPR methodology likely does not do well for these datasets because we are using data far from the boundary to estimate the overall curvature of the mean function, which leads us to underestimating the curvature of the Lee and Ludwig mean functions at the boundary.

Figure 3 also displays $\frac{\hat{\sigma}_{GP}}{\hat{\ell}}$, the ratio of the maximum-likelihood estimates of the covariance parameters, which was computed within a sliding window of a noiseless version of the seven mean functions. Although there is not a one-to-one correspondence between the absolute second derivative and $\frac{\hat{\sigma}_{GP}}{\hat{\ell}}$, their behavior is notably similar, which reinforces the idea that both the variance $\sigma^{2}_{GP}$ and lengthscale $\ell$ play a role in capturing the curvature of the mean function. Furthermore, this connection between the second derivative and the covariance parameters further suggests a similarity between the covariance parameters in our GPR methodology and the bandwidth in the LLR methodology, because the IK bandwidth is estimated as a nonlinear function of the estimated second derivative at the boundary $b$ (Imbens & Kalyanaraman, 2012).

If one did not believe the Stationary Assumption held, one alternative would be to only use data close to the boundary when fitting our GPR model. This is our second method, “GPR (Cut),” whose coverage, mean IL, absolute bias, and RMSE is also displayed in Figure 2. For each of the 1,000 replications of the seven mean functions, we first estimated the IK bandwidth with a rectangular kernel; then, we fit our GPR model within this estimated bandwidth. This procedure improves upon our GPR model for and only for the Lee and Ludwig datasets—the coverage increased to 94.3% and 88.9%, respectively, and the bias improved for the Lee and Ludwig datasets while staying the same for the other datasets. While the coverage also improved for the Cate 1 and Cate 2 datasets, this is likely due to the increase in the interval length. Because results improved only for these two datasets, this further demonstrates that using our GPR model on the whole dataset can be beneficial when $\mu_{T}(x)$ and $\mu_{C}(x)$ are stationary processes; otherwise, it may be preferable to only fit our GPR model to data close to the boundary.

Furthermore, this suggests that our method can be combined with a local randomization perspective for RDDs (e.g., Li et al. 2015): One can first determine the window around the boundary where units are “as-if randomized” by using covariate balance tests such as Cattaneo et al. (2015) and Li et al. (2015) and then use our GPR methodology within this window around the boundary. This approach is similar to Cattaneo et al. (2017), who combined the local randomization perspective with parametric models for estimating the average treatment effect at the boundary.

Overall, GPR performs well compared to LLR and robust LLR. In particular, our GPR method tends to yield better interval length and RMSE properties than LLR and robust LLR, and it also yields better coverage when the underlying mean functions are stationary across the running variable $X$. The issue of undercoverage in LLR methodologies has been relatively unaddressed in the RDD literature, except for robust LLR (Calonico, Cattaneo, & Titiunik, 2014), and so our GPR method can be viewed as the second method to yield promising coverage properties for RDD analyses while also providing a flexible fit to the underlying mean functions $\mu_{T}(x)$ and $\mu_{C}(x)$. When $\mu_{T}(x)$ and $\mu_{C}(x)$ are nonstationary across $X$, our GPR methodology could be extended to incorporate a lengthscale function $\ell(x)$ instead of a single lengthscale $\ell$. However, such a lengthscale function would either need to be specified (see Rasmussen and Williams 2006, Chapter 4, for an example), or estimated via another Gaussian process prior (Plagemann, Kersting, & Burgard, 2008). In a similar vein, there has also been work on using dimension expansion to model nonstationary processes (Bornn et al., 2012). However, in an RDD, we only need a good estimate of the covariance parameters near the boundary, rather than across the entire mean functions $\mu_{T}(x)$ and $\mu_{C}(x)$. More work needs to be done to determine the optimal amount of data to include in our GPR model for estimating these parameters and the average treatment effect at the boundary in the case of nonstationary processes.

Now we compare how LLR, robust LLR, and GPR perform on a real dataset from the National Basketball Association.

5 Empirical Example: The NBA Draft

The National Basketball Association (NBA) draft, held annually, is divided into $2$ rounds, where each NBA teams gets one selection per round to draft a player of their choice. Because players are picked sequentially, there is no reason to believe there is a marked skill difference between the last pick of the first round and first pick of the second round. However, because of the difference in the perceived value of first-round versus second-round picks, as well as differing contract structures between the two rounds, we suspect that first-round picks are treated more favorably and given more playing time than their second-round colleagues, above and beyond what can be explained by differences in skill. As such, we seek to explore if there is a difference between first- and second-round picks in both skill and playing time.

We want to estimate the treatment effect of having a second-round contract instead of a first-round contract on four basketball player outcomes: box plus-minus, win shares played, number of minutes played, and number of games played. The first two are overall measures of player performance (Kubatko, 2009; Myers, 2015), while the latter two are measures of playing time. Our data include the pick number and the four aforementioned outcomes for 1,238 NBA basketball players drafted between 1995 and 2016. Due to anomalies created by the NBA expanding from 29 teams to 30 teams in 2004, as well as some years with teams forfeiting picks, we shifted the pick numbers to ensure that $b=30.5$ marked the discontinuity between the first- and second-round picks in each year. Figure 4 displays the NBA player data grouped by pick number for each of the four outcomes (after aligning pick numbers). Grouping by pick number allows us to understand the average performance and playing time of players drafted at each pick number, which is a standard approach in draft evaluation (Silver, 2014).

Even though the running variable in this case is discrete—which can cause complications in some regression discontinuity analyses (Lee & Card, 2008; Kolesár & Rothe, 2018)—we nonetheless apply LLR, robust LLR, and our GPR method to these data to compare how they perform in practice. As in Section 4, we used the rdrobust R package to implement LLR and robust LLR using the default MSE-optimal bandwidth. For GPR, we used the squared exponential covariance function assuming the Same Covariance Assumption; furthermore, we took the full-Bayesian approach to our GPR methodology and—similar to Section 4—placed independent $\mathcal{N}(0,100^{2})$ priors on the mean function parameters and half-Cauchy$(0,5)$ priors on the covariance parameters.

Figure 5 shows the estimated mean functions and corresponding confidence intervals for LLR and GPR, and Table 1 shows the treatment effect point estimates and corresponding confidence intervals for LLR, robust LLR, and GPR.111Robust LLR yields an inflated confidence interval for the treatment effect specifically, which depends on a bias correction at the boundary. Thus, robust LLR does not yield confidence intervals for the entire mean functions, because the bias correction is boundary-specific. Thus, only LLR and GPR are displayed in Figure 5, while LLR, robust LLR, and GPR are all discussed in Table 1. Furthermore, for robust LLR in Table 1, we report the bias-corrected point estimate given by the rdrobust package. The estimated mean functions for LLR and GPR are quite similar to each other for all outcomes, with GPR yielding slightly wider confidence intervals, as expected. All three methods find the number of games played for second-round picks to be significantly lower than that of first-round picks. Furthermore, LLR and robust LLR find the box plus-minus for second-round picks to be significantly lower than that of first-round picks; meanwhile, GPR finds this difference to be borderline insigificant. Out of these three methods, the results from our GPR method are most in line with previous reports—both qualitative (Barber, 2016) and quantitative (Koenig, 2012)—on the difference between first- and second-round NBA basketball players, which have claimed that there is a difference in attention given to first-round picks (e.g., games played) but not a difference in player ability (e.g., box plus-minus and win shares).

In summary, using our GPR methodology, we find that the treatment effect of being a second-round pick significantly reduces the number of games played and marginally reduces the number of minutes played. This suggests that there is a drop-off in playing time for second-round players relative to their first-round counterparts beyond that explainable by the natural drop-off in playing time between successive picks. Furthermore, we find that there is not a significant difference in player ability between first- and second-round basketball players at the boundary between the 30th and 31st picks that divides the first and second rounds of the NBA draft.

6 Discussion and Conclusion

Local linear regression (LLR) and its variants are currently the most common methodologies for estimating the average treatment effect at the boundary in RDDs. These methods are popular because they are easy to implement and there is a large literature on their theoretical properties. However, recent works have noted that LLR tends to yield confidence intervals that undercover, and new methodologies—namely that of Calonico et al. (2014)—have tried to address this issue. As an alternative to LLR, we proposed a Gaussian Process regression (GPR) methodology that flexibly fits the treatment and control response functions by placing a general prior on the mean response functions. We showed via simulation that our GPR methodology tends to outperform standard LLR and the state-of-the-art methodology of Calonico et al. (2014) in terms of coverage, interval length, and RMSE. Furthermore, we used our GPR methodology on a real-world sharp RDD in the National Basketball Association (NBA) draft and found that GPR yielded results that were more in line with previous reports on the NBA draft than were results from LLR methods. Overall, our methodology addresses the undercoverage issue commonly seen in RDDs without sacrificing too much power to detect treatment effects, thereby adding to the growing literatures on improving inference for RDDs (Calonico et al., 2014, 2018b) and on Bayesian methods for RDDs (Chib & Greenberg, 2014; Chib & Jacobi, 2015; Geneletti et al., 2015; Li et al., 2015).

Our methodology focuses on flexibly fitting the mean treatment and control responses while also improving coverage properties; however, there are many other issues of interest in the RDD literature. For example, we only consider sharp RDDs; Li et al. (2015) and Chib & Jacobi (2015) provide Bayesian methodologies for fuzzy RDDs which could likely be combined with our Gaussian process approach. Furthermore, we focused on RDDs that only have one background covariate—the running variable—but other covariates could be included in our GPR methodology to improve the precision of our treatment effect estimator (e.g., Calonico et al. 2016). See Imbens & Lemieux (2008) for a review of these other RDD concerns.

Furthermore, while our method can incorporate prior knowledge and uncertainty in various parameters that are typically discussed in the RDD literature, it is not necessarily clear when this should be done for our GPR model. For example, Hall & Kang (2001) found that incorporating the uncertainty in the bandwidth—which is somewhat analogous to the covariance parameters in our GPR model—for density estimators via bootstrapping can be inconsequential or even detrimental in some cases. Although our simulation study suggests that it is beneficial to take a full-Bayesian approach and propogate uncertainty in the GPR model parameters when estimating treatment effects in RDDs, more work needs to be done to determine when it is most appropriate to incorporate uncertainty in these parameters. Furthermore, we only focused on the squared exponential covariance function for Gaussian processes, but other covariance functions should be considered. Because estimating the average treatment effect at the boundary in an RDD is fundamentally an extrapolation issue, covariance functions whose purpose is to extrapolate well (e.g., Wilson & Adams 2013) may be particularly suitable for RDDs.

Additionally, although our GPR methodology exhibits arguably better coverage, interval length, and RMSE properties than standard methodologies in the literature, there are ways our methodology could be improved, even for the case of a sharp RDD with only one covariate. Our methodology does not perform well when the mean response functions $\mu_{T}(x)$ and $\mu_{C}(x)$ are nonstationary across the running variable $X$. More work needs to be done to model nonstationary processes within the context of RDDs, such as by estimating a lengthscale function $\ell(x)$ or by determining the optimal amount of data to use when estimating the covariance parameters for our GPR methodology. Relatedly, our simulation results suggest that our GPR methodology can be used in combination with a local randomization framework for RDDs (such as that seen in Li et al. 2015).

Finally, a promising avenue for future research is extending our GPR methodology beyond one-dimensional sharp RDDs. In particular, recent work has explored geographic RDDs, where spatial variation in outcomes must be accounted for when estimating the treatment effect along a geographic boundary (Keele & Titiunik, 2014; Keele et al., 2015). Because GPR has been widely used in spatial statistics (Banerjee et al., 2014; Cressie, 2015), it may be particularly suitable for geographic RDDs. We explore the use of our GPR methodology for geographic RDDs in Rischard et al. (2018).

7 Appendix

7.1 Complete Simulation Results from Section 4

7.2 Assumptions for Posterior Consistency Proofs

Assumptions on the Running Variable $X$ (Assumptions A1)

Assume the control running variable values $\{x_{i}\}_{i=1}^{n_{C}}$ are known elements of $[b_{C},b]$ and the treatment running variable values $\{x_{i}\}_{i=1}^{n_{T}}$ are known elements of $[b,b_{T}]$, for some $b_{C},b_{T}$, and boundary $b$.

Assumptions on the Response Function (Assumptions A2)

Let $C^{\alpha}[b,b_{T}]$ and $C^{\alpha}[b_{C},b]$ be H$\ddot{\text{o}}$lder spaces of $\alpha$-smooth functions $f:[b,b_{T}]\rightarrow\mathbb{R}$ and $f:[b_{C},b]\rightarrow\mathbb{R}$, respectively. Assume $\mu_{T}(x)\in C^{\alpha}[b,b_{T}]$ and $\mu_{C}(x)\in C^{\alpha}[b_{C},b]$, where $\mu_{T}(x)$ and $\mu_{C}(x)$ are the treatment and control response functions. 2. 2.

The treatment and control responses $\{y_{i}\}_{i=1}^{n_{T}}$ and $\{y_{i}\}_{i=1}^{n_{C}}$ have the following relationships with the running variable:

[TABLE]

for mean response functions $\mu_{T}(x)$ and $\mu_{C}(x)$ and independent errors $\{\epsilon_{i}\}_{i=1}^{n_{C}}\sim N(0,\sigma^{2}_{y0})$ and $\{\epsilon_{i}\}_{i=1}^{n_{T}}\sim N(0,\sigma^{2}_{y1})$.

Assumptions on the Noise (Assumptions A3)

The priors on $\sigma_{y0}$ and $\sigma_{y1}$ have support on compact intervals that are subsets of $(0,\infty)$ which contain the true errors $\sigma_{y0}$ and $\sigma_{y1}$, respectively.

Assumptions on the Prior for $\ell$ (Assumptions A4)

The lengthscale $\ell$ has a prior distribution $\kappa$ such that, for positive constants $C_{1},D_{1},C_{2},D_{2}$, nonnegative constants $p,q$, and every sufficiently large $a>0$,

[TABLE]

7.3 Proof of Theorems 1 and 2

Proof of Theorem 1: Assume that the Stationary Assumption holds, the covariance functions $K_{T}(x,x)$ and $K_{C}(x,x)$ are fixed, and Assumptions A1, A2, and A3 hold. Then, according to Theorem 2.2 in van der Vaart & van Zanten (2009), for sufficiently large $M_{T}$ and $M_{C}$,

[TABLE]

for some contraction rates $\epsilon_{n_{T}}$ and $\epsilon_{n_{C}}$, where $h$ is the Hellinger distance. The nature of the contraction rates $\epsilon_{n_{T}}$ and $\epsilon_{n_{C}}$ are discussed in van der Vaart & van Zanten (2009), and depend on the differentiability and smoothness of the true $\mu_{T}^{*}(x)$ and $\mu_{C}^{*}(x)$. Note that the only covariate value for which both of these hold is $x=b$, because by Assumption A1, the intersection of the supports of $\{x_{i}\}_{i=1}^{n_{C}}$ and $\{x_{i}\}_{i=1}^{n_{T}}$ is only the boundary $b$.

Because the Hellinger distance is symmetric and satisfies the triangle inequality,

[TABLE]

where the last inequality holds with posterior probability 1, by (25). Therefore,

[TABLE]

where $M\equiv\sqrt{2}\cdot\text{max}(M_{T},M_{C})$ and $\epsilon_{n}\equiv\sqrt{2}\cdot\text{max}(\epsilon_{n_{T}},\epsilon_{n_{C}})$, so that $M\epsilon_{n}\geq M_{T}\epsilon_{n_{T}}+M_{C}\epsilon_{n_{C}}$. $\blacksquare$

Proof of Theorem 2: Assume that the Stationary Assumption holds, the $\sigma^{2}_{GP}$ parameters in $K_{T}(x,x)$ and $K_{C}(x,x)$ are fixed, and Assumptions A1, A2, A3, and A4 hold. By Theorem 3.1 of van der Vaart & van Zanten (2009), (25) holds, and then the proof of Theorem 2 is identical to that of Theorem 1. $\blacksquare$

7.4 Extending Theorems 1 and 2 to Other Mean and Covariance Functions and Random Variance

Rasmussen & Williams 2006 (Section 2.7) discuss other choices of mean functions besides $m(x)=0$, particularly mean functions of the form $m(x)=\mathbf{h}(x)^{T}\boldsymbol{\beta}$ for some fixed basis functions $\mathbf{h}(x)$. However, Rasmussen & Williams (2006) argue that different choices of $m(x)$ are more for interpretability than predictive accuracy, because $m(x)=0$ does not constrain the posterior mean to be zero, and so different choices of $m(x)$ likely do not affect the consistency results in Section 3. To the best of our knowledge, the literature has focused primarily on the choice $m(x)=0$ for posterior consistency of GPR (e.g., Choi & Schervish 2007 and van der Vaart & van Zanten 2009).

van der Vaart & van Zanten (2009) discuss how their Theorems 2.2 and 3.1 (which correspond to our Theorems 1 and 2, respectively) extend to covariance functions besides the squared exponential. Specifically, their results hold for processes whose spectral measure has subexponential tails; the squared exponential process falls under this class, but van der Vaart & van Zanten (2009) discuss other processes that fall under this class as well, and our Theorems 1 and 2 also hold for those processes.

Finally, the literature has focused on posterior consistency for the case when a prior is placed on the lengthscale $\ell$ but not on the variance $\sigma^{2}_{GP}$, as in van der Vaart & van Zanten (2009). To the best of our knowledge, Choi (2007) is the only work to consider posterior consistency when priors are placed on both $\ell$ and $\sigma^{2}_{GP}$; however, these results only hold for binary Gaussian process regression, and their necessary assumptions are more restrictive than those presented in this paper. We leave posterior consistency of our GPR method when priors are placed on both $\ell$ and $\sigma^{2}_{GP}$ as future work.