Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
60,862 characters · 11 sections · 0 citation commands
Regression Discontinuity Design with Multiple Groups for Heterogeneous Causal Effect Estimation
\thispagestyle{empty}
\pagestyle{plain} \setcounter{page}{1}
Regression discontinuity (RD) design originated in Thistlethwaite and Campbell (1960) that study the effect of the student scholarships on future academic outcomes. In RD design, for evaluation of the intervention of which status is determined by whether a covariate exceed a fixed known threshold or not, subjects with values just below the threshold and those above the threshold are compared, where the intervention status is as good as randomly assigned.
RD design is well applied by empirical researchers to estimate the treatment effect at the target population, similar to other quasi-experimental methods. Applications of RD design are found in various empirical fields in economics such as labor, public, education, and development economics. Detailed literature survey is found in Imbens and Lemieux (2008), and in Lee and Lemieux (2010).
As there have been numerous empirical applications of RD design, more methodological and theoretical extensions are suggested in different directions such as the case of fuzzy discontinuity by Hahn et al. (2001), as for the selection of bandwidth (Ludwig and Miller, 2007; Imbens and Kalyanaraman, 2012; Calonico et al., 2014; Arai and Ichimura, 2018), and for different tests for estimations (Lee, 2008).
Our goal in this paper is to propose a new method for estimation of counterfactual functions and heterogeneous causal effects considering a RD design with multiple groups which have different thresholds. In the standard RD designs, one of the serious limitations is that the intervention effect only at the discontinuity point is evaluable. Angrist and Rokkanen (2015) proposed a method for identification of the causal effects away from the cutoff. However, their approach is that the running variable is assumed to be ignorable if conditional on the other available predictors, and it is different from our attempt to estimate the counterfactual functions themselves. To consider what kind of assumption and estimation method are needed is an important task in this research. In addition, we provide a method for optimization of threshold as an application of our method. Most past studies have considered thresholds as a fixed value and not dealt with threshold itself as an object of study. However, the real interest of researchers should lie not only in evaluation of past interventions under a given threshold but also in an appropriate threshold setting as a support for decision-making for future interventions. Therefore, we develop a method to estimate an optimal threshold in terms of cost effectiveness.
Our methodological development is closely related to the multiple thresholds in RD method. Indeed, the empirical literature utilizes the standard RD design with single threshold. However, it is not uncommon to have multiple thresholds in actual datasets. We often observe multiple thresholds to assign one treatment in a target population. For example, it is often the case that local governments determine the cutoff value for running variables such as test scores, poverty indexes, birth weight, geolocation, and income. When the different administrative districts set each unique threshold of admission test score, it leads multiple thresholds exist in the target population (Lucas and Mbiti, 2014). Similarly, the geographical division often sets own eligible cutoff value for social welfare programs (Crost et al., 2014). In Japan, age limits of the local goverments' programs to make medical expenses for children free vary by the local governments. In this way, countless situations are applicable for multiple thresholds, while RD application is merely concentrated on single threshold method.
There is scarce methodological literature that handles multiple threshold situations. The past literature that deal with multiple thresholds is Papay et al. (2011). Papay et al. (2011) shows how to incorporate multiple dimensions of running variables in the RD design with single dataset, which is different from our model setup with multiple datasets. Literature has also moved to the situations where thresholds or cutoff points were unknown for researchers (Henderson et al., 2014; Porter and Yu, 2015; Chiou et al., 2018). Our method is clearly differed from their works as we assume the situation where value of cutoff is observed from datasets.
This paper is organized as follows. In Section 2 we describe the standard RD design settings as basics of the proposed method. In addition we provides the details of our design that using a special structure that there are multiple groups with different thresholds makes it possible to estimate counterfactuals and causal effects. Section 3 we propose a new AIPW kernel estimator making the best use of the observed data in our design. In Section 4 we investigate the asymptotic properties of the estimator proposed in Section 4 and show its double-robustness. In Section 5 we provide a method to estimate an optimal bandwidth as an application of our method. In Section 6 we report a simulation for studying the properties of the proposed estimator in the finite sample. In Section 7 we summarize this paper and discuss future outlook on this research.
In this section, we briefly summarize framework and theory of the conventional regression discontinuity design. Then we extend the discussion to the case with multiple thresholds and propose our method to estimate the unobservable counterfactuals in the conventional RD designs and heterogeneous causal effects by using them. In this paper consider only the situation where there is just two groups for simplicity. Our discussion and notation are based on Imbens and Lemieux (2008) and modern literature using Rubin Causal Model (RCM) setup with a concept of potential outcomes (Rubin, 1974; Holland, 1986).
As is the usual case with RCM, consider the situation that there are two types of interventions, special intervention (i.e. treatment) and normal intervention (i.e. control), and researchers are interested in the causal effect of the intervention. Corresponding to those two types of interventions, there are two potential outcomes for each unit. Denote by $Y_{ji}$ the potential outcomes of unit $i\in N$, where $N=\{1,...,n\}$ is a set of $n$ units, and the potential outcome for treatment is $Y_{1i}$ and the potential outcome for control is $Y_{0i}$.
Now let the intervention assignment indicator of unit $i$ denote $Z_i\in\{0,1\}$, which is 1 when unit is exposed to treatment and 0 otherwise. The observed outcome variable can be expressed as
In addition, let a finite dimensional vector of pretreatment covariate variables except $X_i$ denote $\bm{W_i} \in \mathbb{R}^m$.
In the setting of RD designs, the type of intervention allocated to unit $i$ is determined by whether a running variable $X$ is above a threshold $c$. RD designs are generally divided into two types, the sharp RD (SRD) design and the fuzzy RD designs depending on how to determine the assignment of intervention. In this study we limit the discssion to the sharp RD design. In the sharp RD design the assignment $Z_i$ is based on a deterministic function of the running variable $X_i$ defined as
Under this function, all the units observing $X_i$ above $c$ are exposed to treatment and the others are exposed to control.
In the sharp RD design, although the running variable $X_i$ does not overlap between the treatment group and the control group, the assignment $Z_i$ is only depending on $X_i$, therefore Missing at random (MAR) (Rubin, 1976), that is,
is satisfied.
Under MAR, if the models of $E(Y_1|X)$ and $E(Y_0|X)$ are parametric, $E(Y_1|X)$ can be extrapolated even below the threshold and $E(Y_0|X)$ also can be extrapolated above the threshold. Therefore, $E[Y_1-Y_0|X=a]$ at any arbitrary point $X=a$ can be estimated and $E[Y_1-Y_0]$ also can be. However, nonparametric regression do not permit extrapolation and only $E(Y_1|X)$ for $X>c$ and $E(Y_0|X)$ for $X<c$ and the difference of those at the discontinuity point ,that is, the local average treatment effect (LATE)
can be estimated. This is the main goal in the RD designs. However only the treatment group can include units who observe $X_i=c$ and the control group cannot, hence the conditional expectation of the observed outcomes $Y_i$ given $X_i$ is discontinuous at $c$. Thus $\tau_{SRD}$ can be regarded as
and obtained by point estimations of the limits from the left and right.
RD design is useful in many practical cases, however it is one of the major limitations that only LATE at the discontinuity point can be estimated and thus the result may lack generalizability (Lee and Lemieux, 2010). This problem is due to the structure that there is no overlap in $X_i$ between the treatment and control groups and the counterfactual cannot be obtained. To solve this problem at least partially, we propose a new method when different datasets with different thresholds are available.
In this paper, to estimate the unobserved potential outcome in the standard RD design (i.e. counterfactual), we consider the RD designs with multiple groups which have different thresholds. We assume the case where the same intervention is provided to several groups (e.g. geographical regions) and those groups have different thresholds from each other on a same running variable and the types of intervention for units are determined by the thresholds of the groups to which they belong. Other basic settings are the same as the case with the standard RD design described in the previous part; there are two types of intervention, or treatment and control, and corresponding to those interventions there are two potential outcomes, and we focus only on the sharp RD design.
In the following, we consider only the case with two groups. Each unit belongs to either of the two groups. Let the group assignment indicator for unit $i \in N$, where $N=\{N_0,N_1\}=\{1,...,n_0,n_0+1,...,n_0+n_1\}$, $N_0=\{1,...,n_0\}$ and $N_1=\{1,...,n_1\}$, be denoted by $D_i\in\{0,1\}$, which takes 0 if $i\in N_0$ and takes 1 if $i\in N_1$. In addition, $c_k(k=0,1; c_0<c_1)$ denotes the thresholds of the two groups, $c_0$ is the one in the group of $N_0$ and $c_1$ is the other. By using the subscript $d_i \in \{0,1\}$ representing the group to which unit $ i $ belongs, the function of intervention assignment is
According to this function, observable outcomes for unit $i$ from $N_0$ are $Y_{0i}$ for $X_i \leq c_0$ and $Y_{1i}$ for $X_i>c_0$, and for unit $i$ from $N_1$, $Y_{0i}$ for $X_i\leq c_1$ and $Y_{1i}$ for $X_i>c_1$ are obsereved. Thus, different potential outcomes are observed depending on the groups for $c_0<X_i<c_1$, while $Y_{0i}$ for $X_i<c_0$ and $Y_{1i}$ for $X_i>c_1$ are commonly observed from both of the two groups, as shown in the Figure(ref).
Now, let the conditional expectation functions given $X_i$ depending on the group assignment be denoted by
This expression allows the regression function to be different by the group assignment, however we are not interested in the individual functions for each group. Our main interest lies in the functions in the target common population:
Especially, between the two thresholds, both of the potential outcomes $Y_1$ and $Y_0$ are observed and thus it should be potentially possible to estimate $g_1(x)$ and $g_0(x)$ for $c_0<X<c_1$, which overlap each other. If we can estimate them, we can also estimate the average treatment effects at arbitrary points between the two thresholds defined as
as shown in the Figure(ref).
Nevertheless what can be estimated from the data is only function ((ref)) and we cannot estimate function ((ref)) directly. The conditional expectation ((ref)) can be rewritten as
and if we knew all factors of the right hand side in equation ((ref)), we could estimate function ((ref)) following equation ((ref)). However, the observed potential outcome is limited as described above, what can estimate directly from the data are only
and the other parts cannot be estimated directly in general, as shown in Figure (ref).
Therefore, whereas $g_0(x)$ for $x<c_0$ and $g_1(x)$ for $x>c_1$ can be estimated according to the equations ((ref)), $g_{00}(x)$ and $g_{11}(x)$ for $c_0<x<c_1$ cannot be estimated and thus we cannot estimate $g_0(x)$ and $g_1(x)$ between the thresholds of most interest. In what follows, we consider what kind of assumption is necessary to realize unbiased estimation of the function $g_j(x)$.
First, consider the most optimistic situation, where the group is randomly assigned for units and the estimated functions are independent of the data assignment. In this case, $E[Y_j|X=x,D=k]=E[Y_j|X=x]$ holds and using only one data from observed group of either $D=0$ or $D=1$ does not generate bias. One of the situations in which this assumption holds is where a type of randomized controlled trial (RCT) can be conducted, where units are randomly distributed to two groups with different thresholds. However, in the field of medicine or social science such as ecnomics, there are not many situations where it is possible to implement random assignment for structural or ethical reasons.
In the following, we investigate the case where the conditions that group assignment is randomly determined and the conditional expectation functions do not depend on the group assignment are not satisfied; that is,
This means that there is a selection bias between the two groups. In this case, a standard approach can cause biased estimates. In order to achieve the unbiased estimator of the conditional expectation functions, we additionally assume ignorability.
This assumption also can be rewritten in another way using Bayes' theorem.
In this form it can be interpreted as meaning that given the covariates $X$ and $W$ the simultaneous distribution of $Y_0$ and $Y_1$ is independent of the group assignment $D$.
Under this assumption, the conditional expectations satisfy
However, when the covariates $W$ is high-dimensional, as in many cases, correct identification of parametric function form is mostly impracticable, furthermore, if using nonparametric regression including the local linear kernel regression, practitioners are faced with the problem known as the Curse of Dimensionality\footnote{The Curse of Dimensionality is the phenomena that the amount of data required for estimation increases exponentially when there are many explanatory variables (Hoshino, 2009). More specifically, let $d$ denote the number of dimension, then asymptotic mean squared error is proportional to $N^{-4/(d+4)}$ (H\"{a}rdle et al., 2004).}. To avoid these problems, we introduce the propensity score.
The propensity score is the concept proposed by Rosenbaum and Rubin(1983) that enables covariate adjustment through a single variable into which the information of multiple covariate variables is aggregated; it is the coarsest one-dimensional balancing score\footnote{A balancing score $b(x)$ is a function of observed covariates $x$ such that the conditional distribution of $x$ given $b(x)$ is independent of assignments $z$; that is,
Balancing scores are not uniquely determined but various functions of $x$. The coarsest balancing score, i.e. the propensity score, is the function of any other balancing scores (Rosenbaum and Rubin, 1983).}. In general propensity score analysis, a selection probability of a missing in the context of missing data analysis or a treatment assignment in the context of causal inference is usually used as a propensity score. In this study, on the other hands, since what determines which of the potential outcome $Y_j$ is the group assignment, the selection probability of $D$ given the covariates $X$ and $W$ is regarded as the propensity score. Under the ignorability assumption ((ref)), we can estimate the conditional expectations as
The specific procedure to estimate as above is described in the next section.
Estimation in the conventional RD designs has been considered as nonparametric estimation problem since the misspecification of the function form may cause bias in estimation of the causal effect (Hahn et al., 2001; Lee and Lemieux, 2010). Therefore we consider nonparametric estimation of $g_j(x)$, in particular, using the local linear regression model taking advantage of the fact that $X$ is one dimensional variable. Note that considering that the purpose of this research is estimation of counterfactual between the two thresholds and estimation of causal effect using it, it is sufficient to estimate even a regression function between thresholds. However, if the estimation target is limited to the interval between the thresholds, the bad boundary behavior of the kernel regression as above occur in the neighborhood of the thresholds. Since data exist outside the thresholds in this design, we use them to improve the stability of estimation; the target of estimation is not limited to the interval between the thresholds.
Now consider a nonparametric regression model $Y_i=g(X_i)+\varepsilon_i$, where $g(x)$ is a unknown smooth function. The local linear estimates of is $g(x)$ formed by minimizing
where $K_h(X_i-x)=K(X_i-x/h)/h$ is the kernel weight with bandwidth $h$ and $\alpha \equiv (\alpha_0(x),\alpha_1(x))^T$, $G(X_i-x)\equiv(1,X_i-x)^T$. The estimated function is $\hat{g}(x)=\hat{\alpha}_0(x)$. When complete data exists, $\bm{\alpha}$ solving the equation ((ref)) gives the correct regression function; however, actually, the presence of missing due to the design in our study make it biased in general.
When focusing on estimate of $E(Y_0|X)$, data of $D=1$ is complete case for $X<c_1$. Consistent estimation of $E(Y_0|X)$ for $X<c_1$ can be implemented using data of $D=1$ and the inversed probability weighted (IPW) method or augmented inversed probability weighted (AIPW) method propsed by Wang et al.(2010), which is more robust than IPW. It is similar for estimate of $E(Y_1|X)$ for $X>c_0$ and data of $D=0$. However, those estimation method ignore the other data ($D=0$ for $g_0(x)$ or $D=1$ for $g_1(x)$) except in estimation of the selection probability model although those data are available. Especially, the observed data of $D=0$ for $X<c_0$ ((c) in Figure (ref)) and $D=1$ for $X>c_1$ ((d) in Figure (ref)) including both the auxiliary variables and even the outcome can be used to estimate in the neighborhood of the thresholds $c_0$ and $c_1$, but the information of those is totally ignored. Those methods are not efficient in this respect. Therefore we propose more efficient method which is capable of exploiting the information from even (c) or (d) in Figure (ref).
We develop a new estimation method for the design of this study based on the AIPW kernel regression proposed by Wang et al. (2010).
As mentioned in the previous section, we consider covariate adjustment using propensity score under the ignorability assumption ((ref)) in order to implement unbiased estimation. Let $\pi_i=Pr(D_i=1|X_i,W_i)$ denote the data selection probability as propensity score and we assume a parametric model:
where $\bm{\gamma}$ is a finite dimensional parameter vector. This model can be specified as logit model or probit model, for example, and we estimate $\hat{\pi}_i=\pi(X_i,W_i;\hat{\bm{\gamma}})$ using $\hat{\bm{\gamma}}$, the maximum likelihood estimate of $\bm{\gamma}$. By weighting the units by the inverse of the estimated $\hat{\pi}_i$ or the inverse of the true selection probability $\pi_i$, if known, we obtain a inversed probability weighted (IPW) estimator.
Denote by $\delta_j(X_i,W_i)$ an arbitrary regression function of $X_i$ and $W_i$. To estimate $\delta_j(X_i,W_i)$ we postulate a parametric model
where $\eta_j$ is a finite dimensional parameter vector. We can estimate $\hat{\delta}_j(X_i,W_i;\hat{\eta}_j)$ by using $\hat{\eta}_j$, the estimate of $\eta_j$ obtained by the standard method such as OLS and by using data satisfying $Z=j$; $\hat{\eta}_0$ is estimated from the part as shown as (a) and (c) in Figure (ref) and $\hat{\eta}_1$ is estimated from (b) and (d).
Now we define the estimating equation for $g_0(\cdot)$ as
where
and for $g_1(\cdot)$ as
where
with $\alpha^j=(\alpha_0^j(x),\alpha_1^j (x))$ solving equation ((ref)) or equation ((ref)) is the local linear estimator of $g_j(x)$, $V_{ji}=V[G(X_i-x)^T \alpha^j;\zeta_j]$ with a known working variance function $V(\cdot, \cdot)$ and an unknown finite dimensional parameter $\zeta_j$. The consistency of the estimation is guaranteed even if $V_j$ is arbitrarily decided under certain conditions (Hoshino, 2009). If we estimate $\zeta_0$ based on the data, we can use the inverse probability weighted moment equations $\sum_{l=1}^n D_l\hat{\pi}_l^ {-1} V_{0l}^{(1)} \left[\left\{Y_l - \hat{\alpha}^0_{0,l} ( \zeta_0 ) \right\}^2 - V \left\{ \hat{\alpha}^0_{0,l} ( \zeta_0 ) , \zeta_0 \right\} \right] = 0$, where $V _ l ^ { ( 1 ) } = \partial V \left\{ \hat { \alpha}^0_{0,l} (\zeta_0 ) ; \zeta_0 \right\} / \partial \zeta_0 $, and $\hat { \alpha } _ { l } ( \zeta_0 ) = \left\{ \hat { \alpha } _ { 0 , l } ( \zeta_0 ) , \hat { \alpha } _ { 1 , l } ( \zeta_0 ) \right\} ^ { T }$ solve ((ref)) with $x = X _ { l } , l = 1 , \dots , n$. We can estimate $\zeta_1$ in a similar way. The estimated conditional expectation function is $\hat{g}_j(x)=\hat{\alpha}_0^j(x)$. The first term of equation ((ref)) and ((ref)) is what constitutes the IPW estimation equation as $\sum U_{IPW,i}^j(\alpha^j)=0$ and the second term $A_i^j(\alpha^j)$ is called an augmented term.
We inevestigate properties of the estimation equations focusing on for $g_0(\cdot)$. These equations allow us to use data of $D=0$ in addition to $D=1$. When $D_i=1$, the first terms in the right hand side in equation ((ref)) and ((ref)) are left and the scond terms are equal to 0, and thus this estimating equation is equal to the one proposed by Wang et al. (2010). When $D_i=0$, the second terms are left and the first terms are equal to 0. For unit $i\in N_0$, $Z_i$ differs depending on either $X_i\leq c_0$ or $X_i>c_0$. If $X_i\leq c_0$, i.e. $Z_i=0$, since complete data including outcomes exists, outcomes and covariates can be included in the estimation as well as $D_i = 1$ in the form changing weight to $1-\hat{\pi}_i$. On the other hand, if $X_i>c_0$, i.e. $Z_i=1$, potential outcome $Y_{0i}$ is regarded as missing but covariates are obtained. In this case, whereas $U_{IPW,i}$ is equal to 0 by $1-Z_i=0$, the augmented term $A_{i}$ is left with weight $-1$. Therefore the information of covariates can be exploited. Now if only data of $D=0$ is used to estimate the parameter $\eta_0$ in $\delta_0(X_i,W_i; \eta_0)$, since the potential outcomes $Y_0$ are obtained only for $X_i\leq c_0$, applying estimated parameters to units satisfying $X_i>c_0$ is an extrapolation and it is not desirable. However in equation ((ref)) units satisfying $D_i=1$ and $X_i>c_0$ are weighted by the selection probability and included in addition to the data of $D_i=0$ and thus it can be interpreted as an interpolation. Figure (ref) shows in what forms units are included in the estimation depending on $D_i$ and $Z_i$.
The estimators solving equation ((ref)) or ((ref)) also have the double-robustness similar to other AIPW estimators including the one proposed by Wang et al.(2010). The estimator is consistent when either of the two following conditions is satisfied (but not necessarily both): (i) the selection probability model is correctly specified, and (ii) the regression function of all covariates is correctly specified. The double-robustness is prooved in the next section.
Appropriate choice of bandwidth is an important issue in kernel regression. The least squares cross validation (LSCV) is one of the most widely used bandwidth selection methods (Li and Racine, 2007). Let $\hat{g}_{0,-i}(X_i)$ and $\hat{g}_{1,-i}(X_i)$ denote the leave-one-out local linear estimator of $g_{0}(X_i)$ and $g_{1}(X_i)$. $\hat{g}_{0,-i}(X_i)$ is the solution in the equation
where
and $\hat{g}_{1,-i}(X_i)$ solves a equation similar to the above. The LSCV method choose the bandwith minimizing a function of bandwidth $h$, denoted as $LSCV_j(h)$, as the optimal bandwidth. $LSCV_0(h)$ and $LSCV_1(h)$ are respectively defined as
and
Therefore the optimal bandwidth is defined as
See Li and Racine (2007) for the mathematical details of the local linear cross validation.
In this section, we describe the asymptotic properties of the estimator proposed in this paper. We can investigate it in a similar way to Wang et al.(2010). Throughout this section we assume the following: (I) $n\rightarrow \infty$, $h\rightarrow 0$, and $nh\rightarrow \infty$: (II) $x$ is in the interior of the support of $X$: (III) the regularity conditions: (i) $g(\cdot)$ and the densitiy function of X, $f_X(\cdot)$ satisfy the smoothness assumptoions of Fan et al. (1996); (ii) the right hand side of the estimating equation are twice continuously differentiable with respect to $\alpha$ at a target point $x$ and second derivatives are uniformly bounded.
The proposed doubly robust (DR) local linear estimator of $g_j(x)$ is $\hat{g}_{j,DR}(x)$ solving equation ((ref)) or ((ref)) and this asymptotic limit is denote by $\tilde{g}_{j,DR}(x)$. The proposed DR kernel estimating equations ((ref)) or ((ref)) should have a sequence of solutions $(\hat{\alpha}^j_{0,DR}(x),\hat{\alpha}^j_{1,DR}(x))$ at $x$ such that as the sample size $n\rightarrow \infty$, and the sequence converges in probability to a vector $(\tilde{\alpha}^j_{0,DR}(x), \tilde{\alpha}^j_{1,DR}(x))$, of which the first component $\tilde{\alpha}^j_{0,DR}(x)$ is denoted by $\tilde{g}_{j,DR}(x)$, and $\tilde{g}_{0,DR}(x)$ satisfies
and $\tilde{g}_{1,DR}(x)$ satisfies
where $\tilde{\pi}=\pi(X_i,W_i; \tilde{\gamma})$ and $\tilde{\gamma}$ is the probability limit of $\hat{\gamma}$, and $\tilde{\delta}_j(X,W)=\delta_j(X,W; \tilde{\eta_j})$ and $\tilde{\eta_j}$ is the probability limit of $\hat{\eta_j}$. Theorem (ref) provides the consistency of the proposed estimator under certain conditions.
Theorem (ref) shows the double-robustness of the proposed estimator as mentioned previously. The proof of Theorem (ref) about $\hat{g}_{0,DR}$ is shown in what follows.
Theorem (ref) about $\hat{g}_{1,DR}$ can be easily proved in a similar way.
Next we invesitgate the asymptotic distribution of the proposed estimator. Theorem (ref) shows the asymptotic bias and variance of the proposed estimator.
Theorem (ref) shows that the asymptotic bias of the proposed estimator is of order $O(h^2)$, and the variance of it is of order $O(1/nh)$, and additionally, it is independent of the working variance $V(\cdot)$ in the proposed DR kernel estimating equations ((ref)). A proof of Theorem (ref) is provided in the Appendix.
The principal aim of this study are expanding the conventional RD design the purpose of which is evaluating the causal effect at the discontinuous point to estimate counterfactual between the two thresholds and to enable evaluation of the causal effect at arbitrary points between the thresholds themselves. In this section, moreover, we propose to estimate optimal thresholds in terms of cost effectiveness as an application of this study. Our position here is to support policy makers' decisions.
In general, it is considered desirable to target as many subjects as possible if special interventions yield better results. However, in practice, special interventions require more costs than regular interventions and the intervention practitioners (e.g. governments or companies) need to bear additional costs. For above reasons, they limit subjects by setting uniform criteria and that is why the RD design is useful in many cases. Considering such background, it is obvious that the question of where to set the threshold to maximize cost performance is one of the most important issues for practitioners. In the web marketing example described above, who pay the expense for the privilege of greater membership or for the coupons are the companies providing such services and it is easy to imagine that they cannot help limiting the target customers due to budgetary reasons. Setting criteria to maximize the return on investment in this example is an important management challenge.
In what follows, we describe how to optimize the threshold using the estimated counterfactuals. Attention should be given to the fact that the following discussion is based on the presupposition that outcomes and cost are measured by the same unit, which is supposed to be money in most cases. Indeed, there are cases where outcomes and costs are variables of different measurement especially in political cases, however this problem has been dealt with in another research area, namely cost benefit analysis. We regard this problem as a issue deviating from the range of this research and do not deal with it here.
We postulate that the optimal thresholds can be estimated by maximization of the function of a threshold $c$ representing the total benefit obtained in the treatment group and the control group minus the additional costs with constraint subject to $c_0<c<c_1$; that is,
where $m(c)$ is a known function of a threshold $c$ that represents the additional cost of treatment and $f_X(x)$ is a probability density function of $X$. The benefits obtained from $X<c_0$ and $X>c_1$ are constant for every $c\in [c_0, c_1]$, thus, practically, we need to consider only maximization of the total benefits and costs across the thresholds. Therefore, the optimal threshold can be defined as
Since it is assumed that the same intervention is performed for all subjects, it is considered reasonable to assume that the additional cost per unit is constant. Therefore, the cost function can be defined as
where $MC(c)$ is the additional cost per unit when threshold is set to $c$. Using this definition and equation ((ref)), the objective function of optimization is
When the intervention providers are beneficiaries at the same time such as the web marketing example mentioned above and $\forall x\in [c_0,c_1]$, $\tau(x)-MC(x)<0$, in other words, $\max\tau(x)<MC(x)$ is satisfied, the objective function is monotonically decreasing and hence the optimal threshold is estimated as $c_{opt}=c_0$. However, this result means that the additional benefit due to the treatment (i.e. the causal effect) is less than the additional cost at any point the intervention does not pay off and implies that the validity of the intervention itself might have to be reviewed from the viewpoint of cost effectiveness.
The practical estimator of the optimal threshold is $\hat{c}_{opt}$ solving equation ((ref)) with $g_j(x)$ replaced by $\hat{g}_j(x)$ estimated in the method proposed in Section 3 and either the true probability density function $f_X(x)$, if known as prior information, or an estimator of it $\hat{f_X}(x)$ estimated by the kernel density estimation, for instance.
In this section, we describe simulation conducted to investigate the properties of the proposed estimator in the finite samples. We evaluate our proposed estimator by comparing it with IPW local linear estimator and the naive local linear estimator. The IPW local linear estimator solves the first terms of equation ((ref)) and ((ref)) $\sum U^j_{IPW,i}=0$ using the data of either $D=0$ or $D=1$. The naive local linear estimator solves equation formed by specification of $\pi$ in the IPW estimating equation to be 1. We generate 100 data sets and estimate using each data set under the following conditions to evaluate the estimators from some viewpoints. First, in order to study the efficiency of the proposed estimator, we performed three types of estimation for each data set using the proposed estimator with all units, the IPW local linear estimator with either of the groups except in the estimation of selection probabilities and the naive local linear estimator with complete case. Note that model specifications here in the IPW and proposed estimator are correct. Next, for the evaluation of the robustness of the estimators, compare the results in the following four cases; (i) the selection probability model in the IPW estimation is incorrect; (ii) the selection probability model in the proposed estimation is incorrect; (iii) the regression model in the proposed estimation is incorrect; (iv) both of the models of $\pi$ and $\delta$ are incorrect. Finally, we examine dependency on the settings of the distribution of the running variable by generating the running variable from either the normal distribution or the log normal distribution. We evaluate the estimation results by comparing mean integrated squared error (MISE) limited to between the two thresholds defined as $\int_{c_0}^{c_1} \{ \hat{g}(x)-g(x)\}f_X(x)dx$. We use the LSCV method to choose the optimal bandwidth as described in Section 4.4.
In what follows, describe the data generating process. The running variable $X$ is generated from a normal distribution with mean 4 and variance $\sigma^2=1.7^2$. We assume that in this simulation the covariates $\boldsymbol{W}$ other than $X$ is 2-dimensional and correlated with $X$ to induce selection bias. Thus we generate $\boldsymbol{W}=(w_1,w_2)^T$ according to a model: $\boldsymbol{W}=\boldsymbol{\eta_0}+\boldsymbol{\eta_1} X+\bm{\xi}$, where $\boldsymbol{\eta_0}$ and $\boldsymbol{\eta_1}$ are $2 \times 1$ parameter vectors and $\boldsymbol{\eta_0}=(-1.5, 2.4)^T$ and $\boldsymbol{\eta_1}=(0.6, 0.4)^T$, $\bm{\xi}=(\xi_1,\xi_2 )^T$ is the disturbance term generated from normal distribution with mean $0$ and $\sigma^2=4$ independently. For the data assignment probability $\pi_i$ we postulate the logit model
where $\gamma_0=0.8$, $\gamma_1=0.5$, $\gamma_2=2$ and $\gamma_3=-0.8$. Then the data assignment indicator $D_i$ is sampled from Bernoulli distribution with probability $\pi_i$. Following $D_i$, the treatment assignment $Z_i$ is determined by the function $Z_i=1\left(X_i>c_{d_i}\right)$, with the lower thresholds $c_0=2$ and the upper thresholds $c_1=6$. Finally we generate the observed outcomes $Y_i=Z_iY_{1i}+(1-Z_i)Y_{0i}$, where
with $(\beta^0_0, \beta^0_1, \beta^0_2, \beta^0_3, \beta^0_4)=(0, 16, -1, 42, 36)$ and $(\beta^1_0, \beta^1_1, \beta^1_2, \beta^1_3, \beta^1_4)=(80, -2, 2, 40, 48)$. Figure (ref) shows a scatter plot of the $(X,Y)$ from one of the generated data sets with the lines that indicate $E(Y_j|X)=E_{W|X}(E(Y_j|X,W))$.
In what follows we report the results when sample size $n=2000$ . Figure (ref) and Table (ref) show the result of the naive, IPW and DR local linear estimators of $g_0(x)$ and $g_1(x)$. Table (ref) summarizes the MISEs of the naive, IPW and DR local linear estimators of $g_0(x)$ and $g_1(x)$ as the performance with correct models. The naive local linear estimates have much larger MISEs than the IPW and AIPW local linear estimates for both of $g_0(x)$ and $g_1(x)$. The DR local linear estimates have smaller MISEs than the IPW local linear estimates. For instance, the DR local linear estimate has approximately 59% gain in MISE efficiency in comparison with the IPW local linear estimate in estimation of $E(Y_0|X)$.
Since it is expected that results for $g_0(x)$ and $g_1(x)$ have similar tendencies from the theory and the results shown in Table (ref), we focus on estimation of $g_0(x)$ in the following simulations. Next consider the case that $\pi$ and/or $\delta$ of the IPW and DR are incorrectly specified as described above. The incorrect model of $\pi$ is specified as the model ((ref)) without the $w_{1}$ term and the incorrect model of $\delta$ is specified as the model ((ref)) without the $X$ squared term. Table (ref) shows the results with incorrectly specified models. The DR estimate with a misspecified $\pi$ has relatively close to the DR estimate with correct models and it is better than the IPW estimate with correct $\pi$. The DR estimate with a misspecified $\delta$ is not as good as the DR estimate with a misspecified $\pi$, however its MISE is still better than the Naive estimate and the IPW estimate with an incorrect $\pi$, and naturally better than the DR estimates when both the model of $\pi$ and $\delta$ are misspecified.
In this paper we proposed a new framework of the regression discontinuity designs for estimation of two conditional expectation functions of potential outcomes, i.e. counterfactuals, between two thresholds by using multiple groups which have difference thresholds. We considered how to realize estimation of them in the two cases with and without selection bias. We showed that we can simply estimate them in the absence of selection bias but cannot generally in the presence of it using the normal estimation method such as the naive local linear regression. In order to estimate consistently and to make the best use of the available data, we proposed the new estimator based on the AIPW kernel estimator with the ignorability assumption. We showed that the proposed estimator has double-robustness and it can exploit the auxiliary information of covariates from even subjects with missing outcomes. In finite samples, the proposed estimator is more efficient compared with the naive local linear estimator and the IPW kernel estimator and has the double-robustness property.
One of the concerns about this study is whether a regression model with a mixture of two data sets as a population is meaningful even if the ignorability ((ref)) is assumed. If we wish to infer the results for a more general population, that is possible when covariates including a running variable are obtained from the more general population and we can assume that which group subjects belong to is determined by the covariates.
In addition, in this paper we have chosen nonparametric regression to estimate conditional expectation functions for some reasons, however parametric regressions are also used in many empirical RD designs. If parametric conditional expectation functions are postulated, counterfactual can be estimated at any point on a running variable by extrapolation with the estimated parameters and covariates from beyond the observation range without using the proposed estimation method.
There is room for the further development of this research. First of all, we should apply the proposed method to real data to confirm its usefulness in empirical cases. As for the theoretical side, in this paper we proposed our method focusing on limited case in some respects. We have considered only the case with two groups, however our method can be extended to cases with three or more groups. Another important topic of future study is an extension to the fuzzy RD design not limited to the sharp RD design.