2.1 Additional Details

We now consider a set of interrelated scaling and hyperparameter selection choices, including the kernel choice. It is not clear whether an optimal solution to these questions exists in a general sense, yet users must make these choices, and some choices can lead to poor performance. We offer guidance that we believe is at least reasoned, practical, and shows good performance to the limits of our testing.

Demeaning and Rescaling: We demean and scale the covariates to have unit variance. This is common for kernel-based machine learning methods, and avoids dependency on unit-of-measure decisions. As noted, we also (globally) demean and scale Y to have unit variance, as in some other kernelized approaches (see, e.g., Hainmueller and Hazlett Reference Hainmueller and Hazlett2014). This is also required for the interpretation of $\sigma ^2$ as 1- $R^2$ noted below.

Bandwidth Selection: For the Gaussian kernel $k(X_i,X_j) = exp(-||X_i-X_j||^2/b)$ , we must choose the bandwidth or “length-scale,” b. There are numerous strategies for doing so. One reasonable approach is to set b equal to (or proportional to) the number of dimensions of X, so that the result does not explode or go to zero as the number of dimensions changes (as in Hainmueller and Hazlett Reference Hainmueller and Hazlett2014). Another recent proposal, offered in Hartman, Hazlett, and Sterbenz (Reference Hartman, Hazlett and Sterbenz2025), chooses the value of b that leads to the highest estimated variance across the (off-diagonal) elements of ${\mathbf {K}}$ . This simply ensures that the columns of ${\mathbf {K}}$ stand to be highly informative rather than having ${\mathbf {K}}$ approach the identity matrix (if b is too small) or a block of ones (if b is too large). This is the approach we take.

Choosing $\ \sigma ^2$ : Having fixed $\sigma _f=1$ and selected the kernel bandwidth based directly on the distribution of ${\mathbf {K}}$ , we are left with only one hyperparameter to tune. As advocated by Rasmussen and Williams (Reference Rasmussen and Williams2006), we tune the residual variance term, $\sigma ^2$ , by maximizing the log marginal likelihood,

(5) $$ \begin{align} \log~p(Y \mid X) = -\frac{1}{2} Y^\top ({\mathbf{K}} + \sigma^2 I)^{-1}Y - \frac{1}{2}(\log |{\mathbf{K}}+\sigma^2 I|) - \frac{n}{2} \log 2\pi. \end{align} $$

In principle, since one could always choose a “wigglier CEF” (even for a fixed kernel bandwidth) together with a smaller $\sigma ^2$ , this quantity might not always be well-identified on the data. While we remain concerned about this in the general case, empirically, this approach appears to perform very well at least over the cases we examine. In particular, we see that for each DGP we attempt, the estimated $\sigma ^2$ values are very close to 1- $R^2$ , as our scaling choice should lead to. This may fail when the DGP does not have the Gaussian, homoskedastic noise on which this likelihood estimate is based.

Combining Kernels: A useful fact about kernels used for such kernel regression tasks is that they can be added together, with the resulting kernel remaining valid (positive semi-definite), and providing access to a space of functions that combines the underlying features of the constituent kernels (see, e.g., Schulz et al. Reference Schulz, Tenenbaum, Duvenaud, Speekenbrink and Gershman2017). For example, consider a kernel ${\mathbf {K}}^{poly}$ that provides access to polynomial models, i.e., ${\mathbf {K}}^{poly}_{i,j}= 1 + \langle x_i, x_j \rangle ^p$ , where p is the desired degree of the polynomial. Such a kernel matrix can be added pointwise to the Gaussian kernel matrix. The resulting kernel matrix allows the GP to represent functions that are both locally smooth (from the Gaussian kernel), and exhibit polynomial growth rather than translation invariance (from the polynomial). Similarly, including a periodic kernel enables the model to capture cyclic structures such as seasonality. We demonstrate this below to explore models for ITS analysis.

