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.

Figure 3 Long description

A line graph with X on the horizontal axis ranging from minus 3 to 3 and Outcome on the vertical axis ranging from minus 6 to 2. A scatter plot of gray dots is concentrated between X equals minus 3 and 1. A solid dark purple line represents the true function, showing a wave-like oscillation.

Between X equals minus 3 and 1 (the training data region):

* The Quadratic L M (dotted gray line) shows a slight downward curve with a narrow gray confidence band.

* The G P (dashed blue line) and K R L S (dashed red line) closely follow the true function's oscillations with narrow, overlapping purple and pink confidence bands.

Between X equals 1 and 3 (the extrapolation region):

* The true function continues to oscillate.

* The Quadratic L M continues its downward slope, with its gray confidence band widening significantly.

* The G P and K R L S models deviate from the true function. The G P blue dashed line and the K R L S red dashed line both trend upward before leveling off.

* The uncertainty bands for G P (light blue) and K R L S (light red) expand dramatically in this region, with the G P band becoming much larger than the others, covering almost the entire vertical range from minus 4 to 2 at X equals 3.

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.

Figure 4 Long description

The figure consists of two rows and three columns. The x-axis for all panels is Covariate, ranging from negative 3 to 3.

Top Row: Outcome vs. Covariate.

- L M (Left): Shows linear trends for Treated (red dots) and Control (blue diamonds). Two intersecting black solid lines represent estimated C E F with narrow gray uncertainty bands. Dotted lines show true C E F.

- G P (Middle): Shows non-linear, U-shaped curves. The gray uncertainty bands widen significantly at the far left and right edges where data overlap is poor.

- B A R T (Right): Displays highly jagged, step-like estimation lines. Uncertainty bands are relatively constant but wider than L M.

Bottom Row: C A T E estimates vs. Covariate.

- A horizontal red dashed line at y equals 3 represents the True Effect.

- L M (Left): A downward sloping straight line. A T E is 1.62 with interval 0.73 to 2.5.

- G P (Middle): A smooth, undulating curve that stays near the true effect in the center but drops and gains massive uncertainty bands at the edges. A T E is 2.49 with interval 0.23 to 4.75.

- B A R T (Right): A jagged, descending line with moderate uncertainty bands throughout. A T E is 2.28 with interval 0.35 to 4.21.

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.

Figure 5 Long description

Top-left line graph: X-axis is Covariate from minus 3 to 3. Y-axis is Pointwise Coverage from 0.00 to 1.00. A horizontal dashed line marks 0.95. L M (solid black) shows a sharp peak at 0.95 near covariate 0 but drops to near 0.10 elsewhere. G P (blue dotted) maintains high coverage near 1.00 across most of the range, dropping at the extremes. BART (green dashed) shows high coverage between minus 2 and 1.5, with steep declines outside that range.

Top-right line graph: X-axis is Covariate. Y-axis is Pointwise Interval Length from 4 to 16. L M is the lowest, forming a shallow U-shape. G P is a wide U-shape, very low between minus 2 and 1 but rising sharply to 16 at the edges. BART is relatively flat around 5 to 6.

Bottom-left boxplot: X-axis is Method (L M, G P, BART). Y-axis is A T E estimates from 0.0 to 7.5. A red dashed line at 3.0 indicates the True Effect. L M has the widest distribution and many outliers. G P and BART are more tightly centered on the red line. Below the plot, a table lists: L M cover 0.554, R M S E 1.317; G P cover 1, R M S E 0.866; BART cover 1, R M S E 0.821.

Bottom-right boxplot: X-axis is Method. Y-axis is Estimated Interval Length (A T E) from 0 to 12. L M has the shortest intervals (median near 2). G P and BART have much larger intervals (median near 5) with several high-value outliers reaching up to 15.

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.

Figure 6 Long description

The figure consists of four panels arranged in a two-by-two grid.

Top-Left Panel: District of Columbia time-series. The x-axis is Month from January 2002 to July 2008. The y-axis is Handgun Background Checks per 100,000. Observed data points (black dots) remain near zero until July 2008, where they spike sharply to nearly 10.0. A vertical dashed line marks 2008-07-01. A red dashed line with a red shaded 95% predictive interval shows the counterfactual prediction staying near zero.

