Dirichlet Process Mixture Models for Regression Discontinuity Designs
Federico Ricciardi Silvia Liverani Gianluca Baio
stat.AP stat.ME





Dirichlet Process Mixture Models for Regression Discontinuity Designs
    Federico Ricciardi^1, Silvia Liverani^2,3, Gianluca Baio^1
        ^1 Department of Statistical Sciences, University College London, London WC1E 6BT, UK

    ^2 School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, UK

^3 The Alan Turing Institute, The British Library, London NW1 2DB, UK 


===================================================================================================================================================================================================================================================================






The Regression Discontinuity Design (RDD) is a quasi-experimental design that estimates the causal effect of a treatment when its assignment is defined by a threshold value for a continuous assignment variable. The RDD assumes that subjects with measurements within a bandwidth around the threshold belong to a common population, so that the threshold can be seen as a randomising device assigning treatment to those falling just above the threshold and withholding it from those who fall just below. 

Bandwidth selection represents a compelling decision for the RDD analysis as the results may be highly sensitive to its choice. A number of methods to select the optimal bandwidth, mainly originating from the econometric literature, have been proposed. However, their use in practice is limited. 

We propose a methodology that, tackling the problem from an applied point of view, consider units' exchangeability, i.e., their similarity with respect to measured covariates, as the main criteria to select subjects for the analysis, irrespectively of their distance from the threshold. We carry out clustering on the sample using a Dirichlet process mixture model to identify balanced and homogeneous clusters. Our proposal exploits the posterior similarity matrix, which contains the pairwise probabilities that two observations are allocated to the same cluster in the MCMC sample. Thus we include in the RDD analysis only those clusters for which we have stronger evidence of exchangeability. 

We illustrate the validity of our methodology with both a simulated experiment and a motivating example on the effect of statins to lower cholesterol level, using UK primary care data.

















§ INTRODUCTION



The Regression Discontinuity Design (RDD) is a quasi-experimental design that estimates the causal effects of a treatment by exploiting the presence of a pre-determined treatment rule (either naturally occurring or regulated by on-going policies). The first publication on RDD was an application in education by <cit.>. Since then this framework has proved to be effective in a wide range of applications in other disciplines, including economics <cit.> and politics <cit.>. 
More recently there has been some interest in the RDD for epidemiology <cit.> and health and primary care applications <cit.>.


The RDD can be applied in any context in which a particular treatment or intervention is administered according to a pre-specified rule linked to a continuous variable, referred to as the `assignment' or `forcing' variable: the treatment is then administered if the units' value for the assignment variable (X) lies above or below a certain threshold (x_0), depending on the nature of the treatment. 
If thresholds are strictly adhered to when assigning treatment, the design is termed sharp, while when this is not the case it is termed fuzzy.

The regression discontinuity design has become of particular interest in the definition of public health policies as it enables the use of routinely collected electronic medical records to evaluate the effects of drugs when these are prescribed according to well-defined decision rules. This is useful as government agencies such as the Food and Drug Administration (FDA) in the USA and the National Institute for Health and Care Excellence (NICE) in the UK are increasingly relying on guidelines for drug prescription in primary care. In fact we will use prescription of statins in the UK as our motivating example, but there is a wide range of potential applications including the prescription of anti-hypertensive drugs when systolic blood pressure exceeds 140mmHg or initiating antiretroviral therapy in patients with HIV-1 when their CD4 count has fallen to 350 cells/mm3 or below.





The RDD can mimic a randomised experiment around the threshold and the treatment effect at the threshold can be obtained averaging the outcomes in `small' bins in its proximity. 
The choice of the `bandwidth' is an important decision for an RDD analysis. because the results are highly sensitive to its choice, especially in all those cases in which the relationship between the assignment variable and the outcome, on both sides of the threshold, deviates from linearity.



In many applied studies <cit.>, a standard strategy adopted to address the bandwidth issue is to produce local linear regression estimates obtained using data within a limited number of bandwidths (often not more than 3 or 4, sometimes defined with the guidance of experts in the field of study). 
Alternatively, more complex approaches can be adopted. 

Historically, these methods find their roots in the econometric literature and have close connection with the non-parametric estimation of the effect for RDD. Their common rationale is that the `optimal' bandwidth must be selected according to some criteria aimed to minimise an error term. The first proposal, by <cit.>, was based on a leave-one-out cross validation (CV) strategy in order to find the estimator minimising the mean integrated square error. Later, <cit.> and <cit.> demonstrated that the CV method was a potential source of bias and that it was not reliable in any case when the design is fuzzy, and hence devised two slightly different minimisation methods based on the asymptotic mean square error. <cit.> give an overview of these approaches.

More recently, Local Randomization (LR) has been proposed by <cit.> and used since in several applied papers <cit.> in an attempt to select a window around the threshold where the units can be seen as part of a randomised experiment. This approach, although motivated by a different intuition, shares a 

common trait with the other approaches outlined above (and further described in Section <ref>): they all aim at finding one bandwidth, having optimal properties under certain criteria and then use it within the RDD framework. As a consequence, they rely on what we named `all-or-nothing' selection mechanism: all units within the bandwidth are considered for the RDD analysis, but none of those outside.

