# Analysis of regression discontinuity designs using censored data

Youngjoo Cho Corresponding author. Email address<EMAIL_ADDRESS>Department of
Mathematical Sciences, The University of Texas at El Paso, El Paso, TX 79912
Chen Hu  Division of Biostatistics and Bioinformatics, Sidney Kimmel
Comprehensive Cancer Center, Johns Hopkins University School of Medicine,
Baltimore, MD 21205  Debashis Ghosh  Department of Biostatistics and
Informatics, University of Colorado Anschutz Medical Campus, Aurora, CO 80045

###### Abstract

In medical settings, treatment assignment may be determined by a clinically
important covariate that predicts patients’ risk of event. There is a class of
methods from the social science literature known as regression discontinuity
(RD) designs that can be used to estimate the treatment effect in this
situation. Under certain assumptions, such an estimand enjoys a causal
interpretation. However, few authors have discussed the use of RD for censored
data. In this paper, we show how to estimate causal effects under the
regression discontinuity design for censored data. The proposed estimation
procedure employs a class of censoring unbiased transformations that includes
inverse probability censored weighting and doubly robust transformation
schemes. Simulation studies demonstrate the utility of the proposed
methodology.

Keywords: Causal effect; Double robustness; Instrumental variable;
Observational studies; Survival analysis.

##  1 Introduction

In observational studies, scientific interest typically focuses on estimation
of a causal estimand. The presence of confounding variables makes its
estimation difficult. To perform causal inference, the analyst typically
relies on several assumptions. One important assumption is the ‘no unmeasured
confounders’ assumption, which implies that treatment assignment is
independent of potential outcome given confounders. This has also been
referred to as the unconfoundedness assumption. However, this assumption is
typically not empirically testable.
A study design that has been used in the social sciences is known as the
regression discontinuity (RD) design. One appealing feature to the RD design
is that the treatment assignment is determined either deterministically or
probabilistically by continuous variable of interest, termed the forcing
variable. It turns out for such a design, the no unmeasured confounders
assumption is not required for inferring causality. The study of RD designs
was initiated by Thistlethwaite and Campbell (1960) and has been developed
further in many subsequent studies. For example, Hahn, Todd and Van der Klaauw
(1999) and Hahn, Todd and Van der Klaauw (2001) prove theoretical results on
RD designs. Ludwig and Miller (2005) and Ludwig and Miller (2007) propose
bandwidth selection procedures for nonparametric estimation with application
to evaluating the effects of funding on educational program. Imbens and
Kalyanaraman (2012). propose an optimal bandwidth selection procedure.
Recently, Calonico, Cattaneo and Titiunik (2014) have developed bias-corrected
nonparametric estimation approaches whose confidence intervals demonstrate
improved coverage relative to those from other RD estimators.
Much of the work in the previous paragraph dealt with the case of uncensored
data. In many biomedical settings, the outcome will be a time to event that is
potentially subject to right censoring. For this problem, there has been
limited work in the area of RD designs. Recently, Bor et al. (2014) and
Moscoe, Bor and Bärnighausen (2015) discuss the use RD designs in medical and
epidemiological studies. In their paper, the goal is to test for the effect of
early versus late treatment initiation to HIV patients on survival.
In practice, RD designs are useful for detection of treatment effect in the
studies with censored data. For example, in the Prostate, Lung, Colorectal and
Ovarian (PLCO) cancer study, the role of prostate cancer screening remains
controversial partially due to the findings of PLCO which failed to show
reduction in prostate cancer-specific survival or overall mortality by a
prostate-specific antigen (PSA) screening strategy (Andriole et al. 2009).
Shoag et al. (2015) consider only the treatment group and propose that a PSA
level $\geq$ 4.0mg/nl is a reasonable cutoff for some outcomes by using RD
design, but they did not consider the censored nature of the outcome variable
in their analysis. It would be of great interest to clinicians for
understanding the screening effect of PSA level 4.0mg/nl on time to death or
to first cancer incidence.
For the use of RD designs in survival analysis, one challenge is the existence
of censoring. Applying standard nonparametric estimation procedures from the
uncensored data setting is not feasible. One way to solve this issue is to use
existing RD estimation procedures with some transformation of the response
that behaves in a manner analogous to the uncensored data case. Fan and
Gijbels (1994) propose local linear regression based on a transformed response
for censored data. Their proposed transformation includes the inverse
probability weighted censoring (IPCW) method, a commonly used method in the
missing data literature to handle censoring. However, this method is
inefficient in that it does not include information of censored observations
in estimation. Rubin and Van der Laan (2007) overcome this difficulty by
proposing a doubly robust transformation of the response which requires
modeling of failure time distribution as well as censoring distribution. This
approach shows promise compared to IPCW methods in a prognostic modeling
(Steingrimsson et al. 2016; Steingrimsson, Diao and Strawderman, 2019), but no
studies have shown its efficiency gains for estimation of parameters in
regression modeling with purpose of inference.
In this paper, we propose a class of estimation procedures in the RD design
for censored data. In Section 2, we review the relevant data structures and
discuss approaches of RD design for uncensored data. Section 3 and 4 describe
the extension of RD designs to censored data as well as laying out the
methodology with attendant asymptotic results. In Section 5, simulation
studies are performed which show the utility of our approach. In Section 6, we
apply our method to the PLCO dataset to test effect of treatment assignment by
PSA. Some discussion concludes Section 7.

##  2 Review of RD design for uncensored data

Before discussing the proposed methodology, we first introduce the RD design
for uncensored data using a potential outcomes framework. Let
$(Y_{*}^{(1)},Y_{*}^{(0)})$ be the potential outcomes under treatment and
control; we use $Z$ to define treatment. We define ${\bf W}$ to be a vector of
forcing variables; for simplicity, we consider only one forcing variable $W$.
As the name suggests, the forcing variable determines the treatment (Imbens
and Lemieux, 2008). Since only one of the potential outcomes is observable,
the observed response is $Y_{*}=ZY_{*}^{(1)}+(1-Z)Y_{*}^{(0)}$. The main
characteristic of the RD design is that the treatment assignment $Z$ depends
on a function of $W$, which can be deterministic or probabilistic. This
corresponds to the sharp and fuzzy RD designs, respectively.
One key assumption in RD designs is that the there is no gaming of the forcing
variable $W$ (McCrary, 2008). This is checked in practice by empirically
plotting the distribution of $W$ and checking to see that there is no
‘clumping’ around the cutoff value of interest. When the forcing variable is
not affected by other variables, we have a ‘locally randomized’ study based on
the cutpoint of the forcing variable from this assumption (Lee and Lemieux,
2010). This ‘local randomization’ is a compelling feature compared to usual
observational studies and allows for establishing causality as in randomized
experiments (Bor et al. 2014).
In the sharp RD design, treatment assignment is decided by a deterministic
function of forcing variable. Let $H^{*}$ be a known discontinuous function.
Then for the sharp RD, $Z=H^{*}(W)$. The main causal effect of interest is the
average treatment effect at the discontinuity point $w_{0}$. By design, if the
value of the forcing variable is greater than or equal to the cutpoint,
$E\\{Y_{*}^{(1)}|W\\}=E(Y_{*}|W)$. Similarly,
$E\\{Y_{*}^{(0)}|W\\}=E(Y_{*}|W)$ if the value of forcing variable is less
than cutpoint. Since our interest focuses on the causal effect at $w_{0}$,
with a continuity assumption for $E\\{Y_{*}^{(1)}|W\\}$ and
$E\\{Y_{*}^{(0)}|W\\}$, we can identify limits around the threshold (Bor et
al. 2014).

$\displaystyle E[Y_{*}^{(1)}-Y_{*}^{(0)}|W=w_{0}]=\lim_{w\downarrow
w_{0}}E(Y_{*}|W=w)-\lim_{w\uparrow w_{0}}E(Y_{*}|W=w).$ (1)

In the fuzzy RD design, treatment assignment is a probabilistic function of
the forcing variable. There exists a jump in the probability of treatment
assignment at the threshold (Imbens and Lemieux, 2008) but it is less than one
and depends on the forcing variable. Let $q_{1}$ and $q_{2}$ be monotone
functions of the forcing variable. For fuzzy RD designs, the probability of
receiving treatment is

$\displaystyle P(Z=1|W)=\begin{cases}q_{1}(W)\ \text{if}\ W<w_{0}\\\ q_{2}(W)\
\text{if}\ W\geq w_{0}.\end{cases}$

By contrast, for sharp RD designs, the probability of receiving treatment
given that forcing variable is greater than certain cutoff is one. Hence the
treatment status and treatment assignment are the same. In fuzzy RD designs,

$\displaystyle\lim_{w\downarrow w_{0}}P(Z=1|W=w)\neq\lim_{w\uparrow
w_{0}}P(Z=1|W=w)$

with the difference of $\lim_{w\downarrow w_{0}}P(Z=1|W=w)$ and
$\lim_{w\uparrow w_{0}}P(Z=1|W=w)$ not being equal to one. Hence, for fuzzy
RD, the treatment assignment is not equivalent to the treatment status. One
can consider the sharp RD as a special case of the fuzzy RD when
$\lim_{w\downarrow w_{0}}P(Z=1|W=w)-\lim_{w\uparrow w_{0}}P(Z=1|W=w)=1$. Now
expression (2.1) is not the average treatment effect in a fuzzy RD design.
However, since the forcing variable determines treatment assignment, it
effectively functions as an instrumental variable. By using the arguments in
Angrist, Imbens and Rubin (1996), we can obtain the so-called complier average
treatment effect. This effect is the main identifiable causal estimand in the
fuzzy RD design and is important in practice because participants in the study
may not comply with the initial treatment assignment. Formally, we have

$\displaystyle\frac{\lim_{w\downarrow w_{0}}E(Y_{*}|W=w)-\lim_{w\uparrow
w_{0}}E(Y_{*}|W=w)}{\lim_{w\downarrow w_{0}}P(Z|W=w)-\lim_{w\uparrow
w_{0}}P(Z|W=w)}=E[Y_{*}^{(1)}-Y_{*}^{(0)}|W=w_{0},\text{subject is a
complier}]$

In this case, the unconfoundedness assumption is not realistic because people
with similar values of the forcing variable away from the cutpoint might
receive different treatment. These two subjects will thus not be comparable
(Imbens and Lemieux, 2008).

##  3 Extension of RD designs to censored data: data structure and assumptions

As in the uncensored data case, we define $(T^{(1)},T^{(0)})$ as potential
outcomes under treatment and control assignments, respectively. In survival
data, $(T^{(1)},T^{(0)})$ are potential times to event for the treatment and
control. Without censoring, time to failure $T$ is only observable for either
of treatment and control group, i.e. $T=ZT^{(1)}+(1-Z)T^{(0)}$. Let
$X=(T,W,Z)$ as full data. Let $C$ be time to censoring. We can only observe
$\tilde{T}=T\wedge C$, $\Delta=I(T\leq C)$. Define $O=(\tilde{T},\Delta,Z,W)$.
The observable data are
$O_{i}=(\tilde{T}_{i},W_{i},Z_{i},\Delta_{i}),i=1,\ldots,n$. We assume that
the observed data are independent and identically distributed. We apply a
logarithm transform to the response. Define

$\displaystyle Y^{(1)}=\log T^{(1)}\quad Y^{(0)}=\log T^{(0)}\quad Y=\log
T\quad\tilde{Y}=Y\wedge\log C$

and $Y_{i}$ and $\tilde{Y}_{i}$ are individual realizations of $Y$ and
$\tilde{Y}$, respectively. Usually, the censoring time is not affected by
treatment assignment and status so that it is reasonable to assume that
$C^{(1)}=C^{(0)}=C$ (Bai, Tsiatis, and O’Brien, 2013), where
$(C^{(1)},C^{(0)})$ are the potential outcomes for censoring under treatment
and control. Let $w_{0}$ be a cutpoint for the forcing variable. We now
discuss necessary assumptions for identification of causal effects in RD
designs with censored data.

  1. 1.

$E(Y^{(1)}|W=w)$ and $E(Y^{(0)}|W=w)$ are continuous in $w$.

  2. 2.

$Y^{(0)},Y^{(1)},W,Z$ are independent with $C$.

  3. 3.

Let $S^{(1)}(t|w)=P(T^{(1)}>t|W=w)$ and $S^{(0)}(t|w)=P(T^{(0)}>t|W=w)$.
$S^{(1)}(t|w)$ and $S^{(0)}(t|w)$ are continuous at $w$ for all $t$. Denote
$S(t|w)=P(T>t|W=w)$.

  4. 4.

Let $G(t)=P(C>t)$. $G(t)$ is continuous for all $t$.

  5. 5.

For the fuzzy RD, $P(Z=1|W=w)$ is continuous at $w$ except for $w=w_{0}$ and
$\lim_{w\downarrow w_{0}}P(Z=1|W=w)\neq\lim_{w\uparrow w_{0}}P(Z=1|W=w)$

  6. 6.

For the fuzzy RD, $Z(w^{*})$ is nondecreasing in $w^{*}$ at $w^{*}=w_{0}$
where $Z(w^{*})$ is potential treatment status given point $w^{*}$ such that
$w^{*}$ is neighborhood of $w_{0}$.

Condition 1 is a smoothness assumption for the mean potential outcome
functions around the cutoff point $w_{0}$. Condition 2 is a typical assumption
in survival analysis and is referred to as random censoring. Conditions 3 and
4 guarantee smoothness for the failure time and censoring distributions.
Condition 5 is a standard assumption for fuzzy RD. Condition 6 states that
potential treatment status is monotonic in the cutoff point (Imbens and
Lemieux, 2008). For survival data, in the sharp RD design, the average causal
effect of treatment given the forcing variable is $E(Y^{(1)}-Y^{(0)}|W=w)$.
Under these assumptions,

$\displaystyle E\\{Y^{(0)}|W=w_{0}\\}=\lim_{w\uparrow
w_{0}}E\\{Y^{(0)}|Z=0,W=w_{0}\\}=\lim_{w\uparrow w_{0}}E(Y|W=w)$
$\displaystyle E\\{Y^{(1)}|W=w_{0}\\}=\lim_{w\downarrow
w_{0}}E\\{Y^{(1)}|Z=1,W=w_{0}\\}=\lim_{w\downarrow w_{0}}E(Y|W=w)$

Then under a sharp RD, the average treatment effect is

$\displaystyle\tau_{SRD}=\lim_{w\downarrow w_{0}}E[Y|W=w]-\lim_{w\uparrow
w_{0}}E[Y|W=w]$

For a fuzzy RD, the average treatment effect is

$\displaystyle\tau_{FRD}=\frac{\lim_{w\downarrow
w_{0}}E(Y|W=w)-\lim_{w\uparrow w_{0}}E(Y|W=w)}{\lim_{w\downarrow
w_{0}}P(Z=1|W=w)-\lim_{w\uparrow w_{0}}P(Z=1|W=w)}$

##  4 Proposed Methodology

###  4.1 Censoring unbiased transformations

In sharp RD with uncensored data, we can use nonparametric estimation directly
with the response variable. However, for censored data, it is very difficult
to use response directly due to censoring. The issue is to find a
transformation $q^{*}$ based on observed data such that
$E\\{q^{*}(O)\\}=E(H(T)|W)$ where $H(\cdot)$ is a nondecreasing function of
$T$. The transformation is called a censoring unbiased transformation (Fan and
Gijbels, 1994; Rubin and Van der Laan, 2007; Steingrimsson, Diao, and
Strawderman, 2019). One approach is to use an inverse probability censoring
weighted (IPCW) method to obtain $E(H(T)|W)$, which is

$\displaystyle H_{IPCW}(T)=\frac{\Delta H(T)}{G(T|W)}.$

It is easy to show that $H_{IPCW}(T)$ is censoring unbiased transformation.
However, this approach requires that the censoring distribution is correctly
specified and this approach yields an inefficient estimator. Rubin and Van der
Laan (2007) propose a doubly robust (henceforth DR) censoring unbiased
transformation for the failure time. Steingrimsson, Diao, and Strawderman
(2019) extend this approach to $H(T)$. Let

$\displaystyle M_{G}(u|W)=I(\tilde{T}\leq
u,\Delta=0)-\int_{0}^{u}I(\tilde{T}\geq s)\lambda_{G}(s|W)ds,$

where $\lambda_{G}(s|W)$ is true hazard function of $G$ given covariate $W$.
The form is

$\displaystyle H_{DR}(T)=\frac{\Delta
H(T)}{G(T|W)}+\int_{0}^{\tilde{T}}\frac{Q_{H(T)}(u,W)}{G(u|W)}dM_{G}(u|W).$

where $Q_{H(T)}(u,W)$ is $E\\{H(T)|T>u,W\\}$. This DR transformation requires
estimation of both censoring and failure time distributions. It is a
combination of an IPCW term and a mean zero martingale transform term. IPCW
term is a special case of DR transformation. This martingale transformation
term utilizes information from censored observations, which yields greater
efficiency than using the IPCW approach only. Since it is also a censoring
unbiased transformation, it guarantees that $E(H_{DR}(T)|W=w)=E(H(T)|W=w)$ for
any $w$ (Steingrimsson, Diao, and Strawderman, 2019). Moreover, if either the
censoring distribution or failure time distribution is correctly specified,
then the estimator from DR transformation based on the estimated $Q$ and $G$
is statistically consistent for $E(H(T)|W)$. Moreover, if models for $Q$ and
$G$ are both correctly specified, then the resulting estimator from the DR
transformation is the most efficient estimator given class of the estimator of
$E(H(T)|W)$ (Steingrimsson, Diao, and Strawderman, 2019).

###  4.2 Asymptotic Theory

The key idea of using the transformation from Section 4.1. is to be able to
use uncensored data techniques in the censored situation by applying the data-
dependent transformation to the observed outcome. As can be seen in the last
section, these transformations depend on the censoring distribution or both
the censoring and failure time distributions. In this section, we discuss RD-
based estimation procedures of the treatment effect and their asymptotic
properties with censored data. Here, we take $H(\cdot)=\log(\cdot)$. and

$\displaystyle Y_{IPCW,i}(O_{i};G)=\frac{\Delta_{i}Y_{i}}{G(T_{i})}\quad
Y_{DR,i}(O_{i};G,S)=\frac{\Delta_{i}Y_{i}}{G(T_{i})}+\int_{0}^{\tilde{T}_{i}}\frac{Q_{Y}(u,W_{i})}{G(u)}dM_{G,i}(u),$

