Hierarchical Gaussian Process Models for Regression Discontinuity/Kink under Sharp and Fuzzy Designs

The regression discontinuity design (Thistlethwaite and Campbell, 1960) and its generalizations are non-experimental methods of causal inference in observational studies. In a typical regression discontinuity design, the assignment of a treatment is determined according to whether a running variable is greater or less than a known cutoff while the outcome of interest is supposed to be a smooth function of the running variable except at the cutoff. The abrupt change in the probability of receiving the treatment can be exploited to infer the local causal effect around the cutoff. For general overviews, see Imbens and Lemieux (2008), Lee and Lemieux (2010), Cattaneo and Escanciano (2017), Cattaneo et al. (2020, forthcoming), and Cattaneo and Titiunik (2022).

We develop in this study a family of estimators for Regression Discontinuity (RD) and Regression Kink (RK) under sharp and fuzzy designs. The proposed estimators are based on Gaussian Process (GP) regression/classification, a probabilistic modeling approach that has gained popularity in machine learning and artificial intelligence community (Rasmussen and Williams, 2006). This nonparametric Bayesian approach places a GP prior on the latent function that underlays the outcome of interest. Modeler’s knowledge regarding the underlying relationship is encoded through the mean and variance functions of the GP and inference is conducted according to Bayes’ rule.

Unlike methods that focus on the conditional mean of the outcome immediately below and above the cutoff (e.g., the local polynomial estimators), the proposed GP methods estimate the underlying functional relationship and rely on the resultant predictive distributions at the cutoff to infer various treatment effects of interest. One key advantage GP methods enjoy over other machine learning methods (such as neural network and random forest) is its automatic Uncertainty Quantification (UQ). Probabilistic quantities such as credibility intervals are readily obtained from the predictive distributions. Another appeal of GP methods is that its gradient remains a Gaussian process, facilitating derivative estimation. This is particularly valuable for the regression kink estimation considered in this study. Furthermore, thanks to the likelihood principle inherent in Bayesian inference, the proposed method infers both RD and RK effects from a common generative model and does not require separate models/configurations for the RD and RK estimation.

The fuzzy RD/RK design arises when treatment assignment differs from actual treatment take-up. The identification of the RD/RK effects under the fuzzy design requires accounting for the change in probability of treatment take-up at the cutoff. Estimating this probability is a classification problem in which the GP method excels. The GP classification is a generalization of GP regression; it models a continuous latent function underlaying the observed discrete outcome. We develop a GP RD estimator for binary outcomes and use it in conjunction with the GP regression to formulate estimators for fuzzy RD/RK effects.

The essence of GP regressions is the covariance function that governs the sample path of a Gaussian process. We use the squared exponential covariance function in our GP regressions. It can be shown to be equivalent to a basis function expansion with an infinite number of radial basis functions and therefore rather expressive. It is, however, stationary such that the covariance between two inputs depends only on their distance. Consequently it may be inadequate if the underlying relationship is nonstationary. To tackle this potential limitation, we further develop a multi-layer hierarchical GP estimator. The first layer utilizes a Bayesian neural network that nonlinearly transforms the inputs into a latent feature space. The second layer utilizes the GP to map the intermediate outputs to the observed responses. This estimator can be interpreted as a deep learning model with a multi-layer architecture. It combines the strength of neural network and GP regressions. Moreover using the GP estimator as the penultimate layer is amenable to automatic uncertainty quantification and derivative estimation, facilitating the estimation and inference of RD/RK effects.

We conduct Monte Carlo simulations to assess the performance of the proposed methods. The results suggest that our estimators provide comparable and sometimes noticeably better performance relative to competing estimators. The hierarchical GP estimators improve upon the one-layer GP estimators substantially, demonstrating the merit of deep learning approach. These overall patterns hold for the RD and RK estimations under both sharp and fuzzy designs, and across different data generating processes and sample sizes. For illustration, we apply the proposed methods to examine the electoral advantage to incumbency in the US House election (Lee, 2008). Lastly we present a further extension to accommodate covariate adjustment.