In this paper, we propose an alternative approach to select the units to be included in a RDD analysis.
Similarly to the LR method, our approach originates from a pragmatic and applied point of view, focusing on units' exchangeability, an attribute rooted in the unconfoundness assumption

that guarantees that a RDD mimics a randomised control trial thanks to the similarity of the units above and below the threshold. 

However, our proposal has a more ambitious goal:

not only do we aim at including units for the RDD analysis based on their mutual similarity and not on their proximity to the threshold, but we also want to overcome the need of an `all-or-nothing' approach shared by all other methods existing in the literature. 


Our novel proposal is motivated by 
the idea that that units can be grouped in an unknown yet finite number of clusters in which the available covariates are balanced among units above and below the threshold. Using a Dirichlet process mixture model (DPMM), we cluster the units using continuous and categorical covariates to account for potential sources of confounding. By quantifying the internal similarity of the clusters obtained, only units belonging to the most homogeneous clusters are then used in the RDD analysis, irrespective of their distance from the threshold. Our proposal aims to a more effective sample selection, as it searches for `signal' in the data in farther regions from the threshold generally overlooked by the currently available bandwidth selection approaches and discards the `noise' from data points closer to the cut-off.






The paper is organised as follows. Section <ref> introduces the RDD and gives details about the Bayesian modelling framework we adopt for the analysis.
Section <ref> gives an overview of the current literature on bandwidth selection for regression discontinuity designs. 
Section <ref> presents the methodological core of the paper, where we discuss the use of clustering based on Dirichlet Process Mixture Models (DPMM) within the RDD framework and Section <ref> addresses the issue on cluster selection for the subsequent RDD analysis. Results on both a simulated experiment and a real dataset on the effect of statins on cholesterol level are given in Section <ref>. Finally a closing discussion is presented in Section <ref>.




§ BAYESIAN INFERENCE FOR THE REGRESSION DISCONTINUITY DESIGN



In this section we introduce the basic framework and notation for the RDD. 
Our work is motivated by an application of the regression discontinuity design to statin prescription in primary care. 

In the past years other works from our broader research group have originated from the same practical application and data, every time exploring a different aspect of the RDD <cit.>.
In the UK, according to guidelines given by the National Institution for Health and Care Excellence (NICE), statins must be prescribed to patients whose 10-year risk score of developing a cardiovascular disease, 
predicted using a logistic regression model with a number of clinical and lifestyle indicators as independent variables, exceeds 20% <cit.>. This threshold has been revised in 2014, lowering it to 10%, but we used pre-2014 data in this work and hence we applied the old cut-off value. 

Using the risk score as our forcing variable X ∈{0,1}, a RDD analysis can assess whether statins treatment (T) can cause a reduction in Low-Density Lipoprotein (LDL) cholesterol (our outcome, Y), evaluated at the threshold set to x_0=0.20. To complete the basic notation, let X^c=(X-x_0) be the centred assignment variable and Z be the binary threshold indicator such that Z=1 if the forcing variable X ≥ x_0 and Z=0 otherwise.
Note that Z coincides with the treatment assignment variable T when the design is sharp, but when RDD is applied to health and medical data it is reasonable to expect the design to be fuzzy, and hence the two variables not to coincide. In our motivating example this can be due both to GPs not adhering to NICE guidelines and to patients failing to take statins although prescribed to do so.  

It is widely known that the threshold indicator Z is a special case of binary Instrumental Variable (IV) <cit.>. For this reason, 
in order for the RDD analysis to be performed, a set of assumptions which can be derived from the IV literature must hold <cit.>. 


While further theoretical and technical aspects of the RDD would add very little to the scope of this paper, being extensively covered in the literature 

<cit.>, 


we make use of the next subsection to provide a more detailed overview of the Bayesian modelling framework we aim to use for the the estimation of the causal effect at the threshold.



 §.§ The causal effect







 
Motivated by our example, where GPs' prescribing behaviour may not adhere to NICE guideline, our primary focus is on fuzzy designs, hence the effect we are interested in is the Local Average Treatment Effect (LATE) at the threshold, defined as 

    = E(Y|Z=1) - E(Y|Z=0)E(T|Z=1) - E(T|Z=0).