where $Q_{Y}(\cdot,\cdot)=E(\log T|T\geq\cdot,\cdot)$ and

$\displaystyle M_{G,i}(u)=I(\tilde{T}_{i}\leq
u,\Delta_{i}=0)-\int_{0}^{u}I(\tilde{T}_{i}\geq s)\lambda_{G}(s)ds.$

We only focus on DR transformation. To obtain limits in the sharp and fuzzy RD
designs, parametric methods are not very attractive because modeling and the
discontinuity point(s) depend on the particular parametric distribution. Hahn,
Todd and Van der Klaauw (2001) show that a causal effect in RD designs is
nonparametrically estimable. Now we apply the local linear regression method
of Fan and Gijbels (1996). Let $K(\cdot)$ be a kernel function and $h$ be
bandwidth. For fuzzy RD, we consider

$\displaystyle
U_{R,DR}^{FRD,Y}(\alpha_{R}^{Y},\beta_{R}^{Y};G,S)=\sum_{i=1}^{n}I(W_{i}\geq
w_{0})\\{Y_{DR,i}(O_{i};G,S)-\alpha_{R}^{(Y)}-\beta_{R}^{(Y)}(W_{i}-w_{0})\\}^{2}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}$
$\displaystyle
U_{L,DR}^{FRD,Y}(\alpha_{L}^{Y},\beta_{L}^{Y};G,S)=\sum_{i=1}^{n}I(W_{i}<w_{0})\\{Y_{DR,i}(O_{i};G,S)-\alpha_{L}^{(Y)}-\beta_{L}^{(Y)}(W_{i}-w_{0})\\}^{2}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}$
$\displaystyle
U_{R}^{FRD,Z}(\alpha_{R}^{Z},\beta_{R}^{Z})=\sum_{i=1}^{n}I(W_{i}\geq
w_{0})\\{Z_{i}-\alpha_{R}^{Z}-\beta_{R}^{Z}(W_{i}-w_{0})\\}^{2}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}$
$\displaystyle
U_{L}^{FRD,Z}(\alpha_{L}^{Z},\beta_{L}^{Z})=\sum_{i=1}^{n}I(W_{i}<w_{0})\\{Z_{i}-\alpha_{L}^{Z}-\beta_{L}^{Z}(W_{i}-w_{0})\\}^{2}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}.$

We can similarly define loss functions for IPCW transformation, say
$U_{R,IPCW}^{FRD,Y}(\alpha_{R}^{Y},\beta_{R}^{Y};G)$ and
$U_{L,IPCW}^{FRD,Y}(\alpha_{L}^{Y},\beta_{L}^{Y};G)$. In sharp RD designs,
since the treatment assignment is deterministic, the loss functions with
transformation responses are only necessary and they are identical to ones in
fuzzy RD design.
Let $\\{\hat{\alpha}_{R,DR}^{FRD,Y}(G,S)$, $\hat{\beta}_{R,DR}^{FRD,Y}(G,S)$,
$\hat{\alpha}_{L,DR}^{FRD,Y}(G,S)$, $\hat{\beta}_{L,DR}^{FRD,Y}(G,S)\\}$ be
estimators using DR transformation in fuzzy RD design, respectively.
Furthermore, we define $\\{\hat{\alpha}_{R}^{FRD,Z}$,
$\hat{\beta}_{R}^{FRD,Z}$, $\hat{\alpha}_{L}^{FRD,Z}$,
$\hat{\beta}_{L}^{FRD,Z}\\}$ for the estimators of modeling treatment
assignment. We can similarly define estimators for sharp RD designs, say
$\\{\hat{\alpha}_{R,DR}^{SRD,Y}(G,S)$, $\hat{\beta}_{R,DR}^{SRD,Y}(G,S)$,
$\hat{\alpha}_{L,DR}^{SRD,Y}(G,S)$, $\hat{\beta}_{L,DR}^{SRD,Y}(G,S)\\}$. Then
we can derive estimators for the fuzzy and sharp RD designs:

$\displaystyle\hat{\tau}_{FRD}^{DR}(G,S)=\frac{\hat{\alpha}_{R,DR}^{FRD,Y}(G)-\hat{\alpha}_{L,DR}^{FRD,Y}(G)}{\hat{\alpha}_{R}^{FRD,Z}-\hat{\alpha}_{L}^{FRD,Z}}$
$\displaystyle\hat{\tau}_{SRD}^{DR}(G,S)=\hat{\alpha}_{R,DR}^{SRD,Y}(G,S)-\hat{\alpha}_{L,DR}^{SRD,Y}(G,S).$

Note that we removed dependence on the bandwidth in the definition of these
estimators; we thus assume that it is given. We will discuss how to estimate
bandwidth formally in Section 4.3. In the next set of results, we show that
under certain conditions, we can prove asymptotic convergence results for the
sharp and fuzzy RD estimators. As can be seen, our estimators depend on $G$
and $S$. Let $G_{0}$ and $S_{0}$ be true distributions of failure and
censoring times, respectively. Let $\hat{G}$ and $\hat{S}$ be estimated
distributions of failure and censoring times, respectively. In the estimation,
we correctly estimate censoring distribution while we may incorrectly estimate
survival distribution. As discussed in the Supplementary Materials, we assume
uniform consistency of $\hat{G}$ to $G_{0}$, and uniform consistency of
$\hat{S}$ to $S^{*}$ where $S^{*}$ is possibly incorrect model of $S$. We will
discuss the estimation of these two distributions in the next subsection. For
these two theoretical results, the IPCW and DR estimators are asymptotically
normal with some bias.

######  Theorem 1.

Assume that conditions (C1)-(C5), (R1)-(R9) hold. By Lemma 1-6 in the
Supplementary Materials,

$\displaystyle
n^{2/5}(\hat{\tau}_{FRD}^{IPCW}(\hat{G})-\tau_{FRD}-\varphi_{FRD})\overset{d}{\longrightarrow}N(0,\Sigma_{FRD}^{IPCW}(G_{0}))$
$\displaystyle
n^{2/5}(\hat{\tau}_{FRD}^{DR}(\hat{G},\hat{S})-\tau_{FRD}-\varphi_{FRD})\overset{d}{\longrightarrow}N(0,\Sigma_{FRD}^{DR}(G_{0},S^{*}))$

where $\varphi_{FRD}$, $\Sigma_{FRD}^{IPCW}(G)$ and $\Sigma_{FRD}^{DR}(G,S)$
are defined in the Supplementary Materials.

For the sharp RD case, the result follows easily from Theorem 1 because there
is no need to model $Z|W$ for fuzzy RD.

######  Corollary 1.

Suppose that conditions (C1)-(C5) and (R1)-(R9) in the Supplementary Materials
hold. By Lemma 1-6 and Theorem 1 in the Supplementary Materials,

$\displaystyle
n^{2/5}(\hat{\tau}_{SRD}^{IPCW}(\hat{G})-\tau_{SRD}-\varphi_{SRD})\overset{d}{\longrightarrow}N(0,\Sigma_{SRD}^{IPCW}(G_{0}))$
$\displaystyle
n^{2/5}(\hat{\tau}_{SRD}^{DR}(\hat{G},\hat{S})-\tau_{SRD}-\varphi_{SRD})\overset{d}{\longrightarrow}N(0,\Sigma_{SRD}^{DR}(G_{0},S^{*}))$

where $\varphi_{SRD}$, $\Sigma_{SRD}^{IPCW}(G)$ and $\Sigma_{SRD}^{DR}(G,S)$
are defined in the Supplementary Materials.

Now we demonstrate an efficiency result similar to that given in
Steingrimsson, Diao, and Strawderman (2019) for a separate problem. For the
sharp RD estimator, due to the ‘local randomization’ result, the DR estimator
with true censoring and failure time distributions is more efficient estimator
than the estimators with IPCW and DR transformations, which only involve true
censoring distribution. For the fuzzy RD estimator, to obtain a similar
result, we should account for correlation between the transformed responses
and $Z$. Suppose that

$\displaystyle\tau^{Y}=\lim_{w\downarrow w_{0}}E(Y|W=w)-\lim_{w\uparrow
w_{0}}E(Y|W=w)\quad\tau^{Z}=\lim_{w\downarrow w_{0}}P(Z=1|W=w)-\lim_{w\uparrow
w_{0}}P(Z=1|W=w)$

We need two conditions to hold: i) $\tau^{Y}$ and $\tau^{Z}$ are positive ii)
$Cov(Y_{DR}(O;G_{0},S_{0}),Z)$ is smaller than $Cov(Y_{IPCW}(O;G_{0}),Z)$ and
$Cov(Y_{DR}(O;G_{0},S),Z)$ at the boundary points. Let

$\displaystyle\eta_{DR}^{+}(w_{0};G,S)=\lim_{w\downarrow
w_{0}}Cov(Y_{DR}(O;G,S),Z|W=w)\quad\eta_{DR}^{-}(w_{0};G,S)=\lim_{w\uparrow
w_{0}}Cov(Y_{DR}(O;G,S),Z|W=w)$
$\displaystyle\eta_{IPCW}^{+}(w_{0};G)=\lim_{w\downarrow
w_{0}}Cov(Y_{IPCW}(O;G),Z|W=w)\quad\eta_{IPCW}^{-}(w_{0};G)=\lim_{w\uparrow
w_{0}}Cov(Y_{IPCW}(O;G),Z|W=w).$

We formally state these in the theorem below:

######  Theorem 2.

Suppose that conditions (C1)-(C5) and (R1)-(R9) in the Supplementary Materials
hold. Let $AVar$ denote asymptotic variance. Then for the sharp RD estimator,

$\displaystyle
AVar(\hat{\tau}_{SRD}^{DR}(G_{0},S_{0}))\leq\min\\{AVar(\hat{\tau}_{SRD}^{IPCW}(G_{0})),AVar(\hat{\tau}_{SRD}^{DR}(G_{0},S))\\}.$

For fuzzy RD estimator, suppose that $\tau^{Y}>0$ and $\tau^{Z}>0$, and if

$\displaystyle\eta_{DR}^{+}(w_{0};G_{0},S_{0})\leq\min\\{\eta_{IPCW}^{+}(w_{0};G_{0}),\eta_{DR}^{+}(w_{0};G_{0},S)\\}$
$\displaystyle\eta_{DR}^{-}(w_{0};G_{0},S_{0})\leq\min\\{\eta_{IPCW}^{-}(w_{0};G_{0}),\eta_{DR}^{-}(w_{0};G_{0},S)\\},$

then

$\displaystyle
AVar(\hat{\tau}_{FRD}^{DR}(G_{0},S_{0}))\leq\min\\{AVar(\hat{\tau}_{FRD}^{IPCW}(G_{0})),AVar(\hat{\tau}_{FRD}^{DR}(G_{0},S))\\}.$

###  4.3 Bandwidth selection

Once we estimate the censoring unbiased transformation with respect to failure
time, we can then apply the existing methodology for RD designs with censored
data. The estimated transformation is

$\displaystyle\hat{Y}_{IPCW,i}(O_{i};\hat{G})=\frac{\Delta_{i}\tilde{Y}_{i}}{\hat{G}(\tilde{T}_{i})}$
$\displaystyle\hat{Y}_{DR,i}(O_{i};\hat{G},\hat{S})=\frac{\Delta_{i}\tilde{Y}_{i}}{\hat{G}(\tilde{T}_{i})}+\int_{0}^{\tilde{T}_{i}}\frac{\hat{Q}_{Y}(u,W_{i})}{\hat{G}(u)}d\hat{M}_{G,i}(u),$

where $\hat{Q}_{Y}(\cdot|\cdot)$ is estimator of $Q_{Y}(\cdot|\cdot)$ and

$\displaystyle\hat{M}_{G,i}(u|W_{i})=I(\tilde{T}_{i}\leq
u,\Delta_{i}=0)-\int_{0}^{u}I(\tilde{T}_{i}\geq s)d\hat{\Lambda}_{G}(s),$

and $\hat{\Lambda}_{G}(\cdot)$ is Nelson-Aalen estimator for $G$. Details of
computation of $\hat{G}$ and $\hat{Q}_{Y}(\cdot,\cdot)$ is found in the
Supplementary Materials.
In standard nonparametric regression, one important issue is bandwidth
selection. Ludwig and Miller (2005) and Ludwig and Miller (2007) propose a
mean squared error (MSE) type cross-validation criterion. Let
$\hat{a}_{L}(W,\xi,L)$ be $\xi$ quantile of the empirical distribution of $W$
using observations $W_{i}<w_{0}$ and let $\hat{a}_{R}(W,1-\xi)$ be $1-\xi$
quantile of the empirical distribution of $W$ using observations $W_{i}\geq
w_{0}$. Moreover, let $\hat{\alpha}_{L}^{(Y)}$ and $\hat{\alpha}_{R}^{(Y)}$ be
estimated parameters for $\alpha_{L}$ and $\alpha_{R}$. Criterion from Ludwig
and Miller (2005) and Ludwig and Miller (2007) for uncensored data is

$\displaystyle\frac{1}{n}\sum_{\hat{a}_{L}(W,\xi)\leq
W_{i}\leq\hat{a}_{R}(W,1-\xi)}(Y_{i}-\hat{\gamma}(W_{i}))^{2},$

where

$\hat{\gamma}(w)=\begin{cases}\hat{\alpha}_{L}(w)&\text{if}\ w<w_{0}\\\
\hat{\alpha}_{R}(w)&\text{if}\ w\geq w_{0}\end{cases}.$

This criterion still works for the unbiased censoring transformation because
it uses a MSE-type criterion and does not depend on variance of $Y$. We now
modify the proposal of Ludwig and Miller (2005) and Ludwig and Miller (2007)
(henceforth LM) for censored data. For the sharp RD estimator, we consider

$\displaystyle(\hat{\alpha}_{R}^{Y}(w;\hat{G},\hat{S}),\hat{\beta}_{R}^{Y}(w;\hat{G},\hat{S}))=\operatorname*{arg\,min}\sum_{i=1}^{n}I(W_{i}\geq
w)\\{Y_{DR,i}(O_{i};\hat{G},\hat{S})-\alpha_{R}^{Y}-\beta_{R}^{Y}(W_{i}-w)\\}^{2}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}$
$\displaystyle(\hat{\alpha}_{L}^{Y}(w;\hat{G},\hat{S}),\hat{\beta}_{L}^{Y}(w;\hat{G},\hat{S}))=\operatorname*{arg\,min}\sum_{i=1}^{n}I(W_{i}<w)\\{Y_{DR,i}(O_{i};\hat{G},\hat{S})-\alpha_{L}^{Y}-\beta_{L}^{Y}(W_{i}-w)\\}^{2}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}.$

Then the LM criterion for sharp RD estimator in censored data is given by

$\displaystyle
CV_{Y_{DR}}(h;\hat{G},\hat{S})=\frac{1}{n}\sum_{\hat{a}_{L}(W,\xi)\leq
W_{i}\leq\hat{a}_{R}(W,1-\xi)}(\hat{Y}_{DR,i}(O_{i};\hat{G},\hat{S})-\hat{\gamma}_{DR}^{Y}(W_{i}))^{2},$

where

$\hat{\gamma}_{DR}^{Y}(w)=\begin{cases}\hat{\alpha}_{L,DR}^{Y}(w;\hat{G},\hat{S})&\text{if}\
w<w_{0}\\\ \hat{\alpha}_{R,DR}^{Y}(w;\hat{G},\hat{S})&\text{if}\ w\geq
w_{0}\end{cases}.$

We then choose
$\hat{h}_{DR}(\hat{G},\hat{S})=\operatorname*{arg\,min}\limits_{h}CV_{Y_{DR}}(h;\hat{G},\hat{S})$.
We can derive a similar quantity for
$\hat{Y}_{IPCW,i}(O_{i};\hat{G}),i=1,\ldots,n$. For the fuzzy RD estimator, we
define

$\displaystyle(\hat{\alpha}_{R}^{Z}(w),\hat{\beta}_{R}^{Z}(w))=\operatorname*{arg\,min}_{\alpha_{R}^{Z},\beta_{R}^{Z}}\sum_{i=1}^{n}I(W_{i}\geq
w)\\{Z_{i}-\alpha_{R}^{Z}-\beta_{R}^{Z}(W_{i}-w)\\}^{2}K\bigg{(}\frac{W_{i}-w}{h}\bigg{)}$
$\displaystyle(\hat{\alpha}_{L}^{Z}(w),\hat{\beta}_{L}^{Z}(w))=\operatorname*{arg\,min}_{\alpha_{L}^{Z},\beta_{L}^{Z}}\sum_{i=1}^{n}I(W_{i}<w)\\{Z_{i}-\alpha_{L}^{Z}-\beta_{L}^{Z}(W_{i}-w)\\}^{2}K\bigg{(}\frac{W_{i}-w}{h}\bigg{)}.$

Then we obtain

$\displaystyle CV_{Z}(h)=\frac{1}{n}\sum_{\hat{a}_{L}(W,\xi)\leq
W_{i}\leq\hat{a}_{R}(W,1-\xi)}(Z_{i}-\hat{\gamma}_{Z}(W_{i}))^{2},$

where

$\hat{\gamma}_{Z}(w)=\begin{cases}\hat{\alpha}_{L}^{Z}(w)&\text{if}\
w<w_{0}\\\ \hat{\alpha}_{R}^{Z}(w)&\text{if}\ w\geq w_{0}\end{cases}.$

Then we obtain $\hat{h}_{Z}=\operatorname*{arg\,min}\limits_{h}CV_{Z}(h)$. A
smaller bandwidth is preferable to reduce the bias of the estimator. Hence for
the fuzzy RD estimator, we consider
$\min\\{\hat{h}_{IPCW}(\hat{G}),\hat{h}_{Z}\\}$ (IPCW) and
$\min\\{\hat{h}_{DR}(\hat{G},\hat{S}),\hat{h}_{Z}\\}$ (DR).

###  4.4 Variance estimation

With the estimation of the censoring and failure time distributions and
bandwidth selection, we obtain

