TY - JOUR
AU - Vansteelandt,, Stijn
AB - Summary Time-to-event analyses are often plagued by both—possibly unmeasured—confounding and competing risks. To deal with the former, the use of instrumental variables (IVs) for effect estimation is rapidly gaining ground. We show how to make use of such variables in competing risk analyses. In particular, we show how to infer the effect of an arbitrary exposure on cause-specific hazard functions under a semi-parametric model that imposes relatively weak restrictions on the observed data distribution. The proposed approach is flexible accommodating exposures and IVs of arbitrary type, and enabling covariate adjustment. It makes use of closed-form estimators that can be recursively calculated, and is shown to perform well in simulation studies. We also demonstrate its use in an application on the effect of mammography screening on the risk of dying from breast cancer. 1. Introduction In most observational studies unobserved confounding cannot be ruled out. This can make the results on exposure effects, as obtained via standard regression methods, questionable. Sometimes, however, it may be possible to estimate an exposure effect without (large sample) bias when an instrumental variable (IV) is available. This is a variable which is (i) associated with the exposure, (ii) has no direct effect on the outcome other than through the exposure, and (iii) whose association with the outcome is not confounded by unmeasured variables (e.g. Hernán and Robins 2006). Condition (i) is empirically verifiable, but conditions (ii) and (iii) are not. IVs estimation of exposure effects is well established for continuous outcomes that obey linear models. One popular technique is two-stage predictor substitution (2SLS) estimation, sometimes also referred to as the 2SPS method (Cai and others, 2011). Here, the exposure variable is regressed on the instrument in the first stage, and then the outcome variable is regressed on the predicted exposure value in the second stage. The regression coefficient of the predicted exposure in the second stage is then interpreted as the exposure effect of interest. Recently there has been a focus on extending these methods to handle also right censored failure time data. Robins and Tsiatis (1991) initiated this work, but although they developed a general estimating equations-based method under structural accelerated failure time models, their proposal suffers from a lack of smoothness of the estimating equations due to the way how censoring is handled (Joffe and others, 2012). Tchetgen Tchetgen and others (2015) developed an easy-to-use two-stage estimation approach under additive hazard models for event times, which works when the exposure obeys a location shift model (see Li and others, 2015 for a related approach). Martinussen and others (2017) generalized these methods by working under a less restrictive semiparametric structural cumulative failure time model, imposing no restrictions on distribution (or type) of instrument or exposure. Their proposal has the further advantage of enabling non-parametric estimation of a possibly time-varying exposure effect. Kjaersgaard and Parner (2016) suggested an alternative approach based on pseudo-observations. Their 2SLS method requires a latent additive model for the target parameter which is not so attractive when focusing on a distribution function. Motivated by an analysis of the HIP-study (Health Insurance Plan of Greater New York, see Joffe 2001), which was designed to assess the potential effect of breast cancer screening, here, we aim at extending the methods of Martinussen and others (2017) to handle competing risk data. The HIP-study comprised approximately 60 000 women, who were randomized into two approximately equally sized groups. About 35% of the women, who were offered screening, refused to participate, resulting in a problem of non-compliance. We planned to correct for this using randomization as an IV. In the first 10 years of follow-up, there were 4221 deaths, but only 340 were deemed due to breast cancer, making competing risks a major issue in these data. Richardson and others (2017) proposed a method that can deal with competing risk data also using an IV approach. Their approach requires the instrument as well as exposure variable to be binary variables. Essentially they generalize the standard IV Wald estimator for survival probabilities to estimate cumulative incidence probabilities. In this way, they estimate the so-called complier average treatment effect. The method we propose puts no restriction on the type of instrument nor on the exposure, and it can also incorporate covariates which is not possible using the Wald type estimator of Richardson and others (2017). Our method is thus much more general. Zheng and others (2017) suggest a two-stage approach where they directly model the subdistribution hazard. This quantity in itself is hard to interpret as it involves a risk set of those alive and those that have died from the competing cause, but it can be seen just as a tool to model the cumulative incidence function, and it is widely used. For further discussion on this, we refer to Andersen and Keiding (2012) and Andersen and others (2012). Zheng and others (2017) require a model for the exposure distribution that involves unmeasured variables. This model is assumed to be linear which is geared toward continuous exposure as it may not be natural to use for a binary exposure. They also need that observed and unobserved covariates are independent. Our method avoid such an assumption and puts no restriction on the exposure distribution. In particular, it can handle binary and continuous exposures. The article is structured as follows. In Section 2, we specify the model and outline the estimation procedure. Section 3 contains large sample results. In Section 4, we study by simulations the practical behavior of the proposed estimator and also analyze the HIP-data. Section 5 contains some closing remarks and technical details are deferred to the Appendix. 2. Model specification and estimation We let |$\tilde T$| denotes the time until one of the two competing events happens and let |$\delta=1,2$| denotes which of the two that takes place. Our aim is to assess the effect of an arbitrary exposure |$X$| on the cause-specific hazard of each of these competing events, by making use of an IV |$G$|⁠. This variable is such that, possibly conditional on measured covariates |$L$|⁠, |$G$| is associated with the exposure |$X$|⁠, but is not associated with the event time |$\tilde T$|⁠, nor the event type |$\delta$|⁠, except because of a possible exposure effect. More formally, let |$(\tilde T^x,\delta^x)$| denote the counterfactual event time and event type that would be observed for given subject if the exposure of that subject were set to |$x$|⁠. We will make the consistency assumption that these coincide with the observed event time |$\tilde T$| and event type |$\delta$| for those subjects who happen to have exposure level |$X=x$|⁠. This notation enables us to be clear about our target of inference, which is the contrast between the counterfactual cause-specific hazard functions $$\begin{equation}\label{estimand}\lambda^j_{T^x}(t|X=x,G,L)-\lambda^j_{T^0}(t|X=x,G,L),\end{equation}$$ (2.1) for |$j=1,2$|⁠, where $$\begin{equation}\label{cause_spec_haz}\lambda_{T^x}^j(t|X,G,L)=\frac{ \frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d} t}P\left(\tilde T^x\le t,\delta^x=j| X,G,L\right)}{P\left(\tilde T^x>t| X,G,L\right)},\quad j=1,2.\end{equation}$$ (2.2) Because |$T^0$| and |$\delta^0$| are unobserved for subjects with non-zero exposure, we will rely on the assumption that |$G$| is an IV for the exposure effect (conditional on |$L$|⁠), in the sense that |$(\tilde{T}_0,\delta^0)$| is conditionally independent of |$G$|⁠, given |$L$| (Hernán and Robins, 2006). This assumption expresses that, if all subjects received zero exposure, both events would have cause-specific hazards (conditional on |$L$|⁠) that would be same at all levels of |$G$|⁠. This would be the case when, as in the causal directed acyclic graph of Figure 1, |$G$| shares no common causes with the event time and type, and does not influence those in the absence of exposure. As in other IVs problems, these IVs assumptions will not generally suffice to identify the contrast (2.1) at all levels of |$X,G,$| and |$L$|⁠. For that reason, as well as for reasons of parsimony, we will assume that the following structural model holds Fig. 1. Open in new tabDownload slide Causal directed acyclic graph. |$G$| is the instrument, |$X$| is the exposure variable, and |$\tilde T$| is the time-to-event and |$\delta$| indicates which of the two competing events that is taking place. The potential unmeasured confounders are denoted by |$U$|⁠, and the observed confounders by |$L$|⁠. Fig. 1. Open in new tabDownload slide Causal directed acyclic graph. |$G$| is the instrument, |$X$| is the exposure variable, and |$\tilde T$| is the time-to-event and |$\delta$| indicates which of the two competing events that is taking place. The potential unmeasured confounders are denoted by |$U$|⁠, and the observed confounders by |$L$|⁠. $$\begin{equation} \label{StrucMod} \lambda^j_{T^x}(t|X=x,G,L)-\lambda^j_{T^0}(t|X=x,G,L)=\beta_j(t)x,\end{equation}$$ (2.3) for all |$t>0$| and for |$j=1,2$|⁠, with |$\beta_j(t)$| an unknown, locally integrable function. Our aim is then to estimate |$B_j(t)=\int_0^t\beta_j(s)\, {\rm d}s$| for all |$t>0$|⁠, |$j=1,2.$| As in Martinussen and others (2017), it can be shown that this model is satisfied when the causal directed acyclic graph of Figure 1 holds and, moreover, $$\lambda^j(t|X,G,L,U)=\lambda_{j}(t)+\beta_{j}(t)X+\psi_j(t,U,L)$$ with the function |$\psi_{j}(.),j=1,2$| left unspecified and the functions |$\lambda_{j}(.)$| and |$\beta_j(t)$| unknown. Model (2.3) excludes interactions of the exposure effect with both G and L. Whereas the model can be extended to allow for interactions between X and L (see Section 5), the absence of interactions between X and G is assumed in order to obtain identification of the exposure effect. Others (e.g. Richardson and others 2017) avoid this no-interaction assumption by restricting themselves to settings with dichotomous instrument and exposure, and focusing on the treatment effect in compliers; here, compliers refer to individuals who would take treatment if assigned to one level of the instrument (say, Level 1), but not the other (say, Level 0). They instead invoke a monotonicity assumption that increasing the level of the instrument cannot decrease the exposure for any individual. Here, we avoid this framework because it does not easily extend to continuous instruments or exposures, and because in our opinion, there is usually greater scientific interest in knowing the exposure effect for the general study population. Throughout, we will allow for the event time |$\tilde T$| to be subject to right-censoring. In that case, we only observe whether or not |$\tilde T$| exceeds a random censoring time |$C$|⁠, i.e. we observe |$D=I(\tilde T\leq C)$|⁠, along with the first time either failure or censoring occurs, i.e. we also observe |$T=\min(\tilde T,C)$| as well as |$\delta$| if |$D=1$|⁠. Let |$\delta=0$| if |$D=0$|⁠. Define also the observed counting processes |$N^j(t)=I(T\leq t,D=1,\delta=j)$|⁠, |$j=1,2$|⁠, and the at risk indicator |$R(t)=I(t\leq T)$|⁠. We assume that the censoring time satisfies the following condition $$\begin{equation}\tag{C}\mbox{$\tilde T\perp\!\!\!\perp C|X,G,L$ and $P(C>t|X,G,L)=P(C>t|L)$}\end{equation}$$ The above condition on the censoring distribution can be relaxed to |$P(C>t|X,G,L,U)=P(C>t|L,U)$| for some variable |$U\perp\!\!\!\perp G|L$|⁠. The following proposition lays the basis of the estimation procedure for |$B_j(t)$|⁠, which we will describe next. Proposition 1 Assume the structural model (2.3) with the assumption that |$G$| is an instrumental variable, conditional on |$L$|⁠, and further that the censoring time satisfies condition (⁠|$\mathcal{C}$|⁠). Then $$\begin{equation}\label{Esteq}E\left[\left\{G-E(G|L)\right\}{\rm e}^{B_1(t)X+B_2(t)X}R(t)\left\{{\rm d}N^j(t)-{\rm d}B_j(t)X\right\}\right]=0,\end{equation}$$ (2.4) for each |$t$|⁠, |$j=1,2$|⁠. Proof. By the independent censoring assumption (⁠|$\mathcal{C}$|⁠) and $$\frac{P\left(\tilde T^0>t|X,G,L\right)}{P\left(\tilde T>t|X,G,L\right)}= {\rm e}^{B_1(t)X+B_2(t)X}$$ it follows, for |$j=1,2$|⁠, that $$\begin{align*}&E\left[\left\{G-E(G|L)\right\} {\rm e}^{B_1(t)X+B_2(t)X}R(t)\left\{{\rm d}N^j(t)-{\rm d}B_j(t)X\right\}\right]\\&=E\left[\left\{G-E(G|L)\right\}{\rm e}^{B_1(t)X+B_2(t)X}I(C>t)I(\tilde T>t)\lambda^j_{\tilde{T}^0}(t|X,G,L){\rm d}t\right]\\&=E\left[P(C>t|L)\left\{G-E(G|L)\right\}P(\tilde T^0>t|X,G,L)\lambda^j_{\tilde{T}^0}(t|X,G,L){\rm d}t\right]\\&=E\left[P(C>t|L)\left\{G-E(G|L)\right\}\frac{{\rm d}}{{\rm d}t}P(\tilde T^0\leq t,\delta^0=j|X,G,L){\rm d}t\right]\\&=\frac{{\rm d}}{{\rm d}t}E\left[P(C>t|L)\left\{G-E(G|L)\right\}P(\tilde T^0\leq t,\delta^0=j|G,L){\rm d}t\right]\\&=0\end{align*}$$ because, for any function |$g_t(L)$|⁠, we have $$E\left[g_t(L)\left\{G-E(G|L)\right\}P(\tilde T^0\leq t,\delta^0=j|G,L){\rm d}t\right]=0 $$ since |$G\perp\!\!\!\perp (\tilde{T}^0,\delta^0)|L$|⁠. This completes the proof. □ The above proposition gives two unbiased estimating functions for each time |$t$|⁠, on the basis of which we can construct a consistent estimator of |$B_j(t)$|⁠. In particular, let |$(T_i,D_i,\delta_i,L_i,G_i,X_i)$|⁠, |$i=1,\ldots, n,$| denote |$n$| independent identically distributed replicates under the structural model (2.3) together with the IVs assumptions. Suppose that the counting processes |$N^j_i(t)=I(T_i\leq t,D_i=1,\delta_i=j)$|⁠, |$i=1,\ldots, n$|⁠, |$j=1,2$|⁠, are observed in the time interval |$[0,\tau ]$|⁠, where |$\tau$| is some finite time point. Solving equation (2.4), with population expectations substituted by sample analogs leads to the recursive estimator |$\hat B_j(t)$| defined by $$\begin{equation}\label{Est}\hat B_j(t,\hat \theta)=\int_0^t \frac{\sum_iG^c_i(\hat\theta){\rm e}^{\{\hat B_1(s-)+\hat B_2(s-)\}X_i}{\rm d}N^j_i(s)}{\sum_iG^c_i(\hat\theta)R_i(s){\rm e}^{ \{\hat B_1(s-)+\hat B_2(s-)\}X_i}X_i},\end{equation}$$ (2.5) where |$G^c_i(\theta)=G_i-E(G_i|L_i;\theta)$|⁠, with |$E(G_i|L_i;\theta)$| a parametric model for |$E(G_i|L_i)$| and |$\hat\theta$| a consistent estimator of |$\theta$| (e.g. a maximum likelihood estimator). Note that |$\hat B_j(t,\hat \theta)$| is a step function that is well defined by setting |$\hat B_j(0,\hat \theta)=0$|⁠, |$j=1,2$|⁠. Furthermore, note that—unlike many other IV estimators—the estimator (2.5) can be evaluated for discrete as well as continuous exposures and instruments. It does not require distributional assumptions for the exposure and does not make assumptions as to how measured covariates relate to the event time. In practice, it may be more appealing to report some function of the cumulative incidence functions. In the following, we focus on the situation where the exposure is binary. One may for instance focus on the relative risk function $$\begin{align}\label{RR}\mbox{RR}(t)\equiv \frac{P(T^0\le t,\delta^0=1| X=1,G,L)}{P(T^1\le t,\delta^1=1| X=1,G,L)}\end{align}$$ (2.6) or the absolute difference function $$\begin{align}\label{AD}\mbox{AD}(t)\equiv P(T^0\le t,\delta^0=1| X=1,G,L)-P(T^1\le t,\delta^1=1| X=1,G,L).\end{align}$$ (2.7) We have $$P(T^1\le t,\delta^1=1| X=1,G,L)=P(T\le t,\delta=1| X=1,G,L),$$ which is expressed via the observed data. Also $$\begin{align*}&P(T^0\le t,\delta^0=1| X=1,G,L)\\=&\int_0^t{\rm e}^{-\Lambda^1_{T^0}(s|X=1,G,L)-\Lambda^2_{T^0}(s|X=1,G,L)}{\rm d}\Lambda^1_{T^0}(s|X=1,G,L)ds\\=&\int_0^t{\rm e}^{-\Lambda^1_{T}(s|X=1,G,L)-\Lambda^2_{T}(s|X=1,G,L)+B_1(s)+B_2(s)}\{{\rm d}\Lambda^1_T(s|X=1,G,L)-{\rm d}B_1(s)-{\rm d}B_2(s)\}\end{align*}$$ By formulating a model for the observed data via the cause-specific hazard functions |$\lambda^j_{T}(s|X,G,L)$|⁠, |$j=1,2$|⁠, and using the estimator (2.5) we may then estimate |$\mbox{RR}(t)$| and |$\mbox{AD}(t)$|⁠. In general, however, specific models for the observed data may not be congenial with the structural model meaning that there may no data generating mechanism such that both models hold, see Robins and Rotnitzky (2004). In the situation without the covariate |$L$|⁠, and where exposure and instrument are both binaries we may estimate the cumulative cause-specific hazard functions for the observed data using non-parametric estimators and hence avoiding the issue of non-congeniality. This is the situation for the HIP-data that we study in Section 4.2. In practice, one can calculate |$\hat \Lambda^j_{T}(s|X=1,G)$|⁠, |$j=1,2$|⁠, using the |${\tt aalen}$|-function of the R-package |${\tt timereg}$|⁠. 3. Large sample properties The following proposition, whose proof is sketched in the Appendix, shows that the estimators |$\hat B_j(t)$|⁠, |$j=1,2$|⁠, are uniformly consistent. It moreover gives the asymptotic distribution of the two estimators. Proposition 2 Under model (2.3) with the assumption (C) and the assumption that |$G$| is an IV, conditional on |$L$|⁠, and given the technical conditions listed in the Appendix, the IV estimators |$\hat B_{j}(t)$|⁠, |$j=1,2$|⁠, are uniformly consistent. Furthermore, |$ W^j_n(t)=n^{1/2}\{\hat B_{j}(t,\hat \theta)-B_j(t)\} $| converges in distribution to a zero-mean Gaussian process with variance |$\Sigma_j(t)$|⁠. A uniformly consistent estimator |$\hat\Sigma_j(t)$| of |$\Sigma_j(t)$| is given below. Let |$\epsilon_i^{B_j}(t,\theta),i=1,...,n$| be the iid zero-mean processes given by expression (A.1) in the Appendix. From the proof in the Appendix, it then follows that |$W^j_n(t)$| is asymptotically equivalent to |$n^{-1/2}\sum_{i=1}^n\epsilon_i^{B_j}(t,\theta)$|⁠. The variance |$\Sigma_j(t)$| of the limit distribution can thus be consistently estimated by $$\begin{equation}\label{Variance}\hat \Sigma_j(t)=n^{-1}\sum_{i=1}^n\{\hat \epsilon_i^{B_j}(t,\hat\theta)\}^2,\end{equation}$$ (3.1) where |$\hat \epsilon_i^{B_j}(t,\hat \theta)$| is obtained from |$\epsilon_i^{B_j}(t,\theta) $| by replacing unknown quantities with their empirical counterparts. These results can be used to construct a pointwise confidence band. In the case without additional covariates |$L$|⁠, and with |$X$| and |$G$| being binaries, |$n^{1/2}\{\hat \Lambda^j_{T}(s|X=1,G)-\Lambda^j_{T}(s|X=1,G)\}$| also allows an iid-representation (Martinussen and Scheike, 2006), and it can furthermore be extracted from the |${\tt aalen}$|-function in R. Now using the functional delta theorem we may then obtain iid-representations of |$n^{1/2}\{\widehat{\mbox{RR}}(t)-\mbox{RR}(t)\}$|⁠, and of |$n^{1/2}\{\widehat{\mbox{AD}}(t)-\mbox{AD}(t)\}$|⁠, that may then be used to also calculate confidence bands. We illustrate this using the HIP-data (see Section 4.2). 4. Numerical results In this section, we investigate the practical behavior of the proposed estimator. In Section 4.1, we study the small sample performance using simulations, and in Section 4.2, we give a worked application using the HIP-data on the potential effect of breast cancer screening on death due to breast cancer. 4.1. Simulation study To investigate the properties of our proposed methods with practical sample sizes, we conducted a simulation experiment, whereby we generated data under the causal Directed Acyclic Graph of Figure 1. We took |$G$| to be binary with |$P(G=1)=0.5$|⁠, and generated |$X$| and |$U$|⁠, given |$G$|⁠, from a normal distribution with |$E(X|G=g)=0.5+\gamma_Gg$|⁠, |$E(U|G=g)=1.5$| and with variance-covariance matrix so that |$Var(X|G)=Var(U|G)=0.25$|⁠, and |$Cov(X,U|G)=-1/6$|⁠. The parameter |$\gamma_G$| determines the size of the correlation between exposure and the IV. Specifically we looked at correlation |$\rho$| equal to 0.3 and 0.5. The two cause-specific hazards were given as $$\begin{align*}\lambda^1(t|X,G,U)&=0.1+0.1U\\\lambda^2(t|X,G,U)&=0.1+0.2X+0.1U\end{align*}$$ so |$\tilde T$| is generated according to the hazard $$\lambda_{\tilde T}(t|X,G,U)=0.2+0.2X+0.2U$$ and failure of Type 1 happens with probability $$\frac{\lambda^1(t|X,G,U)}{\lambda^1(t|X,G,U)+\lambda^2(t|X,G,U)},$$ and likewise with failure of Type 2. It is easily seen that model (2) holds under this model. Twenty percent were potentially censored according to a uniform distribution on (0,3.5), and the rest were censored at |$t=3.5$|⁠, corresponding to the study being closed at this time point, leading to an overall censoring rate of around 17. For this scenario, we considered sample size 1600 when |$\rho=0.3$|⁠, and sample size 1000 when |$\rho=0.5$|⁠. Simulation results are given in Table 1 based on 2000 runs for each configuration, where we report average biases at time points |$t=0.5, 1.5, 2.5$| for |$\hat B_{j}(t)$|⁠. We also report the empirical standard errors as well as estimated standard errors based on formula (3.1) along with coverage probability of 95% pointwise confidence intervals CP(⁠|$\hat B_{j}(t)$|⁠). Biases from the naive Aalen estimator, denoted as |$\tilde B_{j}(t)$| in the table, running Aalen’s additive hazards model (Martinussen and Scheike 2006, Ch. 5) for the two cause-specific hazards using |$X$| and |$G$| as covariates are also given. Table 1. Continuous exposure and binary instrument . . |$\rho=0.3$| . . |$\rho=0.5$| . . n . |$t=0.5$| . |$t=1.5$| . |$t=2.5$| . n . |$t=0.5$| . |$t=1.5$| . |$t=2.5$| . Bias |$\hat B_{1}(t)$| 1600 0.001 0.007 –0.007 1000 –0.000 –0.001 –0.004 SD (⁠|$\hat B_{1}(t)$|⁠) 0.061 0.132 0.233 0.042 0.089 0.149 See (⁠|$\hat B_{1}(t)$|⁠) 0.061 0.130 0.231 0.042 0.089 0.149 95% CP(⁠|$\hat B_{1}(t)$|⁠) 95.4 95.3 96.6 95.0 95.5 95.7 Bias |$\hat B_{2}(t)$| 0.002 –0.001 –0.010 0.000 –0.000 –0.005 SD (⁠|$\hat B_{2}(t)$|⁠) 0.077 0.165 0.300 0.053 0.113 0.198 See (⁠|$\hat B_{2}(t)$|⁠) 0.076 0.165 0.296 0.054 0.117 0.199 95% CP(⁠|$\hat B_{2}(t)$|⁠) 94.8 95.6 96.0 95.2 96.1 96.5 Bias |$\tilde B_{1}(t)$| –0.033 –0.099 –0.167 –0.034 –0.101 –0.168 Bias |$\tilde B_{2}(t)$| –0.034 –0.102 –0.170 –0.035 –0.102 –0.168 . . |$\rho=0.3$| . . |$\rho=0.5$| . . n . |$t=0.5$| . |$t=1.5$| . |$t=2.5$| . n . |$t=0.5$| . |$t=1.5$| . |$t=2.5$| . Bias |$\hat B_{1}(t)$| 1600 0.001 0.007 –0.007 1000 –0.000 –0.001 –0.004 SD (⁠|$\hat B_{1}(t)$|⁠) 0.061 0.132 0.233 0.042 0.089 0.149 See (⁠|$\hat B_{1}(t)$|⁠) 0.061 0.130 0.231 0.042 0.089 0.149 95% CP(⁠|$\hat B_{1}(t)$|⁠) 95.4 95.3 96.6 95.0 95.5 95.7 Bias |$\hat B_{2}(t)$| 0.002 –0.001 –0.010 0.000 –0.000 –0.005 SD (⁠|$\hat B_{2}(t)$|⁠) 0.077 0.165 0.300 0.053 0.113 0.198 See (⁠|$\hat B_{2}(t)$|⁠) 0.076 0.165 0.296 0.054 0.117 0.199 95% CP(⁠|$\hat B_{2}(t)$|⁠) 94.8 95.6 96.0 95.2 96.1 96.5 Bias |$\tilde B_{1}(t)$| –0.033 –0.099 –0.167 –0.034 –0.101 –0.168 Bias |$\tilde B_{2}(t)$| –0.034 –0.102 –0.170 –0.035 –0.102 –0.168 The two causes are given by |$j=1,2$|⁠. Bias of |$\hat B_{j}(t)$|⁠, average estimated standard error, sd(⁠|$\hat B_{j}(t)$|⁠), empirical standard error, see(⁠|$\hat B_{j}(t)$|⁠), and coverage probability of 95% pointwise confidence intervals CP(⁠|$\hat B_{j}(t)$|⁠) based on the instrumental variables estimator, in function of sample size |$n$| and at different strengths |$\rho$| (correlation) of the instrumental variable. Bias of |$\tilde B_{j}(t)$| is the bias of the naive Aalen estimator. Open in new tab Table 1. Continuous exposure and binary instrument . . |$\rho=0.3$| . . |$\rho=0.5$| . . n . |$t=0.5$| . |$t=1.5$| . |$t=2.5$| . n . |$t=0.5$| . |$t=1.5$| . |$t=2.5$| . Bias |$\hat B_{1}(t)$| 1600 0.001 0.007 –0.007 1000 –0.000 –0.001 –0.004 SD (⁠|$\hat B_{1}(t)$|⁠) 0.061 0.132 0.233 0.042 0.089 0.149 See (⁠|$\hat B_{1}(t)$|⁠) 0.061 0.130 0.231 0.042 0.089 0.149 95% CP(⁠|$\hat B_{1}(t)$|⁠) 95.4 95.3 96.6 95.0 95.5 95.7 Bias |$\hat B_{2}(t)$| 0.002 –0.001 –0.010 0.000 –0.000 –0.005 SD (⁠|$\hat B_{2}(t)$|⁠) 0.077 0.165 0.300 0.053 0.113 0.198 See (⁠|$\hat B_{2}(t)$|⁠) 0.076 0.165 0.296 0.054 0.117 0.199 95% CP(⁠|$\hat B_{2}(t)$|⁠) 94.8 95.6 96.0 95.2 96.1 96.5 Bias |$\tilde B_{1}(t)$| –0.033 –0.099 –0.167 –0.034 –0.101 –0.168 Bias |$\tilde B_{2}(t)$| –0.034 –0.102 –0.170 –0.035 –0.102 –0.168 . . |$\rho=0.3$| . . |$\rho=0.5$| . . n . |$t=0.5$| . |$t=1.5$| . |$t=2.5$| . n . |$t=0.5$| . |$t=1.5$| . |$t=2.5$| . Bias |$\hat B_{1}(t)$| 1600 0.001 0.007 –0.007 1000 –0.000 –0.001 –0.004 SD (⁠|$\hat B_{1}(t)$|⁠) 0.061 0.132 0.233 0.042 0.089 0.149 See (⁠|$\hat B_{1}(t)$|⁠) 0.061 0.130 0.231 0.042 0.089 0.149 95% CP(⁠|$\hat B_{1}(t)$|⁠) 95.4 95.3 96.6 95.0 95.5 95.7 Bias |$\hat B_{2}(t)$| 0.002 –0.001 –0.010 0.000 –0.000 –0.005 SD (⁠|$\hat B_{2}(t)$|⁠) 0.077 0.165 0.300 0.053 0.113 0.198 See (⁠|$\hat B_{2}(t)$|⁠) 0.076 0.165 0.296 0.054 0.117 0.199 95% CP(⁠|$\hat B_{2}(t)$|⁠) 94.8 95.6 96.0 95.2 96.1 96.5 Bias |$\tilde B_{1}(t)$| –0.033 –0.099 –0.167 –0.034 –0.101 –0.168 Bias |$\tilde B_{2}(t)$| –0.034 –0.102 –0.170 –0.035 –0.102 –0.168 The two causes are given by |$j=1,2$|⁠. Bias of |$\hat B_{j}(t)$|⁠, average estimated standard error, sd(⁠|$\hat B_{j}(t)$|⁠), empirical standard error, see(⁠|$\hat B_{j}(t)$|⁠), and coverage probability of 95% pointwise confidence intervals CP(⁠|$\hat B_{j}(t)$|⁠) based on the instrumental variables estimator, in function of sample size |$n$| and at different strengths |$\rho$| (correlation) of the instrumental variable. Bias of |$\tilde B_{j}(t)$| is the bias of the naive Aalen estimator. Open in new tab It is seen from Table 1 that the suggested estimators are unbiased and also that the estimated standard errors estimate well the variability resulting in satisfactory coverage probabilities. As expected, the naive estimators |$\tilde B_{j}(t)$|⁠, |$j=1,2$|⁠, are biased. An additional simulation study is reported in the supplementary material available at Biostatistics online. In the supplementary material available at Biostatistics online, we took both the exposure variable and the instrument to be continuous variables. To our knowledge, there are no other available methods to handle such a situation. 4.2. Application to the HIP trial on effectiveness of screening on breast cancer mortality The HIP of Greater New York was a randomized trial of breast cancer screening that began in 1963. About 60 000 women aged 40–60 were randomized into two approximately equally sized groups. Study women were offered the screening examinations consisting of clinical examination, and a mammography. Further three annual examinations were offered in this group. Control women continued to receive their usual medical care. There were 30 565 women in the control group and 30 130 in the screening group of which 9984 (35%) refused to participate (non-compliers). There were large differences between the study women who participated and those who refused (Shapiro, 1977), and therefore, the results from the “as treated” analysis may be doubtful due to potential unobserved confounding. In this setting, the exposure is the treatment and the instrument is the original randomization, which is unrelated to any unobserved confounders. Further, there is no reason to believe that the instrument should have any effect in itself on the risk of dying of any cause, only (potentially) via the treatment. Exposure and instrument were correlated, estimated correlation coefficient is 0.71, being highly significant. We applied the estimator given by (2.5) to these data focusing on the first 10 years of follow-up. This estimator is shown in Figure 2 along with 95% pointwise confidence intervals. The left panel gives results for breast cancer and the right panel for other causes. The intention to treat estimator is also shown (dotted curves). It is seen from Figure 2 that breast cancer screening appears to lower the risk of dying from breast cancer while there is no evidence of an effect of screening on the risk of dying from other causes. The impact of the screening on the risk of dying from breast cancer seems to be slightly more pronounced than what is indicated by the intention to treat estimator. As the specific value of |$B_1(t)$| may be hard to interpret, we suggest also to report the estimated relative risk function |$\widehat{\mbox{RR}}(t)$| and/or the absolute risk function |$\widehat{\mbox{AD}}(t)$|⁠. Figure 3 displays |$\widehat{\mbox{RR}}(t)$| along with 95% pointwise confidence bands. It is seen that for women who received screening on the active arm the risk of dying from breast cancer within 5 (10) years would have been approximately twice (1.5 times) as large had they not received screening. Fig. 2. Open in new tabDownload slide HIP-study. Estimated causal effect of screening, |$\hat B(t)$| along with 95% pointwise confidence bands and the intention to treat estimate (broken curve). Curves are given for the two competing events. Fig. 2. Open in new tabDownload slide HIP-study. Estimated causal effect of screening, |$\hat B(t)$| along with 95% pointwise confidence bands and the intention to treat estimate (broken curve). Curves are given for the two competing events. Fig. 3. Open in new tabDownload slide HIP-study. Estimated relative risk function (breast cancer death), see display (2.6), along with 95% pointwise confidence bands (dotted curves). Fig. 3. Open in new tabDownload slide HIP-study. Estimated relative risk function (breast cancer death), see display (2.6), along with 95% pointwise confidence bands (dotted curves). Figure 4 displays |$\widehat{\mbox{AD}}(t)$| along with 95% pointwise confidence bands. Again it is seen that for women who received screening on the active arm the risk of dying from breast cancer within 5 (10) years would have a larger risk had they not received screening. We also see, however, that the risk of dying of due to breast cancer is low for the women in this study giving a small (but significant) difference of the two cumulative incidence functions. Figure 4 also displays the Wald-type estimator of Richardson and others (2017), dotted curve. Their estimator estimates the causal effect among compliers and is constructed to the situation where both exposure and instrument are binaries. For the HIP-data, those randomized not to receive screening never got access to screening (so-called no contamination), and therefore the Wald-type estimator, in this special case, estimates the causal effect among the treated as is the case for our estimator (2.5). Hence in the special setting in the HIP-study the two estimators estimate the same quantity, and indeed it is seen from Figure 4 that they give almost identical results. Fig. 4. Open in new tabDownload slide HIP-study. Estimated absolute risk function concerning breast cancer death (full line), see display (2.7), along with 95% pointwise confidence bands (broken curves). The Wald-type estimator of Richardson and others (2017) is also shown (dotted curve). Fig. 4. Open in new tabDownload slide HIP-study. Estimated absolute risk function concerning breast cancer death (full line), see display (2.7), along with 95% pointwise confidence bands (broken curves). The Wald-type estimator of Richardson and others (2017) is also shown (dotted curve). 5. Concluding remarks In this article, we have proposed an approach to estimate causal effects in a competing risk setting where there may be unobserved confounding. The proposal is based on the availability of an IV. It can accommodate adjustment for baseline covariates, which is sometimes needed to make the IVs assumptions more plausible. Unlike available IV methods, it makes no restriction on the type, nor the distribution of exposure or instrument. We can in particular deal with a situation where the exposure is continuous and the instrument is categorical, or where both are continuous. Dichotomization of the exposure, which is sometimes considered by simpler proposals, is no valid remedy in such cases as it may entail a violation of the exclusion restriction. A further strength of the approach is that it naturally adjusts for censoring whenever censoring is independent of the event time, exposure and instrument conditional on the confounders (even when unmeasured). Although this assumption could fail, it can be remedied by applying inverse probability of censoring weighting; that is, by redefining |$ {\rm e}^{B_1(t)X+B_2(t)X}R(t) $| in (2.4) to $$\[{\rm e}^{B_1(t)X+B_2(t)X}R(t)\frac{P(C\geq t|\mathcal{F}^N_{t}\vee L)}{P(C\geq t|\mathcal{F}^N_{t}\vee X\vee G\vee L)},\]$$ where |$\mathcal{F}^N_{t}$| denotes the history spanned by the counting processes. Using this modification requires postulating two models for the cause-specific hazard of censoring: one conditional on |$X,G,$| and |$L$|⁠, and one conditional on |$L$| only. However, only the former model must be correctly specified to maintain a consistent estimator. The results in this article may be extended to more general models than (2.3), to also handle interactions with the confounder |$L$|⁠. For instance, suppose that instead of (2.3) the following model holds, for |$j=1,2$|⁠, $$\begin{equation}\label{Model2}\lambda^j_{T^x}(t|X=s,G,L)-\lambda^j_{T^0}(t|X=x,G,L)=\beta_j(t)x+\beta^T_{jXL}(t)xl\end{equation}$$ (5.1) where |$\beta_{jXL}(t)$| is of dimension corresponding to |$L$|⁠, say |$p$|⁠. Let |$B_{\ast}(t)$| denote the integral from 0 to |$t$| of the corresponding |$\beta_{\ast}(t)$|⁠, and define |$B(t)=\{B_1(t),B_2(t),B^T_{1XL}(t),B^T_{2XL}(t)\}$|⁠. Let |$dN(t)=\{ {\rm d}N_1^1(t), {\rm d}N_1^2(t),\ldots , {\rm d}N_n^1(t), {\rm d}N_n^2(t)\}^T$|⁠. We can the write the estimator of |$B(t)$| as $$\hat{B}(t,\theta)=\int_0^t\left [\{J_1(s,\theta)J_2(s,\theta)\}^T\{J_1(s,\theta)J_2(s,\theta)\}\right ]^{-1}[\{J_1(s,\theta)J_2(s,\theta)\}^TJ_1(s,\theta) {\rm d}N(s),$$ where |$J_1(t,\theta)$| is |$1\times 2n$|-matrix with |$i$|th row $$\left\{G_i-E(G_i|L_i)\right\}\exp{\left\{\sum_{j=1}^2\hat{B}_j(t-,\theta)X_i+ \hat{B}^T_{jXL}(t-,\theta)X_iL_i\right\}}R_i(t),$$ and |$J_2(t)$| is |$2n\times (2+2p)$|-matrix consisting of |$n$| blocks of size |$2\times (2+2p)$|⁠, where the |$i$|th block is $$\begin{bmatrix}X_i & 0 & X_iL_i^T & 0 \\0 & X_i & 0& X_iL_i^T\end{bmatrix}.$$ We have shown how to estimate functions of the cumulative incidence functions, |$\mbox{RR}(t)$| and |$\mbox{AD}(t)$|⁠, in the case without covariates |$L$| and with the exposure and instrument both being binaries. In principle, this may be extended to the general situation with covariates |$L$| and without restrictions on the exposure and instrument. However, one then need to specify the cause-specific hazard functions, and such models may not be congenial with the assumed structural model. We plan to investigate this further in future work. 6. Software Simulation script with data is available at https://github.com/TorbenMart/IV_Compet_Risk_2018. Acknowledgments Conflict of Interest: None declared. Funding Torben Martinussen’s work is part of the Dynamical Systems Interdisciplinary Network, University of Copenhagen. Stijn Vansteelandt was supported by IAP research network [grant P07/05] from the Belgian government (Belgian Science Policy). Appendix: Large sample properties Consistency Let |$\mu(L;\theta)=E(G|L;\theta)$| be the conditional mean of the instrument given observed confounders |$L$|⁠, which is function of an unknown finite-dimensional parameter |$\theta$|⁠. In the case of no observed confounders |$\mu(\theta)=\theta=E(G)$| and |$\hat \theta=\overline{G}$|⁠. We assume that |$n^{1/2}(\hat\theta-\theta)=n^{-1/2}\sum_i\epsilon^{\theta}_i+o_p(1)$|⁠, where the |$\epsilon^{\theta}_i$|’s are zero-mean iid variables. We write |${\| {g} \|_\infty} = \sup_{t \in [0,\tau]} |g(t)|$|⁠. Let |$B^\circ_j(t)$| denote the true value of |$B_j(t)$|⁠, and let |$M_j^\circ = {\| {B_j^\circ} \|_\infty}<\infty$|⁠. Conditions: (i) We assume that |$X$| and |$G$| are bounded. (ii) Define |$a(s,h) = E[R(s)XG^c e^{h X}]$|⁠. We assume that there exist |$M > M_j^\circ$|⁠, |$j=1,2$|⁠, and |$\nu>0$| such that |$\inf_{s \in [0,\tau], h \in[-M,M]} a(s,h) \geq 1.01 \nu$|⁠. The quantities |$M_j^\circ$| and |$M$| do not necessarily need to be known. Under these assumptions we may modify the arguments given in Martinussen and others (2017) to also cover the competing risk situation described here. Hence, consistency can be inferred similarly. Asymptotic normality Let |$N(t)$| be the |$n\times 2$| matrix with |$i$|th row |$\{N_i^1(t),N_i^2(t)\}$| and |$X=(X_1,\ldots , X_n)^T$|⁠. For known |$\theta$|⁠, we can write $$\hat B(t,\theta)=\{\hat B_1(t,\theta),\hat B_2(t,\theta)\}=\int_0^tH_{\theta}\{s,\hat B_1(s-,\theta)+\hat B_2(s-,\theta)\}dN(s),$$ where the |$k$|th element of the |$n$|-vector |$H_{\theta}\{t,a\}$| is $$R_k(t)G_k^c(\theta)e^{aX_k}/\sum_{i=1}^nG_i^c(\theta)R_i(t)e^{aX_i}X_i.$$ with |$a=\hat B_1(t-,\theta)+\hat B_2(t-,\theta)$|⁠. Let |$W_n(t,\theta)=n^{1/2}\{\hat B(t,\theta)- B(t)\}$|⁠, |$b=(1,1)^T$|⁠, and let |$\dot{H}$| denote the derivative of |$H$| with respect to its second argument. It is then easy to see that $$\begin{align*}W_n(t,\theta)=&n^{1/2}\int_0^tH\{s, B_1(s-)+B_2(s-)\}\left [{\rm d}N(s)-X{\rm d}B(s)\right ]\\&+\int_0^tV(s-,\theta)\{1+o_p(1)\}b\dot{H}\{s,B_1(s-)+B_2(s-)\}{\rm d}N(s) \end{align*}$$ which is a Volterra-equation, see Andersen and others (1993, p 91). The solution to this equation is given by $$W_n(t,\theta)=\int_0^tn^{1/2}H\{s,B_1(s-)+B_2(s-)\}\left [{\rm d}N(s)-X{\rm d}B_X(s)\right ]\mathcal{F}(s,t)+o_p(1),$$ where $$\mathcal{F}(s,t)=\prod_{(s,t]}\left \{I+b\dot{H}(\cdot,B_1(\cdot )+B_2(\cdot )){\rm d}N(\cdot)\right \}$$ with the latter being a product integral that converges in probability to some limit. This leads to the iid-representation $$W_n(t,\theta)=n^{-1/2}\sum_{i=1}^n\epsilon^B_i(t)$$ with the |$\epsilon^B_i(t)$|’s being zero-mean iid terms. Specifically $$\epsilon^B_i(t)=\int_0^tn^{1/2}\{H\{s,B_1(s-)+B_2(s-)\}\}_i\left [{\rm d}N(s)-X{\rm d}B_X(s)\right ]_i\mathcal{F}(s,t)$$ with |$a_i$| being the |$i$|th element of the vector |$a$|⁠, and |$[A]_i$| being the |$i$|th row of the matrix |$A$|⁠. This together with $$\begin{align*}n^{1/2}\{\hat B(t,\hat\theta)-B(t)\}&=n^{1/2}\{\hatB(t,\theta)-B(t)\}+n^{1/2}\{\hat B(t,\hat\theta)-\hat B(t,\theta)\}\\&=n^{1/2}\{\hat B(t,\theta)-B(t)\}+D_{\theta}(\hatB(t,\theta)^T)_{|\hat\theta}n^{1/2}(\hat\theta-\theta)+o_p(1),\end{align*}$$ where |$D_{\theta}\{\hat B(t,\theta)^T\}$| is the first order derivative of |$\hat B(t,\theta)^T$| w.r.t. |$\theta$| gives an iid-decomposition of |$n^{1/2}\{\hat B(t,\hat\theta)-B(t)\}$|⁠: $$n^{1/2}\{\hatB(t,\hat\theta)-B(t)\}=n^{-1/2}\sum_{i=1}^n\epsilon_i^B(t,\theta)+o_p(1),$$ where $$\begin{equation}\label{Eps}\epsilon_i^B(t,\theta)=\epsilon^B_i(t)+D_{\theta}(\hatB(t,\theta)^T)_{|\theta}\epsilon^{\theta}_i.\end{equation}$$ (A.1) It thus follows that $$n^{1/2}\{\hat B(t,\hat\theta)-B(t)\}$$ converges to a zero-mean Gaussian process with a variance that is consistently estimated by $$n^{-1}\sum_{i=1}^n\hat\epsilon_i^B(t,\hat\theta)^2.$$ The derivative |$D_{\theta}(\hat B(t,\theta))_{|\hat\theta}$| can be calculated recursively as |$\hat B(t,\hat\theta)$| is constant between the observed death times. Denote the jump times by |$\tau_1,\ldots ,\tau_m$|⁠. Hence $$\hat B(\tau_j,\theta)=\hat B(\tau_{j-1},\theta)+d\hat B(\tau_j,\theta)$$ which then also holds for the derivative. Since |$\hat B(0,\theta)=0$| and the derivative of the increment in the first jump time, |$d\hat B(\tau_1,\theta)$|⁠, is easily calculated, then we have a recursive way of calculating the derivatives of |$\hat B(\cdot,\theta)$|⁠. We now argue that the process |$W_n^j(t,\theta)$|⁠, |$j=1,2$|⁠, converges in distribution as a process using arguments similar to what is done in Lin and others (2000, p. 726). It is seen from (2.3) that |$B_j(t)$| can be written as a difference of two monotone functions. Let |$\tilde H_{ijk}(s)$| be the limit in probability of |$\mathcal{F}_{kj}(s,t)H_i(s,B_X(s-))$|⁠. Now, split |$\tilde H_{ijk}(s)$| into its positive and negative parts, |$\tilde H_{ijk}^+(s)$| and |$\tilde H_{ijk}^-(s)$|⁠, and similarly with |$X_i$|⁠, |$X_i^+$| and |$X_i^-$|⁠. Then |$\int_0^t\tilde H_{ijk}(s)[ {\rm d}N_i(s)-X_i{\rm d}B_X(s) ]_{ik}$| can be written as a difference of two monotone functions, and then we follow the arguments of Lin and others (2000) (or use example 2.11.16 of van der Vaart and Wellner, 1996). Convergence in distribution for the process |$W_n^j(t,\hat \theta)$| also holds using the above Taylor expansion. It thus follows that $$n^{1/2}\{\hat B(t,\hat\theta)-B(t)\}$$ converges to a zero-mean Gaussian process. References Andersen P. K. , Borgan O. , Gill R. D. and Keiding N. Statistical Models Based on Counting Processes . Berlin : Springer-Verlag . Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Andersen P. K. , Geskus R. B. , De Witte T. and Putter H. ( 2012 ). Competing risks in epidemiology: possibilities and pitfalls . International Journal Epidemiology , 41 , 861 – 70 . Google Scholar Crossref Search ADS WorldCat Andersen P. K. and Keiding N. ( 2012 ). Interpretability and importance of functionals in competing risks and multistate models . Statistics in Medicine 31 , 1074 – 1088 . Google Scholar Crossref Search ADS PubMed WorldCat Cai B. , Small D. S. and Ten Have T. R. ( 2011 ). Two-stage instrumental variable methods for estimating the causal odds ratio: analysis of bias . Statistics in Medicine , 30 , 1809 – 1824 . Google Scholar Crossref Search ADS PubMed WorldCat Hernán M. A. and Robins J. M. ( 2006 ). Instruments for causal inference: an epidemiologist’s dream? Epidemiology 17 , 360 – 372 . Google Scholar Crossref Search ADS PubMed WorldCat Joffe M. M. ( 2001 ). Administrative and artificial censoring in censored regression models. Statistics In Medicine 20 , 2287 – 2304 . Google Scholar Crossref Search ADS PubMed WorldCat Joffe M. M. , Yang W. P. and Feldman H. ( 2012 ). G-estimation and artificial censoring: problems, challenges, and applications. Biometrics 68 , 275 – 286 . Google Scholar Crossref Search ADS PubMed WorldCat Kjaersgaard M. I. S. and Parner E. T. ( 2016 ). Instrumental variable method for time-to-event data using a pseudo-observation approach. Biometrics 72 , 463 – 472 . Google Scholar Crossref Search ADS PubMed WorldCat Li J. , Fine J. and Brookhart A. ( 2014 ). Instrumental variable additive hazards models. Biometrics , 71 , 122 – 130 . Google Scholar Crossref Search ADS PubMed WorldCat Martinussen T. and Scheike T. H. ( 2006 ). Dynamic Regression Models for Survival Data. New York, NY : Springer . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Martinussen T. , Vansteelandt S. , Tchetgen Tchetgen E. J. and Zucker D. M. ( 2017 ). Instrumental variables estimation of exposure effects on a time-to-event endpoint using structural cumulative survival models . Biometrics . 10.1111/biom.12699 . OpenURL Placeholder Text WorldCat Crossref Richardson A. , Hudgens M. G. , Fine J. and Brookhart A. ( 2017 ). Nonparametric binary instrumental variable analysis of competing risks data . Biostatistics 18 , 48 – 61 . Google Scholar Crossref Search ADS PubMed WorldCat Robins J. M. and Rotnitzky A. ( 2004 ). Estimation of treatment effects in randomised trials with non-compliance and a dichotomous outcome using structural mean models . Biometrika 91 , 763 – 783 . Google Scholar Crossref Search ADS WorldCat Robins J. M. and Tsiatis A. ( 1991 ). Correcting for non-compliance in randomized trials using rank-preserving structural failure time models . Communications in Statistics 20 , 2609 – 2631 . Google Scholar Crossref Search ADS WorldCat Shapiro S. ( 1977 ). Evidence of screening for breast cancer from a randomised trial . Cancer 39 , 2772 – 2782 . Google Scholar Crossref Search ADS PubMed WorldCat Tchetgen Tchetgen E. J. , Walter S. , Vansteelandt S. , Martinussen T. , Glymour M. ( 2015 ). Instrumental variable estimation in a survival context . Epidemiology 26 , 402 – 410 . Google Scholar Crossref Search ADS PubMed WorldCat Zheng C. , Dai R. , Hari P. N. and Zhang M. J. ( 2017 ). Instrumental variable with competing risk model . Statistics in Medicine 10.1002/sim.7205 . OpenURL Placeholder Text WorldCat Crossref © The Author 2018. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Instrumental variables estimation with competing risk data
JF - Biostatistics
DO - 10.1093/biostatistics/kxy039
DA - 2020-01-01
UR - https://www.deepdyve.com/lp/oxford-university-press/instrumental-variables-estimation-with-competing-risk-data-0MeTv0kSPR
SP - 158
VL - 21
IS - 1
DP - DeepDyve
ER -