The LATE numerator is equal to the Average Treatment Effect (ATE).

The denominator, obtained as the difference in the expected treatment probabilities above and below the threshold, scales the ATE to account for the fuzziness of the design. 
In our motivating example, the LATE quantifies the change in LDL cholesterol at the 10-year risk threshold of 20%.




More details about the assumptions that allow the identification of the above effect under a fuzzy observational regime can be found in <cit.>. 




  §.§.§ Models for the ATE

Let the index l ∈{a, b} specify whether a unit's forcing variable value lies above or below the threshold. We decided to model the outcome, i.e., LDL cholesterol, separately for l=a and l=b as

    y_il∼ N(μ_il, σ^2);  
    μ_il = β_0l + β_1l x^c_il,

where x^c_il is the centred distance of variable X from the threshold x_0 for the i-th individual belonging to l. 



In our examples in Section <ref>, both for the simulated scenarios and for the real data analysis, the relatively large sample size reduces the impact on posterior inference of distributional assumptions, especially for σ which is likely dominated by information from observed data. With smaller samples or to ensure further robustness to prior on σ, other models are obviously possible, e.g., by considering an Half-Cauchy distribution <cit.>. 

For the regression parameters, their prior distributions are chosen to reflect plausible LDL cholesterol levels for the observed range of risk scores. Prior specifications are defined as follows













    β_0a    = β_0b + λ;  
    β_0b   ∼ N(3.7; σ^2_0b=0.25);  
    β_1l   ∼ N(0; σ^2_1l=2);   
    σ   ∼(0,5).


To encode in the model some available information from the literature <cit.> about the effect of statins in lowering cholesterol levels, we specify the prior distribution of λ in order to be moderately informative, i.e.,

    λ∼ N(-2,1).




Finally the ATE is calculated as  Δ_β = β_0a- β_0b.


















  §.§.§ Models for the denominator of the LATE

The total number of subjects treated on each side of the threshold is modelled, again separately for l ∈{a, b} as

    ∑_i=1^n_l t_il∼(n_l, π_l),

where n_l is the number of units either above or below the threshold.

Depending on the desired prior structure for (π_b, π_a), we specify two models which, analogously to those in <cit.>, have been named unconstrained and flexible difference model.

This means that for the unconstrained model we use vague Beta distributions

    π_l ∼(1,1),

with l ∈ (a, b).

For the flexible difference model, we impose a mild prior structure acknowledging an actual difference between the treatment probabilities above and below the threshold, defining

    logit(π_a) ∼ N(2, 1)     and     logit(π_b) ∼ N(-2, 1).

These distributions keep the bulk on the prior probability of treatment distributions, above and below the threshold, reasonably separate from one another, limiting the possibility that they result to be similar, while not constraining them to have a fixed difference.

The denominator for the LATE is then given by the difference

    Δ_π = π_a - π_b.






  * * *
  
 
 
 
 






Depending on the chosen model for the denominator we get two different effects, i.e.,

    ^unct = Δ_βΔ_π^unct    and    ^flex = Δ_βΔ_π^flex

for the unconstrained and flexible difference model respectively.





§ A CONCISE REVIEW OF BANDWIDTH SELECTION METHODS



In recent years there has been a surge in the interest of researchers for the choice of the bandwidth, as accounted by <cit.> in their comprehensive review on the topic. In fact the definition of the bandwidth represents a fundamental decision for the RDD as there is both a clear 
link between the size of the bandwidth and the assumption of exchangeability and a trade-off with the precision of the estimates. If the bandwidth is small, units can be reasonably considered more similar to one another. 
If the bandwidth is too large, the converse is true, i.e., units could no longer be considered homogeneous.


In this section, we give an overview of the most prominent methods for neighbourhood selection in the literature.









 §.§ Cross Validation based approach









The first approach found in the literature is based on a Cross Validation procedure as proposed in <cit.>[This is a working paper, later published as peer-reviewed article in a shortened version <cit.>], also discussed 
by <cit.>. 

Let 

    m_h(X_i)=
    α_a + β_a X_i^c,    if X_i ≥ x_0,
    α_b + β_b X_i^c,    if X_i<x_0

be the predicted value, using a bandwidth equal to h, of the outcome Y regressed on the centred assignment variable X_i^c when the i-th unit is left out from the calculation. The Cross Validation criterion is defined as:

    CV_Y,δ(h)= 1N∑_i:q_X,δ,b≤ X_i ≤ q_X, 1-δ,a^N( Y_i - m_h(X_i) )^2.

Here m_h(X_i) is estimated using only observations on one side of X_i

to mimic the fact that RDD estimates are based on regression estimates at the boundary. As a result, equation (<ref>) is an average of boundary prediction errors.