$\displaystyle(\hat{\alpha}_{R,DR,\hat{h}}^{FRD,Y}(\hat{G},\hat{S}),\hat{\beta}_{R,DR,\hat{h}}^{Y}(\hat{G},\hat{S}))=\operatorname*{arg\,min}_{\alpha_{R}^{Y},\beta_{R}^{Y}}U_{R,DR,\hat{h}}^{FRD,Y}(\alpha_{R}^{Y},\beta_{R}^{Y};\hat{G},\hat{S})$
$\displaystyle(\hat{\alpha}_{L,DR,\hat{h}}^{FRD,Y}(\hat{G},\hat{S}),\hat{\beta}_{L,DR,\hat{h}}^{Y}(\hat{G},\hat{S}))=\operatorname*{arg\,min}_{\alpha_{L}^{Y},\beta_{L}^{Y}}U_{L,DR,\hat{h}}^{FRD,Y}(\alpha_{L}^{Y},\beta_{L}^{Y};\hat{G},\hat{S})$
$\displaystyle(\hat{\alpha}_{R,\hat{h}}^{FRD,Z},\hat{\beta}_{R,\hat{h}}^{FRD,Z})=\operatorname*{arg\,min}_{\alpha_{R}^{Z},\beta_{R}^{Z}}U_{R,\hat{h}}^{FRD,Z}(\alpha_{R}^{Z},\beta_{R}^{Z})\quad(\hat{\alpha}_{L}^{Z},\hat{\beta}_{L}^{Z})=\operatorname*{arg\,min}_{\alpha_{L}^{Z},\beta_{L}^{Z}}U_{L,\hat{h}}^{FRD,Z}(\alpha_{L}^{Z},\beta_{L}^{Z}),$

where $U_{R,DR,\hat{h}}^{FRD,Y}$ and $U_{R,\hat{h}}^{FRD,Z}$ correspond to
$U_{R,DR}^{FRD,Y}$ and $U_{R}^{FRD,Z}$ with an estimated bandwidth. We can
define similarly for estimation functions with IPCW transformation. Estimators
for sharp RD designs can be similarly defined. Hence the proposed fuzzy RD and
sharp RD estimators based on $\hat{G}$ and $\hat{S}$ are

$\displaystyle\hat{\tau}_{FRD,\hat{h}}^{DR}(\hat{G},\hat{S})=\frac{\hat{\alpha}_{R,DR,\hat{h}}^{FRD,Y}(\hat{G},\hat{S})-\hat{\alpha}_{L,DR,\hat{h}}^{FRD,Y}(\hat{G},\hat{S})}{\hat{\alpha}_{R,\hat{h}}^{Z}-\hat{\alpha}_{L,\hat{h}}^{Z}}$
$\displaystyle\hat{\tau}_{SRD,\hat{h}}^{DR}(\hat{G},\hat{S})=\hat{\alpha}_{R,DR,\hat{h}}^{SRD,Y}(\hat{G},\hat{S})-\hat{\alpha}_{L,DR,\hat{h}}^{SRD,Y}(\hat{G},\hat{S}).$

where $\\{\hat{\alpha}_{RDR,\hat{h}}^{SRD,Y}(\hat{G},\hat{S})$,
$\hat{\alpha}_{L,DR,\hat{h}}^{SRD,Y}(\hat{G},\hat{S})\\}$ are sharp RD
estimators with estimated bandwidth from DR transformation. For variance
estimation, one may use the methods based on the asymptotic results in Section
4.2. Let $e_{1}=(1,0)^{T}$ and $a(u)=(1,u)^{T}$. Define $g(\cdot)$ to common
density of $W_{i}$. Using the expressions in the Supplementary Materials, the
asymptotic variance of $\hat{\tau}_{SRD}^{DR}(\hat{G},\hat{S})$ can be
expressed as

$\displaystyle\frac{1}{nh}\\{g(w_{0})\\}^{-1}\\{e_{1}^{T}(\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})\mathbf{\Gamma}^{-1}\boldsymbol{\vartheta}\mathbf{\Gamma}^{-1}+\sigma_{DR}^{2-}(w_{0};G_{0},S^{*})\mathbf{\Gamma}^{-1}\boldsymbol{\vartheta}\mathbf{\Gamma}^{-1})e_{1}\\}[1+o_{p}(1)],$

where

$\displaystyle\boldsymbol{\vartheta}=\int_{0}^{\infty}\\{K(u)\\}^{2}a(u)a^{T}(u)\quad\mathbf{\Gamma}=\int_{0}^{\infty}K(u)a(u)a^{T}(u)$
$\displaystyle\sigma_{DR}^{2+}(w_{0};G,S)=\lim_{\epsilon\downarrow
w_{0}}Var(Y_{DR}(O;G,S)|W=w)\quad\sigma_{DR}^{2-}(w_{0};G,S)=\lim_{\epsilon\uparrow
w_{0}}Var(Y_{DR}(O;G,S)|W=w).$

Let $\mathbf{X}_{h}$ be $n\times 2$ matrix with the first column being 1 and
the second column being $\dfrac{W_{i}-w_{0}}{h},i=1,\ldots,n$. Moreover, let
$\mathbf{W}_{h+}$ and $\mathbf{W}_{h-}$ be diagonal matrix with diagonal
elements with $I(W_{1}\geq w_{0})K\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}$ and
$I(W_{1}<w_{0})K\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)},i=1\ldots,n$,
respectively. Define

$\displaystyle\mathbf{\Gamma}_{h+}=\mathbf{X}_{h}^{T}\mathbf{W}_{h+}\mathbf{X}_{h}\quad\mathbf{\Gamma}_{h-}=\mathbf{X}_{h}^{T}\mathbf{W}_{h-}\mathbf{X}_{h}.$

By Calonico, Cattaneo and Titiunik (2014),

$\displaystyle\frac{1}{nh}\\{g(w_{0})\\}^{-1}\\{e_{1}^{T}(\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})\mathbf{\Gamma}^{-1}\boldsymbol{\vartheta}\mathbf{\Gamma}^{-1}+\sigma_{DR}^{2-}(w_{0};G_{0},S^{*})\mathbf{\Gamma}^{-1}\boldsymbol{\vartheta}\mathbf{\Gamma}^{-1})e_{1}\\}[1+o_{p}(1)]$
$\displaystyle=\frac{1}{n}e_{1}^{T}(\boldsymbol{\Gamma}_{h+}^{-1}\boldsymbol{\phi}_{YY+,DR}\boldsymbol{\Gamma}_{h+}^{-1}+\boldsymbol{\Gamma}_{h-}^{-1}\boldsymbol{\phi}_{YY-,DR}\boldsymbol{\Gamma}_{h-}^{-1})e_{1},$

where

$\displaystyle\boldsymbol{\phi}_{YY+,DR}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}\geq
w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}b_{i}b_{i}^{T}\sigma_{DR,(1)}^{2}(W_{i};G_{0},S^{*})$
$\displaystyle\boldsymbol{\phi}_{YY-,DR}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}<w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}b_{i}b_{i}^{T}\sigma_{DR,(0)}^{2}(W_{i};G_{0},S^{*}),$

where $b_{i}=(1,h^{-1}(W_{i}-w_{0}))^{T}$ and

$\displaystyle\sigma_{DR,(1)}^{2}(w;G_{0},S^{*})=Var\bigg{\\{}\frac{\Delta
Y(1)}{G_{0}(T)}+\int_{0}^{\infty}\frac{Q_{Y(1)}(u,W;S^{*})}{G_{0}(u)}dM_{G}(u)\bigg{|}W=w\bigg{\\}}$
$\displaystyle\sigma_{DR,(0)}^{2}(w;G_{0},S^{*})=Var\bigg{\\{}\frac{\Delta
Y(0)}{G_{0}(T)}+\int_{0}^{\infty}\frac{Q_{Y(0)}(u,W;S^{*})}{G_{0}(u)}dM_{G}(u)\bigg{|}W=w\bigg{\\}}.$

where $Q_{Y(k)}(u,W;S^{*})=E_{S^{*}}(Y(k)|W,T\geq u),k=0,1$. We can obtain a
similar result for the IPCW estimator. The first approach is to use plug-in
residuals for the sandwich variance estimator. Let
$\hat{\epsilon}_{Y+,DR,i}=\hat{Y}_{DR,i}(O_{i};\hat{G},\hat{S})-\hat{\alpha}_{R,DR,\hat{h}}^{SRD,Y}(\hat{G},\hat{S})$
and
$\hat{\epsilon}_{Y-,DR,i}=\hat{Y}_{DR,i}(O_{i};\hat{G},\hat{S})-\hat{\alpha}_{L,DR,\hat{h}}^{SRD,Y}(\hat{G},\hat{S})$.
Then one may wish to use

$\displaystyle\hat{\boldsymbol{\phi}}_{YY+,DR,\hat{h}}^{pir}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}\geq
w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}b_{i}b_{i}^{T}\hat{\epsilon}_{Y+,DR,i}^{2}$
$\displaystyle\hat{\boldsymbol{\phi}}_{YY-,DR,\hat{h}}^{pir}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}<w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}b_{i}b_{i}^{T}\hat{\epsilon}_{Y-,DR,i}^{2}.$

The second approach is to use a nonparametric nearest neighborhood (NN)
variance estimator as in Calonico, Cattaneo and Titiunik (2014). This approach
is advantageous in that it does not require nonparametric smoothing and is
robust (Abadie and Imbens, 2006). Mimicking the approach from Abadie and
Imbens (2006) and Calonico, Cattaneo and Titiunik (2014),

$\displaystyle\hat{\boldsymbol{\phi}}_{YY+,DR,\hat{h}}^{rb}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}\geq
w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}b_{i}b_{i}^{T}\hat{\sigma}_{YY+,DR}^{2}(W_{i})$
$\displaystyle\hat{\boldsymbol{\phi}}_{YY-,DR,\hat{h}}^{rb}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}<w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}b_{i}b_{i}^{T}\hat{\sigma}_{YY-,DR}^{2}(W_{i}),$

where

$\displaystyle\hat{\sigma}_{YY+,DR}^{2}(W_{i})=I(W_{i}\geq
w_{0})\frac{K}{K+1}\bigg{(}\hat{Y}_{i}^{DR}(O_{i};\hat{G},\hat{S})-\frac{1}{K}\sum_{k=1}^{K}\hat{Y}_{l_{+,k}(i)}^{DR}(O_{i};\hat{G},\hat{S})\bigg{)}^{2}$
$\displaystyle\hat{\sigma}_{YY-,DR}^{2}(W_{i})=I(W_{i}<w_{0})\frac{K}{K+1}\bigg{(}\hat{Y}_{i}^{trans}(O_{i};\hat{G},\hat{S})-\frac{1}{K}\sum_{k=1}^{K}\hat{Y}_{l_{-,k}(i)}^{DR}(O_{i};\hat{G},\hat{S})\bigg{)}^{2},$

where $\hat{Y}_{l_{+,k}(i)}^{DR}(O_{i};\hat{G},\hat{S})$ is $k$th closest unit
to unit $i$ among $\\{W_{i}:W_{i}\geq w_{0}\\}$ and
$\hat{Y}_{l_{-,k}}^{DR}(O_{i};\hat{G},\hat{S})$ is $k$th closest unit to unit
$i$ among $\\{W_{i}:W_{i}<w_{0}\\}$ for variable
$\hat{Y}_{i}^{DR}(O_{i};\hat{G},\hat{S})$, respectively. Let
$\mathbf{X}_{\hat{h}}$, $\mathbf{W}_{\hat{h}+}$ and $\mathbf{W}_{\hat{h}-}$
are $\mathbf{X}_{h}$, $\mathbf{W}_{h+}$ and $\mathbf{W}_{h-}$ with estimated
bandwidth. Then the variance estimator will be

$\displaystyle\frac{1}{n}e_{1}^{T}(\hat{\boldsymbol{\Gamma}}_{\hat{h}+}^{-1}\hat{\boldsymbol{\phi}}_{YY+,DR,\hat{h}}\hat{\boldsymbol{\Gamma}}_{\hat{h}+}^{-1}+\hat{\boldsymbol{\Gamma}}_{\hat{h}-}^{-1}\hat{\boldsymbol{\phi}}_{YY-,DR,\hat{h}}\hat{\boldsymbol{\Gamma}}_{\hat{h}-}^{-1})e_{1},$

where

$\displaystyle\hat{\boldsymbol{\Gamma}}_{\hat{h}+}=\mathbf{X}_{\hat{h}}^{T}\mathbf{W}_{\hat{h}+}\mathbf{X}_{\hat{h}}\quad\hat{\mathbf{\Gamma}}_{\hat{h}-}=\mathbf{X}_{\hat{h}}^{T}\mathbf{W}_{\hat{h}-}\mathbf{X}_{\hat{h}},$

and $\hat{\boldsymbol{\phi}}_{YY+,DR,\hat{h}}$ is either
$\hat{\boldsymbol{\phi}}_{YY+,DR,\hat{h}}^{pir}$ or
$\hat{\boldsymbol{\phi}}_{YY+,DR,\hat{h}}^{rb}$ and define
$\hat{\boldsymbol{\phi}}_{YY-,DR,\hat{h}}$ similarly. For variance estimation
with the fuzzy RD estimator, by extending the method from the sharp RD design
case, we can calculate $\hat{\boldsymbol{\phi}}_{YY+,DR,\hat{h}}^{pir}$,
$\hat{\boldsymbol{\phi}}_{YY-,DR,\hat{h}}^{pir}$
$\hat{\boldsymbol{\phi}}_{YY+,DR,\hat{h}}^{rb}$,
$\hat{\boldsymbol{\phi}}_{YY-,DR,\hat{h}}^{rb}$. In these calculations,
$\hat{\epsilon}_{Y+,DR,i}$ and $\hat{\epsilon}_{Y-,DR,i}$ are based on
$\hat{\alpha}_{R,DR,\hat{h}}^{FRD,Y}(\hat{G},\hat{S})$ and
$\hat{\alpha}_{L,DR,\hat{h}}^{FRD,Y}(\hat{G},\hat{S})$, resepctively. Now the
covariance term between the transformed response and $Z$ should be reflected
in the estimation. Let
$\hat{\epsilon}_{Z+,i}=Z_{i}-\hat{\alpha}_{R,\hat{h}}^{FRD,Z}$ and
$\hat{\epsilon}_{Z-,i}=Z_{i}-\hat{\alpha}_{L,\hat{h}}^{FRD,Z}$. Then

$\displaystyle\hat{\boldsymbol{\phi}}_{YZ+,DR,\hat{h}}^{rb}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}\geq
w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}b_{i}b_{i}^{T}\hat{\sigma}_{YZ+,DR}^{2}(W_{i})$
$\displaystyle\hat{\boldsymbol{\phi}}_{YZ-,DR,\hat{h}}^{rb}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}<w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}b_{i}b_{i}^{T}\hat{\sigma}_{YZ-,DR}^{2}(W_{i})$
$\displaystyle\hat{\boldsymbol{\phi}}_{ZZ+,\hat{h}}^{rb}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}\geq
w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}b_{i}b_{i}^{T}\hat{\sigma}_{ZZ+}^{2}(W_{i})$
$\displaystyle\hat{\boldsymbol{\phi}}_{ZZ-,\hat{h}}^{rb}=\frac{1}{n}\sum_{i=1}^{n}I(W_{i}<w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}K\bigg{(}\frac{W_{i}-w_{0}}{\hat{h}}\bigg{)}b_{i}b_{i}^{T}\hat{\sigma}_{ZZ-}^{2}(W_{i}),$

where

$\displaystyle\hat{\sigma}_{YZ+,DR}^{2}(W_{i})=I(W_{i}\geq
w_{0})\frac{K}{K+1}\bigg{(}\hat{Y}_{i}^{DR}(O_{i})-\frac{1}{K}\sum_{k=1}^{K}\hat{Y}_{l_{+,k}(i)}^{DR}(O_{i})\bigg{)}\bigg{(}Z_{i}-\frac{1}{K}\sum_{k=1}^{K}Z_{l_{+,k}(i)}\bigg{)}$
$\displaystyle\hat{\sigma}_{YZ-,DR}^{2}(W_{i})=I(W_{i}<w_{0})\frac{K}{K+1}\bigg{(}\hat{Y}_{i}^{DR}(O_{i})-\frac{1}{M}\sum_{k=1}^{K}\hat{Y}_{l_{-,k}(i)}^{DR}(O_{i})\bigg{)}\bigg{(}Z_{i}-\frac{1}{K}\sum_{k=1}^{K}Z_{l_{-,k}(i)}\bigg{)}$
$\displaystyle\hat{\sigma}_{ZZ+}^{2}(W_{i})=I(W_{i}\geq
w_{0})\frac{K}{K+1}\bigg{(}Z_{i}-\frac{1}{K}\sum_{k=1}^{K}Z_{l_{+,k}(i)}\bigg{)}^{2}$
$\displaystyle\hat{\sigma}_{ZZ-}^{2}(W_{i})=I(W_{i}<w_{0})\frac{K}{K+1}\bigg{(}Z_{i}-\frac{1}{K}\sum_{k=1}^{K}Z_{l_{-,k}(i)}\bigg{)}^{2}.$

where $Z_{l_{+,k}(i)}$ and $Z_{l_{-,k}(i)}$ are defined similarly as
$\hat{Y}_{l_{+,k}(i)}^{DR}(O_{i})$ and $\hat{Y}_{l_{-,k}(i)}^{DR}(O_{i})$. For
plug-in approach, in addtion to quantities from sharp RD, we can calculate
$\\{\hat{\boldsymbol{\phi}}_{YZ+,DR,\hat{h}}^{pir},\hat{\boldsymbol{\phi}}_{YZ-,DR,\hat{h}}^{pir},\hat{\boldsymbol{\phi}}_{ZZ+,\hat{h}}^{pir}$,$\hat{\boldsymbol{\phi}}_{ZZ-,\hat{h}}^{pir}$}
similarly by replacing $\\{\hat{\sigma}_{YZ+,DR}^{2}(W_{i})$,
$\hat{\sigma}_{YZ-,DR}^{2}(W_{i}),\hat{\sigma}_{ZZ+}^{2}(W_{i})$,
$\hat{\sigma}_{ZZ-}^{2}(W_{i})\\}$ by
$\\{\hat{\epsilon}_{Y+,DR,i}\hat{\epsilon}_{Z+,i},\hat{\epsilon}_{Y-,DR,i}\hat{\epsilon}_{Z-,i},\hat{\epsilon}_{Z+,i}^{2},\hat{\epsilon}_{Z-,i}^{2}$}
in fuzzy RD.
Let
$\hat{\tau}_{DR}^{Y}(\hat{G},\hat{S})=\hat{\alpha}_{R,DR,\hat{h}}^{FRD,Y}(\hat{G},\hat{S})-\hat{\alpha}_{L,DR,\hat{h}}^{FRD,Y}(\hat{G},\hat{S})$
and
$\hat{\tau}_{Z}=\hat{\alpha}_{R,\hat{h}}^{FRD,Z}-\hat{\alpha}_{L,\hat{h}}^{FRD,Z}$.
Denote $\hat{\tau}_{DR}^{Y}$ to be $\hat{\tau}_{DR}^{Y}(\hat{G},\hat{S})$.
Moreover, $\hat{\boldsymbol{\phi}}_{ZZ+,\hat{h}}$ is either
$\hat{\boldsymbol{\phi}}_{ZZ+,\hat{h}}^{pir}$ or
$\hat{\boldsymbol{\phi}}_{ZZ+,\hat{h}}^{rb}$.
$\hat{\boldsymbol{\phi}}_{ZZ-,\hat{h}}$ is defined similarly. By the plug-in
approach from asymptotic results (see Supplementary Materials),