Stationary vs. Non-Stationary Kernels for Mild vs. Extreme Extrapolation. The next question is whether one is interested in making inferences nearer to observed data (e.g., interpolation, edge estimation, or very mild extrapolation), where we can rely on the covariance between observed points and target points directly, or if one wishes to make inferences well beyond the data’s edge. In the former case, we can employ a “stationary kernel” that operates only on the relative distances between observations. In these settings, we rely on the Gaussian kernel because (i) it directly works with the logic that observations with more similar values should have greater covariance while covariance should drop to zero as observations grow far apart in the covariate space; (ii) it is commonly taken to be the work-horse kernel for many machine learning approaches involving kernels; and (iii) it is an example of a “universal kernel” for the continuous functions, meaning that it can represent any continuous function given sufficient data (Micchelli, Xu, and Zhang Reference Micchelli, Xu and Zhang2006).

By contrast, for some applications, inferences are attempted farther from the edge of the data. Our ITS application below is a leading example, though our poor-overlap case also requires moderate to severe extrapolation. Under stationary kernels like the Gaussian, the predicted CEF will “return to the mean” far from the edge of the data. While this may be reasonable, the uncertainty can reach a maximum (determined by $\hat {\sigma }^2$ ) and stop growing. Without addressing this issue, it would appear that the GP does not extrapolate well (see, e.g., Wilson and Adams Reference Wilson and Adams2013). However, when investigators need to extrapolate farther (relative to the scale of the kernel bandwidth), they may need to entertain functions that show (i) periodic behavior (in the case of time series data) and/or (ii) continued growth, e.g., with some continuity in the first or higher derivative, as in linear or polynomial prediction. The GP can accommodate this by combining kernels, for example, adding a Gaussian kernel with a polynomial kernel of the desired order. In these settings, the GP becomes a device for illustrating the wide range of possible functions (and our uncertainty over them). We illustrate this use in Section 3.2 below.

Prior Mean Functions: As noted above, we set $\mu =0$ . This is because prior to seeing any data, we would have no particular reason to prefer a given functional form, direction of change, or any other feature of the CEF not already encoded by choosing the multivariate normal distribution with covariance given by the kernel function. After conditioning, the posterior CEF can approximate whatever smooth functional form the data support.Footnote ¹ This choice does not restrict the posterior fit, which adapts around the observed data as in Expression (4) (see, e.g., Figure 2). Regularization nevertheless occurs through $\sigma ^2$ , which governs the trade-off between fit and smoothness and plays the same role as the regularization parameter in kernel regression (e.g., $\lambda $ in KRLS; Hainmueller and Hazlett Reference Hainmueller and Hazlett2014).

Taking RD as an example, the GP can represent CEFs that are approximately linear, polynomial, or otherwise smooth on each side of the cutoff. If extrapolation far beyond the observed data is needed, prior expectations about functional form can be encoded through a (non-stationary) kernel rather than a mean function. This distinguishes our strategy from approaches in Branson et al. (Reference Branson, Rischard, Bornn and Miratrix2019), Rischard et al. (Reference Rischard, Branson, Miratrix and Bornn2021), and Ornstein and Duck-Mayr (Reference Ornstein and Duck-Mayr2022).Footnote ²

2.2 Comparison to Conventional Uncertainty Estimation under Extrapolation

For tasks involving extrapolation, a key advantage of the GP is its ability to incorporate what we refer to as “extrapolation uncertainty” in its reported posterior uncertainty intervals. Conventional approaches to uncertainty estimation in parametric models involve selecting a single best-fitting model and using that fitted model alone for extrapolation (as well as to compute residuals that inform uncertainty intervals). This is problematic because, within a sufficiently flexible model space, many model fits can explain the observed data similarly well but will make radically different predictions outside of the observed data. Our uncertainty over which of those fits to choose must be maintained if we wish to sustain appropriate uncertainty over the predictions they make when extrapolated. This “extrapolation uncertainty,” however, is lost once we choose a single fitted model, discarding all others. Contrast this with the GP. While the GP does commit to a specific model space and set of assumptions described above, the process of conditioning on the data and inspecting the posterior is different from that of choosing and extrapolating only a single best-fitting model. The GP’s posterior uncertainty band instead reflects the range of CEFs that remain plausible, to varying degrees, given the observed data. The posterior uncertainty at new points effectively extends these functions beyond the data, encompassing their divergent predictions and thereby incorporating increased uncertainty as we move farther from the conditioning data.