Furthermore q_X,δ,b and q_X,1-δ,a are the δ-th and (1-δ)-th quantiles of the empirical distribution of X for the sub-samples `below' and `above' the threshold, respectively.

<cit.> suggest δ = 0.95 to be appropriate, 
while <cit.> and <cit.> state that δ=0.5 represents a reasonable value, but the choice of an appropriate value varies according to the problem at hand and should be evaluated with care. 


The choice for the bandwidth given by this CV method is then represented by

    h_CV^opt = min_h CV_Y,δ(h).

This criterion leads to the bandwidth choice that minimises an approximation of the Mean Integrated Square Error (MISE): 

    (h)= E [ ∫_x ( m_h(x) - m(x) ) f(x) dx]

where m(x)=E[Y_i|X_i=x] and f(x) is the density of the forcing variable.

In the case of a fuzzy RDD, <cit.> suggest to use the smallest bandwidth selected by two CV criteria applied separately to the outcome and to the treatment:

    h_CV^opt = min( min_h CV_Y,δ(h), min_h CV_T,δ(h) ),

where T denotes the treatment received and the formulation for CV_T,δ(h) is similar to that in (<ref>). 








 §.§ MSE expansion bandwidth selection






Both <cit.> and <cit.> criticise the CV based approach, stating that this criterion relies on fitting the entire regression line between the δ-quantile for the observation on the left and the (1-δ)-quantile for those on the right, so that the result is not optimal for the problem at hand, being the aim of a RDD to estimate the effect at the threshold. 

Let τ be the estimated effect at the threshold for the RDD, the proposal of <cit.> is based on minimising its asymptotic Mean Squared Error (MSE), 
i.e., (τ - τ)^2. Hence the MSE is defined as:

    (h) = E[(τ - τ)^2] = E[((μ_a - μ_a)-(μ_b - μ_b))^2]

where μ_b = lim_x ↑ x_0m_h(x) and μ_a = lim_x ↓ x_0m_h(x), i.e, the two regression estimators for the `true' models on the two sides of the threshold, i.e., μ_b = lim_x ↑ x_0 m(x) and μ_a = lim_x ↓ x_0 m(x).

To overcome some issues arising when trying to minimise the (h) directly, the authors use a first-order approximation around h=0 of the above quantity, which they term Asymptotic Mean Squared Error or (h).  

The optimal bandwidth is therefore:

    h_IK = min_h (h) = C_K (σ^2_a(x_0) + σ^2_b(x_0)f(x_0)(m_a^''(x_0) + m_b^''(x_0))^2)^1/5 N^-1/5

where C_K is a constant value depending on the choice of the kernel function K(·); σ^2_b(x_0) and σ^2_a(x_0) are the left and right limit at the threshold of the variance σ^2(x) = Var(Y_i|X_i=x); f(x) is the density of the forcing variable; m_a^''(x_0) and m_b^''(x_0) are the right and left limits of the second derivative of m(x)=E[Y_i|X_i=x]. 
The authors propose a data-dependent method to estimate h_IK in three steps.



























<cit.> considered that both previous methods produce bandwidths that are too wide, leading to confidence intervals 
with poor asymptotic coverage. The authors prove that correct asymptotic coverage is reached only if the bandwidth can satisfy the bias condition nh_n^5 → 0, a requirement that none of the above mentioned methods can guarantee, 

leading to a first order bias in the distributional approximation. 
As a result, the conventional confidence intervals may substantially over-reject the null hypothesis of no treatment effect.

The authors propose a bias correction to address this problem that is able to improve the performance in finite samples. The final result is a generalisation of h_IK, which we term h_CCT, which allows for higher order polynomial to be used for the inference

and provides more robust confidence interval estimators.










 §.§ Local Randomization


The Local Randomization (LR) approach selects a window around the cutoff in which the randomization assumption is likely to hold <cit.>.

The rationale behind LR is that, because treatment assignment is assumed to be randomised by the threshold inside the window, the distribution of pre-intervention covariates should be the same for treated and untreated units. This observation is directly related
to the non-testable unconfoundness assumption needed for the RDD to infer valid causal estimators. For the RDD framework to be useful, the distribution of these covariates for treated and untreated units should be unaffected by the treatment T within the bandwidth h but should be affected by the treatment outside the window.

To find such desired bandwidth an iterative selection method is implemented. Starting from a arbitrary `small' bandwidth h_1, for each one of the covariates, multiple tests of the null hypothesis of no effect of the treatment on the covariates is conducted and the
minimum p-value taken. 