$\displaystyle\frac{1}{\hat{\tau}_{Z}^{2}}(\hat{V}_{YY+,DR,\hat{h}}+\hat{V}_{YY-,DR,\hat{h}})-\frac{2\hat{\tau}_{DR}^{Y}}{\hat{\tau}_{Z}^{3}}(\hat{V}_{YZ+,DR,\hat{h}}+\hat{V}_{YZ-,DR,\hat{h}})+\frac{(\hat{\tau}_{DR}^{Y})^{2}}{\hat{\tau}_{Z}^{4}}(\hat{V}_{ZZ+,\hat{h}}+\hat{V}_{ZZ-,\hat{h}}),$

where

$\displaystyle\hat{V}_{YY+,DR,\hat{h}}=\frac{1}{n}e_{1}^{T}\hat{\boldsymbol{\Gamma}}_{\hat{h}+}^{-1}\hat{\boldsymbol{\phi}}_{YY+,DR,\hat{h}}\hat{\boldsymbol{\Gamma}}_{\hat{h}+}^{-1}e_{1}\quad\hat{V}_{YY-,DR,\hat{h}}=\frac{1}{n}e_{1}^{T}\hat{\boldsymbol{\Gamma}}_{\hat{h}-}^{-1}\hat{\boldsymbol{\phi}}_{YY-,DR,\hat{h}}\hat{\boldsymbol{\Gamma}}_{\hat{h}-}^{-1}e_{1}$
$\displaystyle\hat{V}_{YZ+,DR,\hat{h}}=\frac{1}{n}e_{1}^{T}\hat{\Gamma}_{\hat{h}+}^{-1}\hat{\boldsymbol{\phi}}_{YZ+,DR,\hat{h}}\hat{\boldsymbol{\Gamma}}_{\hat{h}+}^{-1}e_{1}\quad\hat{V}_{YZ-,DR,\hat{h}}=\frac{1}{n}e_{1}^{T}\hat{\boldsymbol{\Gamma}}_{\hat{h}-}^{-1}\hat{\boldsymbol{\phi}}_{YZ-,DR,\hat{h}}\hat{\boldsymbol{\Gamma}}_{\hat{h}-}^{-1}e_{1}$
$\displaystyle\hat{V}_{ZZ+,\hat{h}}=\frac{1}{n}e_{1}^{T}\hat{\boldsymbol{\Gamma}}_{\hat{h}+}^{-1}\hat{\boldsymbol{\phi}}_{ZZ+,\hat{h}}\hat{\boldsymbol{\Gamma}}_{\hat{h}+}^{-1}e_{1}\quad\hat{V}_{ZZ-}=\frac{1}{n}e_{1}^{T}\hat{\boldsymbol{\Gamma}}_{\hat{h}-}^{-1}\hat{\boldsymbol{\phi}}_{ZZ-,\hat{h}}\hat{\boldsymbol{\Gamma}}_{\hat{h}-}^{-1}e_{1}.$

The nonparametric bootstrap is another method we will consider for standard
error estimation. For implementation of our method, we can use existing
software for uncensored data. The R package `rdboust` (Calonico, Cattaneo and
Titiunik, 2015b) is a powerful tool to perform statistical inference for RD
designs. To implement the proposed methods, one can simply transform the
response by methods in Section 4.3. Then with the transformed quantities as a
new response, we can estimate regression coefficients along with standard
errors.

##  5 Simulation results

We performed simulation studies to evaluate the finite-sample properties of
our proposed estimators. The forcing variable $W$ is generated as a
$Unif(0,1)$ random variable. The error variable is generated as $\epsilon\sim
N(0,0.5)$. Regression coefficients are set to be $\beta_{10}=2$,
$\beta_{20}=1$ and $\beta_{30}=1$. The response is generated from the
following model:

$\displaystyle T=\exp(\beta_{10}+\beta_{20}W+\beta_{30}I(W\geq
0.5)+\epsilon).$

Censoring is generated as a $Unif(0,50)$ random variable that is independent
of $T$ and $W$. Three conditional expectation methods were considered in the
simulation study: Cox model, Log-normal and Log-Logistic model. We use Kaplan-
Meier to estimate $G$. To ensure positivity of $\hat{G}$, we truncate
$\tilde{T}$ by time point $\omega$ where $\omega$ is 95th percentile of
observed time for estimation of $G$ (Steingrimsson et al. 2016). Censoring
distribution is fitted using Kaplan-Meier. Sample sizes are $n=200$ and
$n=400$. The number of bootstraps within each simulation is 50. To select
bandwidth by cross-validation, it is important to select $\xi$ based on the
range of dataset. The value of $\xi$ is 0.5. The amount of observed censoring
is approximately 51% across the simulations. Kernel function is triangular
function, which is

$\displaystyle K(u)=1-|u|.$

Table 1 shows finite-sample properties of the estimator $\hat{\tau}_{SRD}$. In
the columns denoting standard error calculation and coverage, NN, Plug-in and
Boot denote the nearest neighborhood, plug-in residual and bootstrap
approaches, respectively. For coverage, all the calculations are based on the
normal approximation except the Boot_ED, which denotes coverage based on the
empirical distribution of bootstrapped samples. The IPCW approach is more
biased than doubly robust approach. Except for the bootstrap, in general, the
coverage of the estimators satisfies the 95% nominal level. The efficiency
gain of DR approach compared to IPCW approach is noticeable. The performances
of DR approaches across the conditional expectations are very stable. The
results from the DR approach confirms the augmentation theory results from
Tsiatis (2007).
In the simulation study, we consider the fuzzy RD based on a modification of
the simulation setting in Yang (2013). We generate $W\sim Unif(-1,1)$, and let
$V=I(W\geq 0)$. Next, we generate $\kappa\sim N(0,0.25)$ and independent of
all these aforementioned variables. Then treatment variable is defined as
$Z=I(-0.5+V+W+\kappa>0)$. Then we generate $\epsilon\sim N(0,0.25)$ that is
independent of aforementioned variables. Failure time is defined as
$T=\exp(\beta_{10}+\beta_{20}W+\beta_{30}Z+\epsilon)$ where regression
coefficients are set to be $\beta_{10}=2$, $\beta_{20}=1$ and $\beta_{30}=1$.
The censoring variable is generated as $Unif(0,50)$. The average censoring
rate is approximately 39%. In this case, the denominator for the true value is
calculated as

$\displaystyle\lim_{w\downarrow
0}P(Z=1|W=w)=P(\kappa>-0.5)=1-P(\kappa\leq-0.5)=1-\Phi(-2)$
$\displaystyle\lim_{w\uparrow 0}P(Z=1|W=w)=P(\kappa>0.5)=1-P(\kappa\leq
0.5)=1-\Phi(2),$

where $\Phi$ is an inverse function of the standard normal cumulative
distribution function. Hence the denominator should be $\Phi(2)-\Phi(-2)$. For
the numerator,

$\displaystyle\lim_{w\downarrow
0}E\\{\log(T)|W=w\\}=\beta_{10}+\beta_{20}\times 0+\lim_{w\downarrow
0}P(Z=1|W=w)$ $\displaystyle\lim_{w\uparrow
0}E\\{\log(T)|W=w\\}=\beta_{10}+\beta_{20}\times 0+\lim_{w\uparrow
0}P(Z=1|W=w).$

Hence the numerator and denominator are equal so that the average treatment
effect for those who comply with the treatment assignment is 1. We use the
same conditional expectation methods as in the sharp RD case. Table 2 shows
the numerical results for sample sizes $n=250$ and $n=500$. For all
approaches, the bias is greater than those reported in Table 1. This makes
sense because the estimator in fuzzy RD has a denominator which requires
estimation by nonparametric method, which introduces bias. As with the sharp
RD situation, the DR method shows good performance regardless of choice of
conditional expectation. The IPCW method has a larger bias than the DR
methods. The coverage probability tend to perform better in larger sample
sizes. It is interesting to note that the coverage based on the empirical
distribution using the bootstrap works better for fuzzy RD relative to sharp
RD.

|  | Bias | SE | Cover
---|---|---|---|---
|  | NN | Plug-in | Boot | NN | Plug-in | Boot | Boot_ED
$n=200$ | IPCW | 0.074 | 1.422 | 1.390 | 1.542 | 0.900 | 0.894 | 0.916 | 0.912
| Cox | -0.004 | 0.123 | 0.119 | 0.121 | 0.942 | 0.944 | 0.942 | 0.908
| Log-norm | -0.008 | 0.136 | 0.132 | 0.136 | 0.946 | 0.944 | 0.942 | 0.910
| Log-log | -0.008 | 0.136 | 0.132 | 0.139 | 0.942 | 0.942 | 0.946 | 0.912
$n=400$ | IPCW | 0.091 | 1.008 | 0.995 | 1.084 | 0.930 | 0.930 | 0.944 | 0.924
| Cox | -0.003 | 0.086 | 0.085 | 0.084 | 0.940 | 0.932 | 0.938 | 0.898
| Log-norm | -0.004 | 0.095 | 0.093 | 0.094 | 0.922 | 0.920 | 0.936 | 0.892
| Log-log | -0.004 | 0.096 | 0.094 | 0.096 | 0.920 | 0.924 | 0.936 | 0.890
Table 1: Numerical results when sample size $n=200$ and $n=400$ in sharp RD |  | Bias | SE | Cover
---|---|---|---|---
|  | NN | Plug-in | Boot | NN | Plug-in | Boot | Boot_ED
$n=250$ | IPCW | 0.051 | 1.271 | 1.147 | 1.269 | 0.906 | 0.876 | 0.896 | 0.942
| Cox | 0.010 | 0.205 | 0.179 | 0.198 | 0.910 | 0.880 | 0.904 | 0.940
| Log-norm | 0.012 | 0.215 | 0.188 | 0.207 | 0.918 | 0.880 | 0.882 | 0.936
| Log-log | 0.011 | 0.215 | 0.188 | 0.207 | 0.920 | 0.880 | 0.890 | 0.942
$n=500$ | IPCW | 0.136 | 0.916 | 0.831 | 1.363 | 0.926 | 0.913 | 0.950 | 0.952
| Cox | 0.022 | 0.144 | 0.130 | 0.182 | 0.915 | 0.909 | 0.954 | 0.938
| Log-norm | 0.026 | 0.150 | 0.136 | 0.193 | 0.934 | 0.924 | 0.954 | 0.948
| Log-log | 0.026 | 0.150 | 0.136 | 0.194 | 0.932 | 0.920 | 0.954 | 0.946
Table 2: Numerical results for mean response when sample size $n=250$ and
$n=500$ in fuzzy RD

##  6 Real Data Analysis

We now apply the proposed methodology to men who participated the Prostate,
Lung, Colorectal and Ovarian (PLCO) cancer screening trial and randomized to
receive annual prostate-specific antigen (PSA) screening for 6 years. Among
these 76,682 men received annual PSA screening from 1993 to 2001, those with a
PSA of 4.0 ng/ml at anytime was recommended for further workup and biopsy,
e.g., PSA-based screening strategy, for prostate cancer diagnosis. In the
context of RD design, this practice naturally creates a sharp RD design. We
therefore evaluated the role of additional workup and biopsy, as prompted by a
PSA cutoff of 4.0 ng/ml, on patient survival and cancer incidence. To simplify
our discussions, here we focus on the role of PSA-based screening at the time
of study entry among those who were tested for PSA previously. While the role
of PSA-based screening among those who had been exposed to PSA may be also of
interest, its analysis involves methodology that is still being developed and
will not be discussed here.
Although the local randomization property holds for RD design, this property
assumes that treatment assignment is independent of other covariates. To
alleviate the potential impacts due to the associations between PSA-based
screening strategy and covariates at study entry, we conducted propensity
score matching with age, weight, height and BMI at the baseline. Using
matching, we have 2681 people in the two treatment groups, those with PSA
level greater than 4.0 ng/ml and those with PSA level less than or equal to
4.0 ng/ml. In this sample, censoring rates for mortality and cancer incidence
are 78.6% and 61.3%, respectively.
We select DR approach with conditional expectation method on parametric AFT
model. For parametric AFT model, we choose lognormal and loglogistic
distributions. Table 3 shows results from sharp RD. We use nearest neighbor
(NN) approach for calculation of standard error.

| Mortality | First cancer incidence
---|---|---
| Effect | SE | Effect | SE
Log-norm | 0.054 | 0.243 | -0.203 | 0.260
Log-log | 0.043 | 0.195 | -0.193 | 0.242
Table 3: Results from PLCO dataset

Results show that there is no significant screening effect from baseline PSA
threshold level 4.0mg/nL for mortality and cancer incidence. The screening
effect with cutpoint 4.0mg/nL for mortality is slightly positive, which
implies that screening slightly increase survival time although the effect is
not statistically significant. On the other hand, screening decreases the time
to first cancer incidence. We also create a data-driven RD plot proposed by
Calonico, Cattaneo and Titiunik (2015a) in Figure 1. This data-driven plot
reflects the variability in the data by using local sample means with evenly-
spaced bins. Although it is not intuitive to use our transformed response
directly with theory in Calonico, Cattaneo and Titiunik (2015a), these plots
are still useful to capture the variability of the transformed response with
cutoffs and to check causal effect graphically. These two plots also indicate
no treatment effect using a baseline PSA threshold level of 4.0mg/nL.

Figure 1: Data-driven RD plot by Calonico, Cattaneo and Titiunik (2015a) for
mortality (left) and the first cancer incidence (right)

##  7 Conclusion

In this manuscript, we have proposed new estimators of causal effects with
censored data in RD designs. Simulation studies show that the DR approach
yields a more efficient estimator than IPCW. Moreover, the bias of DR method
is smaller than one of DR method. Moreover, our method uses existing software,
so researchers who easily apply this methodology.
In this manuscript, only one forcing variable is considered for analysis.
However, in practice data may contain several forcing variables, and they may
provide additional information for treatment effect. There are two possible
scenarios : (i) the forcing variable is a function of multiple covariates;
(ii) there is one forcing variable correlated with other covariates as shown
in our data analysis. Imbens and Zajonc (2009) and Zajonc (2012) study the
situation of multiple forcing variables. Recently, Calonico et al. (2018)
propose covariate adjustment approach in RD. It is great of interest to
include covariates or consider the composite forcing variable in RD analysis.
Our future work is to propose estimation procedure of treatment effect in RD
adjusting for effects of other covariates. Although we only have one forcing
variable without covariates, it is expected for the forcing variable to have
multiple cutoffs. Cattaneo et al. (2016) discuss the multi-cutoff problem.
This is also an interesting future work.
We have adapted the LM approach for bandwidth selection. Imbens and
Kalyanaraman (2012) propose optimal bandwidth selection based on mean square
error approximation and Calonico, Cattaneo and Titiunik (2014) propose
bandwidth selection which helps bias correction. These have elegant asymptotic
theory for their bandwidth selection proposals. This is currently under
investigation.

Acknowledgements

The authors thank Dr. Jon Steingrimsson, Dr. Liqun Diao and Dr. Robert
Strawderman for sharing computer code related to the IPCW and DR
transformations.

## Supplementary Materials

##  Appendix A Calculation of $\hat{G}$ and $\hat{Q}_{Y}(\cdot,\cdot)$

In our estimation procedure, it is important to be able to model the censoring
distribution and compute the relevant conditional expectation. We follow the
approach from Steingrimsson et al.(2016). The censoring distribution can be
modeled by a Kaplan-Meier estimator. The conditional expectation can be
modeled in several different ways. Mathematically, for a general nondecreasing
function $H$, it is represented as

$\displaystyle E[H(T)|T\geq
u,W]=\frac{\displaystyle\int_{u}^{\infty}H(T)dF(u|W)}{P(T\geq u|W)},$

where $F$ is distribution of the failure time $T$. Various approaches can be
taken to calculate this expectation. One can use tree-based methods to find
observations in the terminal node by the split, and then compute a Kaplan-
Meier estimator based on the observations in the relevant terminal node. When
using nonparametric and semiparametric methods, this integral is changed to
summation at the time point when failure occurs. Note that in the
nonparametric or semiparametric methods, the integral is computed as long as
uncensored observation(s) greater than $u$ exists. If not, we let the
estimator be maximum of observed time points (Steingrimsson et al. 2016). This
kind of adjustment is required because the survival function estimated through
a nonparametric or semiparametric method may not be proper, i.e., for a given
time $t_{0}$, there may exist an interval with survival function being equal
to $S(t_{0})$. Estimation with parametric distributions is preferred because
the parametric approach will always yield a proper survival function. Standard
AFT models with parametric distributions (e.g. lognormal and log-logistic) can
be used for the estimation (Steingrimsson et al. 2016).

##  Appendix B Proof of asymptotic results