In this section, we introduce GP estimators for RD/RK treatment effect under sharp and fuzzy designs. As mentioned above, a GP is characterized by its mean and covariance functions. Usually the zero mean function suffices. However when extrapolating out of sample, predictions of GP regression gravitate towards the mean function. Incorporating an explicit mean function may mitigate the ‘regressing to zero’ tendency associated with a zero mean function in out of sample prediction. This is particularly desirable for RD estimation of local treatment effect at the cutoff, especially if the cutoff is not covered by the sample range.

A common choice of mean function is a polynomial $m(x)=\bm{h}^{T}(x)\beta$, where $\bm{h}(x)=[1,x,\dots,x^{S}]^{T}$ and $\beta=[\beta_{0},\beta_{1},\dots,\beta_{S}]^{T}$. Within the Bayesian paradigm, we can absorb this mean function into a general composite covariance. With a Gaussian prior $\beta\sim\mathcal{N}(0,\Sigma)$, $m(x)\sim\mathcal{N}(0,\bm{h}^{T}(x)\Sigma\bm{h}(x))$. Define

$$k_{\textsf{p}}(x,x^{\prime}|\Sigma)=\bm{h}^{T}(x)\Sigma\bm{h}(x^{\prime})$$

This function turns out to be the so-called polynomial covariance function. One important appeal of GP regression is its expressiveness. One can not only choose from a large collection of covariance functions with varying features, but also construct new ones from the product and/or sum of covariances. In this study, we consider an additively composite covariance that is the sum of a polynomial covariance and a squared exponential covariance:

$$k_{\textsf{p+s}}(x,x^{\prime}|\Sigma,l,\alpha^{2})=k_{\textsf{p}}(x,x^{\prime}|\Sigma)+k_{\textsf{s}}(x,x^{\prime}|l,\alpha^{2})$$

(7)

We can now proceed to RD estimations. We partition the sample into two subsets, denoting the subsample with $x<c$ by $\mathcal{D}_{0}=(\bm{y}_{0},\bm{x}_{0})$ and the rest by $\mathcal{D}_{1}=(\bm{y}_{1},\bm{x}_{1})$, with respective subsample sizes $N_{0}$ and $N_{1}$. We first consider the sharp RD design. For $j=0,1$, we assume that

$$y_{i,j}=f_{j}(x_{i,j})+u_{i,j},i=1,\dots,N_{j}$$

where $u_{i,j}\sim\mathcal{N}(0,\sigma_{j}^{2})$. The latent function $f_{j}$ has a GP prior with zero mean and composite covariance given by (7). We adopt the fully Bayesian approach for inference with independent prior distributions for $\theta_{j}=(l_{j},\alpha_{j}^{2},\Sigma_{j},\sigma_{j}^{2})$.

Denote by $\tilde{\theta}_{n,j},n=1,\dots,\tilde{N}$, an MCMC sample of size $\tilde{N}$ drawn from the posterior distribution $p(\theta_{j}|\mathcal{D}_{j})$. Define $\bm{K}_{x,x^{\prime}}(\theta)=k_{\textsf{p+s}}(x,x^{\prime}|\theta)$. The predictive mean, conditional on $\tilde{\theta}_{n,j}$, at the cutoff is given by

$$\tilde{m}_{n,j}=\bm{K}_{c,\bm{x}_{j}}(\tilde{\theta}_{n,j})(\bm{K}_{\bm{x}_{j},\bm{x}_{j}}(\tilde{\theta}_{n,j})+\tilde{\sigma}_{n,j}^{2}\mathcal{I}_{N_{j}})^{-1}\bm{y}_{j}$$

(8)

with conditional variance

$$\tilde{v}_{n,j}=K_{c,c}(\tilde{\theta}_{n,j})-{\bm{K}}_{c,\bm{x}_{j}}(\tilde{\theta}_{n,j})(\bm{K}_{\bm{x}_{j},\bm{x}_{j}}(\tilde{\theta}_{n,j})+\tilde{\sigma}^{2}_{n,j}\mathcal{I}_{N_{j}})^{-1}{\bm{K}}_{\bm{x}_{j},c}(\tilde{\theta}_{n,j})$$