If the minimum p-value obtained, p_1, is less than some pre-specified level the initial window was too large, hence one should decrease the initial window and start over. Otherwise, if p_1 is greater than the selected significance level, choose a larger window  h_2 ⊇h_1, and go back to calculate a second iteration minimum p-value, p_2. The process continues until the minimum p-value is smaller than the desired level and a final bandwidth h_LR is defined. 





  * * *



The limited literature available and the lack of an unequivocal methodology for the bandwidth selection motivates our work: in the following we develop a more general RDD framework in which the choice of the bandwidth is not required, with positive effect on our results.










































§ DIRICHLET PROCESS MIXTURE MODELS








In this paper, we propose a Dirichlet process mixture model to identify units that are similar (and so will be treated as exchangeable), above and below the threshold. We propose to identify these units by exploiting the characteristics of the clusters obtained with a Dirichlet process mixture model as described in the previous section. 

The Dirichlet process mixture model is a Bayesian nonparametric method for (unsupervised) clustering and applied in a variety of areas, such as retail analysis <cit.>, language processing and classification <cit.>, medical imaging <cit.>, epidemiology <cit.> and genetics <cit.>.

The Dirichlet process was first introduced by <cit.> and is defined as a probability distribution over random probability measures. The distribution of a Dirichlet process is (almost surely) discrete, in that a random sample drawn from a Dirichlet process has a non zero probability that multiple draws will have identical values.  It is this discreteness property which makes the Dirichlet process ideal for clustering, as it avoids the need to determine the number of clusters a priori <cit.>. The basic Dirichlet process mixture model is formulated as follows:

    w_i|θ_i   ∼   p(w_i|θ_i)
    θ_i|G   ∼   G
    
    G   ∼   DP(α,G_0).


The Dirichlet process models the distribution from which data w_1,…, w_n are drawn as a mixture of distributions, p(w_i|θ_i),  where  each  parameter θ_i is  drawn  from  a  mixing distribution G <cit.>.  G_0 is the base distribution, that is  the prior expectation of G, i.e., E[G] =G_0, and the concentration parameter α acts as an inverse variance where larger values of α result in smaller variances. 
Posterior inference from a DPMM utilises Markov chain Monte Carlo (MCMC) posterior simulation and our implementation uses the slice sampling procedure <cit.>. Moreover, due to the nature of the stick-breaking construction of the Dirichlet process <cit.>, label-switching moves are also implemented, to prevent the slice sampler from getting stuck in local modes <cit.>. 

In  this  paper,  the  DPMM  is  implemented to model a mixture of Gaussian and categorical components. 

Let S_i be the latent allocation variable so that if S_i=c then individual i is in cluster 

c ∈{1, C}, then conditional on each cluster c, the likelihood for observable data 𝐃_i = (𝐃^1_i, 𝐃^2_i) is 














    p(𝐃_i|S_i=c,Θ_c) = p(𝐃^1_i|μ^DP_c,Σ_c)p(𝐃^2_i|Φ_c)

where 𝐃^1_i = (D^1_i,1, ..., D^1_i, J_1) is the subset of the J_1 continuous random variables in 𝐃_i and 𝐃^2_i= (D^2_i,1, ..., D^2_i, J_2) is the subset of the J_2 categorical random variables in 𝐃_i. Note that we are assuming independence between continuous and categorical data conditional on the cluster allocations. The cluster specific parameters are given by Θ_c = (μ^DP_c,Σ_c,Φ_c), which are defined in detail below. 

For the continuous random variables, we have 






 

    p(𝐃^1_i|μ_c^DP,Σ_c)=(2π)^-J_1/2|Σ_c|^-1/2exp{-1/2(D^1_i-μ_c^DP)^⊤Σ_c^-1(D^1_i-μ_c^DP)}

and we choose μ_c^DP∼Normal(μ_0^DP,Σ_0) and Σ_c∼InvWishart(R_0,κ_0) (for each c) for our prior model to obtain a conjugate model, permitting Gibbs updates for the parameters μ^DP and Σ. 

For the discrete random variables, we have

    p(𝐃^2_i|Φ_c)=∏_j=1^J_2ϕ_S_i,j,X_i,j.
 
For each individual i, 𝐃^2_i = (D^2_1, …, D^2_J_2) is a vector of J_2 locally independent discrete categorical random variables, where the number of categories for covariate j=1,2,…,J_2 is R_j. 
Then we can write Φ_c=(Φ_c,1,…,Φ_c,J_2) with Φ_c,j = (ϕ_c,j,1, ϕ_c,j,2, …, ϕ_c,j,R_j).
Letting 𝐚=(𝐚_1, 𝐚_2, …, 𝐚_𝐉_2), where 

a_j=(a_j,1,…,a_j,R_j) 

