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.
55,444 characters · 17 sections · 34 citation commands
A nonparametric instrumental approach to endogeneity in competing risks models
{{ Key Words:} Duration Models; Competing risks; Endogeneity; Instrumental variable; Nonseparability; Partial identification.} \\
\setcounter{footnote}{0} \setcounter{equation}{0}
Competing events are events that prevent the statistician from observing the time until the event of interest. A typical example is in biological studies with multiple causes of death. When the competing events are independent from the main outcome, then methods from the survival analysis literature that are designed for random right censoring can be used (Kaplan-Meier, Cox, ...). However, when the causes are dependent, a more careful analysis is necessary and only some features of the model can be identified. As a result competing risks problems can be seen as dependent censoring problems, which may arise if, for instance, some study participants decide not to attend a follow-up interview.
We consider a setting where the researcher is interested in the effect of a treatment on a randomly right-censored duration outcome in the presence of competing risks. The treatment is endogenous, it is not independent of the potential outcomes of the duration. In this case a naive analysis based on data conditional to treatment status could yield biased estimates, if, for example, treated study participants are the ones with the most positive treatment effect. Various approaches have been proposed in the literature to solve the endogeneity issue. Our method is based on an instrumental variable which is sufficiently dependent of the treatment but only affects the outcomes through the treatment. In a typical randomized experiment with binary treatment where there is noncompliance, a natural instrument is the treatment/control group assignment.
This paper focuses on the case where both the treatment and the instrument are categorical. We show that the competing risks model gives rise to a nonparametric nonadditive regression problem. The typical features of interest in competing risks models are functionals of the regression function, which we identify thanks to the instrumental variable. The present framework differs from the usual setting of the nonparametric nonseparable instrumental regression literature (CH, chernozhukov2006instrumental, C, wuthrich2020comparison) because we allow for random right censoring (arising for instance from the end of the observation period) and competing risks (dependent censoring). The regression function cannot be identified for every value of the residual (which corresponds to quantiles of the distribution of the potential outcomes). The cause of identification failure can be censoring, competing risks or both. We single out the points at which the regression function is exactly identified and for the other quantiles we provide partial identification results. We discuss an estimation procedure and assess its performance through simulations. The strategy is applied to the Health Insurance Plan of Greater New York experiment.
When there are no competing risks, this problem has been thoroughly studied (AVdB, AVdB2, BR,chernozhukov2015quantile, frandsen2015treatment, tchetgen2015instrumental, li2015instrumental, chan2016reader, sant2016program, beyhum2021nonparametric) with both parametric and nonparametric methods. In particular, beyhum2021nonparametric also studies a nonparametric instrumental regression problem. The present paper makes two main contributions with respect to the latter works. First, we show how a competing risks model implies a nonparametric regression model. In the absence of competing risks, the duration model can naturally be formulated as a regression model, but in our case the link is far from obvious. Second, we tackle new identification and estimation challenges that arise uniquely because of the presence of competing risks. Note that the mechanism through which competing risks affect identification (non strict monotonicity of the regression function) are different from the effect of random right censoring (non identification of conditional survival functions).
Recent works have tackled the problem in the presence of competing risks. Unlike the present paper, kjaersgaard2016instrumental, zheng2017instrumental, ying2019two, martinussen2020instrumental have semiparametric frameworks. In the case where the treatment and the instrument are binary, richardson2017nonparametric, blanco2019bounds study nonparametric models. Their approaches rely on a monotonicity assumption stating that there are no defiers (see angrist1996identification). This is a restriction on the compliance of study participants to their treatment assignment, which has been criticized (see de2017tolerating). Thanks to this condition, richardson2017nonparametric identify the treatment effects on the cause-specific cumulative incidence functions for the population of compliers, that is agents whose treatment status follows treatment assignment. In blanco2019bounds bounds on the treatment effects for the population of compliers who experience the duration spell of interest are provided, they are valid when the censoring is dependent, which could correspond to a competing risk. In contrast, the present paper makes a rank invariance assumption. It restricts the distribution of potential outcomes, such that ranks (in a sense that is given in Appendix A) between subjects cannot be reversed by the treatment. This condition allows us to identify and estimate the treatment effects on the full population, which is arguably more interesting when one weighs whether or not to make a treatment compulsory. See wuthrich2020comparison for a discussion of the trade-off between the two strategies.
The paper is organized as follows. In Section (ref), we outline the model. Then, identification is discussed in Section (ref). Estimation is studied in Section (ref). In Section (ref), we present simulations. The method is applied to the Health Insurance Plan of Greater New York experiment in Section (ref). Finally, concluding remarks are given in Section (ref).
The objective of this paper is to analyze treatment models when the outcome of the treatment is a duration time that is subject to competing risks and random right censoring. Suppose that there are two latent durations $T_j,j=1,2$. For $j=1,2$, let $U_j$ be a residual which represents the dependence of $T_j$ on unobserved heterogeneity. Without loss of generality, the distribution of $U_j$ is normalized to be unit exponential. The dependence of $U_1$ and $U_2$ is left unrestricted (hence, the independent case is nested by our model). We assume that there exists a continuous and strictly increasing mapping $\psi_j: {\mathbb{R}}_+\mapsto {\mathbb{R}}_+$ such that $$T_j=\psi_j(U_j);\ j=1,2.$$ The objects of interest fail at time $T=\min(T_1,T_2)$ and the cause of this failure is $E=\min\{\operatorname*{arg\,min}_{j=1,2}T_j\}$. Let us now introduce the probability of failure from cause $1$, that is $p={\mathbb{P}}(E=1)$, and the conditional cumulative distribution function of $U_j$ given $E=j$, that is $F_j(u_j)={\mathbb{P}}(U_j\le u_j|E=j)$ for $u_j\in {\mathbb{R}}_+$. We assume that $F_j,j=1,2$ is continuous and strictly increasing on the support of $U_j$ given $E=j$. We define the residual
where $I(\cdot)$ is the indicator function. It can be shown that the distribution of $U$ is uniform on $[0,1]$ (See Lemma A.1 in Appendix A.1). The variable $T$ is generated by $U$ through the following relationship:
Remark that the model does not say anything about what happens at $U=p$. On this event, it is possible that either $E=1$ and $E=2$. Because ${\mathbb{P}}(U=p)=0$, this is innocuous and ensures the symmetry of the model.
Let us now introduce a categorical treatment $Z$ with support $\{z_1,\dots,z_L\}$. We are interested in the potential outcomes (see rubin2005causal) of $(T,E)$ under treatment status $z\in\{z_1,\dots,z_L\}$, which are denoted by $(T(z),E(z))$. We assume that the potential outcomes follow a model of the type (ref), in the sense that for $z\in\{z_1, \dots,z_L\}$ there exists $p_z\in(0,1)$ and continuous and strictly increasing mappings $\varphi_1(z,\cdot):[0,p_z)\mapsto {\mathbb{R}}_+$ and $\varphi_2(z,\cdot) :(p_z,1]\mapsto {\mathbb{R}}_+$ such that
and $T=T(Z),E=E(Z)$. When there are no $Z$, (ref) becomes (ref). Remark that $p_z={\mathbb{P}}(E(z)=1)$ since $U$ is uniform on $[0,1]$, and that the residual $U$ does not vary among treatment statuses, this is the usual rank invariance assumption from the Econometrics literature (see dong2018testing). In Appendix A.2 we argue that this assumption allows for a wide variety of treatment effects.
This paper is concerned with identification and estimation of some features of the model (see Section (ref)) when $Z$ is endogenous and $T=T(Z)$ is randomly right censored. We treat the confounding issue thanks to an instrumental variable (henceforth, IV). Formally, $Z$ and $U$ are dependent but we possess a categorical instrumental variable $W$ with support $\{w_1,\dots,w_K\}$ such that $U$ and $W$ are independent. There also exists a censoring variable $C$ with support in ${\mathbb{R}}_+$ such that we observe $Y=\min(T,C)$ and the censoring indicator $\delta=I(T\le C)$.
Our model can be regarded as a quantile regression model. Indeed, if we introduce the random variable $$T^j(z)=T(z)I(E(z)=j)+\infty I(E(z)\ne j)$$ for $z\in\{z_1,\dots,z_L\}$, then we have $T^j(z)=\varphi^j(z,U)=\varphi^j_z(U),$ where $\varphi^j_z(\cdot):[0,1]\mapsto [0,\infty]$ is such that, for all $u\in[0,1], u\ne p_z$, $$\varphi^1_z(u)=\left\{
\right.;\ \varphi^2_z(u)=\left\{
\right.$$ The quantity $\varphi_z^1(u)$ is the $u$-quantile of the distribution of $T^1(z)$ and, hence, the model $T^1=\varphi^1(Z,U)$ can be seen as a IV model of quantile treatment effects as in \citet{CH}. There are however two distinguishing differences with the usual model of \citet{CH}. First, the variable $T$ is randomly right censored. Second, the mapping $\varphi^1_z(\cdot)$ is not strictly increasing on all of its support. These differences pose additional identification and estimation issues which we tackle in this paper.
Throughout the paper, we suppose that the distribution of $U$ given $\{Z=z,W=w\}$ is continuous, with strictly positive density on $[0,u_{z,w})$, where $$u_{z,w}=\sup\{u\in{\mathbb{R}}_+:\ {\mathbb{P}}(U\le u|Z=z,W=w)<1\}$$ is the upper bound of the support of the distribution of $U$ given $\{Z=z,W=w\}$. Remark that this allows the support of the latter distribution to differ from $[0,1]$. We conclude the presentation of the model $T^1=\varphi^1(Z,U)$ by summarizing the underlying assumptions.
Note that no assumptions on $\varphi^2$ are needed to identifiy $\varphi^1$.
Let us now discuss the features that we seek to recover. It is well-known that it is not possible to identify the distribution of the latent durations in a competing risks model (see tsiatis1975nonidentifiability). However, some characteristics which we introduce below are identified when $Z$ is exogenous (see geskus2020competing). For $z\in\{z_1,\dots,z_L\}$, we introduce the left continuous inverse of the function $\varphi^j_z$, that is $(\varphi^j_z)^{-1}:\ t\in[0,\infty]\mapsto \inf\{u\in [0,1]:\varphi^j_z(u)\ge t\}$. Let us define the cause-specific cumulative incidence function at time $t\in{\mathbb{R}}_+$ for cause $j$ under treatment status $z$,
This equation implies that for $u\in[0,1]$, $\varphi^j_z(u)$ is the $u$-quantile of the subdistribution function $F^j$. Remark that $F_{z}^j(t)$ is a structural function and therefore it is different from the conditional probability ${\mathbb{P}}(T^j\le t|Z=z)={\mathbb{P}}(T\le t, E=j|Z=z)$. We also introduce the subdistribution hazard at time $t\in{\mathbb{R}}_+$ for cause $j$ under treatment status $z$, that is
It is the hazard rate of the subdistribution function $F_{z}^j(t)$, i.e. $h_{z}^j(t)=f_{z}^j(t)/(1-F_{z}^j(t))$, where $f_{z}^j(t)$ is the derivative of $F_{z}^j(t)$ at $t$. Let also the cause-specific hazard at time $t$ for cause $j$ under treatment status $z$ be
We are interested in the effect of the treatment status $z$ on these three quantities, that is the treatment effects. Consider, for instance, the $u$-quantile treatment effect on the subdistribution of the duration until failure from cause $1$, that is $\varphi_1^1(u)-\varphi_0^1(u)$. As is clear from equations (ref), (ref), and (ref), these various treatment effects are identified from knowledge of $\varphi^j, j=1,2$. Therefore, without loss of generality, we focus on identification and estimation of $\varphi^1$ in the remainder of the paper.
Remark that the fact that we assumed that there are only two causes is not restrictive. Indeed, if there are $J\ge 3$ event types, the treatment effects for cause $j^*\in\{1,\dots,J\}$ can be identified from $\varphi^1$ if we label an exit from cause $j^*$ as $j=1$ and an exit from any other cause as $j=2$. Applying this modeling for each cause, that is $J$ times, allows to recover the treatment effects for all causes.
In this section, we discuss identification of $\varphi^1$ when the distribution of the observables $(Y,\delta E, Z,W)$ is known. As usual in nonparametric instrumental regression, the model generates a nonlinear system of equations. Let $p^1=\min_{\ell=1}^Lp_{z_\ell}$. For $u\in[0,p^1)$, $(\varphi^1_{z_\ell}(u))_{\ell=1}^L$ is a solution to the following system of equations in $\theta\in {\mathbb{R}}^L$:
where $S^1(t,z|w) = {\mathbb{P}}(T^1 \ge t,Z=z|W=w)$. This property holds because
where (ref) holds because $\varphi_z^1(\cdot)$ is strictly increasing on $[0,p^1)$. Usually, one would simply make assumptions ensuring that this system has a unique solution. However, the particular context of this paper brings two additional difficulties. First, $S^1(\cdot,z|w)$ may not be identified on all of its support because of censoring. Second, (ref) is only valid for $u\in[0,p^1)$ because $\varphi^1_z(\cdot),z=z_1,\dots,z_L$ are not strictly increasing on all $[0,1]$, which is a direct consequence of the presence of competing risks. These two difficulties restrict the set of values of $u$ for which it is possible to identify $(\varphi_{z_\ell}(u))_{\ell=1}^L$ thanks to (ref). In the next two subsections, we characterize the latter set.
We work throughout the paper with the usual independent censoring assumption:
Let us also define the upper bounds of the support of $C$ conditional on $\{Z=z,W=z\}$ and $\{Z=z\}$, that is
We make the assumption that the upper bound of the censoring time only depends on the treatment status:
This hypothesis simplifies the analysis but could be relaxed. It is likely to hold when the follow-up scheme of the study only depends on the treatment status. Moreover, in this paper we assume that $T$ depends on $W$ only through $Z$, therefore it makes sense to make the same assumption about $C$. Let $$S^1(t|z,w)={\mathbb{P}}(T^1\ge t|Z=z, W=w)={\mathbb{P}}(T\ge t, E=1|Z=z, W=w)$$ be the cause-1-specific survival function conditional to $\{Z=z,W=z\}$. We have $S^1(t|z,w)=1-\int_0^t S(s|z,w)\lambda^1(s|z,w) ds$, where
are the survival function of $T$ conditional on $Z=z, W=w$ and the cause-1-specific hazard rate given $Z=z, W=w$, respectively. The function $S(\cdot|z,w)$ is identified on $[0,c_z]$ by standard arguments from the survival analysis literature. Moreover, remark that, for all $t$ in $[0,c_{z})$, we have
Hence, $\lambda^1(s|z,w)$ is identified on $[0,c_z)$ and so are $S^1(\cdot|z,w)$ and $S^1(t,z|w)=S^1(\cdot|z,w){\mathbb{P}}(Z=z|W=w)$. Next, we have two cases. In order to distinguish them, we introduce the (finite or infinite) upper bound of the support of the distribution of $T$ given $\{E=1, Z=z,W=w\}$: $$t^1_{z,w}=\sup\{t\in{\mathbb{R}}_+:\ {\mathbb{P}}(T\le t|E=1,Z=z,W=w)<1\} .$$ If $c_z\le t^1_{z,w}$, then $S^1(\cdot,z|w)$ is not identified outside $[0,c_z]$. If instead $c_z>t^1_{z,w}$, then the event $\{Y> t^1_{z,w}, \delta E\le 1, Z=z, W=w\}$ would have measure $0$ which would imply ${\mathbb{P}}\left(T> t^1_{z,w}, E=1, Z=z, W=w\right)=0$ leading to $S^1(t^1_{z,w},z|w)=S^1(\infty, z|w)$ and, therefore, $S^1(\cdot,z|w)$ is identified on $[0,\infty]$. Let $$u_C=\sup\{u\in [0,1]:\ \varphi_{z_\ell}^1(u)< c_{z_\ell}\text{ or $\varphi_{z_\ell}^1(u)=\infty$ for all $\ell=1,\dots,L$}\}.$$ The above analysis implies that for $u> u_C$, we do not know the mapping on the left hand side of (ref) at the point $\theta= (\varphi_{z_\ell}(u))_{\ell=1}^L$. As a result, (ref) cannot point identify $ (\varphi_{z_\ell}(u))_{\ell=1}^L$ for $u> u_C$.
The set for which (ref) has the potential to deliver point identification can be further reduced. Let us introduce $$t^1_z=\sup\{t\in{\mathbb{R}}_+:\ {\mathbb{P}}(T\le t|E=1,Z=z)<1\}.$$ The quantity $t^1_z$ is the maximum of $t^1_{z,w}$ over $w\in\{w_1,\dots,w_K\}$. The mapping $S^1(\cdot, z|w)$ is flat on $[t^1_{z,w},\infty]$ by definition of $t^1_{z,w}$. Hence, for all $w\in\{w_1,\dots,w_k\}$, $S^1(\cdot, z|w)$ is constant on $[t^1_{z},\infty]$. Therefore, if $\varphi_z^1(u)$ belongs to $[t^1_{z},\infty]$, then $(\varphi^1_z(u))_{\ell=1}^L$ cannot be identified by the system (ref). Because it is sufficient that this happens for one value of $z$ for identification to break down, it is not possible to identify $(\varphi^1_z(u))_{\ell=1}^L$ with (ref) for all $u\ge u_E$, where $u_E=\min_{\ell=1}^L(\varphi_{z_\ell}^1)^{-1}(t^1_{z_\ell})$. The definition of $t^1_z$ ensures that $\varphi^1_z$ is invertible on $[0,t^1_z)$ and, hence, that $u_E$ is properly defined.
It turns out that $u_E\le p^1$. Indeed, the support of the distribution of $T$ given $\{E=1,Z=z\}$ is equal to the support of the distribution of $\varphi_z^1(U)$ given $\{E=1,Z=z\}$. Since on the event $\{E=1,Z=z\}$ we have $U\le p_z$, we obtain that the support of the distribution of $T$ given $\{E=1,Z=z\}$ is included in the image of $\varphi_z^1(\cdot)$ on $[0,p_z]$. Hence, we have $t^1_z\le \varphi^1_z(p_z-)$, where for a mapping $f:{\mathbb{R}}_+\to {\mathbb{R}}_+$ and $t\in{\mathbb{R}}_+$, $f(t-)$ is the left limit of $f$ at $t$. Therefore $(\varphi^1_z)^{-1}(t^1_z)\le p_z$ because $\varphi^1_z$ is strictly increasing on $[0,p_z]$. This leads to $u_E \le p^1$ as claimed.
Notice that if $t^1_z$ were to be equal to $\varphi^1_z(p_z-)$, we would have that $(\varphi_z^1)^{-1}(t^1_z)=p_z$ and, hence, $u_E=p^1$. The latter happens when the upper bound of the support of $U$ given $Z=z$ is $p_z$. This shows that the fact that the system of equations cannot identify $(\varphi^1_{z_\ell}(u))_{\ell=1}^L$ for $u\in(u_E,p^1)$ is due to the dependence of $U$ and $Z$, that is the selection into treatment.
At the moment, we have that the system (ref) can only identify $(\varphi^1_{z_\ell}(u))_{\ell=1}^L$ for $u\in [0,u_Y)$ with $u_Y= u_E\wedge u_C$ (or $u\in[0,u_Y]$ if $u_C<u_E$). Remark that even if the system (ref) may not have a unique solution at $u_Y$, $(\varphi^1_{z_\ell}(u_Y))_{\ell=1}^L$ can always be identified by continuity as long as identification holds on $[0,u_Y)$. By definition $(\varphi^1_{z_\ell}(u))_{\ell=1}^L\in \prod_{\ell^=1}^L[0,y^1_{z_\ell})$ for these values of $u$, where $y^1_z=t^1_{z_\ell}\wedge c_{z_\ell}$.
Now, we make a strong conditional completeness assumption similar to that in Appendix C of CH which ensures that (ref) has a single solution $(\varphi^1_{z_\ell}(u))_{\ell=1}^L\in \prod_{\ell^=1}^L[0,y^1_{z_\ell})$. We follow the presentation of Appendix A in FFV of this condition. Assume that the density of $(U,Z)$ given $W$ is perturbed in the direction of a function $\Delta \in\mathcal{P}$, where $\mathcal{P}$ is the set of mappings $\Delta:\{z_1,\dots, z_L\}\times{\mathbb{R}}_+\mapsto {\mathbb{R}}_+$ such that $\varphi_z^1(u)+\Delta_z(u)\in[0,y^1_z)$ and $(\varphi_z^1)'(u)+\Delta_z'(u)>0$ for all $u\in[0,u_Y)$. Let $f^1(t,z|w)=-\frac{\partial S^1}{\partial t}(t,z|w)$. For $\mu \in[0,1]$, let us define $$ g_{\mu,\Delta}(u,z|w)=((\varphi_z^1)'(u)+\mu(\Delta_z)'(u))f^1(\varphi_z^1(u)+\mu\Delta_z(u),z|w). $$ The main identification assumption is
It can be interpreted as a strong conditional completeness condition of $Z$ and a uniform random variable $\mu$, given $(W,U)$ for the distribution $g$. In the case where $Z$ and $W$ are binary with support $\{0,1\}$ a simpler identification condition can be given. By Theorem 2 in CH, it suffices to assume
Assuming that $f^1(\cdot,z|w)$ is continuous, this full rank condition implies that the determinant of $G(t)$ is either $>0$ or $<0$ for all $t\in[0,y^1_0)\times [0,y^1_1)$, that is $$\frac{f^1(t_2,1|1)}{f^1(t_1,0|1)}> \frac{f^1(t_2,1|0)}{f^1(t_1,0|0)}$$ or the same inequality with $<$ instead of $>$. Following the semantic of CH, this is a monotone likelihood ratio condition: the instrument increases (or decreases) the probability of being treated ($Z=1$) for all levels of outcomes $t\in[0,y^1_0)\times [0,y^1_1)$. In the case of one-sided noncompliance, where ${\mathbb{P}}(Z=1|W=0)=0$, the condition is trivially satisfied as long as ${\mathbb{P}}(Z=1|W=1)>0$.
To conclude on point identification, we need to be sure that $u_Y$ is identified because otherwise we are not able to know for which value of $u$ the quantity $(\varphi^1_{z_\ell}(u))_{\ell=1}^L$ is identified. Remark that $y^1_z=t^1_{z}\wedge c_{z}$ is identified because it is the upper bound of the support of the distribution of $Y$ given $\{E=1, Z=z\}$. Hence, $u_Y=u_E\wedge u_C=\min_{\ell=1}^L(\varphi_{z_\ell}^1)^{-1}(t^1_{z_\ell}\wedge c_{z_\ell}-)$ is identified (Under Assumption (G), one can solve the system (ref) for increasing values of $u$ until the left limit of $\varphi^1_z(u)$ in $u$ becomes $y^1_z$ for one value of $z\in\{z_1,\dots,z_L\}$). The next theorem summarizes the results of Sections (ref) to (ref).
When $u> u_Y$, it is possible to partially identify $(\varphi^1_{z_\ell}(u))_{\ell=1}^L$. Indeed, (ref) becomes $$\sum_{\ell=1}^L {\mathbb{P}}(\varphi_{z_\ell}^1(U)\ge \varphi^1_{z_\ell}(u),Z=z_\ell|W=w_k)\ge \sum_{\ell=1}^L {\mathbb{P}}(U\ge u,Z=z_\ell|W=w_k)$$ because $\{U\ge u\}\subset\{\varphi_{z}^1(U)\ge \varphi^1_{z}(u)\}$ by (weak) monotonicity of $\varphi_z^1$. As a result, we have $\sum_{\ell=1}^LS^1(\varphi_{z_\ell}^1(u),z_\ell|w_k)\ge u\, \text{for } k=1,\dots,K,$ but because $S^1(\cdot,z|w)$ is only identified on $[0,c_z]$ we cannot directly use these equations. As $S^1(t\wedge c_{z,w},z|w)\ge S^1(t,z|w)$ for all $t\in {\mathbb{R}}$, we leverage instead the following proposition which gives an outer set to the identified set.
This outer set is not sharp because it does not take into account the constraints of continuity and monotonicity of $\varphi^1$. Let us now discuss how to compute this outer set as a finite union of product of intervals. We begin with the case where $Z$ is binary with support $\{0,1\}$. The set can have four types of shapes given in Figure (ref):
for some $\bar\theta_1,\bar \theta_2\ge 0$.
Now, we clarify how to obtain Figure (ref). Take $\theta \in[0,\infty]^2$ outside $[0,y^1_0)\times[0,y^1_1)$. Remark that $S^1(t\wedge c_{z},z|w)$ does not depend on $t$ on $[y^1_{z},\infty)$. Indeed, if $y^1_{z}=t^1_{z}$ (i.e. $t^1_{z}\le c_{z}$), then $t\ge t^1_{z}$ and therefore $S^1(t\wedge c_{z,w},z|w)=S^1(t^1_{z},z|w)$ because $t^1_z$ is larger than the upper bound of the distribution of $T$ given $\{E=1,Z=z,W=w\}$. If instead $y^1_{z}=c_{z}$ (i.e. $c_z\le t^1_z$), we have $S^1(t\wedge c_{z},z|w)=S^1(t^1_{z},z|w)$. As a result, $\theta$ is in the outer set if and only if $(\theta_1+ (y^1_0-\theta_1)I(\theta_1>y^1_0), \theta_2+ (y^1_1-\theta_2)I(\theta_2>y^1_1))^\top$ belongs to the outer set. Therefore, it is enough to study the value of $\theta \mapsto \min_{k=1}^K R_{k,u}(\theta)$ on $[0,y^1_0] \times \{y^1_1\}\cup\{y^1_0\}\times [0,y_1^1]$ to draw the outer set. We begin with $[0,y^1_0)\times \{y^1_1\}$. Because $S^1(\cdot ,0|w)$ is continuously decreasing on $[0,y^1_0)$, the set of vectors $(\theta_1,y^1_1)^\top$ in $[0,y^1_0)\times \{y^1_1\}$ such that $\min_{k=1}^K R_{k,u}(\theta)\ge 0$ is either a segment $[0,\bar{\theta}_1]\times\{y^1_1\}$ where $\bar{\theta}_1\in[0,y^1_1]$ or empty. If this set is not empty, as $S^1(\cdot\wedge c_{0} ,0|w)$ is constant on $[y^1_0,\infty]$, $[0,\bar{\theta}_1]\times[y^1_1,\infty]$ is in the outer set. If rather this set is empty, $[0,y^1_1)\times[y^1_1,\infty]$ is not part of the outer set. Then, implement the approach on $ \{y^1_0\}\times [0,y^1_1)$. Finally, when $\min_{k=1}^K R_{k,u}((y_0^1,y^1_1)^\top)\ge 0$, $[y_0^1,\infty]\times [y^1_1,\infty]$ is in the outer set too.
If $L\ge 3$, one should use a recursive procedure. The outer set in dimension $L$ can be computed as the union of the extrapolations of $L$ outer sets in dimension $L-1$. Indeed, one can begin with fixing the first coordinate $\theta_1$ to $y^1_{z_1}$ and compute the outer set for the other $L-1$ coordinates. If this set is not empty, it should be extrapolated by allowing $\theta_1$ to belong to $[y^1_{z_1},\infty]$. In the same manner, one should use this approach for the $L-1$ other coordinates. The outer set is then the union of the $L$ sets obtained by fixing each coordinate.
We consider estimation with an i.i.d.\ sample of size $n$, $\{Y_i, \delta_i E_i, Z_i,W_i\}_{i=1}^n$. The goal is to estimate $ (\varphi^1_{z_\ell}(u))_{\ell=1}^L$ for $u< u_Y$, that is at quantiles where the regression function is point identified.
We assume that we have consistent estimators $\widehat{S}^1$ of $S^1$, $\widehat{u}_Y$ of $u_Y$ and $\widehat{y}^1_z$ of $y^1_z$ for $z\in \{z_1,\dots,z_L\}$. Choices of these estimators are discussed in Sections (ref) and (ref).
We introduce further notations. For $u\in [0,u_Y)$, let $V(u)$ be a positive definite $K\times K$ weighting matrix. For a $K\times K$ matrix $\bar V$ and a vector $v\in \mathbb{R}^K$, we define $\|v\|_{\bar V}^2=\sqrt{v^{\top}\bar Vv}$. The estimator of $ (\varphi^1_{z_\ell}(u))_{\ell=1}^L$ is defined by
and only the values of $u$ such that $u< \widehat{u}_Y$ should be reported. In practice, we cannot compute the estimator at all $u\in[0,1]$. Instead, we choose a grid $0\le u_1<u_2,<\dots<u_M\le 1$, $M$ values at which we want to estimate $\varphi^1$. Take, for instance, $\{1/M,2/M,\dots,1\}$.
When there are no competing risks, the estimator defined in (ref) is the same as the one in beyhum2021nonparametric, except that in the latter paper the term $1-u$ in (ref) is replaced with $e^{-u}$ because a different normalization of the distribution of $U$ is used. When there are competing risks, the main difference between the estimator of the present paper and the one of beyhum2021nonparametric is the fact that here the function $S^1$ is the survival function of a subdistribution while in beyhum2021nonparametric it is the survival function of the duration of interest. In the present paper this function can be estimated by smoothing the Aalen-Johansen estimator (see Section (ref) below), while in beyhum2021nonparametric it is estimated by smoothing the Kaplan-Meier (Beran) estimator. This does not change the theoretical properties of solutions to (ref), because the Kaplan-Meier estimator of the survival function and the Aalen-Johansen estimator of the survival of the subdistribution of dying from cause $1$ have similar asymptotic properties. Therefore, in the present paper, we skip the theoretical analysis of solutions to (ref) and refer the reader to beyhum2021nonparametric for properties, proofs and further discussion. One major remaining difference between the estimation procedure in beyhum2021nonparametric and the one in the present paper is that, here, we provide a theoretical analysis of how to estimate $u_Y$, see Section (ref) below.
In this subsection, we discuss choices of $\widehat{S}^1$ that work under competing risks and random right censoring. Remark that $S^1(t,z\vert w)=S^1(t\vert z,w)p_{z,w}$ where $p_{z,w}={\mathbb{P}}(Z=z\vert W=w )$. We define the following stochastic processes:
We estimate $S^1(\cdot\vert z,w)$ using the Aalen-Johansen estimator of the cause-specific survival function (see aalen1978empirical,geskus2020competing), that is $$ \widehat{S}_{AJ}^1(t\vert z, w)=1-\sum_{s\le t} \left(\prod_{u\le s}\left[1-\frac{dN_{z,w}(u)}{Y_{z,w}(u)}\right]\right) \frac{dN^1_{z,w}(s)}{Y_{z,w}(s)}.$$ To ensure that (ref) has a unique solution, one should smooth $\widehat{S}^1_{AJ}$. Various techniques are available in the literature, including local polynomials and kernel smoothing. For instance, concerning the latter, if we use a kernel $K$ with a bandwidth $\epsilon$, we obtain $$\widetilde{S}^1(t\vert z,w)=\int \widehat{S}^1_{AJ}(t-s\epsilon\vert z, w)K(s)ds.$$ Our final estimator of $S^1$ is $$\widehat{S}^1(t,z\vert w)= \widetilde{S}^1(t\vert z,w)\widehat{p}_{z,w},$$ where $\widehat{p}_{z,w} =Y_{z,w}/Y_w$.
Let us now introduce estimators of $y^1_z,z=z_1,\dots,z_L$ and $u_Y$. These estimators are new and not discussed in beyhum2021nonparametric, which is why we provide theoretical results. Because $y^1_z$ is the upper bound of the support the distibution of $Y$ given $\{E=1,Z=z\}$, a natural estimator is $$\widehat{y}^1_z=\max\limits_{i\in\{1,\dots,n\}:\ Z_i=z,\ \delta_iE_i=1}Y_i.$$ We have the following result.
The proof is given in Appendix B. In turn, the proposed estimator of $u_Y$ is $\widehat{u}_Y=u_{\widehat{m}_Y}$, where
for some small $\Delta_\ell>0,\ell=1,\dots, L$. The rationale of this estimator is as follows. By definition, $u_Y$ is the lowest value of $u$ such that there exists $\ell \in \{1,\dots,L\}$ for which the left limit of $\varphi^1_{z_\ell}(\cdot)$ in this value is equal to $y^1_{z_\ell}$. If $\widehat{y}^1_z$ is consistent, then $\widehat{y}^1_z-\Delta$ should be close to $y^1_z$. As a result, since $\widehat{\varphi}^1$ is consistent, $\widehat{u}_Y$ is close to $u_Y$. The main role of $\{\Delta_\ell\}_{\ell=1}^L$ is to provide a cushion which ensures that the probability of $\{\widehat{u}_Y>u_Y\}$ is small. This event is not desirable because it would lead the researcher to report results for values of $u$ for which $(\varphi^1_{z_\ell}(u))_{\ell=1}^L$ is not identified. The following proposition formalizes these ideas.
The proof is given in Appendix B. An important remark is that the theoretical analysis in beyhum2021nonparametric concerns the ideal estimator
for all $u< u_Y$ where it is assumed that $y^1_z,z=z_1,\dots,z_L$ and $u_Y$ are known. Our estimator $\widehat{\varphi}^1$ and $\widetilde{\varphi}^1$ coincide on the event $$\mathcal{E}=\left\{u \in [0,1] : (\widetilde{\varphi}^1_{z_\ell}(u))_{\ell=1}^L\in \prod_{\ell=1}^L[0,\widehat{y}_{z_\ell}^1), u< \widehat{u}_Y\wedge u_Y\right\}.$$ For $u<u_Y$, when $(\widetilde{\varphi}^1_{z_\ell}(u))_{\ell=1}^L$ is consistent, we have $\lim\limits_{n\to \infty}{\mathbb{P}}(\mathcal{E})=1$ because $\widehat{y}_{z}^1, z=z_1,\dots,z_L$ and $\widehat{u}_Y$ are consistent. Remark that $\widehat{\varphi}^1=\widetilde{\varphi}^1I(E)+\widehat{\varphi}^1 (1-I(E))$ and $(\widehat{\varphi}^1_{z_\ell})^L_{\ell=1}\in\prod_{\ell=1}^L[0,\widehat{y}_{z_\ell}^1)$ is bounded because $\widehat{y}_{z}^1\le y^1_z, z=z_1,\dots,z_L$ almost surely. Hence, by the continuous mapping theorem and Slutsky's theorem, $\widehat{\varphi}^1$ has the same asymptotic properties as $\widetilde{\varphi}^1$ and the results of beyhum2021nonparametric are not altered.
Let us consider the following data generating processes (henceforth, DGPs). The variable $U$ has a uniform distribution on the interval $[0,1]$ and $W$ is a Bernoulli random variable with parameter $2/3$, independent of $U$. We generate $$Z=I\left(4U+\epsilon-1\ge 0\right)W,$$ where $\epsilon\sim\mathcal{N}(0,1)$ is independent of $(U,W)$. Hence, in this simulation experiment there is one-sided noncompliance as in our empirical application (see Section (ref)). We set $E=I(U>p_Z)+1$, where $p_0=1/2$, $p_1=3/4$. The duration is $$ T=\left\{
\right.;\ $$ $$ \varphi_1(z,u)=\left\{
\right.;\ \quad \varphi_2(z,u)=\left\{
\right. $$ This leads to $$\varphi^1_0(u)=\left\{
\right.;\ \quad \varphi^1_1(u)=\left\{
\right.$$ The treatment increases the probability that $E=1$ and reduces the duration until cause $1$ happens for a given value of $u$. Note that the support of $U|Z$ is $[0,1]$ because of the presence of the random noise $\epsilon$. Therefore, we have $t^1_0= 1$, $t_1^1=3/4$ and $u_E=1/2$. We propose two designs for the censoring variable:
In the first experiment, we have $u_C= 1/3<u_E$, hence identification fails for $u>u_Y=1/3$ because of censoring. In the second exercise, it holds that $u_C=1$ and partial identification arises at $u_Y=1/2$ due to competing risks.
With these DGPs, the probability of treatment given $W=1$ is ${\mathbb{P}}(Z=1|W=1)=73\%$. Under design 1, $30\%$ of the observations are censored, while under design 2, it is only $10\%$ (all these quantities are averages over 1,000,000 Monte Carlo replications). The sample size was set at $n=10,000$. We generated $1,000$ replications of the model under both designs. The estimator (ref) was computed on the grid $u_m=0.01m,\ m=1,\dots,100.$ We used a local polynomial of degree $1$ with an Epanechnikov Kernel to smooth the Aalen-Johansen estimator of the cause-specific survival function. The bandwidth was selected according to the usual rule of thumb for normal densities. To estimate $u_Y$, we used the estimator defined in (ref). The quantity $\Delta_\ell$ was chosen as the optimal bandwidth for normal density estimation with Epanechnikov kernel of the sample $\{Y_i:\ i\in\{1,\dots,n\} \text { such that } Z_i=z_\ell,\ \delta_iE_i=1\}$.
Figures (ref) and (ref) display the histograms of the values of $\widehat{u}_Y$ for design 1 (left) and 2 (right). The estimator $\widehat{u}_Y$ almost always yields a quantile lower than but close to $u_Y$, which avoids reporting results at points at which $\varphi^1$ is not identified. In Figures (ref) and (ref), we present the average value of our estimator of the quantile treatment effect $\widehat{\varphi}^1_1-\widehat{\varphi}^1_0$ for points of the grid $\{u_m\}_{m=1}^M$ smaller than $u_Y$. We also show the average of the naive estimator of the quantile treatment effect, which estimates $\varphi^1_z$ by inverting the Aalen-Johansen estimator of the survival function ${\mathbb{P}}(T\ge t, E=1|Z=z)$. The naive estimator directly compares treated observations to untreated ones. It ignores the endogeneity issue and, hence, is biased. The figures also exhibit the average of the bounds of the $95\%$ confidence intervals of the quantile treatment effects computed with 200 bootstrap draws. The true quantile treatment effects are not reported because they are indistinguishable from our estimator. Finally, Figures (ref) and (ref) show that the coverage of these 95% confidence intervals is almost nominal.
We apply our methodology to the Health Insurance Plan of Greater New York experiment. This clinical trial aimed to evaluate the effect of periodic screening examinations (which aim to detect breast cancer) on breast cancer mortality. The study started in 1963 and lasted until 1986. The experiment follows 60,695 women between 40 and 60 years old. About half (30,565) of the participants were randomized into the control group ($W=0$) while the rest (30,100) were assigned to the intervention group ($W=1$). Members of the intervention arm were offered a treatment consisting of an initial breast examination and mammography and three yearly subsequent screens. 9,984 participants randomized into the intervention group refused the treatment, which corresponds to a rate of non-compliance of 33%. The women in the study group that accepted the treatment exhibited largely different observable characteristics from the ones that refused screening (see shapiro1997periodic). This suggests that the treatment is endogenous. Let $Z$ be the variable which equals $1$ for women who were offered and accepted the treatment and $0$ otherwise.
All participants were followed through three mail surveys, respectively, 5, 10 and 15 years after their entry into the study. The outcome duration of interest $T$ is the time from the initial randomization into the study until death. Censoring arises because some subjects are lost to follow-up before they die. Hence, $C$ is the time between registration into the experiment and the last response to a follow-up survey. We consider two competing risks: deaths from breast cancer ($j=1$) and death from any other cause ($j=2$). In the sample, there are $786$ deaths from breast cancer, $13,798$ deaths from other causes and $46,111$ censored observations.
As in the simulations, the Aalen-Johansen estimators of the survival functions are smoothed using a local polynomial of degree $1$. The bandwidths are selected similarly. The confidence intervals were computed using 200 bootstrap replications. For death by breast cancer, we computed the results on a grid of values of $U$ starting at $u_1=2\times 10^{-4}$ with step $2\times 10^{-4}$. The value of the estimator of $u_Y$ was $0.013$, corresponding to the 1.3% quantile. The quantile treatment effects can only be estimated for such low quantiles because the value of $T^1$ is almost always infinity (most women do not die from breast cancer) and censored (most women do not die in the 15 years following screening). The estimates of the quantile treatment effects along with the 95% confidence intervals are reported in Figure (ref). The quantile treatment effects are insignificant except for some low quantiles. This seems in contrast with findings in shapiro1997periodic, who concluded that the treatment reduced the probability of dying from breast cancer before a certain time, but in the latter paper the endogeneity issue is ignored. For cause $2$ (death by other causes than breast cancer), a grid of values of $U$ starting at $u_1=0.02$ with step $0.02$ was chosen. We found $\widehat{u}_Y=0.12$. However, the estimated quantile treatment effects close to $\widehat{u}_Y$ exhibit a very irregular behavior. Hence, we chose to present in Figure (ref) only the results for quantiles between $0$ and $9.8\%$. As expected, the treatment does not have a significant effect on the subdistribution of the time until of death from another cause than breast cancer.
This paper exhibits the link between competing risks models and nonparametric regression models in the presence of endogeneity. Thanks to this relationship, we are able to formulate our problem in terms of a quantile instrumental regression model. We study identification and estimation of the model. Numerical experiments assess the small sample performance of the method. We show how to revisit an empirical application using this paper's approach.
Many possible directions for future research are interesting. A valuable generalization would allow for continuous or even dynamic treatments, instrumental variables and covariates $X$. With continuous variables, this could be done by replacing the system of equations ((ref)) by a (potentially infinite) set of integral equations $$\int S^1(\varphi^1_z(u),z,x|w)dz=1-u, \text{ for all $x,w,u$ in the support of $X$, $W$ and $U$}.$$ This is an ill-posed problem and regularization of the estimator would be required.
It would also be interesting to identify directly the properties of the marginal distributions of the risks. Following this goal, one could assume that the risks are independent conditional on a set of covariables. Under this condition, the competing events can be treated as independent censoring. Another approach could be to extend models from the dependent censoring literature (such as in basu1978identifiability,emoto1990weibull, deresa2019semiparametric,czado2021) to the case where endogeneity is allowed.
Finally, it should be noted that because the duration $T^1$ has a mass at infinity, all the discussion of this paper can easily be adapted to the estimation of the latency in a cure model with endogenous treatment. In fact a competing risks model can be seen as a cure model where the cure status is known for some observations (see betensky2001nonparametric). A promising project could investigate the identification of the cure fraction (the cause-specific probability of failure in the context of competing risks) when the treatment is endogenous.