Top-Right Panel: District of Columbia Pointwise Effect. The x-axis is Month after Heller (0 to 11). The y-axis is Pointwise Effect. The line starts at zero, peaks at approximately 9.5 at month 4, and then fluctuates between 5.0 and 7.5 for the remainder of the year. A grey confidence band surrounds the line.

Bottom-Left Panel: Vermont time-series. The x-axis is Month from January 2002 to July 2008. The y-axis is Handgun Background Checks per 100,000, ranging from 40 to 160. The data shows high seasonality with peaks and troughs. After the July 2008 vertical dashed line, the observed data (black dots) and counterfactual (red dashed line) both show an upward trend, with the red shaded interval overlapping many observed points.

Bottom-Right Panel: Vermont Pointwise Effect. The x-axis is Month after Heller (0 to 11). The y-axis is Pointwise Effect ranging from -25 to 75. The line fluctuates significantly, starting near 20, dropping below zero at month 2, peaking near 50 at month 7, and ending near 5. The grey confidence band is much wider than the D.C. panel, often crossing the zero baseline.

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 7 Long description

A three-panel line graph showing Handgun Background Checks on the Y-axis (0 to 300) against years on the X-axis (2002 to 2012). A vertical dashed line at late 2008 marks the treatment point. Black dots represent observed data points across all panels.

* Left Panel (Gaussian kernel): Shows a smooth, slightly undulating grey trend line with a narrow confidence interval pre-2008. Post-2008, a red dashed line extrapolates a slight downward curve, while the red shaded confidence interval remains wide and rectangular. The actual data points trend sharply upward, far above the prediction.

* Center Panel (Gaussian plus periodic plus linear kernel): The pre-2008 grey trend line follows the seasonal oscillations of the data. Post-2008, the red dashed line continues this periodic oscillation with a slight linear increase. The red shaded confidence interval is jagged, reflecting the periodic nature, but the actual data points still exceed the predicted trend significantly after 2010.

* Right Panel (Gaussian plus periodic plus quadratic kernel): Similar to the center panel, the pre-2008 trend captures seasonality. Post-2008, the red dashed line and red shaded confidence interval show an accelerating quadratic upward curve. This model's prediction and confidence interval better encompass the rising data points compared to the other two kernels.

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 Long description

The figure consists of two rows of three line graphs. All graphs share a horizontal axis labeled with the Greek letter sigma, ranging from 1.5 to 3.0. A legend indicates red lines for G P r d and black lines for R D robust.

Top row (a) True effect of zero, n = 500:

* Left graph: Coverage Rate on the vertical axis from 0.5 to 1.0. A dashed horizontal line sits at 0.95. The red line stays slightly above 0.95, while the black line stays slightly below it.

* Middle graph: Avg Interval Length on the vertical axis from 2.0 to 3.5. Both red and black lines show a near-identical linear increase from bottom-left to top-right.

* Right graph: R M S E on the vertical axis from 0.5 to 0.9. Both lines show a stepped increase, with the black line consistently higher than the red line.

Bottom row (b) True effect of zero, n = 100:

* Left graph: Coverage Rate on the vertical axis. The red line remains near the 0.95 dashed line, while the black line is lower, around 0.90.

* Middle graph: Avg Interval Length on the vertical axis from 4 to 8. Both lines show a linear increase, but the black line is significantly higher than the red line.

* Right graph: R M S E on the vertical axis from 1.0 to 2.5. Both lines increase with sigma, but the black line has a much steeper slope and higher values than the red line.

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.

Figure 9 Long description

The image consists of two side-by-side line graphs with a shared vertical axis labeled Vote Share ranging from 0.00 to 1.00 and a shared horizontal axis labeled Latent Variable ranging from minus 0.2 to 0.2.

* Left Panel: Titled steepness equals 5. A black line shows a gradual, nearly linear positive increase. It starts at a Vote Share of approximately 0.22 when the Latent Variable is minus 0.25 and rises to approximately 0.78 when the Latent Variable is 0.25.

