{'fraction_non_alphanumeric': 0.06260794473229707, 'fraction_numerical': 0.022428516599501055, 'mean_word_length': 4.207219585311017, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Bayesian Nonparametric Causal Model for Regression Discontinuity Designs

For non-randomized studies, the regression discontinuity design (RDD) can be used to identify and estimate causal effects from a"locally-randomized"subgroup of subjects, under relatively mild conditions. However, current models focus causal inferences on the impact of the treatment (versus non-treatment) variable on the mean of the dependent variable, via linear regression. For RDDs, we propose a flexible Bayesian nonparametric regression model that can provide accurate estimates of causal effects, in terms of the predictive mean, variance, quantile, probability density, distribution function, or any other chosen function of the outcome variable. We illustrate the model through the analysis of two real educational data sets, involving (resp.) a sharp RDD and a fuzzy RDD.


Introduction
A basic objective in scientific research is to infer causal effects from data. Randomized studies are the gold standard of causal inference. In an ideal randomized study, the investigator randomly assigns each subject into one of the treatment conditions, with equal probability, and each subject complies with her/his treatment assignment. Then, treatment subjects are the same as non-treatment subjects, in terms of the distribution of all observed and unobserved pre-treatment covariates, aside from sampling error (e.g., [45]); and the outcome variable is independent of the chosen treatment intervention, conditionally on the treatment variable [15]. Then the causal effect is given by a comparison of the outcome variable under the treatment intervention, against the outcome variable under the non-treatment intervention.
Often, it is necessary to estimate causal effects from a non-randomized, observational study, because a randomized study can be infeasible due to financial, ethical, or time constraints [45]. However, causal inferences from a nonrandomized study is more challenging, because without randomization, treated and non-treated subjects differ almost-surely in terms of the pre-treatment covariates.
The regression discontinuity design (RDD) ( [49], [11]) is a type of non-randomized design, where a continuous-valued assignment variable [36] assigns each subject to the treatment (non-treatment, resp.) condition, whenever her/his observed value of the assignment variable equals or exceeds (is less than, resp.) a fixed cutoff value. Under relatively mild conditions, notably when subjects have imperfect control of the assignment variable, the RDD provides a "locally-randomized experiment," so that "are as good as randomly assigned" for the subgroup of subjects with assignment variable values near the cutoff [35]. Then the causal effect is identifiable for that subgroup. As proven [24], the RDD can empirically produce causal effect estimates that are similar to those estimates of a standard randomized study ( [1], [10], [8], [46], [7], [47]).
The RDD has existed over 50 years, with little initial interest [11]. However, since 1997, more than 74 RDD-based empirical studies have emerged from these fields ( [36], [9], [52], [37]), for at least three reasons ( [17], [36]). First, many nonrandomized studies employ treatment assignment rules that can be easily conceptualized as RDDs. Second, the empirical results of RDDs are intuitive and can be easily conveyed graphically, say, by a plot of the outcomes against the assignment variable. Third, the identification of causal effects in a RDD requires weaker and hence more credible assumptions, compared to the stronger assumptions that are required by other popular causal models, mentioned below. This gives the researcher the flexibility to choose from a range of causal estimation methods.
The other popular causal models for non-randomized studies assume a "potential outcomes" (counterfactual) framework of causal inference (e.g., [42], [43]). This is typically using notation simplified by the Stable Unit Treatment Value As-sumption (SUTVA), which implies no interference between subjects and no versions of treatments [44]. The popular models make further assumptions of unconfoundedness (i.e., treatment and non-treatment outcomes are independent of treatment assignments, conditionally on all pre-treatment covariates) and overlap (i.e., there is a chance to receive either the treatment or the non-treatment, conditionally on any value of the pre-treatment covariates) [28]. These models are defined either by a regression of the outcome variable, on variables of treatment receipt and observed pre-treatment characteristics, and/or involves matching/weighting subjects on the observed pre-treatment variables and/or on propensity scores (e.g., [28]). The regression may also be on a hypothesized set of unobserved pre-treatment covariates, in order to study the sensitivity of causal effect estimates over varying degrees of hidden bias (e.g., [41]), i.e., over changes in the distribution of these covariates. However, it may be argued that for typical non-randomized studies, unconfoundedness and overlap are not credible assumptions (e.g., [28], [35]), and even SUTVA is questionable. We provide a more detailed discussion of these issues in Section 2.4.
For RDDs, the mainstream causal models are either linear, polynomial, or local-linear models that employ a regression of the outcome variable on the assignment variable. Such models aim to provide causal inferences in terms of mean comparisons of treatment outcomes and non-treatment outcomes, and to provide sufficiently-flexible modeling of the regression function ( [28], [36]), in a neighborhood around the cutoff. However, outliers in a linear model may bias causal effect estimates. Furthermore, in many settings, it may also be of interest to base causal inferences on comparisons of additional features of the outcome variable, such as the variance, quantiles (percentiles), and/or the entire probability density function.
To address these modeling issues, we propose a Bayesian nonparametric regression model [33] for causal inference in RDDs. It is an infinite-mixture model that allows for the entire probability density of the outcome variable to change flexibly as a function of covariates. Our model can provide inferences of causal effects in terms of how the treatment variable impacts the mean, variance, a quantile, probability density function (p.d.f.), distribution function, and any other chosen function of the outcome variable. Finally, the accurate estimation of causal effects relies on an appropriate model for the data. Karabatsos and Walker [33] showed that their Bayesian nonparametric regression model tended to have better predictive performance than other parametric and flexible nonparametric regression models of common usage, over many real data sets. Also, our model can be extended to handle causal inferences from a fuzzy RDD [50]. In contrast to a standard "sharp" RDD, a fuzzy RDD involves a study where not all subjects adhere to the treatment assignment rule. This is because, for example, some subjects do not comply with their respective treatment assignments, or because some subjects receive treatments for which they are not eligible. Moreover, the model can be easily extended to handle an assignment variable that has multiple cutoffs, or to handle a multivariate assignment variable that has a vector-valued cutoff.
In Section 2, we review the data assumptions that are required to identify and estimate causal effects from a RDD. Unlike all previous guides to performing causal inference from RDDs (e.g., [30], [36], [9], [52]), we avoid using the potential outcomes approach to causal inference. This approach has questionable foundations, and besides, causal inference can be handled entirely by standard probability theory ( [14], [15]).
In Section 3, we review the current models, and then describe our Bayesian nonparametric model that can estimate causal effects from the various RDDs. In Section 4, we illustrate our model through the analysis of three real educational data sets. These data sets involve (resp.) a sharp RDD, a sharp multivariate RDD, and a fuzzy RDD. Section 5 ends with conclusions.

Identifying Causal Effects in a RDD
A non-randomized study from a RDD involves three variables that are observable from each of a sample of n subjects, indexed by i = 1, . . . , n. They are the outcome variable, Y ; a binary treatment variable T , where T = 1 refers to treatment receipt and T = 0 refers to non-treatment receipt; and a continuous-valued assignment variable R. Each subject i is assigned the treatment whenever R i ≥ r 0 , and is assigned the non-treatment whenever R i < r 0 , given a known fixed cutoff r 0 . The treatment assignment variable is denoted 1 R≥r 0 , with 1 (·) the indicator (0,1) function. An RDD study gives rise to a sample data set, , including derived observations 1 r i ≥r 0 , and possibly observations of p pre-treatment covariates A characterizing assumption of the RDD is that the conditional probability of treatment receipt is discontinuous at r 0 , i.e., where r + 0 = lim r↓r 0 r and r − 0 = lim r↑r 0 r. Equation (1) is named Assumption RD.
There are two types of RDDs. In the classical, sharp RDD ( [49], [11]), the probability function Pr(T = 1 | r) has a discontinuous jump of size 1 at R = r 0 , with point mass probability density: Then treatment receipt is identical to treatment assignment, with T = 1 R≥r 0 .
In the fuzzy RDD [50], Pr(T = 1 | r) has a discontinuous jump that is smaller than 1 at r 0 . Then f (t|r) is not a point mass density. The smaller jump results from imperfect treatment adherence (e.g., treatment non-compliance), where some of the subjects of the given study were either assigned 1 r≥r 0 = 0 but received T = 1, or assigned 1 r≥r 0 = 1 but received T = 0.
For either type of RDD, a typical measure of the causal effect is given by the difference of conditional means (expectations) of Y at R = r 0 : where E(Y | r,t) = y dF(y | r,t), and r + 0 = r − 0 = r 0 when R is continuous. Conditioning on R = r 0 is motivated by the fact that Pr(T = 1 | r) jumps at that point.
In general, for any choice of function H{·} of Y , the causal effect is: Therefore, depending on the choice of function H{·}, causal effects are not only interpretable in terms of the mean of Y (when H{Y } = Y ), but also in terms of comparisons of the variance (H{Y } = {Y − E(Y | r,t)} 2 ), cumulative distribution function (c.d.f.) F(y | r,t) (H{Y } = 1 Y ≤y ), probability density function (p.d.f.) ( f (y | r,t)), survival function 1 − F(y | r,t), and so on. Inverting the c.d.f. obtains F −1 (u | r,t), for u ∈ [0, 1], so then causal effects can also be interpreted in terms of the u th quantile of Y . The outcome Y (or H{Y }) is observable only under the treatment (non-treatment, resp.) when R ≥ r 0 (R < r 0 , resp.), and there are no observations of R = r 0 due to the continuity of R. Therefore, the estimation of a causal effect τ H , from data D n , requires some degree of extrapolation.
Next we describe how a causal effect τ H is identified from the sharp and fuzzy RDDs.

Identification in the Sharp RDD
Given sample data D n from a sharp RDD, a causal effect τ H is inferred from D n using a model for the conditional probability density function (p.d.f.), f (y | r,t), Suppose that there is reason to believe that in the absence of treatment, subjects close to the threshold r 0 are similar. Then the causal effect τ H is identified by assumption RD (1) and that: Equation (5) is a density version of the assumption in [27], who only look at mean shifts and hence assume that We believe it is important to model mean shifts not by having a mean shift model but rather by modeling the density of the observations and then picking out the mean from this. This is the correct approach and hence the usefulness of the nonparametric model. Assuming (5) rather than only (6) allows the treatment effect to exhibit itself in more ways than a mean shift. For example, a variance shift would also be informative, even in the absence of a mean shift. The model we employ is quite general and allows many aspects of treatment effect to be explored by studying any differences between density estimate either side of the cut-off at r 0 . However, we can obviously estimate the key mean shift, having modeled the density functions either side of the cut-off point r 0 , simply by estimating the means of the two density functions.
The assumptions of RD (1) and continuity (5) are sufficient for estimating the treatment effect, when R is sufficient for the outcome, in the sense that no other variables effect Y at R = r 0 . Authors such as [36] further elaborate on the continuity assumption, given the existence of observed and unobserved pre-treatment covariates of the subjects, collectively labeled as W , which could also influence the outcome Y . Their elaborations are as follows.
If subjects have imprecise control of R at r 0 , i.e., then: (a) local randomization at r 0 holds, i.e., f (w | r) is continuous in r at r 0 for all W = w; (b) the identifying, continuity conditions (5, 6) hold ( [27], [35], [36]); and (c) the marginal density f (r) is continuous in r at r 0 , to provide an empirically-testable consequence of the imprecise control assumption (see [39]). Local randomization is a consequence of imprecise control because the continuity of f (r | w) in r implies that the right hand side of the Bayes' rule equation, , the distribution of W is the same for treatment subjects and for non-treatment subjects, having R values located near r 0 . Hence, treatments are "as good as randomly assigned" for this local subgroup of subjects (e.g., [35]).
Under local randomization, (W, R) represent all the information which could influence the outcome Y . So, if is continuous in r at r = r 0 , then any non-zero value for can only be caused by the difference in treatment. In particular, the continuity of f (w | r) at r 0 implies that (9) does not depend on w. Therefore, in the remainder of the paper we will assume that R is sufficient for Y , for the sharp RDD.

Conditional Independence Properties of the Sharp RDD
The sharp RDD can also be characterized by the extended conditional independence framework of causal inference, extending the ideas from [15], Sections 6.2, 7. Again, we consider the random variables (R, T,Y ), and make assumptions RD (1) and imprecise control at r 0 (7), implying that R is sufficient covariate and that W is ignorable. Also, for this framework, we introduce a non-random regime parameter Ψ T ∈ { / 0, 0, 1}. An "idle" or "observational" regime is indicated by Ψ T = / 0, when the distribution of (Y, T, R) "arises naturally." An intervention regime is indicated by setting Ψ T = 0 or 1, and such an intervention modifies the idle distribution, as will be described below.
Under a sharp RDD, the joint p.d.f. of (R, T,Y ) is: This p.d.f. admits the property Y ⊥ ⊥ Ψ T | R = r 0 , T , where ⊥ ⊥ denotes conditional independence [13]. This is a causal property, because it says that the distribution of Y is unaffected by the choice of interventional regime Ψ T , conditionally on (R = r 0 , T ) [15]. We therefore focus on the conditional p.d.f.: If the intervention parameter take on a null value Ψ T = / 0, then the joint distribution of the random variables (T,Y ) arises naturally, and then (11) is the point mass density of equation (2). In contrast, an intervention that sets , and then changes the p.d.f. (11) to: where || denotes "conditioning by intervention" [34]. The equality in (12) holds by virtue of the causal property. Then for a general choice of function H{·}, the causal effect is given by a comparison of E(H{Y } | r 0 ,t), for t = 0, 1, including the p.d.f.s f (y | r 0 ,t).

Identification in the Fuzzy RDD
In the fuzzy RDD, the probability function Pr(T = 1 | r) has a discontinuous jump that is smaller than 1 at r 0 , meaning that R does not determine T . Then 1 R≥r 0 and T are distinct variables since the event 1 R≥r 0 = T is possible; and then T and Y may both depend unobserved confounding variables, collectively labelled as U .
Given these considerations, we may extend the joint p.d.f. (10) to: But for the purposes of making causal inferences from real data, we cannot condition on (R, T,U ), because U is unobserved.
However, for the fuzzy RDD, the assumptions of RD (1) and continuity (5) at r 0 (incl. imprecise control) together imply that the causal effect can be represented by: for a general choice of function H{·} of Y . This follows from the equality: and the fact that the term (15b) is zero under the continuity assumption (5) [27].
In the ratio (14), the numerator is the Intention to Treat (ITT) effect. The denominator is a measure of treatment adherence, which decreases as noncompliance increases. In the sharp RDD, where 1 R≥r 0 = T , the denominator is 1 (i.e., perfect adherence), and then the ITT effect coincides with the causal effect τ H .

On Unconfoundedness and Overlap Assumptions
Arguably, the assumptions of RD (1) and imprecise control at r 0 (7) (implying the continuity condition at r 0 (5)), which together identify a causal effect τ H in either a sharp or a fuzzy RDD, are weaker and hence more credible ( [35], p. 679), compared to the unconfoundedness and overlap assumptions of some popular causal models that employ the potential outcomes framework of causal inference (e.g., [42], [28]). The potential outcomes framework employs notation such as (Y (1),Y (0)), denoting the outcomes to treatment (T = 1) and to non-treatment (T = 0), at a given time point. The use of this specific notation requires SUTVA, which states that the 2 n potential outcomes Y i (t n ), t n ∈ {0, 1} n are reducible to two potential outcomes, (Y i (1),Y i (0)) [44], for every subject i of the n subjects of a given study. Unconfoundedness requires (Y (1),Y (0)) ⊥ ⊥ T | w, and overlap requires 0 < Pr(T = 1 | w) < 1, for all w. If the assumptions of unconfoundedness and overlap hold, then the causal effect of T on Y is identified for all values of the observed and unobserved pre-treatment covariates, w [28]. However, the assumption that unconfoundedness and overlap hold, for all w, is questionable [28], and moreover uncounfoundedness is not empirically testable. Besides, it may be argued that it is virtually impossible for an investigator of a non-randomized study to have precise prior knowledge of, and then collect data on, all the important pre-treatment covariates.
Albeit, the RDD can identify causal effects only for the subgroup of subjects with R values near r 0 , pertaining to only a subset of all possible values of the pre-treatment covariates, w. Thus, the RDD trades off generalizability, for more credible assumptions in causal effect identification. Also the RDD requires a sufficiently-large number of observations of R near r 0 , to allow for statistically precise estimation of a causal effect τ H . Finally, because our review of the RDD avoids the potential outcomes framework, we do not address SUTVA. This assumption does not make sense outside of this framework [14].

Estimating Causal Effects In a RDD
Here, we review the current (published) models for estimating causal effects from sample data D n of an RDD. Then we propose our Bayesian nonparametric model for causal inference for RDDs. Afterwards, we discuss how any of these models can be extended to handle RDDs with either multiple cutoffs, or with a multivariate assignment variable.
Similarly, β x (x) = ∑ p k=1 ∑ q k l=1 β xkl x l k is a sum of polynomials for p observed pre-treatment covariates x, a subset of the covariates that comprise W . Often in practice, x is excluded from the model (16) (i.e., β x (x) = 0), because R is a sufficient covariate under assumptions RD (1) and continuity (5), making W (and hence x) ignorable. Then, excluding the covariates x from the model has little effect on the bias of the least-squares estimate τ of the causal effect (3). However including these covariates in the model may improve estimation efficiency [30].
Under the classical linear model (16), outliers of y i can bias estimates of the causal effect τ, even for observations r i far from r 0 [11]. Therefore, local linear models [18] were proposed for RDDs, because they can provide more outlierrobust estimation of τ, by placing more emphasis on observations (y i , r i ) with values r i near r 0 . A weighted local linear model (16) assigns smooth kernels K({r i − r 0 }/h) > 0 to observations (y i , r i ) with chosen bandwidth h > 0. Kernels that have been considered in the RDD literature include the edge kernel (K(u) = (1 − |u|)1 |u|≤1 ) and the uniform kernel (K(u) = 1 |u|≤1/2 ) ( [27], [30], [29]). The parameter of the local linear model can be estimated by generalized least-squares. However, the choice of bandwidth parameter can impact the bias of the causal estimate τ. For the edge and uniform kernels, [29] describe methods for selecting the optimal bandwidths for RDDs, having large-sample justifications. Another local approach is to specify a model (16) with terms β m (r) (m = 1, 2) each constructed by a spline basis [16]. For example, the truncated linear spline basis, β m (r) = β m1 r + ∑ L m l=1 β ml max(0, r − r ml ), with L k knots r ml including r 0 . The weighted local linear model was extended to handle quantile regression [19]. This approach entails independent regressions of multiple indicator vari- However, this model suffers from the "quantile crossing" problem [6], and relies on an ad-hoc resorting correction that has only large sample justifications. The quantile crossing problem refers to the possibility of incoherent quantile estimates F −1 (u | r, 1 r≥r 0 ) > F −1 (u ′ | r, 1 r≥r 0 ) where u < u ′ , given some (r, 1 r≥r 0 ).
Any of the linear models can be extended to provide an estimate τ H of a causal effect for a fuzzy RDD. One approach is to take the ratio estimate τ H = τ (F) H = τ HY / τ HT , for (14), where τ HY is the coefficient estimate of 1 R≥r 0 in a linear regression of the outcome Y (e.g., H{Y }) on (R, 1 R≥r 0 ), and τ HT is the coefficient estimate of 1 R≥r 0 in a linear regression T on (R, 1 R≥r 0 ) ( [30], [9], [19]). Again, estimation efficiency can be improved by also conditioning on covariates x. The standard error (SE) of the estimate τ (F) H can be approximated by SE( τ HY )/ τ HT [9], with the assumption that T (r) is non-increasing in r at r 0 , under the potential outcomes framework [30].
Finally, for assessing the fit of a linear model to data, we may inspect standardized residuals z i = (y i − η i )/ σ of the observations y i (or H{y i }), given leastsquares estimates η i = ((1, r i , 1 r i ≥r 0 , x ⊺ i ) ⊺ β , σ). If | z i | > 2, then y i can be judged as an outlier.

Bayesian Nonparametric Model
For the sharp RDD, our Bayesian nonparametric model [33] is defined by: where the mixture weights ω j (η(r), σ(r)) sum to 1 at each value of r, and T = 1 R≥r 0 in the sharp RDD. Also, the terms (17a) and (17b) may be deconstructed via the generation of a latent indicator variable Z ∼ N(η, σ 2 ), and then taking (17) allows the entire probability density of the outcome variable Y to change flexibly as a function of covariates. The parameter σ (r) measures the multimodality of f (y | r, 1 r≥r 0 ) [33]. Specifically, as σ (r) → ∞, the density f (y | r, 1 r≥r 0 ) becomes more multimodal, with weights ω j (η(r), σ(r)) converging to a discrete uniform distribution; and as σ (r) → 0, the density f (y | r, 1 r≥r 0 ) becomes more unimodal, and "local," with f (y | r, 1 r≥r 0 ) ≈ n(y i | µ j , σ 2 j ) and ω j (η(r), σ(r)) ≈ 1 if j − 1 < η ≤ j. Furthermore, the model has a discontinuity at r 0 due to the presence of the term 1 r≥r 0 in both (17c) and (17d). The effect, controlled by the coefficients (β , λ ), is to reallocate the weights either side of r 0 , resulting in different densities either side of this value. Obviously, there is a discontinuity if and only if either of the coefficients (λ 2 , β 2 ) is non-zero. The normal prior N(β , λ | 0,vI) consists of a prior variance parameter v, which controls for both the prior support for the range of the mixture density component indices j = 0, ±1, ±2, . . . (via the parameter β ), and for the range of the level of multimodality in f (y | r, 1 r≥r 0 ). As v → ∞, a wider range of component densities and multimodality is supported; and as v → 0, f (y | r, 1 r≥r 0 ) becomes a normal density.
When prior information is limited about the model parameters, we may attempt to specify non-informative priors, for example, by choosing µ 0 = 0, σ 2 0 → ∞, a 0 → 0, b 0 → 0, and v = 10 5 , and by choosing b σ µ according to prior knowledge about range of the Y variance. For instance, if Y is known to have a variance of 1, then b σ µ = 5 provides a vague prior choice. For such choices of prior parameters, the Bayesian model (17), over 22 real data sets, demonstrated very good predictive accuracy, and better predictive accuracy compared to many other regression models, and compared to the Bayesian model under different choices of prior [33].
The model (17) has infinite-dimensional parameter, updates the prior π(ζ ) to a posterior density, given by , with Π(ζ ) (and Π(ζ | D n ), resp.) the c.d.f. of π(ζ ) (of π(ζ | D n )). Also, let F(y | r, a; ζ ) be the c.d.f. of f (y | r, 1 r≥r 0 ; ζ ). Then the posterior predictive density, f n (y | r, a), and the conditional posterior predictive expectation (E n ) and variance (V n ) of the outcome H{Y } are given (resp.) by: Depending on the choice of function H{·}, the posterior mean E n and variance V n of the conditional expectation E(Y | r, a), variance V(Y | r, a), c.d.f. F(y | r, a) at a point y, are given (resp.) by E n (E(Y | r, a)) and V n (E(Y | r, a)); E n (V(Y | r, a)) and V n (V(Y | r, a)); and E n (F(y | r, a)) = F n (y | r, a) and V n (F(y | r, a)). For assessing the fit of the Bayesian model to data, a standardized residual for each observation y i may be computed by If |z i | > 2, then y i can be judged as an outlier. Alternatively, we may consider a version of our model that performs covariate selection, via the specification of spike-and-slab priors [21] for (β , λ ). Such priors may be defined by prior densities say, for large v 1 = 1000 and small v 0 = .1 [33]. Then γ j = 1 (γ j = 0, resp.) means that the given covariate is included (excluded, resp.) from the model. In this version of the model, γ would be added to the full parameter vector ζ , and Bayesian inference of the model would proceed in the same way as described above. A covariate can then be judged to be a "significant predictor" when the posterior probability is Pr[γ j = 1| D n ] ≥ .5 [5].

Estimating Causal Effects With the Bayesian Model
For either version of our Bayesian model, the posterior estimates of the causal effect of T on H{Y }, conditionally on R = r 0 , are given as follows, under the assumptions RD (1) and continuity at r 0 (5).
For the sharp RDD, the estimate τ H of the causal effect is given by  , 1), we may judge whether τ H is significantly different from zero, by using a P-P plot [51] to check for non-overlap of the 95% posterior credible intervals (F n(.025) (y | r, a), F n(.975) (y | r, a)) at u, for a = 0, 1, and over a wide range of y ∈ Y . Here, F n(.025) (y | r, a) (F n(.975) (y | r, a), resp.) denotes the posterior 2.5%ile (97.5%ile, resp.) of F(y | r, a).
For the fuzzy RDD, the causal effect estimate τ H , in terms of the ratio estimator (14), may be obtained by two independent regressions. The first involves estimating the numerator using our regression model (17), and the second involves a regression of T on (R, 1 R≥r 0 ) to estimate the denominator, via the posterior predictive expectations: Our model (17) can be extended to binary regression, by replacing the model likelihood (17a) with This provides flexible modeling of the inverse link function by a covariate dependent, infinite mixture of normal c.d.f.s [33]. Since the posterior distributions of these two models are independent, the joint posterior of both models can either be estimated separately or jointly. Then for the fuzzy RDD, the estimate τ H of the causal effect, in terms of the ratio (14), is given by: The general estimator (19) provides a computationally-fast, first-order Taylor approximation to the posterior expectation E n (τ Using Markov Chain Monte Carlo (MCMC), Gibbs sampling methods, along with a slice sampling step for σ µ , can be used to estimate all of the aforementioned posterior quantities [33]. We use Rao-Blackwell (RB) methods to estimate all the posterior linear functionals, such as E n (H{Y } | r, a), V n (H{Y } | r, a),
Alternatively, it is possible to consider a model with covariates 1 R≥r 0 including 1 R≥r 0d . This approach adds D − 1 extra linear predictors (covariates) into the model, and thus this may improve the estimation efficiency of τ H . For the Bayesian model, it is possible to assign spike-and-slab priors to the coefficients of 1 R≥r 0 , in order to provide a type of change-point regression model.

Extensions to Multivariate RDDs
A RDD may involve a multidimensional continuous assignment variable, denoted by R = (R 1 , . . ., R D ) ⊺ , with vector cutoff r 0 = (r 01 , . . . , r 0D ) ⊺ (for d = 1, . . . , D). For such a design, it is possible to reduce R to a scalar assignment variable, R, using one of at least two approaches. First, the centering approach [32] defines a new assignment variable by R = min d (R d − r 0d ). The second approach defines a new assignment variable R by a thin-plate spline [25], with R = || R − r 0 || log || R − r 0 ||, and with || · || the Euclidean norm. If either approach is used, then all the same approaches to identifying and estimating the causal effect τ H are applicable (Sections 2, 3.1-3.2) to the new scalar assignment variable R and cutoff r 0 = 0.

Illustrative Applications
The Bayesian nonparametric model was illustrated through the analysis of three data sets, using menu-driven software that was developed by the first author ( [32], [31]). The first two data sets were collected from four Chicago University schools of education, which established a new curriculum that aims to train and graduate teachers to improve Chicago public schools. Each of these data sets involved a sharp RDD, using (resp.) a scalar-valued and a multivariate assignment variable. The third data set, obtained from [4], involves a fuzzy RDD, from a study of the effect of class size on student achievement [3]. For each of the three data sets, it seems reasonable to make the assumptions of RD (1) and continuity (5) at r 0 (i.e., imprecise control at r 0 ), in order to identify the causal effects of treatment on the outcome, conditionally on r 0 (see Section 2).
For all three data sets, the Bayesian nonparametric model assumed the same vague priors that were mentioned in Section 3.2. The model under these priors provided good predictive accuracy for each of these three data sets, and better predictive accuracy compared to the model under different choices of these priors. Also, the model under the vague priors yielded similar estimates of causal effects, compared to the model with the spike-and-slab priors. Hence, in this section, we will focus on presenting the results of the Bayesian model under the vague priors.
All posterior estimates of this model, reported in the next three subsections, are based on 40K MCMC samples that were obtained from every 5th iterate of a run of 200K MCMC sampling iterations, after discarding the first 2K burn-in samples. This provided accurate posterior estimates according to standard convergence assessments [22]. Specifically, univariate trace plots displayed good mixing of model parameters and posterior predictive samples, and all posterior predictive estimates obtained 95% MC confidence intervals with half-width sizes near .01. Also, for comparative purposes, we also analyzed each of the three data sets using different versions of the standard linear model (16), under least-squares estimation. They including the locally-weighted linear models with optimal bandwidth selection [29].

Learning Math Teaching: Time Series Data
For the first data set, the aim is to estimate the effect of the new teacher education curriculum on math teaching ability, among n = 347 undergraduate teacher education students attending one of four Chicago universities. This data set involves a sharp RDD, an interrupted time-series design [12] using an assignment variable of time, ranging from fall semester 2007 through spring semester 2013. The new curriculum (treatment) was instituted in Fall 2010 (the cutoff, r 0 ), and the old teacher curriculum (non-treatment) was active before then. The outcome variable (Y ) is the number-correct score on the 25-item Learning Math for Teaching (LMT) test [38]. Each of the students completed the LMT test (89.9% female; 135 and 212 students under the old and new curriculum), after finishing a course on teaching algebra. Among them, the average LMT score was 12.9 (s.d. For the model, R-squared was .99, and nearly all the standardized residuals ranged between −1 and 1 over the 347 observations, with one residual slightly exceeding 2. Figure 1 presents the model's posterior predictive density estimate of the LMT outcome, for the new curriculum (treatment) and for the old curriculum (non-treatment), at Fall 2010. The new curriculum, compared to the old, increased the LMT scores, by shifting the density of LMT scores to the right. This shift corresponds to an increase in the mean (from −.17 to −.13), the 10%ile (−2.01 to −1.97), and the 25%ile (−1.43 to −1.31), but these increases were not statistically significant from zero according to 95% credible intervals of the predictive mean and of the posterior c.d.f. estimates. Also, each density presents two modes (clusters) of students.
Most U.S. schools of education base their undergraduate admissions decisions on the ability to pass basic skills tests. We investigate the causal effect of basic skills on teacher performance, as previously studied in the literature on teacher education (e.g., [23]).
Here, the assignment variable is a 4-variate random variable, defined by Illinois Test of Basic Skills subtest scores, in reading, language, math and writing. Each subtest has a minimum passing score of 240. The dependent variable (Y ) is total score on the 50-item Haberman Teacher Pre-screener assessment. A score in the 40-50 range indicates a very effective teacher. This assessment has a testretest reliability of .93, has a 95% accuracy rate in predicting which teachers will stay and succeed in the teaching profession, and is used by many schools to assess applicants of teaching positions [26]. Over the 205 students, the Haberman scores Using the Bayesian model, we analyzed the data set to estimate the causal effect of passing the reading basic skills exam (treatment), versus not passing (nontreatment), on students' ability to teach in urban schools. The model included the Haberman z-score as the outcome (dependent) variable (Y ), and included as covariates the assignment variable (R), B240d10, and an indicator of passing all four basic skills subtests BasicPass = 1 B240d10≥0 .
For the data, the Bayesian model obtained an R-squared of .98, with standardized fit residuals ranging from −1.7 to 1.3 over the 205 observations. Thus the model had no outliers. Figure 2 presents the model's posterior predictive density estimates, for the treatment versus the non-treatment (resp.). A more detailed inspection revealed that passing the basic skills reading test causally increased the Haberman z-score, in terms of the mean (from −.10 to .09), 25%ile (−1.11 to −1.06), median (−.19 to .05), 75%ile (.74 to 1.16), 90%ile (1.44 to 2.36), and 95%ile (2.18 to 2.82), and causally decreased in terms of the 5%ile (−1.69 to −3.66). Also, the treatment density and the non-treatment density each showed 2 modes (clusters) of students.
------ Figure 2 here ------We also analyzed the data using several linear regression models. We considered strictly linear and polynomial models, and linear/polynomial interaction models, up to order 5. Among these models, the linear interaction and strictly cubic models attained the best AIC scores, and (resp.) yielded estimated causal effects of 1.22 (2-sided p-value = .00) and 1.06 (p = .02). These models had between 6 and 7 outliers. The local linear models with edge and uniform kernels (resp.), yielded estimated causal effects of 1.11 (p = .01), and 1.21 (p = .02), and 6 and 2 outliers.

Maimonides' Data: Fuzzy RDD
The twelfth-century rabbinic scholar Maimonides proposed a rule that specifies a maximum class size of 40, under the belief that smaller class sizes promotes higher student achievement. Specifically, for a given class c in school s, the rule assigns average class size (Psize sc ) as a function of beginning-of-the-year school enrollment (e s ), according to the prediction equation Psize sc = e s /floor[((e s − 1)/40) + 1]. The rule (equation) assigns students of a school into a single classroom, when the school's enrollment is less than 41, assigns students into two classrooms of average size 20.5 when school enrollment reaches 41; assigns students into three classrooms of average size 27 when enrollment reaches 81, and so on. The cutoff number 20.5 distinguishes between small and large classes.
Here, we study the effect of class size on average class verbal achievement, through the analysis of data on 4th grade students who each attended one of 2,056 classes in Israeli public schools during 1991. These schools used Maimonides' rule to allocate students into classrooms. Demographic statistics are reported in [3] (three other classes were not analyzed because they had missing achievement data). For the Bayesian model, the dependent variable (Y ) is average class verbal score (avgverb), which we transformed to z-scores with mean 0 and variance 1. The covariates include the assignment variable (R), defined by the rule-predicted class size centered at the cutoff 20.5 (i.e., Psize205 = Psize − 20.5), and include the indicator of large (vs. small) class assignment, Plarge = 1 Psize≥20.5 . Now, while Maimonides' rule may assign a given class to be a large (small, resp.), the class could become small (large, resp.). For example, one school in the data set had an enrollment of 41, leading to some students receiving a large class of 21, and other students receiving a small class of 20. Therefore, the data arise from a fuzzy RDD, and for the data analysis, we also consider a variable defined by the indicator of large class receipt, large = 1 classize≥20.5 . We also fit the Bayesian model, with the treatment (T ) variable, large, as the dependent variable, and with covariates Psize205 and Plarge.
For the avgverb (Y ) dependent variable, the Bayesian model obtained an Rsquared of .88, with standardized fit residuals ranging from −1.1 to 1.3 over the 2, 056 observations. Thus the model had no outliers. For the treatment (T ) dependent variable, large, the Bayesian model had no outliers, and estimated .93 as the denominator of the causal effect estimator (14). Figure 3 presents the model's posterior predictive density estimates, for the treatment versus the non-treatment. The left panel of the figure presents the density estimates (expectations) of avgverb, in terms of the numerator of the effect estimator (14), obtained from the Bayesian model with avgverb as the dependent variable (Y ). The right panel presents the density estimates divided by .93. It was found that passing the basic skills reading test causally increased class size, in terms of the 5%ile (−2.45 to −2.19), 10%ile (−1.81 to −1.66), 25%ile (−.94 to −.88), and causally decreased in terms of the 75%ile (.90 to .79), 90%ile (1.48 to 1.37), and 95%ile (1.76 to 1.66). Each of these estimates are based on taking the predictive quantile estimates of the Bayesian model for avgverb (Y ), and dividing them by .93.
------ Figure 3 here ------Various linear models were also fit to the data. For the linear regression modeling of avgverb (Y ) on (Psize205, Plarge), and for the linear regression modeling of large (T ) on (Psize205,Plarge), we considered strictly linear and polynomial models, and linear and polynomial interaction models, up to order 9. For the avgverb dependent variable, the septic (order 7) interaction model yielded the best AIC value, however, it also yielded a large least-squares estimate τ Y = 3.33 (SE = 1.59) and 75 outliers. For the large (T ) dependent variable, the linear interaction model (order 1) attained the best AIC value, had no outliers, and yielded τ T = .94 (SE = .02). Then the estimate of the causal effect of large (versus small) class is τ = τ Y / τ T = 3.33/.94 = 3.54, with approximate SE = 1.59/.94 = 1.69. The local linear models with edge and uniform kernels (resp.), for the avgverb (Y ) dependent variable, yielded τ Y estimates of .82 (SE = .23), and .39 (SE = .27), and 144 and 31 outliers (resp.); and for the large (T ) dependent variable, yielded τ T estimates of .80 (SE = .05), and .84 (SE = .05), with no outliers. Then the local edge kernel model yielded an estimate of the causal effect of large (versus small) class as τ = τ Y / τ T = .82/.80 = 1.03, and the local linear uniform model yielded a causal effect estimate of τ = τ Y / τ T = .39/.84 = .46, with approximate standard errors (resp.) .23/.80 = .29 and .27/.84 = .32. For all the linear models, the average causal effects varied widely, and there were many outliers for Y .