In this section, we provide proofs for the main asymptotic results in the
paper. Since the proofs for IPCW and DR are similar, we mainly prove the
result for DR. Moreover, results for IPCW depending on the censoring
distribution only are similar. The proofs extend the arguments of Hahn, Todd
and Van der Klaauw (1999).
Recall that $O=(\tilde{T},\Delta,W,Z)$. The observed data is given by
$\mathcal{O}=\\{O_{i}\\}_{i=1}^{n}$ where
$O_{i}=(\tilde{T}_{i},\Delta_{i},W_{i},Z_{i})$. Let $G$ and $S$ be survival
probabilities of censoring and failure time based on possibly incorrect
models. Let $G_{0}$ and $S_{0}$ be survival probabilities of censoring and
failure time based on the true models for failure time and censoring,
respectively. We assume that $G_{0}(u)$ and $S_{0}(u|w)$ are continuous and
nonincreasing functions in $u$ for each $w$ with $0\leq
G_{0}(u),S_{0}(u|w)\leq 1$. Furthermore, we assume that $G(u)$ and $S(u|w)$
are right-continuous and nonincreasing in $u$ for each $w$ with $0\leq
G(u),S(u|w)\leq 1$ and $G_{0}(0)=G(0)=S_{0}(0|w)=S(0|w)=1$. Define
$F_{0}(u|w)=1-S_{0}(u|w)$, $\bar{G}_{0}(u)=1-G_{0}(u)$, $F(u|w)=1-S(u|w)$ and
$\bar{G}(u)=1-G(u)$. Now we express $Q_{Y}(\cdot,\cdot)$ by
$Q_{Y}(\cdot,\cdot,S)$ to indicate dependence of conditional expectation on
failure time distribution. We first assume the following regularity
conditions, similar to those in Steingrimsson, Diao, and Strawderman (2019):

  * (C1)

$I_{1}=\displaystyle\int_{0}^{\infty}\log(u)\frac{G_{0}(u)}{G(u-)}dF_{0}(u|w)<\infty$

  * (C2)

For $a>0$,

$\displaystyle
D_{1}(a)=\displaystyle\int_{0}^{a}\frac{S_{0}(u|w)}{S(u|w)}\frac{d\bar{G}_{0}(u)}{G(u-)}<\infty\quad
D_{2}(a)=\displaystyle\int_{0}^{a}\frac{G_{0}(u)S_{0}(u|w)}{G(u)S(u|w)}\frac{d\bar{G}(u)}{G(u-)}<\infty$

  * (C3)

$I_{2}=\displaystyle\int_{0}^{\infty}\log(u)[D_{1}(a-)-D_{2}(a-)]dF(u|w)<\infty$

  * (C4)

$\displaystyle\int_{0}^{\infty}\frac{[\log(u)]^{2}}{G_{0}(u)}dF_{0}(u|w)<\infty$

  * (C5)

$D_{3}(a)=\displaystyle\int_{0}^{a}\frac{Q_{Y}(u,w,S)}{\\{G_{0}(u)\\}^{2}}d\bar{G}_{0}(u)<\infty$
for each $a>0$.

  * (C6)

$\hat{G}$ is uniformly consistent to $G_{0}$.

  * (C7)

$\hat{S}$ is uniformly consistent to $S^{*}$ where $S^{*}$ is possibly
incorrect model of $S$.

As Steingrimsson, Diao, and Strawderman (2019) have shown, from conditions
(C1)-(C5), we can prove that

$\displaystyle
E(Y_{DR}(O;G_{0},S)|W)=E(Y_{DR}(O;G,S_{0})|W)=E(Y_{DR}(O;G_{0},S_{0})|W)=E(Y|W)=\mu(W).$

This is necessary for proving asymptotic normality of
$\hat{\tau}_{FRD}^{IPCW}(\hat{G})$, $\hat{\tau}_{FRD}^{DR}(\hat{G},\hat{S})$,
$\hat{\tau}_{SRD}^{IPCW}(\hat{G})$, $\hat{\tau}_{SRD}^{DR}(\hat{G},\hat{S})$.
Let $p(w)=E(Z|W=w)$ and define

$\displaystyle\mu^{+}(w)=\lim_{w\downarrow
w_{0}}E(Y|W=w)\quad\mu^{-}(w)=\lim_{w\uparrow w_{0}}E(Y|W=w)$ $\displaystyle
p^{+}(w)=\lim_{w\downarrow w_{0}}P(Z=1|W=w)\quad p^{-}(w)=\lim_{w\uparrow
w_{0}}P(Z=1|W=w).$

Now we need conditions similar to Hahn, Todd and Van der Klaauw (1999). Define

$\displaystyle Y_{IPCW}(O;G)=\frac{\Delta Y}{G(T)}\quad
Y^{DR}(O;G,S)=\frac{\Delta
Y}{G(T)}+\int_{0}^{\tilde{T}}\frac{Q_{Y}(u,W,S)}{G(u)}dM_{G}(u)$
$\displaystyle
Y_{IPCW^{*}}(O;G)=Y_{IPCW}(O;G)-\mu^{+}(w_{0})-\mu^{\prime+}(w_{0})(W-w_{0})$
$\displaystyle
Y_{DR^{*}}(O;G,S)=Y^{DR}(O;G,S)-\mu^{+}(w_{0})-\mu^{\prime+}(w_{0})(W-w_{0})$
$\displaystyle Z^{*}=Z-p^{+}(w_{0})-p^{\prime+}(w_{0})(W-w_{0})$
$\displaystyle L_{ih}^{+}=I(W_{i}\geq
w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}\quad
L_{ih}^{-}=I(W_{i}<w_{0})K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}.$

We can define empirical quantities of
$Y_{IPCW}(O;G),Y_{IPCW^{*}}(O;G),Y_{DR}(O;G,S),Y_{DR^{*}}(O;G,S),Z^{*}$.
Denote the empirical quantities by
$Y_{IPCW,i}(O_{i};G),Y_{IPCW^{*},i}(O_{i};G),Y_{DR,i}(O_{i};G,S),Y_{DR^{*},i}(O_{i};G,S),Z_{i}^{*}$.
We further define

$\displaystyle\sigma_{DR}^{2}(w;G,S)=Var(Y_{DR}(O;G,S)|W=w)\quad\sigma_{DR}^{2+}(w_{0};G,S)=\lim_{\epsilon\downarrow
w_{0}}Var(Y_{DR}(O;G,S)|W=w)$
$\displaystyle\sigma_{DR}^{2-}(w_{0};G,S)=\lim_{\epsilon\uparrow
w_{0}}Var(Y_{DR}(O;G,S)|W=w)\quad\eta_{DR}(w;G,S)=Cov(Y_{DR}(O;G,S),Z|W=w)$
$\displaystyle\eta_{DR}^{+}(w_{0};G,S)=\lim_{w\downarrow
w_{0}}Cov(Y_{DR}(O;G,S),Z|W=w)\quad\eta_{DR}^{-}(w_{0};G,S)=\lim_{w\uparrow
w_{0}}Cov(Y_{DR}(O;G,S),Z|W=w).$

We can make similar definitions for
$\sigma_{IPCW}^{2}(G,S),\sigma_{IPCW}^{2+}(w_{0};G,S),\sigma_{IPCW}^{2-}(w_{0};G,S),\eta_{IPCW}(w_{0};G,S),\\\
\eta_{IPCW}^{+}(w_{0};G,S)$ and $\eta_{IPCW}^{-}(w_{0};G,S)$. Define matrices

$\displaystyle\mathbf{X}_{h}=\begin{pmatrix}1&\dfrac{W_{1}-w_{0}}{h}\\\
1&\dfrac{W_{2}-w_{0}}{h}\\\ \vdots&\vdots\\\
1&\dfrac{W_{n}-w_{0}}{h}\end{pmatrix}$
$\displaystyle\mathbf{W}_{h}^{+}=\begin{pmatrix}I(W_{1}\geq
w_{0})K\bigg{(}\dfrac{W_{1}-w_{0}}{h}\bigg{)}&0&0&\ldots&0\\\ 0&I(W_{2}\geq
w_{0})K\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}&0&\ldots&0\\\
\vdots&\vdots&\vdots&\ddots&0\\\ 0&0&\ldots&0&I(W_{n}\geq
w_{0})K\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}\end{pmatrix}.$

  * (R1)

For $W\neq w_{0}$, let $\mu(z)$ and $p(z)$ be twice continuously
differentiable functions. Let $\mu^{\prime}(w)$ and $\mu^{\prime\prime}(w)$ be
the first and second derivatives of $\mu(w)$ and similarly for $p^{\prime}(w)$
and $p^{\prime\prime}(w)$. Let $\mu^{\prime+}(w)$ and $\mu^{\prime\prime+}(w)$
be the first and second derivatives of $\mu^{+}(w)$, and $p^{\prime+}(w)$ and
$p^{\prime\prime+}(w)$ are the first and second derivatives of $p(w)$. Define
$\mu^{\prime-}(w)$ and $\mu^{\prime\prime-}(w)$ to be the first and second
derivative of $\mu^{-}(w)$ and $p^{\prime-}(w)$ and $p^{\prime\prime-}(w)$ are
first and second derivative of $p^{-}(w)$. Assume there exists $B>0$ such that
$|\mu^{+}(w)|,|\mu^{\prime+}(w)|,|\mu^{\prime\prime+}(w)|$ and
$|p^{+}(w)|,|p^{\prime+}(w)|,|p^{\prime\prime+}(w)|$ are uniformly bounded on
$(w_{0},w_{0}+B]$. Similarly,
$|\mu^{-}(w)|,|\mu^{\prime-}(w)|,|\mu^{\prime\prime-}(w)|$ and
$|p^{+}(w)|,|p^{\prime+}(w)|,|p^{\prime\prime+}(w)|$ are uniformly bounded on
$[w_{0}-B,w_{0})$.

  * (R2)

Assume that
$\mu^{+}(w_{0}),\mu^{\prime+}(w_{0}),\mu^{\prime\prime+}(w_{0}),\mu^{-}(w_{0}).\mu^{\prime-}(w_{0}),\mu^{\prime\prime-}(w_{0})$,
$p^{+}(w_{0}),p^{\prime+}(w_{0}),p^{\prime\prime+}(w_{0}),\\\
p^{-}(w_{0}),p^{\prime-}(w_{0})$ and $p^{\prime\prime-}(w_{0})$ are finite.

  * (R3)

Let $g(w)$ be the common density of $W_{i}$. Assume that $g(w)$ is continuous
and bounded away from zero in a neighborhood of $w_{0}$.

  * (R4)

$\sigma_{IPCW}^{2}(w;G_{0}),\sigma_{DR}^{2}(w;G_{0},S^{*})$ and
$\eta_{IPCW}(w;G_{0}),\eta_{DR}(w;G_{0},S^{*})$ are uniformly bounded in a
neighborhood of $w_{0}$.

  * (R5)

Assume that
$\sigma_{IPCW}^{2+}(w_{0};G_{0}),\sigma_{DR}^{2+}(w_{0};G_{0},S^{*}),\sigma_{IPCW}^{2-}(w_{0};G_{0}),\sigma_{DR}^{2-}(w_{0};G_{0},S^{*})$
and $\eta_{IPCW}^{+}(w_{0};G_{0}),\eta_{DR}^{+}(w_{0};G_{0},S^{*})$,
$\eta_{IPCW}^{-}(w_{0};G_{0}),\eta_{DR}^{-}(w_{0};G_{0},S^{*})$ are finite.

  * (R6)

$\lim\limits_{W_{i}\uparrow
w_{0}}E\bigg{[}|Y_{IPCW,i}(O_{i};G_{0})-\mu(W_{i})|^{r}\bigg{|}W_{i}\bigg{]}$
and $\lim\limits_{W_{i}\uparrow
w_{0}}E\bigg{[}|Y_{DR,i}(O_{i};G_{0},S^{*})-\mu(W_{i})|^{r}\bigg{|}W_{i}\bigg{]},r=1,2,3$
are finite. We assume similarly when $W_{i}\downarrow w_{0}$.

  * (R7)

$K$ is continuous and symmetric. Moreover, support of $K$ is compact and for
any $u$, $K(u)\geq 0$.

  * (R8)

The bandwidth satisfies $h\sim n^{-1/5}$ where $\sim$ indicates
“asymptotically equivalent”.

  * (R9)

Let $V_{n}=o_{p}(1)$. Then

$\displaystyle
E\bigg{[}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{j_{1}}(L_{ih}^{+})^{j_{2}}V_{n}\bigg{]}=O(1)\quad
j_{1}=0,\ldots,6\quad j_{2}=1,2,3.$

Conditions (R1) - (R8) are standard ones for proving asymptotic normality of
sharp and fuzzy RD estimators without censoring. (R9) is a required assumption
for making $o_{p}(1)$ term from consistency of $\hat{G}$ and $\hat{S}$.
Modified versions of Lemma 1- Lemma 6 in Hahn, Todd and Van der Klaauw (1999)
are required. Let
$W_{ih}^{r}=\bigg{(}\dfrac{W_{i}-w_{0}}{h_{n}}\bigg{)}^{r},r=0,1$. Consider

$\displaystyle\mathbf{A}_{IPCW^{*},i,h}^{+}(G)=\begin{pmatrix}W_{ih}^{0}Y_{IPCW^{*},i}(O_{i};G)L_{ih}^{+}\\\
W_{ih}^{1}Y_{IPCW^{*},i}(O_{i};G)L_{ih}^{+}\end{pmatrix}\quad\mathbf{A}_{DR^{*},i,h}^{+}(G,S)=\begin{pmatrix}W_{ih}^{0}Y_{DR^{*},i}(O_{i};G,S)L_{ih}^{+}\\\
W_{ih}^{1}Y_{DR^{*},i}(O_{i};G,S)L_{ih}^{+})^{T}\end{pmatrix}$
$\displaystyle\mathbf{A}_{Z^{*},i,h}^{+}=\begin{pmatrix}W_{ih}^{0}Z_{i}^{*}L_{ih}^{+}\\\
W_{ih}^{1}Z_{i}^{*}L_{ih}^{+}\end{pmatrix}.$

Since the proof will be similar for $W<w_{0}$, we will only provide proofs for
the case $W\geq w_{0}$. In the following lemmas, we assume that conditions
(C1) - (C5) and (R1) - (R8) hold. Lemma 1 in Hahn, Todd and Van der Klaauw
(1999) does not involve any modeling of failure time and censoring, so it
clearly holds in our case.

######  Lemma 1.

$\displaystyle\frac{1}{nh}(\mathbf{X}_{h}^{T}\mathbf{W}_{h}^{+}\mathbf{X}_{h})\rightarrow
g(w_{0})\boldsymbol{\Omega}.$

where $\boldsymbol{\Omega}=\begin{pmatrix}\gamma_{0}&\gamma_{1}\\\
\gamma_{1}&\gamma_{2}\end{pmatrix}$ and
$\gamma_{j}=\displaystyle\int_{0}^{\infty}u^{j}K(u)du,j=0,1,2$.

###### Proof.

See Hahn, Todd and Van der Klaauw (1999). ∎

######  Lemma 2.

Two terms

$\displaystyle
E\bigg{\\{}\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{IPCW^{*},i,h}^{+}(O_{i};G_{0})\bigg{\\}}\quad
E\bigg{\\{}\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{DR^{*},i,h}^{+}(O_{i};G_{0},S_{0})\bigg{\\}}$

converge to
$\frac{1}{2}g(z_{0})\mu^{\prime\prime+}(w_{0})h^{2}(\boldsymbol{\delta}+o(1))$
and
$E\bigg{\\{}\dfrac{1}{nh}\displaystyle\sum_{i=1}^{n}\mathbf{A}_{Z^{*},i,h}^{+}\bigg{\\}}$
converges to
$\frac{1}{2}g(z_{0})p^{\prime\prime+}(w_{0})h^{2}(\boldsymbol{\delta}+o(1))$
where

$\displaystyle\boldsymbol{\delta}=\begin{pmatrix}\delta_{0}\\\
\delta_{1}\end{pmatrix}=\begin{pmatrix}\int_{0}^{\infty}u^{2}K(u)du\\\
\int_{0}^{\infty}u^{3}K(u)du\end{pmatrix}.$

###### Proof.

We look at individual terms in
$\mathbf{A}_{IPCW^{*},i,h}^{+}(\mathcal{O};G_{0})$ and
$\mathbf{A}_{DR^{*},i,h}^{+}(\mathcal{O};G_{0},S)$. Let $q=0,1,2$ and for any
model for $S$,

$\displaystyle
U_{IPCW}^{q+}(\mathcal{O};G_{0})=\frac{1}{nh}\sum_{i=1}^{n}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}Y_{IPCW^{*},i}(O_{i};G_{0})L_{ih}^{+}$
$\displaystyle
U_{DR}^{q+}(\mathcal{O};G_{0},S)=\frac{1}{nh}\sum_{i=1}^{n}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}Y_{DR^{*},i}(O_{i};G_{0},S)L_{ih}^{+}.$

Since
$E(Y_{IPCW,i}(G_{0})|W_{i})=E(Y_{DR,i}(G_{0},S)|W_{i})=E(Y_{DR,i}(G,S_{0})|W_{i})=E(Y_{DR,i}(G_{0},S_{0})|W_{i})=\mu(W_{i})$,

$\displaystyle\\{U_{IPCW}^{q+}(\mathcal{O};G_{0})\\}=E\\{U_{DR}^{q}(\mathcal{O};G_{0},S)\\}$
$\displaystyle=\frac{1}{h}E\bigg{[}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}\bigg{(}\frac{1}{2}\mu^{\prime\prime+}(w_{0})\cdot(W_{i}-w_{0})^{2}+\zeta^{Y}(W_{i})\bigg{)}L_{ih}^{+}\bigg{]},$

where

$\displaystyle\zeta^{Y}(w)=\mu(w)-\mu^{+}(w_{0})-\mu^{\prime+}(w_{0})(w-w_{0})-\frac{1}{2}\mu^{\prime\prime+}(w_{0})\cdot(w-w_{0})^{2}.$

By condition (C6) and (C7),

$\displaystyle
E(U_{DR}^{q+}(\mathcal{O};\hat{G},\hat{S}))=\frac{1}{h}E\bigg{[}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}\bigg{(}\frac{1}{2}\mu^{\prime\prime+}(w_{0})\cdot(W_{i}-w_{0})^{2}+\zeta^{Y}(W_{i})+o_{p}(1)+o_{p}(1)\bigg{)}L_{ih}^{+}\bigg{]}$

By Taylor series expansion, $\sup\limits_{w_{0}<w\leq
w_{0}+Bh}|\zeta^{Y}(w)|=o(h^{2})$. Hence by the arguments of Hahn, Todd and
Van der Klaauw (1999),