* Right Panel: Titled steepness equals 20. A black line shows a sharp S-shaped sigmoidal curve. The curve begins near 0.00 at a Latent Variable of minus 0.25, remains flat until minus 0.1, then rises steeply through the center point of 0.50 at a Latent Variable of 0.0, and plateaus near 1.00 as the Latent Variable reaches 0.25.

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 Long description

The grid consists of two columns for steepness equals 5 and steepness equals 20, and two rows for Coverage Rate and Interval Size. The horizontal axis for all panels is R super 2 with values 0.1, 0.3, 0.5, and 0.7. A legend at the bottom identifies four methods. G P r d as dark red circles, gp causal as red triangles, gp causaltrim as orange squares, and rd robust as black plus signs.

* Top Row: Coverage Rate. The vertical axis ranges from 0.80 to 1.00 with a dotted horizontal reference line at 0.95. In both steepness panels, G P r d, gp causal, and gp causaltrim maintain high coverage near or above 0.95. The rd robust method consistently shows lower coverage, fluctuating between 0.90 and 0.95.

* Bottom Row: Interval Size. The vertical axis ranges from 0 to 30. In both steepness panels, all methods show a non-linear decrease in interval size as R super 2 increases. The rd robust method consistently has the largest interval size, appearing as the highest point in each cluster, while G P r d, gp causal, and gp causaltrim are tightly grouped at lower values.

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.

Figure 11 Long description

The x-axis is labeled Cutoff Points (True: 0.5) with values from 0.3 to 0.7. The y-axis is labeled Treatment Effect estimates with values from minus 20 to 40. A dashed horizontal line indicates the zero baseline.

At each cutoff point, three vertical confidence interval bars are clustered together: red (G P rd) on the left, blue (gp causaltrim) in the middle, and black (rdrobust) on the right.

* At the true cutoff of 0.5, all three estimates show a significant positive treatment effect, with point estimates clustered near 45 and confidence intervals well above zero.

* At placebo cutoffs (0.3, 0.35, 0.4, 0.45, 0.55, 0.6, 0.65, and 0.7), the estimates are concentrated near the zero baseline.

* Most placebo confidence intervals cross the zero line, indicating no significant effect, though some variation exists at 0.65 where estimates are slightly positive and at 0.3 where they are slightly negative.

* The red G P rd bars generally show wider confidence intervals compared to the blue and black bars across the placebo points.

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.

Figure 12 Long description

A statistical boxplot chart with a horizontal x-axis labeled Cutoff Points True X equals 0.5 and a vertical y-axis labeled Treatment Effect Estimates for Cutoff Points ranging from negative 100 to 100. A dashed horizontal line sits at zero. A legend at the top left identifies three models: G P r d in red, G P r d causaltrim in blue, and R D robust in black.

Data is grouped into five clusters along the x-axis at points 0.3, 0.4, 0.5, 0.6, and 0.7. Each cluster contains three boxplots corresponding to the models.

* At 0.3 and 0.4: All three models are centered near zero with narrow distributions.

* At 0.5 (the true cutoff): All models show a positive shift, with medians near 50.

* At 0.6 and 0.7: Medians return toward zero, but R D robust (black) shows significantly higher variance and more outliers compared to the G P r d variants.

Below the plot is a table aligned with the x-axis clusters showing cover and r m s e values for each model:

* At 0.3: G P r d (0.99, 8.79), G P r d causaltrim (0.98, 6.1), R D robust (0.92, 18.29).

* At 0.4: G P r d (0.99, 8.01), G P r d causaltrim (0.98, 5.86), R D robust (0.94, 15.09).

* At 0.5: G P r d (0.99, 8.94), G P r d causaltrim (0.98, 6.42), R D robust (0.9, 16.92).

* At 0.6: G P r d (0.98, 11.06), G P r d causaltrim (0.97, 8.22), R D robust (0.92, 20.85).

* At 0.7: G P r d (0.96, 12.5), G P r d causaltrim (0.97, 8.05), R D robust (0.85, 23.58).

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 ¹⁰