The predictive mean and variance are then calculated as

$$\hat{\mu}_{j}=\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}\tilde{m}_{n,j},\quad\hat{v}_{j}=\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}(\tilde{m}_{n,j}-\hat{\mu}_{j})^{2}+\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}\tilde{v}_{n,j}$$

where the first term of $\hat{v}_{j}$ is the variance of conditional predictive mean and the second is the mean of conditional predictive variance. Finally the RD treatment effect is estimated by

$$\hat{\tau}_{\textsf{SRD}}=\hat{\mu}_{1}-\hat{\mu}_{0}$$

with posterior variance

$$\hat{v}_{\textsf{SRD}}=\hat{v}_{1}+\hat{v}_{0}$$

Next we consider the estimation of RK effect. Given a binary treatment, the denominator of the RK estimator (2) equals one. Define $\dot{K}_{x,x^{\prime}}(\theta)=\dot{k}_{\textsf{p+s}}(x,x^{\prime}|\theta)$ and $\ddot{K}_{x,x^{\prime}}(\theta)=\ddot{k}_{\textsf{p+s}}(x,x^{\prime}|\theta)$. The predictive derivative, conditional on $\tilde{\theta}_{n,j}$, at the cutoff is given by

$$\tilde{\dot{m}}_{n,j}=\dot{\bm{K}}_{c,\bm{x}_{j}}(\tilde{\theta}_{n,j})(\bm{K}_{\bm{x}_{j},\bm{x}_{j}}(\tilde{\theta}_{n,j})+\tilde{\sigma}_{n,j}^{2}\mathcal{I}_{N_{j}})^{-1}\bm{y}_{j}$$

with conditional variance

$$\tilde{\dot{v}}_{n,j}=\ddot{K}_{c,c}(\tilde{\theta}_{n,j})-\dot{\bm{K}}_{c,\bm{x}_{j}}(\tilde{\theta}_{n,j})(\bm{K}_{\bm{x}_{j},\bm{x}_{j}}(\tilde{\theta}_{n,j})+\tilde{\sigma}^{2}_{n,j}\mathcal{I}_{N_{j}})^{-1}\dot{\bm{K}}_{\bm{x}_{j},c}(\tilde{\theta}_{n,j})$$

Formulae for $\dot{\bm{K}}_{c,\bm{x}_{j}}$ and $\ddot{\bm{K}}_{c,c}$ are provided in Appendix A. The predictive mean and variance of the derivative at the cutoff are calculated as

$$\hat{\dot{\mu}}_{j}=\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}\tilde{\dot{m}}_{n,j},\quad\hat{\dot{v}}_{j}=\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}(\tilde{\dot{m}}_{n,j}-\hat{\dot{\mu}}_{j})^{2}+\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}\tilde{\dot{v}}_{n,j}$$

The RK treatment effect is then estimated by

$$\hat{\tau}_{\textsf{SRK}}=\hat{\dot{\mu}}_{1}-\hat{\dot{\mu}}_{0}$$

with posterior variance

$$\hat{v}_{\textsf{SRK}}=\hat{\dot{v}}_{1}+\hat{\dot{v}}_{0}$$

Under a fuzzy RD/RK design, the actual treatment take-up may differ from the treatment assignment. Hahn et al. (2001) show that the fuzzy RD design is closely connected to an instrumental variable problem with a discrete instrument, and that the treatment effect of interest is a version of the local average treatment effect or complier average treatment effect. As in the previous section, we start with the RD effect. The estimation under fuzzy RD design entails estimating the ratio of RD effects on the response and treatment take-up probability respectively. The former is given by the sharp RD estimator of the previous section, while the latter is based on the GP regression for binary responses, or GP classification.