and adopting conjugate Dirichlet priors Φ_c,j∼Dirichlet(a_j), each Φ_c,j can be updated directly using Gibbs iterations. 






As each iteration of the MCMC Gibbs sampler provides an estimate of the cluster labels, Partitioning Around Medoids (PAM) was used to obtain an overall estimate of the optimal number of clusters <cit.>. As the number of clusters varies between iterations, the proposed method uses the posterior 
similarity matrix 𝐏. 
The  best  clustering  is selected by maximising an associated clustering score <cit.>.
The Dirichlet process mixture model described above is available in the R package PReMiuM <cit.>.




















































































§ CLUSTER RANKING AND SELECTION 



Once we have identified units that are similar to one another using a Dirichlet process mixture model, above and below the threshold, we must identify the most suitable clusters for the RDD analysis. We propose to identify clusters that are balanced and homogeneous. These concepts have been extensively exploited in several branches of statistics, most notably by the Propensity Score Weighting literature <cit.>, where overlap in covariates between treatment groups is a desired feature to estimate average treatment effects for sub-populations defined according to the propensity score.
 


A cluster is balanced when it has enough units on both sides of the threshold. As many small clusters are usually fully above, or below, the threshold, it is important to ensure that we consider balanced clusters for the RDD analysis. 
Therefore, we define a value π^Z for the proportion of units with assignment value above the threshold, i.e., X_i ≥ x_0 for i ∈ K_c.

We call π^Z_c the proportion of units in cluster c with Z_i =1, i.e., for which the assignment variable is greater than the threshold x_0.
We the empirically set a constant parameter ζ





 , deeming a cluster balanced if the proportion π^Z_c falls within an acceptable range, i.e., 1ζ≤π^Z_c≤ζ-1ζ. We then discard unbalanced clusters, leaving us with C' ≤ C clusters. 

A cluster is homogeneous (or compact) when the observations within it are very similar to one another. 

However, modelling with a mixture model does not always result in clusters of similar observations. For example, a Gaussian mixture model with a fully flexible covariance matrix may incur in large within-cluster dissimilarities compared to a model in which covariance matrices are assumed to be equal or spherical: observations that are modelled well by a common probability distribution are not necessarily close. For example, in the case of a 2-dimensional Gaussian distribution with a high correlation, the maximum distance between the further observations can be significant. Generally, the mixture model does not come with implicit conditions that ensure the separation of clusters <cit.>. Therefore, we employ the Dirichlet process mixture model to exploit its flexibility, but we must take a close look to the homogeneity of each cluster.

We propose to rank clusters based on their homogeneity. The concept of homogeneity is widely explored in the clustering literature <cit.> and relies on the idea that if properly identified, units in a cluster must have a cohesive structure. The most straightforward way to formalise that all objects within a cluster should be similar to each other is the average within-cluster distance, a commonly used index for cluster internal validation <cit.>. We employ a version of this within-cluster index based on the posterior  
similarity matrix 𝐏 obtained post-processing the output of the Dirichlet process mixture model. The values in 𝐏 are the pairwise probabilities that two observations are allocated to the same clusters in the MCMC sample. 
As such, adapting the definition of dissimilarity from <cit.>, 

we can define a 

similarity function s: 𝒱^2 ⟼ℝ^+_0 so that s(v_1, v_2) = s(v_2, v_1)≥ 0 and s(v_1, v_1) = 1, 
where v_1 and v_2 are elements from the space 𝒱, i.e., the objects needed to be clustered. 





This similarity function can be used to compute the within-cluster homogeneity. 




For c' = 1, …, C' 

let 𝒦'={K_1, …, K_C'} with K_c'⊆𝒱 be the clustering set where the unbalanced clusters have been removed so that n_c' = | K_c'| is the number of units in cluster c'. 