We illustrate this in Figure 3. First, consider a linear model (LM) fitted by OLS. The classical variance estimate for the coefficient under homoskedastic and independent (spherical) errors will be $\widehat {\mathbb {V}}(\hat {\beta } \mid {\mathbf {X}},Y) = ({\mathbf {X}}^{\top } {\mathbf {X}})^{-1}\hat {\sigma }^2 I$ . The value of $\hat {\sigma }^2$ is proportional to $\sum _i (Y_i - {\mathbf {X}}_i^{\top }\hat {\beta })^2$ , that is, the sum of the squared fitted residuals.Footnote ³ Consequently, the estimated variance of the predicted value of the CEF at some point ${\mathbf {X}}_j$ will be $\widehat {\mathbb {V}}(\hat {Y}({\mathbf {X}}_j)) = \hat {\mathbb {V}}({\mathbf {X}}_j^{\top }\hat {\beta }) = {\mathbf {X}}_j^{\top } \hat {\mathbb {V}}(\hat {\beta }) {\mathbf {X}}_j$ . While this quantity increases for predictions at ${\mathbf {X}}_j$ points farther from the mean of the data, this is blind to a point’s distance from supporting data. The quadratic LM line and uncertainty envelope in Figure 3 illustrate this in the case of a fitted quadratic model.

Figure 3 Comparison of uncertainty for conditional expectation functions given by GP (blue), KRLS (red), and a quadratic polynomial (gray). Models are fitted on X data between $-$ 3 and 1, then extrapolated over X between 1 and 3. The solid line shows the (randomly drawn) true function, while dashed/dotted lines indicate each model’s estimated CEF. Bands indicate 95% CIs for the CEFs.

This problem is not one of relying on a restrictive model space. We can demonstrate this by comparing GP and KRLS, employing the same kernel (Figure 3). Since these models have identical model spaces for the CEF, their CEFs match very closely, differing only due to hyperparameter tuning. The concern, however, is the uncertainty intervals. KRLS employs the conventional fitted-function variance estimation. This results in a collapse at points farther from the observed data, because the model is structured such that for any one model, we can be sure the CEF it predicts returns to the mean of Y as we leave the support of the data. However, this is the opposite of what we would want, which is to see increased uncertainty about the value of Y as we move further from the support of the data. As Figure 3 shows, the predictive uncertainty from the GP model provides this, widening as we move toward and then beyond the edge of the data.

3.1 Comparing Groups with Poor Overlap

In cross-sectional comparisons, treated and control groups must often be compared as if they had the same covariate distributions. For example, we have an observed confounder X, and assuming the absence of unobserved confounders (i.e., under selection on observables or conditional ignorability), we wish to make comparisons between treated and control units adjusting for X.

Adjustment strategies of all kinds are vulnerable to extrapolation bias when there are regions of X containing only treated or only control units. Here, treated units have values of X between $-$ 3 and 1, whereas control units have values between $-$ 2 and 3. Thus, a model for the treatment outcome must be extrapolated in areas with few or no treated units ( $X> 1$ ) and a model for the control outcome must be extrapolated in areas with few or no control units ( $X < -2)$ . Treatment effect estimates—computed by comparing these two models for each unit (i.e., g-computation, regression imputation, etc.)—consequently suffer in the regions of poor overlap insofar as these model fits fail to predict well in those regions.

For each simulation iteration, we draw a random function to serve as the CEF of Y given X, with a constant treatment effect of 3. To favor non-GP methods, we draw the CEF as a quadratic function $f(x) = \alpha _0 + \alpha _1 x + \alpha _2 x^2$ with the coefficients independently drawn from a standard normal distribution. Following this, “observed” samples are drawn from each of these CEFs ( $N=500$ ) with the addition of noise, distributed normally and with variance calibrated to ensure an overall true $R^2$ of 0.5.Footnote ⁵