$\displaystyle\frac{1}{h}E\bigg{[}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}\bigg{(}\frac{1}{2}\mu^{\prime\prime+}(w_{0})(W_{i}-w_{0})^{2}+o_{p}(1)\bigg{)}L_{ih}^{+}\bigg{]}=$
$\displaystyle=\frac{1}{2h}\mu^{\prime\prime+}(w_{0})\int_{w_{0}}^{\infty}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}(W_{i}-w_{0})^{2}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}g(w)dw+o(h^{2})$
$\displaystyle=\delta_{q}g(w_{0})h^{2}+o(h^{2}).$

We can apply the same arguments for
$\mathbf{A}_{DR^{*},i,h}(\mathcal{O};\hat{G},\hat{S})$

Define
$U_{Z}^{q}=\dfrac{1}{nh}\displaystyle\sum_{i=1}^{n}\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}^{q}Z_{i}^{*}L_{ih}^{+}$
and
$\zeta^{Z}(w)=p(w)-p^{+}(w_{0})-p^{\prime+}(w_{0})(w-w_{0})-\frac{1}{2}p^{\prime\prime+}(w_{0})\cdot(w-w_{0})^{2}$.
Clearly, $\sup\limits_{w_{0}<w\leq w_{0}+Bh}|\zeta^{Z}(w)|=o(h^{2})$. Hence

$\displaystyle
E(U_{Z}^{q})=E\bigg{(}\frac{1}{nh}\displaystyle\sum_{i=1}^{n}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}Z_{i}^{*}L_{ih}^{+}\bigg{)}$
$\displaystyle=\frac{1}{h}E\bigg{[}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}\bigg{(}\frac{1}{2}p^{\prime\prime+}(w_{0})(W_{i}-w_{0})^{2}+\zeta^{Z}(W_{i})\bigg{)}L_{ih}^{+}\bigg{]}$
$\displaystyle=\frac{1}{2h}p^{\prime\prime+}(w_{0})\int_{w_{0}}^{\infty}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{k}(W_{i}-w_{0})^{2}K\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}g(w)dw+o(h^{2})=\delta_{q}g(w_{0})h^{2}+o(h^{2}).$

∎

######  Lemma 3.

Let

$\displaystyle\mathbf{A}_{IPCW^{*},h}^{+}(W_{i};G)$
$\displaystyle=E(\mathbf{A}_{IPCW^{*},i,h}^{+}(G)|W_{i})\quad\mathbf{A}_{DR^{*},h}^{+}(W_{i};G,S)=E(\mathbf{A}_{DR^{*},i,h}^{+}(G,S)|W_{i})$
$\displaystyle\mathbf{A}_{Z^{*},h}^{+}(W_{i})$
$\displaystyle=E(\mathbf{A}_{Z^{*},i,h}^{+}|W_{i})\quad\mathbf{A}_{IPCW^{*},h}^{-}(W_{i};G)=E(\mathbf{A}_{IPCW^{*},i,h}^{-}(G)|W_{i})$
$\displaystyle\mathbf{A}_{DR^{*},h}^{-}(W_{i};G,S)$
$\displaystyle=E(\mathbf{A}_{DR^{*},i,h}^{-}(G,S)|W_{i})\quad\mathbf{A}_{Z^{*},h}^{-}(W_{i})=E(\mathbf{A}_{Z^{*},i,h}^{-}|W_{i})$

Then

$\displaystyle\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{IPCW^{*},h}^{+}(W_{i};\hat{G})=E\bigg{\\{}\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{IPCW^{*},i,h}^{+}(G_{0})\bigg{\\}}+o_{p}(h^{2})$
$\displaystyle\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{DR^{*},h}^{+}(W_{i};\hat{G},\hat{S})=E\bigg{\\{}\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{DR^{*},i,h}^{+}(G_{0},S^{*})\bigg{\\}}+o_{p}(h^{2})$
$\displaystyle\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{Z^{*},h}^{+}(W_{i})=E\bigg{(}\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{Z^{*},i,h}^{+}\bigg{)}+o_{p}(h^{2}).$

###### Proof.

$\displaystyle\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{DR^{*},h}^{+}(W_{i};\hat{G},\hat{S})=\frac{1}{nh}\sum_{i=1}^{n}\boldsymbol{X}_{i,h}L_{i,h}^{+}\bigg{(}\frac{1}{2}\mu^{\prime\prime+}(w_{0})(W_{i}-w_{0})^{2}+\zeta^{Y}(W_{i})+E_{n}^{*}\bigg{)}$

where $E_{n}^{*}=o_{p}(1)$ and $\mathbf{X}_{i,h}=\begin{pmatrix}1\\\
\frac{W_{i}-w_{0}}{h}\end{pmatrix}$. Then using similar calculations to those
in Hahn, Todd and Van der Klaauw (1999),

$\displaystyle
Var\bigg{\\{}\frac{1}{nh}\sum_{i=1}^{n}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}L_{ih}^{+}\bigg{(}\frac{1}{2}\mu^{\prime\prime+}(w_{0})(W_{i}-w_{0})^{2}+\zeta^{Y}(W_{i})+E_{n}^{*}\bigg{)}\bigg{\\}}$
$\displaystyle=\frac{1}{nh}\bigg{[}E\bigg{\\{}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}L_{ih}^{+}\bigg{(}\frac{1}{2}\mu^{\prime\prime+}(w_{0})(W_{i}-w_{0})^{2}+\zeta^{Y}(W_{i})+E_{n}^{*}\bigg{)}\bigg{\\}}^{2}$
$\displaystyle-\bigg{[}E\bigg{\\{}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}L_{ih}^{+}\bigg{(}\frac{1}{2}\mu^{\prime\prime+}(w_{0})(W_{i}-w_{0})^{2}+\zeta^{Y}(W_{i})+E_{n}^{*}\bigg{)}\bigg{\\}}\bigg{]}^{2}\bigg{]}.$

Then we can choose constant $\varrho$ such that

$\displaystyle\frac{1}{nh}\bigg{[}E\bigg{\\{}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}L_{ih}^{+}\bigg{(}\frac{1}{2}\mu^{\prime\prime+}(w_{0})(W_{i}-w_{0})^{2}+\zeta^{Y}(W_{i})+E_{n}^{*}\bigg{)}\bigg{\\}}^{2}$
$\displaystyle-\bigg{[}E\bigg{\\{}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{q}L_{ih}^{+}\bigg{(}\frac{1}{2}\mu^{\prime\prime+}(w_{0})(W_{i}-w_{0})^{2}+\zeta^{Y}(W_{i})+E_{n}^{*}\bigg{)}\bigg{\\}}\bigg{]}^{2}\bigg{]}$
$\displaystyle\leq\frac{1}{nh}\varrho
E\bigg{\\{}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{2q}(L_{ih}^{+})^{2}\bigg{(}(W_{i}-w_{0})^{4}+\\{\zeta^{Y}(W_{i})\\}^{2}\bigg{)}\bigg{\\}}^{2}.$

From Hahn, Todd and Van der Klaauw (1999), we have that

$\displaystyle\frac{1}{nh}\varrho
E\bigg{\\{}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{2q}(L_{ih}^{+})^{2}\bigg{(}(W_{i}-w_{0})^{4}+\\{\zeta^{Y}(W_{i})\\}^{2}\bigg{)}\bigg{\\}}^{2}=o(h^{2}).$

Hence

$\displaystyle\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{DR^{*},h}^{+}(W_{i};\hat{G},\hat{S})=E\bigg{\\{}\frac{1}{nh}\sum_{i=1}^{n}\mathbf{A}_{DR^{*},i,h}^{+}(G_{0},S^{*})\bigg{\\}}+o_{p}(h^{2}).$

Similar arguments can be applied to
$\dfrac{1}{nh}\displaystyle\sum_{i=1}^{n}\mathbf{A}_{IPCW^{*},h}^{+}(W_{i};\hat{G})$
and $\dfrac{1}{nh}\displaystyle\sum_{i=1}^{n}\mathbf{A}_{Z^{*},h}^{+}(W_{i})$.
∎

######  Lemma 4.

Let $v_{q}=\int_{0}^{\infty}u^{q}\\{K(u)\\}^{2},q=0,1,2$ and

$\displaystyle\mathbf{\bar{A}}_{IPCW^{*},h}^{+}(\mathcal{O};\hat{G})=\frac{1}{nh}\sum_{i=1}^{n}\\{\mathbf{A}_{IPCW^{*},i,h}^{+}(O_{i};\hat{G})-\mathbf{A}_{IPCW^{*},h}^{+}(W_{i};\hat{G})\\}$
$\displaystyle\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};\hat{G},\hat{S})=\frac{1}{nh}\sum_{i=1}^{n}\\{\mathbf{A}_{DR^{*},i,h}^{+}(O_{i};\hat{G},\hat{S})-\mathbf{A}_{DR^{*},h}^{+}(W_{i};\hat{G},\hat{S})\\}$
$\displaystyle\mathbf{\bar{A}}_{Z^{*},h}^{+}=\frac{1}{nh}\sum_{i=1}^{n}\\{\mathbf{A}_{Z^{*},i,h}^{+}-\mathbf{A}_{Z^{*},h}^{+}(W_{i})\\}.$

Then

$\displaystyle
Var\\{\mathbf{\bar{A}}_{IPCW^{*},h}^{+}(\mathcal{O};G_{0})\\}=\frac{1}{nh}\sigma_{IPCW}^{2+}(w_{0};G_{0})g(w_{0})\mathcal{V}\quad
Var\\{\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};G_{0},S)\\}=\frac{1}{nh}\sigma_{DR}^{2+}(w_{0};G_{0},S)g(w_{0})\mathcal{V}$
$\displaystyle
Var(\mathbf{\bar{A}}_{Z^{*},h}^{+})=\frac{1}{nh}p^{+}(w_{0})\\{1-p^{+}(w_{0})\\}g(w_{0})\mathcal{V}$
$\displaystyle
Cov(\mathbf{\bar{A}}_{IPCW^{*},h}^{+}(\mathcal{O};G_{0}),\mathbf{\bar{A}}_{Z^{*},h}^{+})=\frac{1}{nh}\eta_{IPCW}^{+}(w_{0};G_{0})g(w_{0})\mathcal{V}$
$\displaystyle
Cov(\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};G_{0},S),\mathbf{\bar{A}}_{Z^{*},h}^{+})=\frac{1}{nh}\eta_{DR}^{+}(w_{0};G_{0},S)g(w_{0})\mathcal{V},$

where

$\mathcal{V}=\begin{pmatrix}v_{0}+o(1)&v_{1}+o(1)\\\
v_{1}+o(1)&v_{2}+o(1)\end{pmatrix}.$

###### Proof.

Let $q=0,1,2$ and

$\displaystyle\mathbf{A}_{i,IPCW}^{ctd,+}(O_{i};G)=\begin{pmatrix}W_{ih}^{0}L_{ih}^{+}(Y_{IPCW,i}(O_{i};G)-\mu(W_{i}))\\\
W_{ih}^{1}L_{ih}^{+}(Y_{IPCW,i}(O_{i};G)-\mu(W_{i}))\end{pmatrix}$
$\displaystyle\mathbf{A}_{i,DR}^{ctd,+}(O_{i};G,S)=\begin{pmatrix}W_{ih}^{0}L_{ih}^{+}(Y_{DR,i}(O_{i};G,S)-\mu(W_{i}))\\\
W_{ih}^{1}L_{ih}^{+}(Y_{DR,i}(O_{i};G,S)-\mu(W_{i}))\end{pmatrix}.$

By expansions,

$\displaystyle\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};\hat{G},\hat{S})=\frac{1}{nh}\sum_{i=1}^{n}\\{\mathbf{A}_{DR^{*},i,h}^{+}(O_{i};\hat{G},\hat{S})-\mathbf{A}_{DR^{*},h}^{+}(W_{i};\hat{G},\hat{S})\\}$
$\displaystyle=\frac{1}{nh}\sum_{i=1}^{n}\\{\mathbf{A}_{DR^{*},i,h}^{+}(O_{i};\hat{G},\hat{S})-\mathbf{A}_{DR^{*},h}^{+}(W_{i};G_{0},S^{*})+\mathbf{A}_{DR^{*},h}^{+}(W_{i};G_{0},S^{*})-\mathbf{A}_{DR^{*},h}^{+}(W_{i};\hat{G},\hat{S})\\}$
$\displaystyle=\frac{1}{nh}\sum_{i=1}^{n}\mathbf{X}_{i,h}L_{ih}^{+}\\{Y_{DR^{*},i}(O_{i};\hat{G},\hat{S})-E(Y_{DR^{*},i}(O_{i};G_{0},S^{*})|W_{i})$
$\displaystyle+E(Y_{DR^{*},i}(O_{i};G_{0},S^{*})|W_{i})-E(Y_{DR^{*},i}(O_{i};\hat{G},\hat{S})|W_{i})\\}$
$\displaystyle=\frac{1}{nh}\sum_{i=1}^{n}\mathbf{X}_{i,h}L_{ih}^{+}\\{Y_{DR^{*},i}(O_{i};\hat{G},\hat{S})-E(Y_{DR^{*},i}(O_{i};G_{0},S^{*})|W_{i})$
$\displaystyle+E(Y_{DR^{*},i}(O_{i};G_{0},S^{*})|W_{i})-E(Y_{DR^{*},i}(O_{i};\hat{G},\hat{S})|W_{i})\\}$
(1)

By consistency of $\hat{G}$ and $\hat{S}$,

$\displaystyle\eqref{eq:eqlem4}=\frac{1}{nh}\sum_{i=1}^{n}\\{\mathbf{X}_{i,h}L_{ih}^{+}\\{Y_{DR^{*},i}(O_{i};\hat{G},\hat{S})-E(Y_{DR^{*},i}(O_{i};G_{0},S^{*})|W_{i})+E_{n}^{*}\\}$
$\displaystyle=\frac{1}{nh}\sum_{i=1}^{n}\\{\mathbf{A}_{i,DR}^{ctd,+}(O_{i};G_{0},S^{*})+\mathbf{X}_{i,h}L_{ih}^{+}E_{n}^{*}\\}$

By looking at variances of the individual terms and and assumption (R9),

$\displaystyle
Var\bigg{\\{}\frac{1}{nh}\sum_{i=1}^{n}\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};\hat{G},\hat{S})\bigg{\\}}=Var\bigg{\\{}\frac{1}{nh}\sum_{i=1}^{n}\\{\mathbf{A}_{i,DR}^{ctd,+}(O_{i};G_{0},S^{*})+\mathbf{X}_{i,h}L_{ih}^{+}E_{n}^{*}\\}\bigg{\\}}$
$\displaystyle=\frac{1}{(nh)^{2}}\bigg{[}\int_{w_{0}}^{\infty}\bigg{(}\frac{w-w_{0}}{h}\bigg{)}^{2q}\bigg{\\{}I(W_{i}\geq
w_{0})K\bigg{(}\frac{w-w_{0}}{h}\bigg{)}\bigg{\\}}^{2}\sigma_{DR}^{2+}(w;G_{0},S^{*})g(w)dw+O(1)+O(1)\bigg{]}$
$\displaystyle=\frac{1}{nh}\\{\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})v_{q}+o(1)+o(1)\\}=\frac{1}{nh}\\{\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})v_{q}+o(1)\\}.$

Note that the covariance term between
$\mathbf{A}_{i,DR}^{ctd,+}(O_{i};G_{0},S^{*})$ and
$\mathbf{X}_{i,h}L_{ih}^{+}E_{n}^{*}$ is finite, hence it is $O(1)$.
Similarly,

$\displaystyle
Var\bigg{\\{}\frac{1}{nh}\sum_{i=1}^{n}\mathbf{\bar{A}}_{IPCW^{*},h}^{+}(\mathcal{O};\hat{G})\bigg{\\}}=\frac{1}{nh}\\{\sigma_{IPCW}^{2+}(w_{0};G_{0})v_{q}+o(1)\\}.$

Derivation of the variance of $\mathbf{\bar{A}}_{Z^{*},h}^{+}$ is similar.
Next, for cross-covariance quantities, we only prove for individual terms of
$Cov(\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};\hat{G},\hat{S}),\mathbf{\bar{A}}_{Z^{*},h}^{+})$.
Note that

$\displaystyle\frac{1}{nh^{2}}\int_{w_{0}}^{\infty}\bigg{(}\frac{w-w_{0}}{h}\bigg{)}^{2q}\bigg{\\{}I(W_{i}\geq
w_{0})K\bigg{(}\frac{w-w_{0}}{h}\bigg{)}\bigg{\\}}^{2}\eta_{DR}(w)g(w)dw$
$\displaystyle=\frac{1}{nh}\\{\eta_{DR}^{+}(w_{0};G_{0},S^{*})v_{q}+o(1)\\}.$

Hence the result follows. ∎

######  Lemma 5.

$\displaystyle\sqrt{nh}\begin{pmatrix}\mathbf{\bar{A}}_{IPCW^{*},h}^{+}(\mathcal{O};G_{0})\\\
\mathbf{\bar{A}}_{Z^{*},h}^{+}\end{pmatrix}\overset{d}{\longrightarrow}\\{g(w_{0})\\}\mathcal{N}\bigg{(}0,\begin{pmatrix}\sigma_{IPCW}^{2+}(w_{0};G_{0})\mathcal{V}&\eta_{IPCW}^{+}(w_{0};G_{0})\mathcal{V}\\\
\eta_{IPCW}^{+}(w_{0};G_{0})\mathcal{V}&p^{+}(w_{0})(1-p^{+}(w_{0}))\mathcal{V}\end{pmatrix}\bigg{)}$
$\displaystyle\sqrt{nh}\begin{pmatrix}\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};G_{0},S)\\\
\mathbf{\bar{A}}_{Z^{*},h}^{+}\end{pmatrix}\overset{d}{\longrightarrow}\\{g(w_{0})\\}\mathcal{N}\bigg{(}0,\begin{pmatrix}\sigma_{DR}^{2+}(w_{0};G_{0},S)\mathcal{V}&\eta_{DR}^{+}(w_{0};G_{0},S)\mathcal{V}\\\
\eta_{DR}^{+}(w_{0};G_{0},S)\mathcal{V}&p^{+}(w_{0})(1-p^{+}(w_{0}))\mathcal{V}\end{pmatrix}\bigg{)}.$