Let p_l,v be the elements of the similarity matrix 𝐏. For each cluster this within-cluster homogeneity index can be calculated as:




    I_c' = 2n_c'(n_c'-1)∑_l=1^n_c'∑_v ≤ l^n_c' p_l,v.


A lower within-cluster index is an indicator of a more homogeneous cluster, with 0 being the minimum value for I_c'. We exploit this measure of homogeneity to rank the clusters from the least homogeneous to the most homogeneous. We relabel the index as I_(c') for c'=1,…,C' such that I_(1) < I_(2) < … < I_(C'). 



Among the balanced clusters, we propose to use homogeneity to select the clusters to include in our model. We propose the following four criteria. 




  * We include clusters until the relative difference between the homogeneity for the c'-th and (c'+1)-th ordered clusters is within a 10% margin, that is, all ordered clusters from 1 to c' such that 

    I_(c'+1)-I_(c')/I_(c') < 0.10

for c'=1,…, C'. We refer to this criteria as inc10.


  * We include the first quartile of the balanced clusters, that is, all clusters c' with

    I_(⌈ h ⌉)     h ≤ C'/4.





We refer to this method as c25.


  * We include clusters starting from the most homogeneous until the sample includes at least half of the units from the entire cohort, that is, all clusters c' with c'=1,…, C' such that 

    ∑_c'=1^C'-1 n_(c')< N/2     ∑_c'=1^C' n_(c')≥ N/2

where n_(c') is the cardinality of the c'-th cluster, ordered according to the homogeneity index I_(s).

We refer to this criteria as n50.


  * We named to this final criteria as n25 as it is similar to n50, but only considering 
one quarter of the units from the entire cohort, that is, all clusters c' with c'=1,…, C' such that 

    ∑_c'=1^C'-1 n_(c')< N/4     ∑_c'=1^C' n_(c')≥ N/4.




The four strategies detailed above define four (possibly) different sub-samples of the partition obtained applying a Dirichlet process mixture model as in Section <ref>. RDD analysis, as detailed in Section <ref>, is hence performed for each of the sub-samples of units irrespective of their distance from the threshold (x_0).




§ APPLICATIONS AND RESULTS








We make use of our methodology for an application to primary care prescription: according to the guidelines given by the National Institute for Health and Care Excellence (NICE) between 2008 and 2014, statins should have been prescribed in the UK to patients with 10-year cardiovascular disease (CVD) risk scores, calculated via the so called Framingham Risk Score <cit.>, in excess of 20%. To illustrate our methodology and check its performance, we use 
statins prescriptions data from The Health Improvement Network (THIN - ) a large primary care database that provides anonymised longitudinal general practice data on patients' diagnostic and prescribing records from more than 500 general practices across the UK. The database is broadly representative of the UK population <cit.>.

In the following Sections we will present results obtained using our methodology both on a realistically simulated dataset (Section 
<ref>) and on a subset of data from THIN patients (Section <ref>). The simulated experiment makes a formal comparison between different methods (i.e., our DPMM clustering based approach and other relevant bandwidth selection criteria), while the real-data application showcases 
how our methodology can be useful in practice. 
In both cases, the values of three key covariates are used to cluster units with our DPMM: age, systolic blood pressure and high-density lipoprotein (HDL) cholesterol. 
With the same data, we have obtained results of RDD analyses using established bandwidth selection methods: those based on MSE (i.e., IK and CCT) and Local Randomization (LR) as detailed in Sections <ref> and <ref> as well as two arbitrarily selected windows, i.e., bandwidth of width 0.05 and 0.1 on each side of the threshold. Appropriate functions from R packages ,  and  are used to estimate h_IK, h_CCT and h_LR respectively.


















































































































































































 §.§ Simulated example




For this example we followed the same approach as <cit.> and used simulated data originated from the THIN database (details about the simulation algorithm can be found on the supplementary material of that paper). In particular data are obtained under a simulation scenario in which the risk score is a strong instrument for the treatment, the treatment effect size is equal to -2 and there is low level confounding. 
Both statins treatment status and the LDL cholesterol outcome are simulated to mimic realistic values.  


We have simulated 100 datasets with N=5,720 units




 and for each of them, separately, we clustered the units using the DPMM approach.
Then we selected the most homogeneous clusters based on the four criteria detailed in Section <ref>. The range for acceptable assignment probability for each cluster is 110≤π^Z_c ≤910, i.e., ζ=9. 




This is done to account for the fact that in most of the clusters the assignment probabilities are not very well balanced (i.e., to prevent a too drastic exclusion of ineligible clusters). 
Finally we performed RDD Bayesian analysis and combined the results to obtain ^unct and ^flex. 














Figure <ref> 





 gives a visual representation of how units are selected according to different bandwidth methods compared with the our DPMM framework combined with the c25 criteria: solid red dots and blue diamonds represent the selected units out of the whole initial sample, represented using void grey markers. Vertical lines show the bandwidths selected with some of the methods described in Section <ref>. Note that LR bandwidths are not shown to avoid confusion, as they are too close to the threshold. 







Table <ref>





  and Figure <ref> show the results of these scenarios. 
It is worth noticing that flexible and unconstrained estimators give very similar results. 
Among the four cluster selection strategies we propose, c25, n50 and n25 all show a reduced bias than those obtained using other established methods - i.e. CCT, IK, LR and arbitrarily-selected fixed-width bandwidths - while results for strategy inc10 are considerably less reliable.
Precision of all estimators is comparable for all strategies but inc10, which shows wider credible intervals. 

























































































 §.§ Real data - Statins prescription in the UK.



In this second example, we considered a subset of patients from THIN: male individuals aged from 50 to 70 who had not previously received a statin prescription nor suffered from  a CVD event and for whom the Framingham risk score was recorded by the GP during the time between 1 January 2007 and 31 December 2008. 
We further restricted the analysis to non-diabetic and non-smoking patients, so that the total number of units is 1386.



Figure <ref> shows why we believe a RDD is appropriate for the data at hand, clearly highlighting a discontinuity at the threshold for both the LDL level and the probability of statins prescription. Also we can appreciate the substantial fuzziness of the data, so that the LATE estimators are appropriate in this setting. 










For the clustering selection process, the range for acceptable assignment probability for each cluster is set to 110≤π^Z_c ≤910, i.e., ζ=9. Table <ref> and Figure <ref> show the results. Obviously there is no real value to compare the results of the estimators with, but there are a few aspects of interest nonetheless. All our DPMM estimators, including inc10, produce similar results, irrespective of which cluster selection method is used, with n50 and n25 strategies both producing more precise estimates. It is also interesting to note how, in this case, the LR method produces very wide credible intervals for both ^unct and ^flex, as this method is not able to pick a large enough subset of similar unit, being constrained to limit the search within nested windows. Results from both MSE based and arbitrary-selected bandwidths appear substantially different: CCT estimators are less precise that the DPMM based ones, and only h=0.10 produces results similar, in median, to those obtained applying our DPMM and cluster selection.




 

































   



§ CONCLUSIONS



We have proposed a novel, data-driven approach to deal with the 
bandwidth selection issue for the regression discontinuity design from a different perspective than those adopted in  the available literature.
Our approach originates from the idea that what matters the most in a regression discontinuity design is the exchangeability of the units included in the analysis, i.e., their homogeneity with respect to know observable characteristics. Following this rationale, it is reasonable to believe that subgroups of units might  share common characteristics irrespective of their distance from the threshold. This, we believe, represents the 

most appealing aspect of this framework: 
instead of relying on the `all-or-nothing' approach, which is implicit with any of the 
currently available bandwidth selection methodologies, we propose a tool that is capable of using all the available information from all the individuals showing homogeneous covariates and balanced forcing variable. 

Furthermore, when compared with the Local Randomization method which is similar in principle to ours, our DPMM clustering approach has the merit of tackling exchangeability more directly: while the former tests the null hypothesis of no effect of the treatment on each observed confounder separately in a univariate way, the latter, relying on clustering methods, evaluates the homogeneity of the considered covariates in a joint, more comprehensive approach.




The results of the RDD analysis using our DPMM clustering framework, especially in combination with c25 cluster selection strategy, compared favourably in terms of bias with those obtained following other bandwidth selection approaches, i.e., CCT and IK (methods that are specifically designed to minimise the bias of the causal estimator ), LR and with the arbitrarily-selected fixed-width bandwidths. 


We are aware of the limitations to our approach. In particular we acknowledge the issue that, due to the complexity of the DPMM which involves the estimation of a latent clustering structure, our analysis is more time consuming than those based on other bandwidth selection methods, an issue that is amplified as the number of clustering covariates increases. Due to label switching and the lack of a specific parameter to target, it is also hard to assess Bayesian DPMM convergence with the usual MCMC diagnostics. We remain convinced that these limitations are a reasonable price to pay in order to be able to overcome the `all-or-nothing' bandwidth approach.
A further limitation is the fact that our method relies on the availability of observed data or known confounders, although this issue is not exclusive of our approach as it is shared with Local Randomization method as well.

As a final remark, we think it is useful to note that we are not advocating the indiscriminate use of our methodology in any given RDD analysis. Expert assessment of any application and a proper evaluation of the plausibility of the RDD assumptions must always constitute the ground for subsequent analyses. Also, availability of covariates data and their role as potential confounders must be assessed beforehand. 
A degree of subjectivity remains in the choice of parameter ζ, for which an assessment of the balance of the forcing variable has been proposed as a way to deal with clusters with unbalanced representation on both sides of the threshold, but the magnitude of the reasonably allowed unbalance represents an application-specific feature and it is left for the practitioner to be determined.


Rather than being a `one-size-fits-all' tool, our proposal offers an alternative approach to identify the units to be included in the RDD analysis in a more targeted way than the bandwidth selection methods currently available. Thoughtful use of our proposed DPMM clustering framework can prove valuable in all those RDD applications where exchangeability is regarded as a key feature and where traditional methods do not offer viable solutions to tackle it.





































§ ACKNOWLEDGEMENTS

This research has been funded by a UK MRC grant MR/K014838/1. Approval for this study was obtained from the Scientific Review Committee of the THIN in August 2014.

plainnat