Figure 4 presents a single example of a random function drawn from this space (dotted lines), along with the simulated data under treatment (red) and control (blue). The three models tested fit the data similarly well within the region of common support ( $-2 < X < 1$ ). In the regions with poor overlap ( $X< -2$ and $X>1$ ), however, the result is highly sensitive to the choice of modeling approach. Notably, the uncertainty estimates for the LM and BART show insufficient uncertainty to accommodate this model dependency under poor overlap.

Figure 4 Uncertainty quantification of CATE under varying overlap. Top row: Illustration from one draw of the simulation setting, with a good overlap region ( $-2<X<1$ ) and no overlap elsewhere. For treated units (red) and control units (blue), the true CEF (dotted lines) and estimated CEF (black solid lines) with corresponding uncertainty bands are drawn. As shown at both ends of the covariate value, GP allows for growing uncertainty bands adaptive to the degree of overlap. The upper bound of the confidence interval in extrapolation depends on the variability in the fitted data. Bottom row: The CATE estimates with 95% confidence bands. The wider bands of GP in poor-overlap regions propagate to higher uncertainty in CATE estimation. The red dashed line represents the true effect size.

The bottom row of Figure 4 shows what these estimated CEFs and uncertainty intervals imply for conditional average treatment effects (CATEs) and the average treatment effect (ATE). Since the true CATE is constant, the apparently linear change in the treatment effect as a function of X for the LM (left) is erroneous. This is a byproduct of the poor overlap, which led to fitted models with very different slopes for treated versus control groups. While the problem is less pronounced for BART, its uncertainty estimates still fail to reflect the distance from relevant observed data. Its uncertainty over the CATE remains nearly constant across areas of good and poor overlap, producing widely insufficient coverage in the areas of weaker overlap. In contrast, GP’s uncertainty estimates adaptively reflect model dependence due to data sparsity—narrower where common support is better, wider where it is weaker.

Figure 5 Pointwise (top) and average (bottom) performance of GP model compared to LM and BART. Top row: The graph on the left displays the pointwise coverage rate of the ATE estimators by the three models. On the right, the lengths of corresponding 95% confidence intervals are shown over 500 simulations ( $n=500$ ). Bottom row: Boxplots represent the distributions of the ATE estimates (left) and average interval lengths (right). The true treatment effect is denoted by the red dashed line. The coverage rates and RMSEs (of ATE) for each method are also shown.

We repeat this process across 500 iterations, drawing different random functions each time with sample sizes of 500 as described above. Figure 5 illustrates the behaviors of each method across these iterations. The coverage rate for the CATE at each possible value of X (top left) is problematic for the LM at all values of X. BART generally overcovers in areas of common support but radically undercovers in non-overlap regions, dropping below 50% on one side and 25% on the other. The GP achieves approximately nominal coverage where common support is good. In non-overlap areas, it sustains good coverage initially on one side, but still transitions to undercoverage, dropping as low as roughly 75%. While this undercoverage is much less severe than the others, it remains imperfect. This imperfection serves as a warning of the limitations of the GP approach. Specifically, the example is an adversarial one with common support failing for half of the range of X, while the underlying DGPs are quadratic functions that can diverge radically beyond the edges of their training data. In practice, we would recommend GP with a non-stationary kernel when this degree of extrapolation is required (see Section 3.2).

The interval lengths required to achieve this (top right) are also of interest. The drastic undercoverage of LM results from inappropriately short intervals. BART maintains almost constant interval lengths throughout, which are excessively long for the region of good common support. GP exhibits adaptive behavior, maintaining narrow intervals in regions of good overlap and wider intervals in areas of poor overlap.

Finally, while conditional uncertainty and coverage are of potential interest for many queries and demonstrate the behavior of interest here, investigators may be interested primarily in average effects. Even here, however, the poor behavior of LM under weak overlap leads to an RMSE 1.52 times larger than that of GP, together with severe undercoverage (55%).

In Section A.2 of the Supplementary Material, we expand these findings with additional simulations across various function spaces and signal-to-noise ratios ( $R^2$ ).

3.2 Interrupted Time Series with GP

3.2.1 Background