Under the fuzzy RD design, the probability of treatment take-up $p(D=1)$ is assumed to be a smooth function of the running variable, except for a jump at the cutoff. Following the convention of the machine learning literature, we set $D=-1$ for no treatment take-up and $D=1$ otherwise. Similarly to the SRD estimator, we model $p(D=1)$ separately for the two subsamples $\mathcal{D}_{0}=(\bm{D}_{0},\bm{x}_{0})$ and $\mathcal{D}_{1}=(\bm{D}_{1},\bm{x}_{1})$ defined according to whether $x\geq c$. For $j=0,1$ we assume that

$$p(D_{i,j}=1|x_{i,j})=\psi(\gamma_{j}+D_{i,j}f_{j}(x_{i,j}))$$

where $\psi$ is a sigmoid function such as the logistic or probit function, $\gamma_{j}$ is an offset coefficient, and $f_{j}$ is a latent function with a GP prior.

The conditional marginal likelihood is given by

	$\displaystyle\log p(\mathcal{D}_{j}|\theta_{j})=$	$\displaystyle\sum_{i=1}^{N_{j}}\log\psi(\gamma_{j}+D_{i,j}f_{j}(x_{i,j}))-\frac{1}{2}f_{j}(\bm{x}_{j})^{T}\bm{K}_{\bm{x}_{j},\bm{x}_{j}}^{-1}f_{j}(\bm{x}_{j})$
		$\displaystyle-\frac{1}{2}\log|\bm{K}_{\bm{x}_{j},\bm{x}_{j}}|-\frac{N_{j}}{2}\log 2\pi$

The hyperparameter $\theta_{j}$ includes $\gamma_{j}$ for the observation model and those for the covariance of the latent function. For notational simplicity, the dependence of the covariance and latent function on the hyperparameters is suppressed.

Since the marginal likelihood $p(D_{j})=\int p(D_{j}|\theta_{j})p(\theta_{j})d\theta_{j}$ is not tractable, we use MCMC method for inference. Let $\tilde{\theta}_{n,j},n=1,\dots,\tilde{N}$, be an MCMC sample of hyperparameters from the posterior distribution $p(\theta_{j}|\mathcal{D}_{j})$. For each $\tilde{\theta}_{n,j}$, we generate an $N_{j}$-dimensional vector of latent values according to

$$\tilde{\bm{f}}_{n,j}\sim\mathcal{N}(0,\bm{K}_{\bm{x}_{j},\bm{x}_{j}}(\tilde{\theta}_{n,j}))$$

The predictive latent value at the cutoff is then given by

$$\tilde{q}_{n,j}=\bm{K}_{c,\bm{x}_{j}}(\tilde{\theta}_{n,j})\bm{K}^{-1}_{\bm{x}_{j},\bm{x}_{j}}(\tilde{\theta}_{n,j})\tilde{\bm{f}}_{n,j}$$

(9)

To ensure numerical stability, a small positive constant is typically added to the diagonal of $\bm{K}_{\bm{x}_{j},\bm{x}_{j}}$ before taking its inverse. The corresponding predictive probability at the cutoff is then estimated by $\hat{p}_{j}=\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}\tilde{p}_{n,j}$, where $\tilde{p}_{n,j}=\psi(\tilde{\gamma}_{n,j}+\tilde{q}_{n,j})$. Its variance is estimated by $\hat{v}_{p,j}=\frac{1}{\tilde{N}_{j}}\sum_{n=1}^{\tilde{N}_{j}}(\tilde{p}_{n,j}-\hat{p}_{j})^{2}$. The RD effect of treatment assignment on treatment take-up probability is calculated as

$$\hat{\tau}_{\mathsf{SRDP}}=\hat{p}_{1}-\hat{p}_{0}$$

with posterior variance

$$\hat{v}_{\mathsf{SRDP}}=\hat{v}_{p,0}+\hat{v}_{p,1}$$

One can proceed to estimate the fuzzy RD effect with the ratio $\hat{\tau}_{\mathsf{SRD}}/\hat{\tau}_{\mathsf{SRDP}}$. This simple ratio estimator, however, is known to be biased and can be improved via a Jackknife refinement (Scott and Wu, 1981; Shao and Tu, 2012). Denote by $\hat{\tau}^{(-n)}_{\mathsf{SRD}}$ the sharp RD estimator on the response calculated with the $n$-th MCMC sample left out, and $\hat{\tau}^{(-n)}_{\mathsf{SRDP}}$ its counterpart on the treatment take-up probability. The Jackknife FRD estimator is given by

