Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
80,947 characters · 13 sections · 56 citation commands
Regression discontinuity design with right-censored survival data
The {{regression discontinuity design}} is a widely used technique for causal inference in economics, political science, and sociology (see van2008regression, imbens2008regression, lee2010regression, and cattaneo2022regression for recent reviews). In the {{regression discontinuity design}}, possible confounders are {`}controlled for{'} by exploiting that units are assigned to treatment based on whether their value of an observed covariate is above or below some known cut-off, the idea being that units with values of this observed covariate just above the known cut-off are similar to the subjects with values just below the cut-off.
For situations where the observed outcomes are noncensored random variables, the inference theory, building on that of local polynomial regression, is well developed (see e.g., fan1996local, hahn2001identification, and calonico2014robust,calonico2014supplement, and the reviews above). In this paper we extend the {{regression discontinuity design}} to the survival analysis setting. In particular, we study the {{regression discontinuity design}} applied to right-censored survival data in an intensity based counting process framework (see ABGK93, and ABG08,aalen2010history). The survival analysis models we study are all instances of the broad class of multiplicative intensity models, meaning that the modelling and inference revolves around the hazard rate function.
The setup is as follows: On a probability space $(\Omega,\mathscr{H},{\rm Pr})$, let $Z \in \mathbb{R}$ be an observed covariate; $X\in \{0,1\}$ a treatment indicator; $U$ a vector of unobserved confounders; $\widetilde{T}^0 \geq 0$ the potential outcome for the non-treated, and $\widetilde{T}^1 \geq 0$ the potential outcome under treatment (for potential outcomes theory, see, for example, holland1986statistics, morgan2015counterfactuals, imbens2015causal, and imbens2020potential\footnote{imbens2020potential is an excellent and very readable article comparing the potential outcomes theory to the directed acyclic graph theory associated with Judea Pearl. See pearl2009causality for the canonical treatise on DAGs, and pearl2018why or pearl2016causal for more accessible accounts.}). Since $\widetilde{T}^0$ and $\widetilde{T}^1$ take values on the positive half of the real line, we refer to these two random variables as {\it potential lifetimes}. That we are in a regression discontinuity setting, means that treatment is determined by the value of an observed covariate. We take $X = I\{Z \geq z_0\}$ for some known cut-off $z_0$, where $I\{A\}$ is the indicator function of the event $A$.\footnote{What we describe here is the {\it sharp} regression discontinuity design. In contrast, a {\it fuzzy} regression discontinuity design occurs when the probability of receiving treatment does not jump from zero to one at the cut-off, rather $z \mapsto {\rm Pr}(X = 1\mid Z = z)$ has a point of discontinuity at the cut-off. This paper deals only with the sharp regression discontinuity design.} The two potential lifetimes $\widetilde{T}^0$ and $\widetilde{T}^1$ are assumed to stem from distributions with hazard rate functions $\alpha_0(t,Z,U)$ and $\alpha_1(t,Z,U)$, respectively, meaning that
In analogy with the average treatment effect (the {{\sc ate}}) and conditional or local average treatment effects ({{\sc cate}}s or {{\sc late}}s) (see, e.g., imbens2015causal), the causal estimands of the present paper are all defined in terms $\alpha_0(t,Z,U)$ and $\alpha_1(t,Z,U)$ when the confounder $U$ is averaged out. Specifically, with
the main estimand in this paper is
evaluated in the cut-off $z = z_0$. The estimands $\theta(t,z_0)$ and $\Theta(t,z_0)$ are, we contend, natural survival analysis counterparts of the classical estimand in the standard {{regression discontinuity design}}, namely the average treatment effect at the cut-off, see, for example, Eq. (2.1) in imbens2008regression.
The paper proceed as follows. In Section (ref) we provide a general presentation of the potential outcomes framework as it applies to right-censored survival data studied in a counting process framework. In Section (ref), we narrow in on the {{regression discontinuity design}} and present a lemma that is key to making the {{regression discontinuity design}} feasible in the present setting. Section (ref) contains a discussion of various causal estimands, and assumptions under which they may be identified. In Section (ref) we lay out the assumptions made about the model generating the data, and present a special case of the covariate-localised Aalen additive hazards estimator. Section (ref) contains large-sample theory for the general version of the estimator, with results on variance estimation and bias correction in Sections (ref) and (ref), respectively.
Let $(\widetilde{T}^0,\widetilde{T}^1,X,Z,U,C)$ be random variables on the probability space $(\Omega,\mathscr{H},{\rm Pr})$. Here, $\widetilde{T}^0$ and $\widetilde{T}^1$ are the potential lifetimes corresponding to whether a subject is not treated or treated, respectively; $X \in \{0,1\}$ is an indicator of treatment; $Z$ is a vector of observed covariates; $U$ is a vector of unobserved possible confounders;\footnote{A confounder is a covariate that is associated with the outcome {\it and} with the treatment variable of interest. In the multiplicative intensity models studied here, this means that a covariate, to be a confounder, must have an effect on the hazard rate {\it and} be correlated with the treatment variable. If a covariate only satisfy one of these two conditions, then it is not a confounder.} and $C \geq 0$ is a censoring variable. We use $g$ as an index for untreated ($g = 0$) and treated ($g = 1$), and write, for example, $\widetilde{T}^g$ to indicate one of the two potential lifetimes. The units under study are observed over the time interval $[0,\tau]$, where $\tau < \infty$. Assume that, conditionally on $(Z,U)$, the potential lifetimes stem from distributions with hazard rate functions $\alpha_{0}(t,Z,U)$ and $\alpha_{1}(t,Z,U)$.\footnote{Here we model the unconfoundedness assumption $(\widetilde{T}^0,\widetilde{T}^1)\perp \!\!\! \perp X \mid (Z,U)$ directly. See yadlowsky2018bounds. In the {{regression discontinuity design}} this assumption is trivially satisfied imbens2008regression.}
What distinguishes the potential outcomes theory when right-censoring is present from when there is no censoring, is that with no censoring only one of the two potential lifetimes is observed for the same unit, while for right-censored data {\it at most} one of the potential lifetimes is actually observed for the same unit. This means that when there is no censoring we observe
while, when the data are right-censored we only observe $\widetilde{T}$ if it is smaller than the censoring time $C$, that is
where $T^0 = \min(\widetilde{T}^0,C)$ and $T^1 = \min(\widetilde{T}^1,C)$ are the possibly right-censored potential lifetimes. Throughout the paper, we assume that the two potential lifetimes are conditionally independent given the covariates, which we express by
and we work under the assumption of {\it random right-censoring} (see ABGK93), which means that $(\widetilde{T}^0,\widetilde{T}^1)$ is independent of the censoring time $C$ given $(X,Z,U)$, with symbols
Two potential lifetimes in turn leads to two indicators of noncensoring,
and to two potential counting and potential at-risk processes
Let $\mathcal{X} = \sigma(X,Z,U)$, define the filtrations $\mathscr{E}_{t}^g = \sigma(\{N^g(s),Y^g(s)\}_{s\leq t})$ for $g = 0,1$, and set
It is assumed that $X$, $Z$, and $U$ are realised at time zero, meaning that $\mathscr{G}_0^g = \mathcal{X}$. With respect to $\mathscr{G}_t^0$ and $\mathscr{G}_t^1$ we have, using the assumptions in (ref) and (ref) (see ABGK93), two {`}potential{'} local square integrable martingales $M^0$ and $M^1$, respectively, given by
This ends our description of our basic modelling assumptions. The reader familiar with the counting process approach to, and martingale methods in, survival analysis will see that the above is nothing more than what one gets the when taking the potential outcomes framework seriously and applying it to the standard intensity based counting process approach to survival analysis. In fact, the presentation so far is just a slight notational and semantic reformulation of the theory presented in Chapter III.2 of ABGK93. It is important to notice that the modelling undertaken up to this point has taken place in the two potential worlds, so to speak, culminating in the {`}potential world{'} martingales of (ref). We now proceed to the consequences of this model for quantities of this world, that is, the observed quantities. Denote the observed counting process
and the observed at-risk process
Let $\mathcal{E}_{t} = \sigma(\{N(s),Y(s)\}_{s\leq t})$ be the filtration generated by the observable counting and at-risk processes, and let $\mathcal{X}^{\rm obs} = \sigma(X,Z)$ be the {$\sigma$-algebra} generated by the observed covariates. Define the filtrations
where, as above, $\mathcal{G}_{0} = \mathcal{X}$ and $\mathcal{F}_0 = \mathcal{X}^{\rm obs}$. The next lemma says that the intensity process of the observed counting process $N$ with respect to $\mathcal{G}_t$ takes the form it it intuitively should take.
An important, though rather intuitive, thing to note in the preceding lemma is that it does {\it not} say that $M^g$ is a $\mathcal{G}_t$-martingale. Instead, the lemma says that $I_{X = g}M^g$ is a $\mathcal{G}_t$-martingale. This is intuitive because $\mathcal{G}_t$ is generated by the observables $N$ and $Y$ in (ref) and (ref), respectively, and if $X=0$, for example, then the path of $N^1$ is certainly not observable.
In the next section we specialise the potential outcomes model for right-censored survival data to the regression discontinuity setting.
We retain the definitions and assumptions from the previous section, with the following two exceptions: To conform with the {{regression discontinuity design}}, we require $Z \in \mathbb{R}$ and set $X = I\{Z \geq z_0\}$ for some known cut-off $z_0$. Henceforth, we often refer to $Z$ as the {\it forcing variable}. Notice also that $\sigma(X)$ is in included in $\sigma(Z)$, with the consequence that $\mathcal{X} = \sigma(Z,U)$ and $\mathcal{X}^{\rm obs} = \sigma(Z)$, and thus $\mathcal{G}_t = \mathcal{E}_t \vee \sigma(Z,U)$ and $\mathcal{F}_t = \mathcal{E}_t \vee \sigma(Z)$ (compare with (ref)).
In order to make what follows clear, we now take a short detour via the standard regression discontinuity design, as presented, for example, in imbens2008regression or calonico2014robust. Let $y^0$ and $y^1$ be two real valued potential outcomes,\footnote{We use lowercase letters to distinguish these two potential outcomes from the potential at-risk processes introduced in Section (ref).} $Z \in \mathbb{R}$ the forcing variable, and $X = I\{Z \geq z_0\}$ the treatment indicator. The observed outcome is $y = Xy^1 + (1 - X)y^0$, and the estimand of interest is $\tau_{\rm srd} = {\rm E}\,( y^1 - y^0 \mid Z = z_0)$, called the average treatment effect at the cut-off (or threshold). To estimate $\tau_{\rm srd}$, the limits $\lim_{z \uparrow z_0}{\rm E}\,( y \mid Z = z)$ and $\lim_{z \downarrow z_0}{\rm E}\,( y \mid Z = z)$ are estimated using local polynomial regressions to the left and to the right of the cut-off, respectively. Since the roles of the hazard rate functions $\alpha_0(t,Z,U)$ and $\alpha_1(t,Z,U)$ in multiplicative intensity models are analogoues to those of the conditional expectations ${\rm E}\,(y^0 \mid Z)$ and ${\rm E}\,(y^1 \mid Z)$ in the standard {{regression discontinuity design}}, we would like to work with hazard rate functions that only depend on the forcing variable (and time). More to the point, if ${\rm E}\,\{\alpha_g(t,Z,U) \mid \sigma(Z)\}$ was the hazard rate function of $I_{X = g} N^g$ with respect to the filtration of observables, then the standard local polyonimial regression theory would be straighforward to mimic. This is not quite the case, but nearly, in a sense made precise by the next lemma.
Our application of this lemma is when $\xi_t$ is one of the hazard rate functions $\alpha_{0}(t,Z,U)$ or $\alpha_{1}(t,Z,U)$. The importance of this lemma derives from the fact that $Y(t){\rm E}\,(\xi_t \mid \mathcal{F}_{t-})$ is, at first sight, a complicated function as it may depend on the forcing variable $Z$ as well as the paths of $N(t)$ and $Y(t)$. But an implication of Lemma (ref), since $Y(t)$ is left-continuous, is that
and since $\mathcal{Z}_t = \{\{T \geq t\},\{T < t\},\emptyset,\Omega\}\vee \sigma(Z)$, the conditional expectation $Y(t) {\rm E}\,(\xi_t \mid \mathcal{Z}_t)$ is a function of the forcing variable only (see, for example, Lemma 1.13, p. 7 and the discussion on p. 106 in kallenberg2002foundations). This motivates defining the functions,\footnote{ That ${\rm E}\,(\xi_t \mid \mathcal{Z}_t) = {\rm E}\,\{I_{T \geq t}\xi_t \mid \sigma(Z)\}/{\rm Pr}\{T \geq t\mid \sigma(Z)\} + {\rm E}\,\{I_{T < t}\xi_t \mid \sigma(Z)\}/{\rm Pr}\{T < t\mid \sigma(Z)\}$ almost surely, can be proved following the steps in Ex. 34.4 of billingsley1995probability.}
Notice that if $\alpha_g(t,Z,U) \leq g(Z,U)$ for some $g$ such that ${\rm E}\,g(Z,U) < \infty$, then we can pass the derivative under the intergal sign in $\partial/\partial t\,{\rm E}\,\{S_g(t,Z,U) \mid \sigma(Z)\} = {\rm E}\,\{\partial/\partial tS_g(t,Z,U) \mid \sigma(Z)\}$,\footnote{The dominated convergence theorem extends to conditional expectations, see, e.g. cohen2015stochastic} and consequently we have the relation $I_{X = g}\bar{\alpha}_g(t,z) = - I_{X = g}(\partial \bar{S}_g(t,z)/\partial t)/ \bar{S}_g(t,z)$. We summarise the above in a lemma that is used repeatedly in the remainder of the paper.
If the difference $\alpha_{1}(t,z,u) - \alpha_{0}(t,z,u)$ is functionally independent of $u$, meaning that $\alpha_{1}(t,Z,Y) - \alpha_{0}(t,Z,U)$ is a $\sigma(Z)$-measurable random variable, then
where $\bar{\alpha}_g(t,z)$ for $g = 0,1$ are defined in (ref) and $\theta(t,z)$ is defined in (ref). The functional independence assumption just introduced is crucial for whether $\theta(t,z_0)$ is identifiable or not. We state the functional independence assumption here for easy reference, note, however, that it is not assumed throughout.
If this assumption is dropped, we are only able to estimate the average treatment effect at the cut-off {\it among those at-risk}, a quantity we denote $\theta_{\rm risk}(t,z_0)$, it is
with cumulative $\Theta_{\rm risk}(t,z) = \int_0^t \theta_{\rm risk}(s,z)\,{\rm d} s$.\footnote{By the definition in (ref) we mean that $\theta_{\rm risk}(t,Z)$ is the $\sigma(Z)$-measurable function such that $ Y(t){\rm E}\,\{\alpha_{1}(t,Z,U) - \alpha_{0}(t,Z,U)\mid \mathcal{Z}_t\} = Y(t)\theta_{\rm risk}(t,Z)$ almost surely.} The assumption of the difference $\alpha_1(t,Z,U) - \alpha_0(t,Z,U)$ being functionally independent of the confounder is crucial for whether we are estimating $\theta(t,z_0)$ or the at-risk version $\theta_{\rm risk}(t,z_0)$, but is otherwise immaterial to the theory developed in the subsequent sections.
In this section, we first elaborate on the assumptions made about what are in statistics and econometrics jargon, respectively, called the true model or the data generating process. Subsequently, in Section (ref) we introduce a special case of our estimator, and provide some theory for this special case. The general large-sample theory is deferred to Section (ref).
Let $(\widetilde{T}_i^0,\widetilde{T}_i^1,Z_i,U_i,C_i),\, i = 1,\ldots,n$ be independent replicates of $(\widetilde{T}^0,\widetilde{T}^1,Z,U,C)$, where this latter is as described in Sections (ref) and (ref). In particular, the forcing variable $Z \in \mathbb{R}$ and $X = I\{Z \geq z_0\}$ for a known cut-off $z_0$. This entails that $(N_1,Y_1),\ldots,(N_n,Y_n)$ are independent replicates of $(N,Y)$ with $N(t) = XN^1(t) + (1 - X)N^0(t)$ and $Y(t) = XY^1(t) + (1 - X)Y^0(t)$ as defined in Section (ref). As above, we assume that the processes $N_i$ and $Y_i$ are observed over the finite time interval $[0,\tau]$. The filtrations $\mathcal{G}_t$ and $\mathcal{F}_t$ are now $\mathcal{G}_t = \mathcal{E}_t \vee \sigma(Z_i,U_i,\,i = 1,\ldots,n)$ and $\mathcal{F}_t = \mathcal{E}_t \vee \sigma(Z_i,\,i = 1,\ldots,n)$, with $\mathcal{E}_t^g = \sigma(\{(N_i(s),Y_i(s)\}_{s \leq t},\,i = 1,\ldots,n)$. For $g = 0,1$, denote $M_1^g,\ldots,M_n^g$ the independent replicates of the martingale $M^g$ in (ref). These are orthogonal local square integrable martingales with respect to the filtration $\mathscr{G}_t^g$. Similarly, $I_{X_1 = g}M_1^{g}(t),\ldots,I_{X_n = g}M_n^{g}(t)$ are orthogonal local square integrable martingales with respect to the filtration $\mathcal{G}_t$ (see Lemma (ref)), and $I_{X_1 = g}\bar{M}_1^{g}(t),\ldots,I_{X_n = g}\bar{M}_n^{g}(t)$ are orthogonal local square integrable martingales with respect to the filtration $\mathcal{F}_t$ (see Lemma (ref)).
The survival functions associated with the hazard rates $\alpha_{0}(t,Z,U)$ and $\alpha_{1}(t,Z,U)$ are denoted $S_{g}(t,Z,U) = \exp\{- \int_0^t \alpha_{g}(s,Z,U)\,{\rm d} s\}$ for $g = 0,1$, and we set
Let $y_g(t,z)$ to be the conditional expectations of the at-risk process $Y^g(t)$ given $Z = z$. Using the assumption in (ref), these functions are
where $H(t)$ is the distribution function of the censoring variable $C$. Without further mention, the following is assumed throughout the paper
Assumption (ref) gives the martingale representation in (ref), and is the key to Lemma (ref), and thereby also to Lemma (ref). Assumption (ref) is standard in survival analysis, and is equivalent to requiring absolute continuity of the survival functions. The assumption also entails that the conditional hazards $\bar{\alpha}_0(t,z)$ and $\bar{\alpha}_1(t,z)$ are continuous in $t$ for all $z$. Assumption (ref) is needed to ensure that $y_g(t,z)$ is bounded below (a fact that is, for example, used in the proof of Theorem (ref)). Assumption (ref) allows for the estimation of the parameters of interest using a regression discontinuity design. For a discussion of the analogue of this assumption in the standard regression discontinuity design, see hahn2001identification, and imbens2008regression. In the next section a special case of our main estimator is presented.
The estimator we propose for $\Theta(t,z_0)$ is a weighted version of the Aalen additive hazards estimator (Aalen80,Aalen89,Aalen93, ABG08, and ABGK93).\footnote{Throughout this section we assume, for notational convenience, that Assumption (ref) holds. If this assumption does not hold, then all the theory of this section translates directly to the estimation of $\Theta_{\rm risk}(t,z_0)$.} The idea is to fit local polynomial regression models in the vicinity of the cut-off $z_0$. Thus, {`}local{'} here refers to an interval on the real line on which the forcing variable $Z$ takes its values, and does {\it not} refer to the time axis. Specifically, the conditional expectations $\bar{\alpha}_{0}(t,z)$ and $\bar{\alpha}_{1}(t,z)$, as defined in (ref), are approximated on intervals to the left and to the right of $z_0$ by $p$th order local polynomials. Since, for $g = 0,1$, the conditional expectation $\bar{\alpha}_{g}(t,z)$ may be approximated by
where $\bar{\alpha}_{g}^{(\nu)}(t,z)$ is the $\nu$th derivative of $\bar{\alpha}_{g}(t,z)$ with respect to $z$, the local polynomial regression version of the Aalen additive hazards estimator that we introduce, is an estimator of $\int_0^t \bar{\alpha}_{g}^{(\nu)}(s,z)/\nu!\,{\rm d} s$ for $\nu = 0,1,\ldots,p$. We start by presenting the local linear estimator, and then move on to general results for $p$th order local polynomial estimators in Section (ref).
For a bandwidth $h > 0$ and a kernel function $K_h(u) = K(u/h)/h$ with $K(u) = k(-u)I\{u < 0\} + k(u)I\{u \geq 0\}$,\footnote{In principle, different kernels could be used on either side of the cut-off. For simplicity, we employ the same kernel on both sides as this does not affect the theory.} where $k$ is a function satisfying Assumption (ref), the local linear estimator $(\widehat{B}_{g,1},\widehat{B}_{g,1}^{(1)})$ is, for $g = 0,1$, given by
with
This is seen to be a version the Aalen additive hazards estimator estimator with kernel weights on the forcing variable. An estimator for $\Theta(t,z_0)$ is then given by
For $t$ fixed, this estimator is simply the difference between the intercepts of two local linear regression to the left and to the right of the cut-off. In other words, it is the analogue of the most common estimator of the average treatment effect at the cut-off in the standard regression discontinuity design, see, for example, the estimator $\widehat{\tau}_{\rm srd}$ in imbens2008regression.
Under the assumptions of Lemma (ref) below, the bias of the estimator $\widehat{\Theta}(t,h)$ can be described by
with
Here $\bar{\alpha}_{g}^{(2)}(s,z),\,g = 0,1$ are the unknown second derivates of the conditional expectations defined in (ref); while $e_{1,0} = (1,0)^{{\rm t}}$, and $\Gamma_1$, and $\vartheta_{1,2}$ are quantities that only depend on the chosen kernel, and need not be estimated from the data. They are
Here and elsewhere in the paper, the notation is borrowed from calonico2014robust,calonico2014supplement. If the bandwidth minimising the mean squared error is chosen, namely $h_n = cn^{-1/5}$ for some constant $c > 0$, then, conditionally on the filtration $\mathcal{F}_t$ (here Lemma (ref) is invoked), we have process convergence in the space of c{\`a}dl{\`a}g functions on $[0,\tau]$,
as $n \to \infty$, where $\bar{M}_{0,1}$ and $\bar{M}_{1,1}$ are independent bivariate Gaussian martingales with variation processes
The general version of this result is the content of Corollary (ref). A consistent estimator of $\langle \bar{M}_{g,1},\bar{M}_{g,1}\rangle_t$ is introduced in Section (ref).
To avoid or to get rid of the bias term that appears on the right hand side of (ref), two approaches are discussed in this paper. First, one can choose a bandwidth $h_n$ such that $nh_n^{5}\to 0$, meaning that the bandwidth must be {`}smaller{'} than the mean squared error optimal one. Second, the bias term may be removed by subtracting off a consistent estimate of the bias term. The standard approch is the following (see, e.g., fan1996local): Suppose that ${\rm Bias}_n(t,b_n)$ is consistent for $\int_0^t {\rm bias}(s)\,{\rm d} s$ (uniformly in $t$) as $n \to \infty$ and the so-called pilot bandwidth $b_n \to 0$, then (see Corollary (ref))
as $n \to \infty$, provided, among other things, that $h_n/b_n \to 0$. This latter condition ensures that the variability of the bias correction estimation disappears, meaning that $\bar{M}_{0,1}$ and $ \bar{M}_{1,1}$ are the Gaussian martingales from (ref) (the variation process does not change). As pointed out in the influential paper calonico2014robust, $h_n/b_n$ is never zero in finite samples, and therefore, the variability associated with the bias correction ought to be accounted for in the limiting distribution. In the present paper, results of this type are presented in Section (ref).
In this section we first consider the general $p$th order local polynomial regression estimator of the $\nu$th derivative function $\int_0^t \bar{\alpha}_g^{(\nu)}(s,z_0)\,{\rm d} s$, and derive a representation for the bias of this estimator. Next, in Section (ref), we present two central limit theorems for this estimator. Throughout this section, Assumption (ref) (the functional independence assumption) is, for notational convenience, assumed to hold. If this assumption does not hold, all subsequent results are true with $\theta(t,z_0)$ replaced by $\theta_{\rm risk}(t,z_0)$ (see the discussion in Section (ref)). Moreover, $g = 0,1$ is the index used to indicate the potential outcomes, or statistics depending on these, for non-treated and treated, respectively. Since the theory we develop is the same on both sides of the cut-off, a result concerning a quantity with subscript $g$ means that it applies for both $g = 0$ and $g = 1$. Most of the proofs of the claims made in the present section are deferred to the appendices.
For an integer $p \geq 1$, let $r_p(x) = (1,x,x^2,\ldots,x^p)^{{\rm t}}$, and define the estimator $\widehat{B}_{g,p}$ by
with
and
If $J_{n,h}(t) = 0$, we take ${\rm d} \widehat{B}_{g,p}(t,h) = 0$ for $g = 0,1$. Let $H_p(h) = {\rm diag}(1,h^{-1},\ldots,h^{-p})$. Using that $H_p(h)r_p(z) = r_p(z/h)$ and $r_p(z) = H_{p}(h)^{-1}r_p(z/h)$, the estimator above can be expressed as
with
and, since $G_{g,p,n}(t,h)$ is positive definite if and only if $\Gamma_{g,p,n}(t,h)$ is positive definite, $J_{n,h}(t) = I\{\text{$\Gamma_{0,p,n}(t,h)$ and $\Gamma_{1,p,n}(t,h)$ are positive definite}\}$. Let $e_{p,\nu}$ be the $(p+1)$-dimensional column vector with its $(\nu+1)$th element equal to $1$, and all other elements equal to zero, for example, $e_{1,0} = (1,0)^{{\rm t}}$, $e_{2,1} = (0,1,0)^{{\rm t}}$, $e_{3,2} = (0,0,1,0)^{{\rm t}}$, and so on. Denote $\bar{\alpha}_{g}^{(\nu)}(t,z)$ the $\nu$th partial derivative of $\bar{\alpha}_{g}(t,z)$ with respect to $z$, so $\bar{\alpha}_{g}^{(0)}(t,z) = \bar{\alpha}_{g}(t,z)$. With this notation,
are estimators for $\int_0^t \bar{\alpha}_0^{(\nu)}(s,z_0)\,{\rm d} s$ and $\int_0^t \bar{\alpha}_1^{(\nu)}(s,z_0)\,{\rm d} s$, respectively. An estimator for the $\nu$th derivative of $\Theta(t,z) = \int_0^t \theta(s,z)\,{\rm d} s$ with respect to $z$ and evaluated in $z_0$ (see (ref)), based on a $p$th order polynomial regression estimator, is then
Before we state the bias lemma, we also need the following quantities,
and
The matrices $\Gamma_p$ and $\Psi_p$ are of dimension $(p+1)\times(p+1)$, while $\vartheta_{p,q}$ are $(p+1)$-dimensional column vectors. It will be assumed throughout that the kernel $K$ is chosen so that $\Gamma_p$, for all relevant $p$, is positive definite. This entails that the $\Gamma_{p}^{-1}$ is also positive definite, and that $\lVert\Gamma_{p}\rVert$ and $\lVert\Gamma_{p}^{-1}\rVert$ are finite, with $\lVert\cdot\rVert$ here denoting the matrix norm. The $\mathcal{F}_t$-compensator of $I_{X= g}N^g(t)$, namely $I_{X=g}\int_0^tY^g(s)\bar{\alpha}_g(s,Z)\,{\rm d} s$, is absolutely continuous. This entails that the $\mathcal{F}_t$-compensator of $\widehat{B}_{g,p}^{(\nu)}(t,h)$ is also absolutely continuous, the derivative of this compensator therefore exists, and we denote these derivatives ${\rm E}\,({\rm d} \widehat{B}_{g,p}^{(\nu)}(t,h) \mid \mathcal{F}_{t-})/{\rm d} t$, and ${\rm E}\,({\rm d} \widehat{A}_{g,p}^{(\nu)}(t,h) \mid \mathcal{F}_{t-})/{\rm d} t = e_{p,\nu}^{{\rm t}}\nu!{\rm E}\,({\rm d} \widehat{B}_{g,p}^{(\nu)}(t,h) \mid \mathcal{F}_{t-})/{\rm d} t$. Introduce the column vectors of partial derivatives
and the vector valued functions
We can now state the bias lemma.
The central limit theorems coming up concern the sequences
properly normalised. Since $\widehat{\Theta}_{p}^{(\nu)}(t,h) = \widehat{A}_{1,p}^{(\nu)}(t,h) - \widehat{A}_{0,p}^{(\nu)}(t,h)$, and the theory for both sides of the cut-off are the same, we consider only $\widehat{A}_{g,p}^{(\nu)}(t,h) - \int_0^t \bar{\alpha}_g^{(\nu)}(s,z_0)\, {\rm d} s$ for a generic $g = 0,1$. These sequences are
and can be decomposed in different ways depending on which filtration we analyse it with respect to. From Lemma (ref) we have that $I_{X=g}M_i^g(t) = I_{X=g}\{N_i^g(t) - \int_0^tY_i^g(s)\alpha_g(s,Z_i,U_i)\,{\rm d} s\}$ are martingales with respect to the filtration $\mathcal{G}_t = \mathcal{E}_t \vee \mathcal{X}$, and from Lemma (ref) that $\bar{M}_i^g(t) = I_{X_i = g}\{N_{i}^g(t) - \int_0^t Y_i^g(s) \bar{\alpha}(s,Z_i)\,{\rm d} s\}$ are martingales with respect to the filtration $\mathcal{F}_{t} = \mathcal{E}_t \vee \mathcal{X}^{\rm obs}$ of observables. Recall that $\mathcal{X}^{\rm obs} \subset \mathcal{X}$, and that the confounders are not measurable with respect to $\mathcal{F}_t$ (see Sections (ref) and (ref) for the definitions). This, in turn, leads to two decompositions of the sequence in (ref), and, consequently, to two different central limit theorems for this sequence. Write,
where ${\rm Bias}_{g,p,n}(t,h)$ is, using Lemma (ref),
The first term on the right in (ref) can be written
where $\bar{M}_{g,p,n}(t,h)$ is a martingale with respect to the filtration $\mathcal{F}_t$, namely
where $\bar{M}_i^g$ are the $\mathcal{F}_t$-martingales of Lemma (ref).
From the above theorem it is seen that the predictable quadratic variation of $\bar{M}_{g,p,n}$ is $\langle \bar{M}_{p,n}(\cdot,h),\bar{M}_{p,n}(\cdot,h)\rangle_t = O_p((nh)^{-1})$, so, in particular
Moreover, from Lemma (ref) the bias is seen to be
This shows that the mean squared error optimal bandwidth is $h = c n^{1/(2p+3)}$ for some constant $c > 0$. Combining the bias lemma with Theorem (ref), we get the following corollary.
The two results above, Theorem (ref) and Corollary (ref), pertain to the sequences of $\mathcal{F}_t$-martingales $H_p(h)\bar{M}_{g,p,n} = \widehat{B}_{g,p}(t,h) - \int_{0}^t {\rm E}\,\{{\rm d}\widehat{B}_{g,p}(s,h)\mid \mathcal{F}_{s-}\}$. Considering these martingales essentially means that, in the central limit theorem, we average out the confounders but condition on the forcing variable.\footnote{There is a parallel here to the difference between the observed information and the Fisher information in likelihood inference for regression models, where, in the latter, the covariates are averaged out.} We now turn to a central limit theorem relative to the filtration $\mathcal{G}_t = \mathcal{E}_t \vee \mathcal{X}$, that is, the filtration with respect to which the confounder is measurable. Because the $\mathcal{F}_t$-martingales $\bar{M}_{g,p,n}$ analysed above are {\it not} martingales with respect to $\mathcal{G}_t$, these results are slightly more involved. Consider the decomposition
where ${\rm Bias}_{g,p,n}(t,h)$ is as defined in (ref); $M_{g,p,n}(t,h)$ is the $\mathcal{G}_t$-martingale
see Lemma (ref), and $L_{g,p,n}(t,h)$ is
where $Q_{g,p,n}(s,h)$ is the average of i.i.d. random variables given by
with $\Delta_{g,i}(t)$ being the difference between true hazard and the conditional hazards defined in (ref), that is
It turns out that the $\mathcal{G}_t$-martingales $M_{g,p,n}$ have the same limiting variance as the $\mathcal{F}_t$-martingales $\bar{M}_{g,p,n}$, however, the sequence $L_{g,p,n}$ is of the same order as the martingales, and thus contribute to the asymptotic variance of the estimator sequence. One might therefore view the variance stemming from $L_{g,p,n}$ as extra variance induced by not averaging out the confounder $U$. For the next theorem, in addition to Assumptions (ref)--(ref), we impose a Lipschitz condition and a boundedness condition on the hazard rate functions of the potential lifetimes, as well as on the conditonal hazard rate functions. These assumptions are likely much stronger than necessary.
For the results in the remainder of the paper, we use the $\mathcal{F}$-central limit theorem in Theorem (ref). There are two reasons for this. First, estimating the limiting variance of the martingale is more straightforward than estimating the variances and covariances appearing in Theorem (ref). Second, the common approach in the regression discontinuity literature is to condition on the forcing variable, and average out all other covariates (see the discussion in the Section (ref)).
The probability limit of the process $nh\langle \bar{M}_{g,p,n}(\cdot,h),\bar{M}_{g,p,n}(\cdot,h)\rangle_t$ is $\langle \bar{M}_{g,p},\bar{M}_{g,p}\rangle_t$, for which an expression is given in (ref). Consistent estimators for the variance processes $e_{p,\nu}^{{\rm t}}H_p(h_n)\langle \bar{M}_{g,p},\bar{M}_{g,p}\rangle_tH_p(h_n)e_{p,\nu}$ can be developed along the lines of the standard variance estimator for the variance of the estimator of the cumulative regression coefficients in the Aalen additive hazards model (see, for example, hjort2021partly). One such estimator is
This means that a consistent estimator of the limiting variance of the sequence $(nh^{2\nu + 1})^{1/2}\{\widehat{A}_{g,p}^{(\nu)}(t,h) - \int_0^{t}\bar{\alpha}_g^{(\nu)}(s,z_0)\,{\rm d} s \}$ of Corollary (ref) is $h^{2\nu} (\nu!)^2 e_{p,\nu}^{{\rm t}} V_{g,p,n}(t,h)e_{p,\nu}$. In particular, with $\Phi(z)$ the standard normal cumulative distribution function,
are approximate $(1 - \alpha)100$ percent pointwise confidence intervals provided $nh^{2p+3} \to 0$. If $nh^{2p+3} \to c$ for $c>0$, these confidence intervals are not valid due to the bias term appearing in the limit in Corollary (ref). The topic of the next section is how this may be fixed.
In this section we follow the conventional bias correction approach in local polynomial regression (see, e.g., fan1996local) and study the estimators given by
where
Here $b$ is a so-called pilot bandwidth sequence, typically larger than $h$ (because it is harder to estimate the $(p+1)$th derivative than the $\nu$th derivative when $\nu \leq p$, as it is here). In (ref) the second term on the right is an estimator of the bias term appearing in the limiting distribution of Corollary (ref). Under certain conditions on the bandwidth sequences $h$ and $b$ made precise below, subtracting off a bias estimate removes the asymptotic bias term from the limiting distribution in Corollary (ref), even when $nh^{2p+3}$ tends to a positive constant. In particular, the mean squared error optimal bandwidth $h \propto n^{-1/(2p+3)}$ may be employed, and the limiting martingale processes are the same as those given in said corollary.
As pointed out by calonico2014robust, this large-sample approximation relies on the condition $h/b \to 0$ as $h,b \to 0$, which makes the variability of the bias correction estimate disappear. That is, provided $h/b \to 0$, the estimator $\widehat{\Theta}_{g,q}^{(p+1)}(t,b)$ in (ref) does not contribute to the limiting variance of $\widehat{\Theta}_{g,p,q}^{(\nu),{\rm bc}}(t,h,b)$. Since $h/b$ is never zero in finite samples, the idea of calonico2014robust is to remove this requirement, and instead let $h/b \to \rho > 0$ as $h,b\to 0$, and thereby get limiting distributions of the estimator sequence where the variability of the bias estimator is accounted for. These ideas are formalised in Corollary (ref) below.
By Lemma (ref) and using the decomposition in (ref) twice, we obtain
From Lemma (ref)(iii) in the appendix, this approximation holds with probability one when $n,h,b$ are so that $J_{n,h}(t) = 1$ and $J_{n,b}(t) = 1$ for all $t \in [0,\tau]$, thus ensuring that $(nh^{2\nu+1})^{1/2}\int_0^t \{J_{n,h}(s) - 1\}\bar{\alpha}_g^{(\nu)}(s,z_0)\,{\rm d} s\ + (nh)^{1/2}h^{p+1}e_{p,\nu}\nu! \kappa_p \int_0^t \{J_{n,b}(s) - 1\}\bar{\alpha}_g^{(p+1)}(s,z_0)\,{\rm d} s = 0$ for all $t$, almost surely.
The ideas underlying designs such as the {{regression discontinuity design}}, the difference in difference design, and the intstrumental variable design are not bound to any particular estimand, model, or type of data. The statistical theory for these designs, however, are much more developed for the type of estimands, models, and data often encountered in economics, than for estimands, models, and data typically encountered in other fields of application. This paper is an attempt at taking one of these designs from the estimation of a conditional average treatment effect ({\sc cate}) based on uncensored data, to the estimation of a {{\sc cate}}-like object, namely the difference of two cumulative hazards, based on right-censored survival data.
A few directions the results of the present paper can be extended in are: First, the estimator developed in this paper can be used as a building block in the estimation of $\theta(t,z_0)$ using kernel smoothing techniques similar to those introduced by ramlau1983smoothing (research in this direction is underway). Second, the results of this paper is limited to the sharp {{regression discontinuity design}}, and ought to be extended to the fuzzy {{regression discontinuity design}}. Third, the estimator of this paper can be used to test whether the parameter of interest in a Cox regression model with unobserved confounders equals zero or not, but it can not provide asymptotically unbiased estimates of this parameter. Whether the {{regression discontinuity design}} can be used to identify the parameter of interest in proportional hazards models under confounding, can be studied.