The ITS design is used to assess the effect of events or shocks experienced universally after a specific time point (Bernal, Cummins, and Gasparrini Reference Bernal, Cummins and Gasparrini2017; Box and Jenkins Reference Box and Jenkins1976; Box and Tiao Reference Box and Tiao1975. We observe a time series of non-treatment outcomes, $Y_t(0)$ prior to the event of interest. The event occurs at time $t=T$ and all $Y_t$ outcomes measured after that time are treatment outcomes, $Y_t(1)$ . We train a model on pre-treatment data then must extrapolate it to the post-treatment era. These extrapolated predictions of $Y_t(0)$ are then compared to actual post-treatment outcomes ( $Y_t(1)$ ) at any $t \geq T$ .

The GP approach is potentially valuable here for several reasons. First, it handles autocorrelation in the outcomes naturally. Second, the covariance function (kernel) can be designed to accommodate not only the “smoothness” we expect but also secular trends that will continue on over time, and periodic trends such as seasonality. Third and most relevant to our discussion, its estimates of uncertainty take into account the distance (in time) to the observed data and the consequent model dependency. We note that regardless of estimation approach, making causal claims from ITS requires additional demanding identification assumptions (see Section A.4 of the Supplementary Material).

Figure 6 Left column: Monthly rate of handgun background checks per 100,000 population in D.C. (top) and Vermont (bottom), seven years before and one year after the Heller ruling. The red envelope shows the 95% predictive interval for the GP-estimated non-treatment outcome, post-treatment. Right column: Point estimates and confidence band for the differences between the observed outcome (treated) and the predicted (non-treatment) counterfactual at each month. The kernel choice combines (adds) the linear, periodic, and Gaussian kernels.

3.3 Regression Discontinuity

The RD design is a widely used tool in the social sciences for cases when decisions about treatment turn on a sharp cutoff in some variable. While this approach rests on relatively credible identification assumptions, it poses estimation challenges as it relies on point estimates with uncertainty intervals precisely where the underlying model(s) are at the very edge of their support.

When fitting models below/above the cutoff in RD, decisions must be made regarding the bias-variance tradeoff: how much should we seek to reduce variance by using observations farther from the cutoff, given that including such data can bias our estimate? Even when optimizing this choice by some criterion, the consequent bias must be dealt with when estimating uncertainty. Stommes, Aronow, and Sävje (Reference Stommes, Aronow and Sävje2023) review varied approaches to dealing with this bias and its implications for inference. Among such approaches, Calonico et al. (Reference Calonico, Cattaneo and Titiunik2014) provide a now standard optimal bandwidth selection technique, bias-adjustment of the point estimate, and bias-adjusted inference, which we refer to by the name of the widely used software implementation, rdrobust (Calonico, Cattaneo, and Titiunik Reference Calonico, Cattaneo and Titiunik2015; Calonico et al. Reference Calonico, Cattaneo, Farrell and Titiunik2017).

The GP approach is attractive for RD, first, because it offers a flexible function space that makes only weak assumptions about smoothness (when using the Gaussian kernel, for example). Second, and more importantly, it provides a principled approach to uncertainty quantification for edge prediction or mild extrapolation. In cases with considerable data near the threshold, an optimally chosen parametric or semi-parametric approach, like those employed by rdrobust, can establish reliable point estimates and uncertainty intervals at the edge of the data, as demonstrated by Calonico et al. (Reference Calonico, Cattaneo and Titiunik2014). However, with smaller samples or sparse data near the cutoff, multiple models may fit the observed data similarly well while yielding different predictions at the boundary. The GP approach provides uncertainty estimates that retain this uncertainty over the CEF, allowing inference without relying on large-sample asymptotics.

3.3.1 RD Case 1: Total Random Simulation

We first consider the simple “total random case,” in which the running variable (X) and the outcome (Y) are both drawn independently at random, with zero treatment effect.Footnote ⁷ For each of 500 simulations, we draw a sample of size 500, with $X_i \sim \mathcal {N}(0, 1)$ , and $Y_i \sim \mathcal {N}(0, \sigma ^2)$ . The value of $\sigma $ is varied from near zero up to 3, reflecting various possible scales of the outcome variable relative to the forcing variable. We use a cutoff of 0 for X to determine treatment eligibility status. We employ rdrobust with software defaults and utilize the bias-adjusted point estimates with the “robust” confidence intervals. To compute the GP estimate (GPrd), we first use the GP function from gpss to estimate a model for Y given X for the treated units ( $X> 0$ ), and another for the control units ( $X \leq 0$ ). From this, we obtain the point estimate and standard error for the prediction of the expected Y at precisely $X=0$ from each model. The treatment effect is estimated as the difference between these point estimates, with standard error given by the square root of the sum of the variances of the two predicted outcomes.

Figure 7 GP results with three different kernels for Vermont. The post-treatment period is extended to four years to better illustrate the behavior under extreme extrapolation. Left: Gaussian kernel; Center: “Gaussian + periodic + linear” kernel; Right: “Gaussian + periodic + quadratic” kernel.

Figure 8 The coverage rate, average length of 95% confidence interval, and RMSE of GPrd (red) and rdrobust with “robust” options (black) in the total random setting. The ratio of the variance of Y to X is given by $\sigma $ on the horizontal axis. Top row: True effect size is zero, 500 observations. Bottom row: Effect size is zero, sample size is reduced to 100. Results with a true effect size of three look identical.

Figure 8 shows simulation results in this setting. Looking first at the top row, we see that GPrd has a higher coverage rate across all levels of $\sigma $ , at 97%–98%, while rdrobust has a coverage rate between 93% and 96%. The two approaches have essentially identical interval sizes, which appropriately scale linearly with $\sigma $ (middle column). GPrd also shows much smaller RMSE across the 500 simulations. This RMSE benefit is amplified in smaller samples (bottom row). We examine and offer an explanation for this behavior in Section 3.3.3.

3.3.2 RD Case 2: Latent Variable Confounding Simulation

The total random case is a useful starting point, but may not produce sufficient risk of bias because there is no relationship between X and $Y(d)$ on either side of the cutoff such that incorporating information farther from the cutoff risks adding bias to the estimate. To better simulate this concern, we consider a latent confounding formulation. In each of 500 replications, we generate a random sample of size 200 according to the following:

1. 1. Draw a latent variable: $\mu \sim \mathcal {U}(-0.25, 0.25)$ .

2. 2. Create our running variable, X, loosely resembling vote share in an electoral RD as a sigmoidal function of $\mu $ , with $x_i = 1/(1 + exp(-s*\mu _i))$ . The parameter s represents a varying degree of steepness of this sigmoid, such that larger values produce a steeper sigmoid, increasing the risk of bias due to misspecification (see Figure 9).

3. 3. The treatment D represents “winning the election,” coded as 1 if the vote share is greater than 0.5, and 0 otherwise.

4. 4. Generate the outcome $Y_i = 1.5 \mu _i + \tau D + \varepsilon ,$ where $\tau =3$ and $\varepsilon \sim \mathcal {N}(0, \sigma ^2)$ . The $\sigma $ will be varied by simulation setting to explore a range of realistic signal-to-noise ratios ( $R^2$ ).

Figure 9 Relationship between simulated latent variable $\mu $ (horizontal axis) and vote share X (vertical axis) at two steepness (s) parameters used in the simulation.

Next, purely data-driven RD specification risks using distant observations to estimate $\mathbb {E}[Y|X]$ at the cutoff, compromising the motivating logic for using RD. To address this concern, we consider options that set the bandwidth or trim the sample (or both) based on outside causal assumptions:

• • gp_causal: The investigator proposes a kernel bandwidth based on what they are willing to believe about how close units may covary in their outcomes. For example, it might be reasonable to propose that two observations whose vote shares are 2% apart may be highly correlated, e.g., at 0.9. This implies a surprisingly small value for b at 0.005.Footnote ⁸

• • gp_causaltrim: After choosing the causaltrim bandwidth, regions that are effectively irrelevant to estimation at the cutoff can be removed. This aids in transparency, but is also helpful because the value of $\sigma ^2$ is a single global parameter, so trimming the data first ensures it is determined over only the range of the relevant range of the sample. In our simulation, we remove observations with vote shares below 0.4 or above 0.6, on the premise that users would likely argue that data outside this range are of no value to modeling what happens at the cutoff.

Figure 10 RD simulation results with latent variable (effect size = 3). Each point represents the coverage rate, or average length of 95% confidence interval of GPrd, gp_causal, gp_causaltrim, and rdrobust across 500 iterations of the latent variable simulations ( $n = 200$ ). The horizontal axis represents the $R^2$ .

Figure 10 summarizes the results over 500 iterations at each choice for the steepness (s) of the latent relationship between $\mu $ and Y and at four different levels of $R^2$ . All the approaches show similar interval sizes, with the exception of rdrobust, which yields wider intervals when steepness equals 20. While rdrobust tends to slightly undercover the true effect, the GP-based approaches exhibit slight overcoverage, except for gp_causaltrim at a steepness of 20.

3.3.3 RD Case 3: Empirical Application with Benchmarking by Placebo Cutoffs

Finally, we study the performance of GPrd on real data, in a case where the “true” effect is knowable for benchmarking purposes. To do so, we consider the data from Lee, Moretti, and Butler (Reference Lee, Moretti and Butler2004), in which the forcing variable is Democrats’ (two-party) vote share in U.S. House Elections (1948–1990), and the outcome of interest is a measure of how liberal each representative’s vote record is assessed to be, called the ADA score. While we cannot know the true effect size, it would be arguably zero when using any value other than 0.5 as a (placebo) cutoff.Footnote ⁹

We consider three estimators: GPrd, gp_causaltrim, and rdrobust. Figure 11 shows results with the full data. The choice of estimator matters little to the point estimates at the true cutoff (0.5), though GPrd and gp_causaltrim are more conservative in their uncertainty estimates. Looking at the eight placebo cutoffs, the differences are relatively small and while rdrobust suggests a statistically significant estimate at the cutoff of 0.35, GPrd and gp_causaltrim very nearly do as well. At the cutoff of 0.65, the rdrobust interval again excludes zero, while the intervals for the GP approaches do not, but the actual difference is again relatively small.

Figure 11 RDD estimates using close election data with placebo cutoff points and different GP RD Specifications. The 95% confidence interval of GPrd (red), gp_causaltrim (blue), and rdrobust (black) across eight different placebo cutoffs and the true cutoff point at 0.5 are presented.

However, the Lee et al. (Reference Lee, Moretti and Butler2004) data initially contain 13,577 observations without missingness on the running variable and outcome. While each of the placebo cutoff analyses uses only a subset of these points, there are still thousands of observations available to each of these estimates. The similarity of the GP and rdrobust approaches in this context fits the prediction that they should be similar with ample data. We also wish to know whether GP provides a suitable approach to RD for smaller sample sizes. We thus use the same data as the basis for a second analysis in which we examine much smaller sub-samples. Specifically, we:

• • fix a (placebo) cutoff point: 0.3, 0.4, 0.5, 0.6, 0.7 (0.5 is the true cutoff);

• • limit the data to just 0.1 above and below the cutoff in question (e.g., for cut = 0.4, data with $0.3<X<0.5$ are used) to maintain symmetry and comparability across cutoff estimates;

• • sub-sample 200 observations (without replacement) from the remaining eligible observations;

• • estimate the treatment effect and its uncertainty using various models.

Figure 12 RDD placebo cutoff subsampling simulation. Boxplots show the distribution of treatment effect estimates from GPrd (red), gp_causaltrim (blue), and rdrobust (black) across four different placebos and true cutoff point at 0.5.

We repeat this 500 times. Figure 12 shows the results. For coverage and RMSE calculations, we use the “full sample” estimates above as if it were the long-run target to reveal the change in behavior due solely to sample size. In general, the GP approaches are similar and perform well. Coverage rates vary from 96% to 99% across the GP methods, while the gp_causaltrim approach shows greater efficiency. As this approach makes stricter (and arguably transparency- and credibility-improving) choices, it is fortunate to also find it has desirable behavior. The rdrobust approach is not designed for the small sample setting, but we include it here to illustrate its sensitivity to sample size for those who might use it in smaller sample settings regardless. It shows occasional erratic estimates (omitted from Figure 12 with RMSE values roughly twice those of GPrd.Footnote ¹⁰