$$\hat{\tau}_{\mathsf{FRD}}=\tilde{N}\frac{\hat{\tau}_{\mathsf{SRD}}}{\hat{\tau}_{\mathsf{SRDP}}}-\frac{\tilde{N}}{\tilde{N}-1}\sum_{n=1}^{\tilde{N}}\frac{\hat{\tau}^{(-n)}_{\mathsf{SRD}}}{\hat{\tau}^{(-n)}_{\mathsf{SRDP}}}$$

with variance

$$\hat{v}_{\mathsf{FRD}}=\frac{1}{\tilde{N}}\left(\frac{\hat{v}_{\mathsf{SRD}}}{\hat{\tau}^{2}_{\mathsf{SRDP}}}+\frac{\hat{\tau}^{2}_{\mathsf{SRD}}{\hat{v}_{\mathsf{SRDP}}}}{\hat{\tau}^{4}_{\mathsf{SRDP}}}\right)$$

The estimator for the fuzzy RK effect is constructed analogously as follows:

$$\hat{\tau}_{\mathsf{FRK}}=\tilde{N}\frac{\hat{\tau}_{\mathsf{SRK}}}{\hat{\tau}_{\mathsf{SRDP}}}-\frac{\tilde{N}}{\tilde{N}-1}\sum_{n=1}^{\tilde{N}}\frac{\hat{\tau}^{(-n)}_{\mathsf{SRK}}}{\hat{\tau}^{(-n)}_{\mathsf{SRDP}}}$$

For fuzzy RD/RK design with a binary treatment, the FRK estimator shares with FRD the same ‘denominator’ estimator $\hat{\tau}_{\mathsf{SRDP}}$. Thus its variance is estimated by

$$\hat{v}_{\mathsf{FRK}}=\frac{1}{\tilde{N}}\left(\frac{\hat{v}_{\mathsf{SRK}}}{\hat{\tau}^{2}_{\mathsf{SRDP}}}+\frac{\hat{\tau}^{2}_{\mathsf{SRK}}{\hat{v}_{\mathsf{SRDP}}}}{\hat{\tau}^{4}_{\mathsf{SRDP}}}\right)$$

To a large degree, the expressiveness of a GP model is determined by its covariance. The estimators introduced in the preceding sections employ the squared exponential covariance function. In spite of its flexibility, it is stationary and can be inadequate if the underlying curve is non-stationary with rapid variations in its slope/curvature. There are two general ways to introduce non-stationarity to GP regressions (Rasmussen and Williams, 2006). One is to make the parameters of the covariance function input dependent. However it remains a challenge to maintain the positive-definiteness of the covariance with varying parameters. Another possibility is to introduce a non-linear transformation of the input $x$, say $g(x)$, and then construct a GP model in the $g(x)$ space.

We adopt the second approach and consider nonlinear transformations that are smooth and bounded, such as the logit, probit and tanh transformations. These sigmoid functions are widely used in neural network (NN) estimations, often referred to as activation functions. They have bounded derivatives and therefore provide desirable numerical stability, especially for RK estimations.

Given an activation function $g$, we construct our first stage mapping $g(x|\theta_{g})=g(\lambda_{0}+\lambda_{1}x)$, where $\theta_{g}=(\lambda_{0},\lambda_{1})$. We then apply the various GP models proposed in the preceding sections to $g(x|\theta_{g})$. For instance, the Gaussian observation model becomes $y_{i}=f(g(x_{i}|\theta_{g}))+u_{i}$. A Bayesian approach is adopted for the first stage as well. Thus the hyperparameters for this hierarchical GP model consist of $\theta=(\theta_{g},\theta_{f},\theta_{y})$ with prior distribution $p(\theta_{g})p(\theta_{f})p(\theta_{y})$. Figure 1 presents a graphical representation of the hierarchical GP model. This model can be interpreted as a four-layer deep learning model. The first and last layers are the inputs and outputs. The second layer is the NN transformation $g(x)$ of the input $x$ and the third layer constructs the GP latent variable based on $g(x)$. Unlike the other layers, the third layer does not assume (conditional) independence and instead models all units jointly as a Gaussian process.