###### Proof.

For any $\mathbf{b}\in\mathbb{R}^{4}$, we will show that

$\displaystyle\mathbf{a}^{T}\begin{pmatrix}\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};\hat{G},\hat{S})\\\
\mathbf{\bar{A}}_{Z^{*},h}^{+}\end{pmatrix}\overset{d}{\longrightarrow}\mathbf{a}^{T}\mathcal{X},$

where $\mathcal{X}$ has multivariate normal distribution with mean 0 and
covariance matrix

$\displaystyle\begin{pmatrix}\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})\mathcal{V}&\eta_{DR}^{+}(w_{0};G_{0},S^{*})\mathcal{V}\\\
\eta_{IPCW}^{+}(w_{0};G_{0})\mathcal{V}&p^{+}(w_{0})(1-p^{+}(w_{0}))\mathcal{V}\end{pmatrix}.$

As in Hahn, Todd and Van der Klaauw (1999), by using the Lyapounov condition
with only considering individual terms, we would like to show that

$\displaystyle\frac{1}{\sqrt{nh}}\frac{1}{h}E\bigg{[}(L_{ih}^{+})^{3}\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}^{3q}|Y_{DR,i}(O_{i};\hat{G},\hat{S})-\mu(W_{i})|^{3}\bigg{]}=o(1)$

By expansion,

$\displaystyle\frac{1}{h}E\bigg{[}(L_{ih}^{+})^{3}\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}^{3q}|Y_{DR,i}(O_{i};\hat{G},\hat{S})-\mu(W_{i})|^{3}\bigg{]}$
$\displaystyle=\frac{1}{h}E\bigg{[}(L_{ih}^{+})^{3}\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}^{3q}|Y_{DR,i}(O_{i};\hat{G},\hat{S})-Y_{DR,i}(O_{i};G_{0},S^{*})+Y_{DR,i}(O_{i};G_{0},S^{*})-\mu(W_{i})|^{3}\bigg{]}$
$\displaystyle=\frac{1}{h}E\bigg{[}(L_{ih}^{+})^{3}\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}^{3q}\\{|Y_{DR,i}(O_{i};\hat{G},\hat{S})-Y_{DR,i}(O_{i};G_{0},S^{*})|^{3}$
$\displaystyle+3\\{Y_{DR,i}(O_{i};\hat{G},\hat{S})-Y_{DR,i}(O_{i};G_{0},S^{*})\\}^{2}|Y_{DR,i}(O_{i};G_{0},S^{*})-\mu(W_{i})|$
$\displaystyle+3|Y_{DR,i}(O_{i};\hat{G},\hat{S})-Y_{DR,i}(O_{i};G_{0},S^{*})|\\{Y_{DR,i}(O_{i};G_{0},S^{*})-\mu(W_{i})\\}^{2}$
$\displaystyle+|Y_{DR,i}(O_{i};G_{0},S^{*})-\mu(W_{i})|^{3}\\}\bigg{]}$ (2)

By consistency of $\hat{G}$ and $\hat{S}$, and by assumptions (R6) and (R9),

$\displaystyle\eqref{eq:eqlem5}=\int_{0}^{\infty}u^{3q}\\{K(u)\\}^{3}\iota(w_{0}+uh)g(w_{0}+uh)du+O(1).$

where $\iota(w)=E\\{|Y_{DR,i}(O_{i};G_{0},S^{*})-\mu(W_{i})|^{3}|W_{i}=w\\}$.
Hence

$\displaystyle\frac{1}{\sqrt{nh}}\frac{1}{h}E\bigg{[}(L_{ih}^{+})^{3}\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}^{3q}|Y_{IPCW,i}(O_{i};G_{0})-\mu(W_{i})|^{3}\bigg{]}=o(1).$

Similarly,
$\dfrac{1}{\sqrt{nh}}\dfrac{1}{h}E\bigg{[}(L_{ih}^{+})^{3}\bigg{(}\dfrac{W_{i}-w_{0}}{h}\bigg{)}^{3q}|Z_{i}-p(W_{i})|^{3}\bigg{]}=o(1)$.
Hence the conditional asymptotic normality of
$\begin{pmatrix}\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};G_{0},S^{*})\\\
\mathbf{\bar{A}}_{Z^{*},h}^{+}\end{pmatrix}$ follows by the Cramér-Wold
theorem. Similarly,

$\displaystyle\frac{1}{\sqrt{nh}}\frac{1}{h}E\bigg{[}(L_{ih}^{+})^{3}\bigg{(}\frac{W_{i}-w_{0}}{h}\bigg{)}^{3q}|Y_{IPCW,i}(\mathcal{O};\hat{G})-\mu(W_{i})|^{3}\bigg{]}=o(1).$

Hence the result for
$\begin{pmatrix}\mathbf{\bar{A}}_{IPCW^{*},h}^{+}(\mathcal{O};\hat{G})\\\
\mathbf{\bar{A}}_{Z^{*},h}^{+}\end{pmatrix}$ follows. ∎

######  Lemma 6.

$\displaystyle\frac{1}{\sqrt{nh}}\sum_{i=1}^{n}\begin{pmatrix}\mathbf{A}_{IPCW^{*},i,h}^{+}(O_{i};\hat{G})\\\
\mathbf{A}_{Z^{*},i,h}^{+}\end{pmatrix}-\frac{1}{2}n^{1/2}h^{5/2}g(w_{0})\begin{pmatrix}\mu^{\prime\prime+}(w_{0})\boldsymbol{\delta}\\\
p^{\prime\prime+}(w_{0})\boldsymbol{\delta}\end{pmatrix}$
$\displaystyle\overset{d}{\longrightarrow}\\{g(w_{0})\\}^{1/2}\mathcal{N}\bigg{(}0,\begin{pmatrix}\sigma_{IPCW}^{2+}(w_{0};G_{0})\mathcal{V}&\eta_{IPCW}^{+}(w_{0};G_{0})\mathcal{V}\\\
\eta_{IPCW}^{+}(w_{0};G_{0})\mathcal{V}&p^{+}(w_{0})(1-p^{+}(w_{0}))\mathcal{V}\end{pmatrix}\bigg{)}$
$\displaystyle\frac{1}{\sqrt{nh}}\sum_{i=1}^{n}\begin{pmatrix}\mathbf{A}_{DR^{*},i,h}^{+}(O_{i};\hat{G},\hat{S})\\\
\mathbf{A}_{Z^{*},i,h}^{+}\end{pmatrix}-\frac{1}{2}n^{1/2}h^{5/2}g(w_{0})\begin{pmatrix}\mu^{\prime\prime+}(w_{0})\boldsymbol{\delta}\\\
p^{\prime\prime+}(w_{0})\boldsymbol{\delta}\end{pmatrix}$
$\displaystyle\overset{d}{\longrightarrow}\\{g(w_{0})\\}^{1/2}\mathcal{N}\bigg{(}0,\begin{pmatrix}\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})\mathcal{V}&\eta_{DR}^{+}(w_{0};G_{0},S^{*})\mathcal{V}\\\
\eta_{DR}^{+}(w_{0};G_{0},S)\mathcal{V}&p^{+}(w_{0})(1-p^{+}(w_{0}))\mathcal{V}\end{pmatrix}\bigg{)}.$

###### Proof.

Note that

$\displaystyle\frac{1}{\sqrt{nh}}\sum_{i=1}^{n}\begin{pmatrix}\mathbf{A}_{DR^{*},i,h}^{+}(O_{i};\hat{G},\hat{S})\\\
\mathbf{A}_{Z^{*},i,h}^{+}\end{pmatrix}=\sqrt{nh}\begin{pmatrix}\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};\hat{G},\hat{S})\\\
\mathbf{\bar{A}}_{Z^{*},h}^{+}\end{pmatrix}+\frac{1}{\sqrt{nh}}\sum_{i=1}^{n}\begin{pmatrix}\mathbf{A}_{DR^{*},h}^{+}(W_{i};\hat{G},\hat{S})\\\
\mathbf{A}_{Z^{*},h}^{+}(W_{i})\end{pmatrix}.$

By Lemma 1 and 2,

$\displaystyle\frac{1}{\sqrt{nh}}\sum_{i=1}^{n}\begin{pmatrix}\mathbf{A}_{DR^{*},h}^{+}(W_{i};\hat{G},\hat{S})\\\
\mathbf{A}_{Z^{*},h}^{+}(W_{i})\\\
\end{pmatrix}=\begin{pmatrix}\frac{1}{2}g(w_{0})(nh)^{1/2}\mu^{\prime\prime+}(w_{0})h^{2}(\boldsymbol{\delta}+o(1))+(nh)^{1/2}o_{p}(h^{2})\\\
\frac{1}{2}g(w_{0})(nh)^{1/2}p^{\prime\prime+}(w_{0})h^{2}(\boldsymbol{\delta}+o(1))+(nh)^{1/2}o_{p}(h^{2})\end{pmatrix}.$

Since $h\sim n^{-1/5}$, we obtain

$\displaystyle\frac{1}{\sqrt{nh}}\sum_{i=1}^{n}\begin{pmatrix}\mathbf{A}_{DR^{*},h}^{+}(W_{i};\hat{G},\hat{S})\\\
\mathbf{A}_{Z^{*},h}^{+}(W_{i})\end{pmatrix}\overset{p}{\longrightarrow}\begin{pmatrix}\frac{1}{2}g(w_{0})(nh)^{1/2}\mu^{\prime\prime+}(w_{0})h^{2}\boldsymbol{\delta}\\\
\frac{1}{2}g(w_{0})(nh)^{1/2}p^{\prime\prime+}(w_{0})h^{2}\boldsymbol{\delta}\end{pmatrix},$

and by Lemma 4,

$\displaystyle\sqrt{nh}\begin{pmatrix}\mathbf{\bar{A}}_{DR^{*},h}^{+}(\mathcal{O};\hat{G},\hat{S})\\\
\mathbf{\bar{A}}_{Z^{*},h}^{+}\end{pmatrix}\overset{d}{\longrightarrow}\\{g(w_{0})\\}\mathcal{N}\bigg{(}0,\begin{pmatrix}\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})\mathcal{V}&\eta_{DR}^{+}(w_{0};G_{0},S^{*})\mathcal{V}\\\
\eta_{DR}^{+}(w_{0};G_{0},S^{*})\mathcal{V}&p^{+}(w_{0})(1-p^{+}(w_{0}))\mathcal{V}\end{pmatrix}\bigg{)}.$

Then by Slutsky’s theorem, the result holds for the doubly robust estimator.
Similarly, we can prove for
$\dfrac{1}{nh}\displaystyle\sum_{i=1}^{n}\mathbf{A}_{IPCW^{*},i,h}^{+}(O_{i};\hat{G})$.
∎

For all lemmas, we can derive similar results for $W<w_{0}$. By combining the
results from the cases $W\geq w_{0}$ and $W<w_{0}$, we are ready to prove the
main theorems. Recall that

$\displaystyle\tau^{Y}=\lim_{w\downarrow w_{0}}E(Y|W=w)-\lim_{w\uparrow
w_{0}}E(Y|W=w)$ $\displaystyle\tau^{Z}=\lim_{w\downarrow
w_{0}}P(Z=1|W=w)-\lim_{w\uparrow w_{0}}P(Z=1|W=w)$
$\displaystyle\hat{\tau}_{FRD}^{IPCW}(G)=\frac{\hat{\alpha}_{R,IPCW}^{FRD,Y}(G)-\hat{\alpha}_{L,IPCW}^{FRD,Y}(G)}{\hat{\alpha}_{R}^{FRD,Z}-\hat{\alpha}_{L}^{FRD,Z}}\quad\hat{\tau}_{FRD}^{DR}(G,S)=\frac{\hat{\alpha}_{R,DR}^{FRD,Y}(G,S)-\hat{\alpha}_{L,DR}^{FRD,Y}(G,S)}{\hat{\alpha}_{R}^{FRD,Z}-\hat{\alpha}_{L}^{FRD,Z}}$
$\displaystyle\hat{\tau}_{SRD}^{IPCW}(G)=\hat{\alpha}_{R,IPCW}^{SRD,Y}(G)-\hat{\alpha}_{L,IPCW}^{SRD,Y}(G)\quad\hat{\tau}_{SRD}^{DR}(G,S)=\hat{\alpha}_{R,DR}^{SRD,Y}(G,S)-\hat{\alpha}_{L,DR}^{SRD,Y}(G,S).$

Let

$\displaystyle\rho^{+}=\frac{\displaystyle\bigg{(}\int_{0}^{\infty}u^{2}K(u)du\bigg{)}^{2}-\bigg{(}\int_{0}^{\infty}u^{3}K(u)du\bigg{)}\bigg{(}\int_{0}^{\infty}uK(u)du\bigg{)}}{2\displaystyle\bigg{(}\int_{0}^{\infty}u^{2}K(u)du\bigg{(}\int_{0}^{\infty}K(u)du\bigg{)}-\bigg{(}\int_{0}^{\infty}uK(u)du\bigg{)}^{2}\bigg{)}}$
$\displaystyle\rho^{-}=\frac{\displaystyle\bigg{(}\int_{-\infty}^{0}u^{2}K(u)du\bigg{)}^{2}-\bigg{(}\int_{-\infty}^{0}u^{3}K(u)du\bigg{)}\bigg{(}\int_{-\infty}^{0}uK(u)du\bigg{)}}{2\displaystyle\bigg{(}\int_{-\infty}^{0}u^{2}K(u)du\bigg{(}\int_{-\infty}^{0}K(u)du\bigg{)}-\bigg{(}\int_{-\infty}^{0}uK(u)du\bigg{)}^{2}\bigg{)}}$
$\displaystyle\upsilon^{+}=\frac{\displaystyle\int_{0}^{\infty}\bigg{\\{}\bigg{(}\int_{0}^{\infty}s^{2}K(s)ds\bigg{)}-\bigg{(}\int_{0}^{\infty}sK(s)ds\bigg{)}u\bigg{\\}}^{2}\\{K(u)\\}^{2}du}{\displaystyle
g(w_{0})\bigg{\\{}\bigg{(}\int_{0}^{\infty}u^{2}K(u)du\bigg{)}\bigg{(}\int_{0}^{\infty}K(u)du\bigg{)}-\bigg{(}\int_{0}^{\infty}uK(u)du\bigg{)}^{2}\bigg{\\}}^{2}}$
$\displaystyle\upsilon^{-}=\frac{\displaystyle\int_{-\infty}^{0}\bigg{\\{}\bigg{(}\int_{-\infty}^{0}s^{2}K(s)ds\bigg{)}-\bigg{(}\int_{-\infty}^{0}sK(s)ds\bigg{)}u\bigg{\\}}^{2}\\{K(u)\\}^{2}du}{\displaystyle
g(w_{0})\bigg{\\{}\bigg{(}\int_{-\infty}^{0}u^{2}K(u)du\bigg{)}\bigg{(}\int_{-\infty}^{0}K(u)du\bigg{)}-\bigg{(}\int_{-\infty}^{0}uK(u)du\bigg{)}^{2}\bigg{\\}}^{2}}.$

Moreover, let

$\displaystyle\varphi_{FRD}=\frac{1}{\tau^{Z}}(\rho^{+}\mu^{\prime\prime+}(w_{0})-\rho^{-}\mu^{\prime\prime-}(w_{0}))-\frac{\tau^{Y}}{(\tau^{Z})^{2}}(\rho^{+}p^{\prime\prime+}(w_{0})-\rho^{-}p^{\prime\prime-}(w_{0}))$
$\displaystyle\Sigma_{FRD}^{IPCW}(G)=\frac{1}{\tau^{Z}}(\upsilon^{+}\sigma_{IPCW}^{2+}(w_{0};G)+\upsilon^{-}\sigma_{IPCW}^{2-}(w_{0};G))-2\frac{\tau^{Y}}{(\tau^{Z})^{3}}(\upsilon^{+}\eta_{IPCW}^{+}(w_{0};G)+$
$\displaystyle\upsilon^{-}\eta_{IPCW}^{-}(w_{0};G))+\frac{(\tau^{Y})^{2}}{(\tau^{Z})^{4}}\\{\upsilon^{+}p^{+}(w_{0})(1-p^{+}(w_{0}))+\upsilon^{-}p^{-}(w_{0})(1-p^{-}(w_{0}))\\}$
$\displaystyle\Sigma_{FRD}^{DR}(G,S)=\frac{1}{\tau^{Z}}(\upsilon^{+}\sigma_{DR}^{2+}(w_{0};G,S)+\upsilon^{-}\sigma_{DR}^{2-}(w_{0};G,S))-2\frac{\tau^{Y}}{(\tau^{Z})^{3}}(\upsilon^{+}\eta_{DR}^{+}(w_{0};G,S)+$
$\displaystyle\upsilon^{-}\eta_{DR}^{-}(w_{0};G,S))+\frac{(\tau^{Y})^{2}}{(\tau^{Z})^{4}}\\{\upsilon^{+}p^{+}(w_{0})(1-p^{+}(w_{0}))+\upsilon^{-}p^{-}(w_{0})(1-p^{-}(w_{0}))\\}.$

and

$\displaystyle\varphi_{SRD}=\rho^{+}\mu^{\prime\prime+}(w_{0})-\rho^{-}\mu^{\prime\prime-}(w_{0})$
$\displaystyle\Sigma_{SRD}^{IPCW}(G)=\upsilon^{+}\sigma_{IPCW}^{2+}(w_{0};G)+\upsilon^{-}\sigma_{IPCW}^{2-}(w_{0};G)$
$\displaystyle\Sigma_{SRD}^{DR}(G,S)=\upsilon^{+}\sigma_{DR}^{2+}(w_{0};G,S)+\upsilon^{-}\sigma_{DR}^{2-}(w_{0};G,S).$

###### Proof of Theorem 1.

Now we use a matrix to express terms. As before, we only prove for doubly
robust estimator. Let
$\mathbf{Y}_{DR}(\mathcal{O};\hat{G},\hat{S})=\\{Y_{DR,i}(O_{i};\hat{G},\hat{S})\\}_{i=1}^{n}$
and $\mathbf{Z}=\\{Z_{i}\\}_{i=1}^{n}$. Note that