Transformation of variables has had a long tradition in statistics. For instance, the power transformation and Box-Cox transformation have been customarily applied to non-Gaussian variables to reduce non-Gaussianity. In nonparametric estimations, the inputs are often transformed into features (for instance splines, wavelets in series estimation and perceptrons in one-layer neural network) and the subsequent curve fitting is conducted on the feature space. Two key innovations distinguish deep learning methods from these conventional methods (Murphy, 2012). First, in contrast to a ‘shallow’ model with a potentially large number of features, deep learning methods employ a multi-layer architecture, each layer nonlinearlly transforming its input into a slightly more abstract and composite representation. Second, many conventional methods that operate on the feature space can be characterized as un-supervised ‘feature engineering’ wherein the construction of features and the subsequent fitting are disconnected. In contrast, deep learning methods jointly train all layers.

The most popular deep learning models are multilayer neural networks (Goodfellow et al., 2016). The proposed hierarchical model can be characterized a hybrid estimator that features an NN layer followed by a GP layer. It is not uncommon that deep NN models feature thousands or even millions of tuning parameters and therefore can be ‘data hungry’. Given the sometimes small samples for RD estimations, we opt for the more economic hybrid model with two hidden layers. We prefer the GP approach to the NN as the penultimate layer for several reasons. First, it is probabilistic and offers automatic uncertainty quantification. Second, it is amenable to derivative estimations. Third, it is more expressive than an NN with a fixed number of nodes. As a matter of fact, the squared exponential covariance function can be obtained as the limit of basis function expansion with an infinite number of radial basis functions (see e.g., Chapter 4 of Rasmussen and Williams (2006)).

The marginal likelihood of the hierarchical GP model, conditional on $\theta$, is given by

$$\log p(\mathcal{D}|\theta)=-\frac{1}{2}\bm{y}^{T}(\bm{K}_{g(\bm{x}),g(\bm{x})}+\sigma^{2}\mathcal{I}_{N})^{-1}\bm{y}-\frac{1}{2}\log|\bm{K}_{g(\bm{x}),g(\bm{x})}+\sigma^{2}\mathcal{I}_{N}|-\frac{N}{2}\log 2\pi$$

and the unconditional marginal likelihood $p(\mathcal{D})=\int p(\mathcal{D}|\theta)p(\theta_{g})p(\theta_{f})p(\theta_{y})d\theta_{g}d\theta_{f}d\theta_{y}$. Let $\tilde{\theta}_{n,j}=(\tilde{\theta}_{g,n,j},\tilde{\theta}_{f,n,j},\tilde{\theta}_{y,n,j}),n=1,\dots,\tilde{N}$, be an MCMC sample generated from the posterior distribution $p(\theta_{j}|\mathcal{D}_{j})$. Define ${\bm{K}}_{g(x),g(x^{\prime})}(\theta)=k_{\mathsf{p+s}}(g(\bm{x}|\theta_{g}),g(\bm{x}^{\prime}|\theta_{g})|\theta_{f})$. The predictive mean, conditional on $\tilde{\theta}_{n,j}$, at the cutoff is given by

$$\tilde{m}_{n,j}(g)=\bm{K}_{g(c),g(\bm{x}_{j})}(\tilde{\theta}_{n,j})(\bm{K}_{g(\bm{x}_{j}),g(\bm{x}_{j})}(\tilde{\theta}_{n,j})+\tilde{\sigma}_{n,j}^{2}\mathcal{I}_{N_{j}})^{-1}\bm{y}_{j}$$

(10)