$\displaystyle\sum_{i=1}^{n}\begin{pmatrix}\mathbf{A}_{DR^{*},i,h}^{+}(O_{i};\hat{G},\hat{S})\\\
\mathbf{A}_{Z^{*},i,h}^{+}\end{pmatrix}=\begin{pmatrix}\mathbf{X}_{h}^{T}\mathbf{W}_{h}^{+}\mathbf{Y}_{DR}(\mathcal{O};\hat{G},\hat{S})\\\
\mathbf{X}_{h}^{T}\mathbf{W}_{h}^{+}\mathbf{Z}\end{pmatrix}$

Similar calculation as Hahn, Todd and Van der Klaauw (1999),

$\displaystyle
n^{2/5}\begin{pmatrix}\hat{\alpha}_{R,DR}^{FRD,Y}(\hat{G},\hat{S})-\mu^{+}(w_{0})\\\
\hat{\beta}_{R,DR}^{FRD,Y}(\hat{G},\hat{S})-\mu^{\prime+}(w_{0})\\\
\hat{\alpha}_{R}^{FRD,Z}-p^{+}(w_{0})\\\
\hat{\beta}_{R}^{FRD,Z}-p^{\prime+}(w_{0})\end{pmatrix}=\boldsymbol{\Psi}\begin{pmatrix}(\mathbf{X}_{h}^{T}\mathbf{W}_{h}^{+}\mathbf{X}_{h})^{-1}\mathbf{X}_{h}^{T}\mathbf{W}_{h}^{+}\mathbf{Y}_{DR}(\mathcal{O};\hat{G},\hat{S})\\\
(\mathbf{X}_{h}^{T}\mathbf{W}_{h}^{+}\mathbf{X}_{h})^{-1}\mathbf{X}_{h}^{T}\mathbf{W}_{h}^{+}\mathbf{Z}\end{pmatrix},$

where $\boldsymbol{\Psi}=\begin{pmatrix}1&0&0&0\\\ 0&h^{-1}&0&0\\\ 0&0&1&0\\\
0&0&0&h^{-1}\end{pmatrix}$. Then by Lemma 6 with Lemma 1,

$\displaystyle
n^{2/5}\begin{pmatrix}\hat{\alpha}_{R,DR}^{FRD,Y}(\hat{G},\hat{S})-\mu^{+}(w_{0})\\\
\hat{\beta}_{R,DR}^{FRD,Y}(\hat{G},\hat{S})-\mu^{\prime+}(w_{0})\\\
\hat{\alpha}_{R}^{FRD,Z}-p^{+}(w_{0})\\\
\hat{\beta}_{R}^{FRD,Z}-p^{\prime+}(w_{0})\end{pmatrix}-\frac{1}{2}\boldsymbol{\Psi}^{-1}\begin{pmatrix}\boldsymbol{\Omega}^{-1}&0\\\
0&\boldsymbol{\Omega}^{-1}\end{pmatrix}\begin{pmatrix}\mu^{\prime\prime+}(w_{0})\boldsymbol{\delta}\\\
p^{\prime\prime+}(w_{0})\boldsymbol{\delta}\end{pmatrix}\overset{d}{\longrightarrow}$
$\displaystyle\mathcal{N}\bigg{(}0,g(w_{0})^{-1}\boldsymbol{\Psi}^{-1}\begin{pmatrix}\boldsymbol{\Omega}^{-1}&0\\\
0&\boldsymbol{\Omega}^{-1}\end{pmatrix}\begin{pmatrix}\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})\mathcal{V}&\eta_{DR}^{+}(w_{0};G_{0},S^{*})\mathcal{V}\\\
\eta_{DR}^{+}(w_{0};G_{0},S^{*})\mathcal{V}&p^{+}(w_{0})(1-p^{+}(w_{0}))\mathcal{V}\end{pmatrix}\begin{pmatrix}\boldsymbol{\Omega}^{-1}&0\\\
0&\boldsymbol{\Omega}^{-1}\end{pmatrix}\boldsymbol{\Psi}^{-1}\bigg{)}.$

Then

$\displaystyle
n^{2/5}\bigg{(}\begin{pmatrix}\hat{\alpha}_{R,DR}^{FRD,Y}(\hat{G},\hat{S})-\mu^{+}(w_{0})\\\
\hat{\alpha}_{R}^{FRD,Z}-p^{+}(w_{0})\end{pmatrix}-\rho^{+}\begin{pmatrix}\mu^{\prime\prime+}(w_{0})\\\
p^{\prime\prime+}(w_{0})\end{pmatrix}\bigg{)}\overset{d}{\longrightarrow}\mathcal{N}\bigg{(}0,\upsilon^{+}\begin{pmatrix}\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})&\eta_{DR}^{+}(w_{0};G_{0},S^{*})\\\
\eta_{DR}^{+}(w_{0};G_{0},S^{*})&p^{+}(w_{0})(1-p^{+}(w_{0}))\end{pmatrix}\bigg{)}.$

Similarly,

$\displaystyle
n^{2/5}\bigg{(}\begin{pmatrix}\hat{\alpha}_{L,DR}^{FRD,Y}(\hat{G},\hat{S})-\mu^{-}(w_{0})\\\
\hat{\alpha}_{L}^{FRD,Z}-p^{-}(w_{0})\\\
\end{pmatrix}-\rho^{-}\begin{pmatrix}\mu^{\prime\prime-}(w_{0})\\\
p^{\prime\prime-}(w_{0})\end{pmatrix}\bigg{)}\overset{d}{\longrightarrow}\mathcal{N}\bigg{(}0,\upsilon^{-}\begin{pmatrix}\sigma_{DR}^{2-}(w_{0};G_{0},S^{*})&\eta_{DR}^{-}(w_{0};G_{0},S^{*})\\\
\eta_{DR}^{-}(w_{0};G_{0},S^{*})&p^{-}(w_{0})(1-p^{-}(w_{0}))\end{pmatrix}\bigg{)}.$

Then

$\displaystyle
n^{2/5}\bigg{(}\begin{pmatrix}\hat{\alpha}_{R,DR}^{FRD,Y}(\hat{G},\hat{S})-\hat{\alpha}_{L,DR}^{FRD,Y}(\hat{G},\hat{S})-(\mu^{+}(w_{0})-\mu^{-}(w_{0}))\\\
\hat{\alpha}_{R}^{FRD,Z}-\hat{\alpha}_{L}^{FRD,Z}-(p^{+}(w_{0})-p^{-}(w_{0}))\\\
\end{pmatrix}\bigg{)}-\begin{pmatrix}\rho^{+}\mu^{\prime\prime+}(w_{0})-\rho^{-}\mu^{\prime\prime-}(w_{0})\\\
\rho^{+}p^{\prime\prime+}(w_{0})-\rho^{-}p^{\prime\prime-}(w_{0})\end{pmatrix}\overset{d}{\longrightarrow}\mathcal{N}(0,\Xi),$

where

$\displaystyle\Xi=\upsilon^{+}\begin{pmatrix}\sigma_{DR}^{2+}(w_{0};G_{0},S^{*})&\eta_{DR}^{+}(w_{0};S^{*})\\\
\eta_{DR}^{+}(w_{0};G_{0},S^{*})&p^{+}(w_{0})(1-p^{+}(w_{0}))\end{pmatrix}+\upsilon^{-}\begin{pmatrix}\sigma_{DR}^{2-}(w_{0};G_{0},S^{*})&\eta_{DR}^{-}(w_{0};G_{0},S^{*})\\\
\eta_{DR}^{-}(w_{0};G_{0},S^{*})&p^{-}(w_{0})(1-p^{-}(w_{0}))\end{pmatrix}.$

By the delta method, the result follows. The proof for
$\hat{\tau}_{FRD}^{IPCW}(\hat{G})$ is similar. ∎

From Theorem 1, we can derive Corollary 1 because the corollary is a special
case of Theorem 1. Next, we prove Theorem 2.

###### Proof of Theorem 2.

By Steingrimsson, Diao, and Strawderman (2019) and Suzukawa (2004), variance
formulae for IPCW and DR estimator are given by

$\displaystyle\sigma_{IPCW}^{2}(w;G_{0})=Var(Y_{IPCW}(O;G_{0})|W=w)=\int_{0}^{\infty}\frac{Y^{2}}{G_{0}(u)}dF_{0}(u|w)-\\{\mu(w)\\}^{2}$
$\displaystyle\sigma_{DR}^{2}(w;G_{0},S)=Var(Y_{DR}(O;G_{0},S)|W=w)=\int_{0}^{\infty}\frac{Y^{2}}{G_{0}(u)}dF_{0}(u|w)-$
$\displaystyle\int_{0}^{\infty}\frac{S_{0}(u|w)\\{Q_{Y}(u,w,S)(Q_{Y}(u,w,S)-2Q_{Y}(u,w,S_{0}))\\}}{\\{G_{0}(u)\\}^{2}}d\bar{G}_{0}(u)-\\{\mu(w)\\}^{2}$
$\displaystyle\sigma_{DR}^{2}(w;G_{0},S_{0})=Var(Y_{DR}(O;G_{0},S_{0})|W=w)=\int_{0}^{\infty}\frac{Y^{2}}{G_{0}(u)}dF_{0}(u|w)-$
$\displaystyle\int_{0}^{\infty}\frac{S_{0}(u|w)\\{Q(u,w,S_{0})\\}^{2}}{\\{G_{0}(u)\\}^{2}}d\bar{G}_{0}(u)-\\{\mu(w)\\}^{2}.$

Hence for any $w$, we have
$\sigma_{DR}^{2}(w,G_{0},S_{0})\leq\min\\{\sigma_{IPCW}^{2}(w;G_{0}),\sigma_{DR}^{2}(w;G_{0},S)\\}$.
Clearly,

$\displaystyle\sigma_{DR}^{2+}(w_{0};G_{0},S_{0})\leq\min\\{\sigma_{IPCW}^{2+}(w_{0};G_{0}),\sigma_{DR}^{2+}(w_{0};G_{0},S)\\}$
$\displaystyle\sigma_{DR}^{2-}(w_{0};G_{0},S_{0})\leq\min\\{\sigma_{IPCW}^{2-}(w_{0};G_{0}),\sigma_{DR}^{2-}(w_{0};G_{0},S)\\}.$

Since the $\upsilon^{+}$ and $\upsilon^{-}$ are positive, the result holds.
For the fuzzy RD estimator,

$\displaystyle\Sigma_{SRD}^{DR}(G_{0},S_{0})-\Sigma_{SRD}^{IPCW}(G_{0})=\frac{1}{\tau^{Z}}[\upsilon^{+}\\{\sigma_{DR}^{2+}(w_{0};G_{0},S_{0})-\sigma_{IPCW}^{2+}(w_{0};G_{0})\\}$
$\displaystyle+\upsilon^{-}\\{\sigma_{DR}^{2-}(w_{0};G_{0},S_{0})-\sigma_{IPCW}^{2-}(w_{0};G_{0})\\}]-2\frac{\tau^{Y}}{(\tau^{Z})^{3}}[\upsilon^{+}\\{\eta_{DR}^{+}(w_{0};G_{0},S_{0})-\eta_{IPCW}^{+}(w_{0};G_{0})\\}$
$\displaystyle+\upsilon^{-}\\{\eta_{DR}^{-}(w_{0};G_{0},S_{0})-\eta_{IPCW}^{-}(w_{0};G_{0})\\}].$

Since calculation at sharp RD case,
$\sigma_{DR}^{2+}(w_{0};G_{0},S_{0})\leq\sigma_{IPCW}^{2+}(w_{0};G_{0})$ and
$\sigma_{DR}^{2-}(w_{0};G_{0},S_{0})\leq\sigma_{IPCW}^{2-}(w_{0};G_{0})$.
Moreover, from conditions
$\eta_{DR}^{+}(w_{0};G_{0},S_{0})\leq\eta_{IPCW}^{+}(w_{0};G_{0})$ and
$\eta_{DR}^{-}(w_{0};G_{0},S_{0})\leq\eta_{IPCW}^{-}(w_{0};G_{0})$, with
$\tau^{Y}>0$ and $\tau^{Z}>0$, we have
$\Sigma_{FRD}^{DR}(G_{0},S_{0})\leq\Sigma_{FRD}^{IPCW}(G_{0})$. Similarly,
$\Sigma_{FRD}^{DR}(G_{0},S_{0})\leq\Sigma_{FRD}^{DR}(G_{0},S)$. Hence

$\displaystyle
AVar(\hat{\tau}_{FRD}^{DR}(G_{0},S_{0}))\leq\min\\{AVar(\hat{\tau}_{FRD}^{IPCW}(G_{0})),AVar(\hat{\tau}_{FRD}^{DR}(G_{0},S))\\}.$

∎

## References

  * 1 Abadie, A. and Imbens, G. W. (2006). Large sample properties of matching estimators for average treatment effects. Econometrica 74, 235-267.
  * 2 Andriole, G. L., Crawford, E. D., Grubb III, R. L., Buys, S. S., Chia, D., Church, T. R., Fouad, M. N., Gelmann, E. P., Kvale, P. A., Reding, D. J. and Weissfeld, J. L. (2009). Mortality results from a randomized prostate-cancer screening trial. New England Journal of Medicine 360, 1310-1319.
  * 3 Angrist, J. D., Imbens, G. W. and Rubin, D. B. (1996). Identification of causal effects using instrumental variables. Journal of the American statistical Association, 91, 444-455.
  * 4 Bai, X., Tsiatis, A. A. and O’Brien, S. M. (2013). Doubly-robust estimators of treatment-specific survival distributions in observational studies with stratified sampling. Biometrics 69, 830-839.
  * 5 Bor, J., Moscoe, E., Mutevedzi, P., Newell, M. L. and Bärnighausen, T. (2014). Regression discontinuity designs in epidemiology: causal inference without randomized trials. Epidemiology 25, 729-737.
  * 6 Calonico, S., Cattaneo, M. D. and Titiunik, R. (2014). Robust nonparametric confidence intervals for regression- discontinuity designs. Econometrica 82, 2295-2326.
  * 7 Calonico, S., Cattaneo, M. D. and Titiunik, R. (2015a). Optimal data-driven regression discontinuity plots. Journal of the American Statistical Association, 110, 1753-1769.
  * 8 Calonico, S., Cattaneo, M. D. and Titiunik, R. (2015b). rdrobust: An R package for robust nonparametric inference in regression-discontinuity designs. R Journal 7, 38-51.
  * 9 Calonico, S., Cattaneo, M. D., Farrell, M. H. and Titiunik, R. (2018). Regression discontinuity designs using covariates. Review of Economics and Statistics 1-10.
  * 10 Cattaneo, M. D., Titiunik, R., Vazquez-Bare, G. and Keele, L. (2016). Interpreting regression discontinuity designs with multiple cutoffs. The Journal of Politics 78(4), 1229-1248.
  * 11 Fan, J. and Gijbels, I. (1994). Censored regression: local linear approximations and their applications. Journal of the American Statistical Association 89, 560-570.
  * 12 Fan, J. and Gijbels, I. (1996). Local polynomial modelling and its applications: monographs on statistics and applied probability 66. CRC Press.
  * 13 Hahn, J., Todd, P. and Van der Klaauw, W. (1999). Evaluating the effect of an antidiscrimination law using a regression-discontinuity design (No. w7131). Technical Report, National bureau of economic research.
  * 14 Hahn, J., Todd, P. and Van der Klaauw, W. (2001). Identification and estimation of treatment effects with a regression-discontinuity design. Econometrica 69, 201-209.
  * 15 Imbens, G. and Kalyanaraman, K. (2012). Optimal bandwidth choice for the regression discontinuity estimator. The Review of economic studies, 79, 933-959.
  * 16 Imbens, G. and Zajonc, T. (2009). Regression discontinuity design with vector-argument assignment rules. Unpublished paper.
  * 17 Imbens, G. W. and Lemieux, T. (2008). Regression discontinuity designs: A guide to practice. J. Econometrics 142, 615-635.
  * 18 Lee, D. S. and Lemieux, T. (2010). Regression discontinuity designs in economics. Journal of economic literature 48, 281-355.
  * 19 Ludwig, J. and Miller, D. L. (2005). Does head start improve children’s life chances? evidence from a regression discontinuity design (No. w11702). Technical report, National Bureau of Economic Research.
  * 20 Ludwig, J. and Miller, D. L. (2007). Does Head Start improve children’s life chances? Evidence from a regression discontinuity design. The Quarterly journal of economics, 122, 159-208.
  * 21 McCrary, J. (2008). Manipulation of the running variable in the regression discontinuity design: A density test. Journal of econometrics, 142, 698-714.
  * 22 Moscoe, E., Bor, J. and Bärnighausen, T. (2015). Regression discontinuity designs are underutilized in medicine, epidemiology, and public health: a review of current and best practice. Journal of clinical epidemiology, 68, 132-143.
  * 23 Rubin, D. and Van der Laan, M. J. (2007). A doubly robust censoring unbiased transformation. Int. J. Biostat. 3, 1-19.
  * 24 Shoag, J., Halpern, J., Eisner, B., Lee, R., Mittal, S., Barbieri, C. E. and Shoag, D. (2015). Efficacy of prostate-specific antigen screening: Use of regression discontinuity in the PLCO cancer screening trial. JAMA oncology, 1, 984-986.
  * 25 Steingrimsson, J. A., Diao, L., Molinaro, A. M. and Strawderman, R. L. (2016). Doubly robust survival trees. Statistics in medicine, 35, 3595-3612.
  * 26 Steingrimsson, J. A., Diao, L. and Strawderman, R. L. (2019). Censoring unbiased regression trees and ensembles. Journal of the American Statistical Association, 114, 370-383.
  * 27 Suzukawa, A. (2004). Unbiased estimation of functionals under random censorship. Journal of the Japan Statistical Society, 34, 153-172.
  * 28 Thistlethwaite, D. L. and Campbell, D. T. (1960). Regression-discontinuity analysis: An alternative to the ex post facto experiment. Journal of Educational Psychology 51, 309-317.
  * 29 Tsiatis, A. (2007). Semiparametric theory and missing data. Springer Science & Business Media.
  * 30 Yang, M. (2013). Treatment effect analyses through orthogonality conditions implied by a fuzzy regression discontinuity design, with two empirical studies. Dep. Econ. Rauch Bus. Cent. Lehigh Univ., Bethlehem, PA 308, 1-61.
  * 31 Zajonc, T. (2012). Essays on Causal Inference for Public Policy. PhD thesis, Harvard University.