with conditional variance

	$\displaystyle\tilde{v}_{n,j}(g)=$	$\displaystyle K_{g(c),g(c)}(\tilde{\theta}_{n,j})$
		$\displaystyle-{\bm{K}}_{g(c),g(\bm{x}_{j})}(\tilde{\theta}_{n,j})(\bm{K}_{g(\bm{x}_{j}),g(\bm{x}_{j})}(\tilde{\theta}_{n,j})+\tilde{\sigma}^{2}_{n,j}\mathcal{I}_{N_{j}})^{-1}{\bm{K}}_{g(\bm{x}_{j}),g(c)}(\tilde{\theta}_{n,j})$

The predictive mean and variance at the cutoff are then calculated as

$$\hat{\mu}_{j}(g)=\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}\tilde{m}_{n,j}(g),\quad\hat{v}_{j}(g)=\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}(\tilde{m}_{n,j}(g)-\hat{\mu}_{j}(g))^{2}+\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}\tilde{v}_{n,j}(g)$$

Finally the treatment effect is estimated by

$$\hat{\tau}_{\textsf{SRD}}(g)=\hat{\mu}_{1}(g)-\hat{\mu}_{0}(g)$$

with variance

$$\hat{v}_{\textsf{SRD}}(g)=\hat{v}_{1}(g)+\hat{v}_{0}(g)$$

Next consider the SRK effect. Let $\dot{g}(x|\theta)=\partial g(x|\theta_{g})/\partial x$. Define $\dot{K}_{g(x),g(x^{\prime})}(\theta)=\dot{k}_{\textsf{p+s}}(g(x|\theta_{g}),g(x^{\prime}|\theta_{g})|\theta_{f})\dot{g}(x|\theta_{g})$ and $\ddot{K}_{g(x),g(x^{\prime})}(\theta)=\ddot{k}_{\textsf{p+s}}(g(x|\theta_{g}),g(x^{\prime}|\theta_{g})|\theta_{f})\dot{g}(x|\theta_{g})\dot{g}(x^{\prime}|\theta_{g})$. The predictive derivative, conditional on $\tilde{\theta}_{n,j}$, at the cutoff is given by

$$\tilde{\dot{m}}_{n,j}(g)=\dot{\bm{K}}_{g(c),g(\bm{x}_{j})}(\tilde{\theta}_{n,j})(\bm{K}_{g(\bm{x}_{j}),g(\bm{x}_{j})}(\tilde{\theta}_{n,j})+\tilde{\sigma}_{n,j}^{2}\mathcal{I}_{N_{j}})^{-1}\bm{y}_{j}$$

with conditional variance

	$\displaystyle\tilde{\dot{v}}_{n,j}(g)=$	$\displaystyle\ddot{K}_{g(c),g(c)}(\tilde{\theta}_{n,j})$
		$\displaystyle-\dot{\bm{K}}_{g(c),g(\bm{x}_{j})}(\tilde{\theta}_{n,j})(\bm{K}_{g(\bm{x}_{j}),g(\bm{x}_{j})}(\tilde{\theta}_{n,j})+\tilde{\sigma}^{2}_{n,j}\mathcal{I}_{N_{j}})^{-1}\dot{\bm{K}}_{g(\bm{x}_{j}),g(c)}(\tilde{\theta}_{n,j})$

The predictive mean derivative is then calculated as

$$\hat{\dot{\mu}}_{j}(g)=\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}\tilde{\dot{m}}_{n,j}(g)$$

with predictive variance

$$\hat{\dot{v}}_{j}(g)=\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}(\tilde{\dot{m}}_{n,j}(g)-\hat{\dot{\mu}}_{j}(g))^{2}+\frac{1}{\tilde{N}}\sum_{n=1}^{\tilde{N}}\tilde{\dot{v}}_{n,j}(g)$$

The SRK treatment effect is estimated by

$$\hat{\tau}_{\textsf{SRK}}(g)=\hat{\dot{\mu}}_{1}(g)-\hat{\dot{\mu}}_{0}(g)$$

with posterior variance

$$\hat{v}_{\textsf{SRK}}(g)=\hat{\dot{v}}_{1}(g)+\hat{\dot{v}}_{0}(g)$$

The hierarchical GP models under the fuzzy RD and RK designs are constructed analogously. For brevity, the details are omitted.