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.
86,286 characters · 21 sections · 90 citation commands
Estimation of the complier causal hazard ratio under dependent censoring
Acknowledgments: The authors thank the National Cancer Institute for access to NCI's data collected by the HIP Breast Cancer Screening Trial (HIPB). The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government department EWI. Moreover, the authors thank Ilias Willems and Sofia Guglielmini for their advice on high performance computing.
Funding: G. Crommen is funded by a PhD fellowship from the Research Foundation - Flanders (grant number 11PKA24N). J. Beyhum acknowledges support from the FWO (Research Foundation Flanders) through the project G018725N. I. Van Keilegom acknowledges funding from the FWO and F.R.S. - FNRS (Excellence of Science programme, project ASTeRISK, grant no. 40007517), and from the FWO (senior research projects fundamental research, grant no. G047524N).
When estimating the causal effect of a binary treatment variable $Z$ on a right-censored duration outcome $T$, the presence of unobserved heterogeneity poses a significant challenge. Unobserved heterogeneity refers to unmeasured confounders that influence both the treatment $Z$ and the duration outcome $T$, such that $Z$ is endogenous. Consequently, the causal effect of $Z$ on $T$ cannot be identified from the conditional distribution of $T$ given $(Z,X^{\top})$, where $X$ represents a vector of observed exogenous covariates. To address the issue of endogeneity, instrumental variable (IV) methods are commonly used LATEimbensangrist1994, angrist1995two, angrist1996identification, ABADIE2003231. These methods use external sources of variation (instruments) that influence the treatment $Z$, but are independent of the unobserved confounders affecting both $Z$ and $T$. However, in the context of censored duration outcomes, most IV approaches assume that the censoring time $C$ is (conditionally) independent of the survival time $T$. This assumption simplifies the analysis but is often unrealistic in practice and leads to biased estimates moeschberger1984consequences, emura2016gene. In their seminal work, kaplan1958nonparametric even stated that: “In either case it is usually assumed in this paper that the lifetime (age at death) is independent of the potential loss time; in practice, this assumption deserves careful scrutiny". Therefore, the main objective of this paper is to propose an IV method that can identify the causal effect of an endogenous binary treatment variable on a dependently censored duration outcome.
Our motivating example is the National Job Training Partnership Act (JTPA) study, which was designed to evaluate the efficacy of publicly funded job training programs on a range of outcomes, including unemployment duration and earnings. The unemployment duration data have been analyzed under the independent censoring assumption by frandsen2015treatment and beyhum2024instrumental among others. In the JTPA study, individuals were randomly assigned to either a treatment group, which was eligible for job training services, or a control group, which was ineligible for 18 months. It is important to note that participants were not obliged to adhere to the treatment that they were assigned, allowing for some movement between groups. The issue of unobserved heterogeneity arises because participants may switch between the treatment and control groups in a way that is related to their time until employment but not measured (e.g. motivation to find employment). In this context, the randomized treatment assignment serves as a natural IV. The assumption of (conditional) independent censoring would be violated in cases where a participant's decision to drop out of the study is influenced by their employment status. For the stratum of fathers who reported having no job at the time of randomization, crommen2024instrumental found a significant negative dependence between the unemployment duration and censoring time, conditional on the measured covariates. Moreover, frandsen2019testing found that the independent censoring assumption is often questionable for unemployment duration data. Therefore, it seems likely that the assumption of (conditional) independent censoring is not satisfied for these data.
Instrumental variable methods have been developed and applied across econometrics and related fields to address the issue of endogeneity. Seminal contributions by robins1989analysis, manski1990nonparametric and balke1997bounds demonstrate that under the standard assumptions of random assignment, exclusion restrictions and monotonicity, the average treatment effect (ATE) is not point-identified. Therefore, LATEimbensangrist1994 and angrist1996identification proposed focusing on specific subpopulations rather than attempting to identify effects for the entire population. In particular, they demonstrated that the IV framework could identify the local average treatment effect (LATE), also referred to as the complier average causal effect, for a subgroup of individuals called the compliers. The identification strategy is based on a monotonicity assumption, which rules out the existence of another subgroup called the defiers. Compliers are defined as those whose treatment status is directly influenced by the instrument, making the treatment assignment effectively act as a randomized experiment within this subgroup. Subsequent work by ABADIE2003231 introduced a new class of IV estimators designed for linear and nonlinear treatment response models with covariates, providing additional flexibility. To estimate population-level causal effects, researchers can impose additional parametric assumptions efron1991compliance or replace the monotonicity assumption with the rank invariance assumption as introduced by chernozhukov2005iv.
Most instrumental variable methods for right-censored duration outcomes assume that the duration time $T$ is (conditionally) independent of the censoring time $C$. This assumption is known as independent censoring, or conditional independent censoring if it is assumed conditional on the measured covariates. Some non-parametric approaches are given by frandsen2015treatment, SantAnnaPedroH.C2016PEwR and BeyhumJad2021NIRW, while other approaches such as BijwaardGovertE2005Cfsc, TchetgenTchetgenEricJ2015IVEi, ChernozhukovVictor2015Qrwc, LiJialiang2015Ivah, kianian2021causal, beyhum2024instrumental and tedesco2023instrumentalvariableestimationproportional are semiparametric. Moreover, van2020nonparametric and beyhum2024instrumentaldynamic consider the problem of identifying and estimating the causal effect of a dynamic treatment. The problem of endogeneity has also been discussed in a competing risks framework by RichardsonAmy2017Nbiv, ZhengCheng2017Ivwc, martinussen2020instrumental and beyhum2023nonparametric among others.
Since the (conditional) independent censoring assumption is often unrealistic and empirically untestable, models that allow for dependent censoring have been developed. Among the most popular approaches are those based on copulas, where a first line of research focuses on fully known copulas ZHENGMING1995Eoms, Rivest2001AMA, BraekersRoel2005Acef, HuangXuelin2008RSAw, SujicaAleksandar2018Tcei. A copula being fully known means that the association parameter specifying the dependence between $T$ and $C$ is known. Since the assumption of a fully known copula is usually equally unrealistic in practice as the independent censoring assumption, czado2021dependent introduced a new approach that does not require the association parameter of the copula to be specified. However, this approach requires the marginals for $T$ and $C$ to be fully parametric to identify the association parameter. This approach was extended in various ways by DeresaNegeraWakgari2020Fpmf, DERESA2020106879diftypecens and semiparderesa2020 among others. More recently, deresa2021copulacox showed that (under some additional assumptions regarding the copula function and the covariates) the model remains identified when the conditional marginal distribution of $T$ follows a semi-parametric proportional hazards model. The additional assumption on the covariates requires that, given 2 distinct continuous covariates $X_1$ and $X_2$, the conditional distribution of $T$ does not depend on $X_1$ and the conditional distribution of $C$ does not depend on $X_2$\footnote{Note that this assumption is not explicitly mentioned by deresa2021copulacox, but is necessary for their identifiability proof.}. This type of covariate restriction was also used by deresa2024semiparametric to allow for both $T$ and $C$ to follow a semi-parametric transformation model, and by hiabu2025identifiability to allow for non-parametric marginals. A non-parametric approach that does not rely on this type of covariate restriction has been proposed by Lo_Wilke_2017, but they only identify the sign of covariate effects on the marginals. Note that all of these approaches assume that the treatment is exogenous.
The body of literature addressing instrumental variable methods in the context of dependent censoring remains relatively limited. KhanShakeeb2009Ioec examined an endogenously censored regression model. However, their framework imposes a restrictive support condition (Assumption IV2, page 110) concerning the relationship between the instruments and covariates. Moreover, blanco2020bounds investigated treatment effects on duration outcomes in the presence of censoring, selection, and noncompliance. Instead of providing point estimates, their work derives bounds for the causal effect. More recently, ying2024proximal introduced a proximal survival analysis framework, which draws on ideas from proximal causal inference. This approach addresses dependent censoring but requires the analyst to classify observed covariates into three specific categories, which can be challenging in practice. Furthermore, crommen2024instrumental proposed a fully parametric model that identifies the causal effect of an endogenous treatment variable on a potentially dependently censored event time. Their method employs a control function approach to handle the endogeneity issue, and the possible dependent censoring is taken into account by relying on the strong parametric assumption of bivariate normal error terms of the joint regression model for the logarithms of $T$ and $C$. To make the normality assumption more realistic, rutten2024flexible extended this work to a competing risks framework where different power transformations are applied to each risk.
We propose an IV method that allows for dependent censoring, enabling the identification and estimation of causal effects in settings where the (conditional) independent censoring assumption may not hold. More precisely, we are interested in the complier causal hazard ratio (CCHR), which is the ratio of compliers' hazard function when the treatment is received versus when it is not. As mentioned before, the compliers are a latent subgroup who adhere to their assigned treatment status. However, a major challenge is that the compliers are a latent subgroup, meaning we cannot directly identify if an individual belongs to the compliers from the observed data. This is because, for any given individual, we only observe their treatment choice under their assigned treatment group. To overcome this issue, we rely on a monotonicity assumption to leverage the weighting scheme from Theorem 3.1 by ABADIE2003231. This theorem establishes a connection between the conditional moment of any measurable real function of the observed data for the compliers and the unconditional moment. The monotonicity assumption rules out the existence of defiers, a latent subgroup for whom the treatment choice is always different from the assigned treatment. Note that this assumption is automatically satisfied when there is only one-sided noncompliance. To account for the potential dependence between $T$ and $C$, conditional on the measured covariates, it is assumed on the stratum of compliers that $T$ follows a semiparametric proportional hazards model, $C$ follows a fully parametric model and their joint distribution is modeled by a parametric copula. In addition, we allow for an independent right censoring time $A$ (e.g. administrative censoring) such that only the minimum of $T,C$ and $A$ is observed. This approach enables us to achieve the following contributions:
In Section (ref), we introduce the potential outcomes framework and specify the model. The identification is discussed in Section (ref) and Section (ref) describes the two-step estimation procedure. Section (ref) establishes the asymptotic properties of the proposed estimator. Simulation results and an empirical application regarding the effect of job training programs on unemployment duration are described in Sections (ref) and (ref), respectively. The proofs, technical details, extra simulation results and a second data application can be found in the Supplementary Material. The R code can be found at \url{https://github.com/GillesCrommen/CCHR}.
Let $T$ and $C$ be the duration and right censoring time respectively. We will allow for the variables $T$ and $C$ to be dependent on each other, even after conditioning on the measured covariates. The covariates that influence both $T$ and $C$ are given by $(Z,X^{\top})$, where $Z$ and $X$ are of dimension 1 and $m$ respectively. Note that $X$ does not include an intercept and that $Z$ represents an endogenous binary treatment variable for which a binary instrumental variable, denoted by $W$, exists. In addition, we allow for another right censoring time $A$, which is independent of $T,C$ given the measured covariates and the measured covariates themselves. A common example of observing $A$ is when a subject is administratively censored. Because $T, C$ and $A$ censor each other, only one of them is observed through the follow-up time $Y = \text{min}\{T, C, A\}$ and the censoring indicators $\Delta_1 = \mathbbm{1}(Y = T)$ and $\Delta_2 = \mathbbm{1}(Y = C)$, where $\mathbbm{1}(\cdot)$ is the indicator function.
To state the necessary assumptions, we will use the potential outcome framework as described by rubin1974estimating, LATEimbensangrist1994 and ABADIE2003231 among others. Firstly, let $Z_{W=w}$ denote the potential treatment selection under instrument $W=w$, where $w \in \{0,1\}$, such that the observed treatment is given by $Z=W Z_{W=1} + (1-W) Z_{W=0}$. Further, let $Y_{Z=z}$ denote the potential outcome given $Z=z$, where $z \in \{0,1\}$, such that the observed outcome is given by $Y=Z Y_{Z=1} + (1-Z) Y_{Z=0}$. Using this notation, we can classify subjects into four latent subgroups: compliers ($Z_{W=1} > Z_{W=0}$), always takers ($Z_{W=1} = Z_{W=0} = 1$), never takers ($Z_{W=1} = Z_{W=0} = 0$) and defiers ($Z_{W=1} < Z_{W=0}$). We will restrict our attention to the compliers, a subgroup for whom the treatment choice acts as randomization. For ease of notation, let $G$ denote a latent variable that indicates the subgroup (i.e. $G=g$ where $g \in \{co,at,nt,df\}$ indicates compliers, always takers, never takers and defiers respectively). We will make the following standard IV assumptions for almost all values of $X$:
Assumption (ref) states that both the potential outcomes and the potential treatments are independent of the instrument given $X$. Note that this assumption is the combination of two commonly made requirements in an IV framework. Firstly, $W$ follows a random assignment conditional on $X$ such that we can identify the causal effect of the instrument on the treatment, which is commonly known as an exogeneity assumption. Secondly, the potential outcomes are not directly affected by the instrument, which is commonly known as an exclusion restriction. Therefore, the instrument only affects the outcome through variation in the treatment selection. Note that this assumption is made on both the potential outcomes of the duration and censoring time, which is a consequence of allowing for some dependence between $T$ and $C$. Assumption (ref) is a standard monotonicity assumption that rules out the existence of defiers. This assumption is automatically satisfied when individuals in the control group ($W=0$) cannot access the treatment. Finally, the first part of Assumption (ref) requires non-trivial assignment of the instrument, while the second part is a standard first-stage assumption. Note that Assumption (ref) implies that the fraction of compliers is non-zero.
We will further assume that:
Assumptions (ref)-(ref) are standard assumptions in the survival analysis literature and Assumption (ref) is necessary for the identification of the model. Lastly, Assumption (ref) is required to ensure that the denominator of the estimator for the baseline cumulative hazard function in Section (ref) is non-zero. Note that because of Assumption (ref), the endpoint of $\mathcal{Y}$, where $\mathcal{Y} \subseteq \mathbb{R}_{> 0}$ is the support of $Y$, should be greater than $\Bar{\tau}$. It is important to emphasize that the following model specifications are only assumed for the subgroup of compliers.
We start by noticing that Assumption (ref) and $W=Z$ for the compliers imply that
such that, in the population of the compliers, the causal effect of $Z$ on $T$ can be identified from $F_{T \mid Z,X,G}(t\mid z,x,co)$. Following deresa2021copulacox, a proportional hazards model (see cox1972regression for details) is assumed for the conditional distribution of $T$ given $(Z,X^{\top})$ for the subgroup of compliers: $$ F_{T \mid Z,X,G}(t\mid z,x,co) =1-\exp\{-\Lambda(t\mid G=co) \exp(z \alpha + x^\top \beta)\}, $$ where $\mu = (\alpha,\beta^\top) \in \mathcal{M} \subset \mathbb{R}^{m+1} $ is a vector of regression coefficients and $\Lambda(t \mid G=co)$, with $\Lambda(0 \mid G= co)=0$, is an unknown strictly increasing differentiable baseline cumulative hazard function. Further, let $\lambda(t\mid G=co)=\frac{ \mathop{}\!d \Lambda(t\mid G=co)}{\mathop{}\!d t}$ be the baseline hazard function for the compliers and $$ f_{T \mid Z,X,G}(t\mid z,x,co)=\lambda(t\mid G=co)\exp(z \alpha + x^\top \beta)\exp\left\{-\Lambda(t\mid G=co) \exp(z \alpha + x^\top \beta)\right\}, $$ the conditional density function of $T$ given $(Z,X^{\top})$ for the compliers. This allows us to define the conditional hazard function for the compliers as $$ \lambda(t\mid Z=z,X=x,G=co) = \frac{f_{T \mid Z,X,G}(t\mid z,x,co)}{1-F_{T \mid Z,X,G}(t\mid z,x,co)}= \lambda(t\mid G=co)\exp(z \alpha+ x^\top \beta). $$ From this it is clear that $$ \alpha =\log\lambda(t\mid Z=1,X=x,G=co)-\log\lambda(t\mid Z=0,X=x,G=co), $$ which shows that $\exp(\alpha)$ has a causal interpretation as the ratio of compliers' hazard function when the treatment is received versus when the treatment is not received. This means that when $\alpha$ is significantly greater (smaller) than zero, the treatment reduces (increases) the survival time of the compliers. We will refer to the causal effect of interest, $\exp(\alpha)$, as the complier causal hazard ratio (CCHR).
For the conditional distribution of the censoring time $C$ of the compliers, given $(Z,X^{\top})$, a fully parametric model is assumed: $$ F_{C \mid Z,X,G}( \cdot \mid z,x,co) \in \big\{ F_{C \mid Z,X,G}( \cdot \mid z,x,co ; \eta) : \eta \in \mathcal{H} \big\},$$ for a finite-dimensional parameter space $\mathcal{H}$, where $F_{C \mid Z,X,G}(c\mid z,x,co)=\mathbb{P}(C\leq c\mid Z=z,X=x,G=co)$ and $f_{C \mid Z,X,G}(c\mid z,x,co)$ the corresponding conditional density for the compliers. This parametric model will be necessary to ensure the identifiability of the model. To account for the possible dependence between $T$ and $C$, given $(Z,X^{\top})$ for the compliers, we propose to use a bivariate copula $\mathcal{C}$. The joint distribution function of $T$ and $C$, conditional on $(Z,X^{\top})$ for the compliers, can therefore be modeled as $$ \mathbb{P}(T \leq t, C\leq c\mid Z=z,X=x,G=co)=\mathcal{C}\left(F_{T \mid Z,X,G}(t\mid z,x,co),F_{C \mid Z,X,G}(c\mid z,x,co)\right), $$ with $t,c > 0$ and for some copula $\mathcal{C}$. The copula $\mathcal{C}$ is defined as a bivariate distribution function with uniform marginals over the unit interval. Note that Assumption (ref) implies that the copula function is unique sklar1959fonctions. We will assume that $\mathcal{C}$ is a parametric copula, meaning that $$ \mathcal{C} \in \{ \mathcal{C}_\xi : \xi \in \Xi \}, $$ for some parameter space $\Xi$. Moreover, for $u,v \in [0,1]$ let $$ \zeta_{1}(u,v)=\frac{\partial}{\partial u}\mathcal{C}(u,v), \quad \zeta_{2}(u,v)=\frac{\partial}{\partial v}\mathcal{C}(u,v). $$ It can be shown that (see czado2021dependent):
The distribution of $A$, denoted by $F_A(a) = \mathbb{P}(A \leq a)$ with $f_A(a)$ the corresponding density, is left completely unspecified.
Further, let $F_{Y,\Delta_1,\Delta_{2} \mid Z,X,G}(y,\delta_1,\delta_2 \mid z,x,co)$ denote the sub-distribution function of $(Y,\Delta_1,\Delta_2)$ given $(Z,X^{\top})$ for the compliers, that is $$ F_{Y,\Delta_1,\Delta_{2} \mid Z,X,G}(y,\delta_1,\delta_2 \mid z,x,co) = \mathbb{P}(Y \leq y,\Delta_1 = \delta_1,\Delta_2 = \delta_2 \mid Z=z, X=x,G=co), $$ with the corresponding sub-density $f_{Y,\Delta_1,\Delta_{2} \mid Z,X,G}(y,\delta_1,\delta_2 \mid z,x,co)$, for which it easily can be shown that:
and $$ f_{Y,\Delta_1,\Delta_{2} \mid Z,X,G}(y,0,0 \mid z,x,co)= S(y \mid z,x,co) \times f_A(y), $$ where
Lastly, we have that
Note that when we want to stress the dependence of these functions on their respective parameters, we will add them to their notation (e.g. $S_{\mu,\Lambda,\eta,\xi}$ instead of $S$).
In this section, we will discuss the identifiability of our model. Note that we have two identifiability issues. Firstly, we cannot identify if a particular individual belongs to the subgroup of compliers since only one of $Z_{W=0}$ and $Z_{W=1}$ is observed. Secondly, we want to show that two different parameter vectors result in two different joint distributions of $S=(Y,\Delta_1,\Delta_2,Z,X^\top)^\top$ for the subgroup of compliers.
We will start by restating Theorem 3.1 by ABADIE2003231. This theorem shows that even though we cannot identify if a particular individual belongs to the subgroup of compliers, we can identify expectations for compliers.
Because of this theorem, we have that any statistical characteristic that is defined in terms of moments of $S$ is identified for the compliers.
Let $\mathcal{S} = [0,\Bar{\tau}] \times \{0,1\}\times \{0,1\}\times \{0,1\} \times \mathcal{X}$ be the support of $S$, where $\mathcal{X} \subset \mathbb{R}^{m}$ is the support of $X$. Note that we have restricted $\mathcal{Y}$ to $[0,\Bar{\tau}]$ as this is the relevant support for $Y$. The problem of identifiability under dependent censoring has been studied by deresa2022copula and czado2021dependent for a fully parametric model and by deresa2021copulacox for a semi-parametric model. To show the identifiability of the model for the subgroup of compliers, we will need the following conditions:
It has already been shown by czado2021dependent that Condition (ref) is satisfied for the parametric families of log-normal, log-logistic, log-Student-t and Weibull densities. More recently, it was shown by delhelle2024copula that this condition is also satisfied for the parametric family of Gamma densities. Furthermore, let the strict Clayton copula be given by $$ \mathcal{C}_{\xi}(u,v) = (u^{-\xi} + v^{-\xi} -1)^{-1/\xi}, \quad \xi > 0, $$ and the rotated copulas by
where the superscript indicates the angle of rotation. For ease of notation, let Clayton($\gamma$) be the rotated strict Clayton copula with $\gamma$ the angle of rotation. We now have the following result, for which the proof can be found in Section A of the Supplementary Material.
Although the conditions for the Gumbel and Gaussian copula given in Lemma (ref) seem different from the ones given in Lemma 3.1 by deresa2021copulacox, they are equivalent. The advantage of presenting the conditions in this form is that they are much easier to interpret. Further, it is to be noted that the non-rotated Clayton copula is not mentioned by Lemma (ref). This is because Condition (ref) can never hold for the Clayton copula, as it would require that $\lim_{y \to 0}F_{T \mid Z,X,G}(y\mid z,x,co)/F_{C \mid Z,X,G}(y\mid z,x,co) \text{ converges to } 0 \text{ and } \infty \text{ at the same time. }$ However, this does not necessarily mean that the model is not identified, since these conditions are only sufficient for identification. We now have the following identifiability theorem:
The proof of Theorem (ref) can be found in Section A of the Supplementary Material. It is to be noted that Conditions (ref)-(ref) are similar to the ones given by deresa2021copulacox, but differ as Condition (ref) is about the quotient of $\zeta_{2,\xi_1}$ and $\zeta_{2,\xi_2}$, instead of the respective copula densities, and we have added Condition (ref). By using these conditions, we avoid having to assume that, given 2 distinct continuous covariates $X_1$ and $X_2$, the conditional distribution of $T$ does not depend on $X_1$ and the conditional distribution of $C$ does not depend on $X_2$. This covariate restriction is not explicitly mentioned by deresa2021copulacox, but is necessary for their identification proof. We have verified in Lemma (ref) that the copula functions meeting the identifiability requirements of deresa2021copulacox also satisfy our proposed conditions, thereby relaxing the identifiability result by removing the untestable covariate restriction. Moreover, Lemma (ref) verifies the identifiability conditions for additional copula functions to allow for more modeling flexibility.
The observed data consists of $n$ i.i.d. realizations of $O = (Y,\Delta_1,\Delta_2,Z,X^\top,W)$, which we denote by $O_i = \{Y_i,\Delta_{1,i},\Delta_{2,i},Z_i,X_i^\top,W_i\}_{i=1,...,n}$. If we let $\theta = (\mu,\eta,\xi)$, the profile pseudo-likelihood function for the compliers can be given by
Note that the distribution of $A$ can be omitted from the likelihood function due to Assumption (ref). However, maximizing this likelihood directly is impossible as $G$ is a latent variable. Moreover, even if we were to observe $G$, maximizing this likelihood still poses a challenge as $\Lambda$ is an unknown function. To solve the problem of not observing $G$, we can use Theorem 3.1 by ABADIE2003231 as described in Section (ref). However, it is important to note that $\omega$ being negative when $Z \neq W$ can cause non-convexity issues to arise during the maximization of the weighted log-likelihood. Moreover, when we define our estimator for $\Lambda$, it will be clear that the negativity of $\omega$ can cause our estimated baseline cumulative hazard for the compliers to be non-monotonic. Therefore, we follow the approach by abadie2002instrumental and replace $\omega$ by $\operatorname{\mathbb{E}}[\omega \mid S]$, where $\operatorname{\mathbb{E}}[\omega \mid S]$ is given by $$ \kappa^*(S) = \operatorname{\mathbb{E}}[\omega \mid S]=\mathbb{P}(G = co \mid S=s) = 1-\frac{Z\mathbb{P}(W=0\mid S)}{\mathbb{P}\big(W = 0 \mid X\big)}-\frac{(1-Z)\mathbb{P}(W = 1\mid S) }{\mathbb{P}\big(W = 1 \mid X\big)}, $$ with $S=(Y,\Delta_1,\Delta_2,Z,X^\top)^\top$. It is shown in Section B of the Supplementary Material that $\kappa^*(S)$ is indeed the probability of being a complier conditional on $S$, which is therefore bounded between $0$ and $1$. In what follows, we start by proposing an estimator for the conditional probability of an individual being a complier given $S$. Secondly, we will replace the unknown function $\Lambda$ with an estimator $\hat{\Lambda}_{\theta}$ for given values of the finite-dimensional parameter $\theta$. This way, we can estimate $\theta$ by solving the score equations that are derived from the pseudo-profile likelihood function $L(\theta,\hat{\Lambda}_{\theta})$.
We propose a non-parametric estimator of $\kappa^*(\cdot)$ as described by wei2021estimation. Firstly, let $R=(Y,X^\top),$ $ \pi(X)=\mathbb{P}(W=1\mid X),$ $ \nu(S)= \mathbb{P}(W=1\mid S)$ and $$ \nu_{j,l,k}(R)= \mathbb{P}(W=1\mid R,\Delta_1 = j,\Delta_2 = l,Z=k), \quad \text{with} \quad j,l,k=0,1. $$ From this, it follows that $$ \nu(S_i)=\sum^1_{j=0} \sum^1_{l=0} \sum^1_{k=0}\mathbbm{1}(\Delta_{1,i} = j,\Delta_{2,i}=l, Z_i=k)\nu_{j,l,k}(R_i). $$ Further, let $K_{h_{1,n}}(x) = h_{1,n}^{-m}K_1(\frac{x}{h_{1,n}})$ and $K_{h_{2,n}}(r)= h_{2,n}^{-(m+1)}K_2(\frac{r}{h_{2,n}})$ with $K_1$ and $K_2$ (multiplicative) kernel functions. Note that we have assumed that all the covariates in $X$ are continuous (otherwise we stratify on the discrete ones) such that we can estimate $\pi(x)$ and $\nu_{j,l,k}(r)$ by $$ \hat{\pi}(x) = \frac{\sum_{i=1}^n K_{h_{1,n}}(x-X_i)W_i}{\sum_{i=1}^n K_{h_{1,n}}(x-X_i)}$$ and $$\hat{\nu}_{j,l,k}(r) = \frac{\sum_{i=1}^n \mathbbm{1}(\Delta_{1,i} = j,\Delta_{2,i}=l, Z_i=k)K_{h_{2,n}}(r-R_i)W_i}{\sum_{i=1}^n \mathbbm{1}(\Delta_{1,i} = j,\Delta_{2,i}=l, Z_i=k)K_{h_{2,n}}(r-R_i)}. $$ From this, we estimate $\nu(S_i)$ by $\hat{\nu}(S_i)=\sum^1_{j=0} \sum^1_{l=0} \sum^1_{k=0}\mathbbm{1}(\Delta_{1,i} = j,\Delta_{2,i}=l, Z_i=k)\hat{\nu}_{j,l,k}(R_i)$. A non-parametric estimator of $\kappa(S_i)$ is then given by $$ \hat{\kappa}(S_i)= 1-\frac{Z_i\big(1-\hat{\nu}(S_i)\big)}{1-\hat{\pi}(X_i)}-\frac{(1-Z_i)\hat{\nu}(S_i) }{\hat{\pi}(X_i)}. $$ Note that $\hat{\kappa}(S_i)$ is a probability and should therefore be bounded between 0 and 1. Therefore, it is proposed to replace $\hat{\kappa}(S_i)$ by $\Tilde{\kappa}(S_i) = \min(\max(\hat{\kappa}(S_i), a_{l,n}),a_{u,n})$, where $a_{l,n}$ and $a_{u,n}$ are positive sequences that, as $n$ increases, approach 0 and 1 respectively.
The non-parametric estimator for $\Lambda(\cdot\mid G=co)$ will be constructed using martingale ideas from Theorem 1.3.1 by fleming2011counting. Firstly, let $I_i(y)=\mathbbm{1}(Y_i \leq y, \Delta_{1,i}=1)$ and $\tilde{I}_i(y)=\mathbbm{1}(Y_i \geq y)$. Following a similar martingale construction as Rivest2001AMA, we have that $$ \mathbbm{M}_i(y)=I_i(y) - \int^y_0 \tilde{I}_i(u) \lambda^{\#} (u\mid Z_i,X_i,co)\mathop{}\!d u, $$ are right-continuous martingales with respect to the $ \sigma$-algebra $$\mathcal{F}^{co,i}_y=\sigma\left\{I_i(u), \tilde{I}_i(u), Z_i,X_i, G = co: 0 < u < y \leq \bar{\tau} \right\},$$ where $\bar{\tau}$ is the maximum follow-up time defined in Assumption (ref) and $\lambda^{\#} (u\mid z,x,co)$ the conditional crude hazard rate, that is, $$ \lambda^{\#} (y\mid z,x,co) = \frac{-\frac{\partial}{\partial u}\mathbbm{P}(T \geq u, C \geq y \mid Z = z, X = x, G = co)\mid_{u=y}}{\mathbbm{P}(T \geq y, C \geq y \mid Z = z, X = x, G = co)}, $$ for all $(z,x^{\top})$. Using the model specification from Section (ref), it can easily be shown that
with
and where $\theta^* = (\mu^*,\eta^*,\xi^*), $ $\Lambda^*$ and $\lambda^*$ are the true values of $\theta,\Lambda$ and $\lambda$ respectively. If we let $I_i(y-)=\lim_{u \text{ } \uparrow \text{ } y}I_i(u)$, it follows that $\mathop{}\!d I_i(y) = I_i(y)-I_i(y-)$ is a binary random variable with conditional probability $\tilde{I}_i(y)\exp\left(\psi_i(\theta^*,\Lambda^*(y \mid G=co))\right)\lambda^*(y \mid G=co)$ of being one given $\mathcal{F}^{co,i}_{y}$. Motivated by $\mathbbm{M}_i(y)$ being a right-continuous martingale, it follows that
with $0 < y \leq \bar{\tau}$ and $\Lambda^*(0\mid G=co)=0.$ However, this equation cannot be used for estimation since it depends on a latent variable $G$. Using Theorem 3.1 from ABADIE2003231 that we described before, we can rewrite this as $$ \operatorname{\mathbb{E}} \bigg[ \kappa^*(S_i) \times \Big\{ \mathop{}\!d I_i(y)-\tilde{I}_i(y)\exp\left(\psi_i\left(\theta^*,\Lambda^*(y \mid G=co)\right)\right)\mathop{}\!d \Lambda^*(y \mid G = co) \Big\} \mid \mathcal{F}^i_{y} \bigg] = 0, $$ with $\kappa^*$ the true value of $\kappa$ and $\mathcal{F}^i_{y}=\sigma\left\{I_i(u), \tilde{I}_i(u), Z_i,X_i: 0 < u < y \leq \bar{\tau}\right\}.$ In this way, the expectation no longer depends on the unobserved variable $G$. Therefore, for a given $\theta$, we could estimate $\Lambda^*(y \mid G=co)$ by solving the following functional estimating equation:
with $0 < y \leq \Bar{\tau}$ and $\Lambda(0\mid G=co)=0.$ It is clear that $\mathop{}\!d I_{i}(y)$ will be equal to 1 at the observed event times $(Y_i, \Delta_{1,i}=1)$ and 0 everywhere else. Combining this with the fact that $\Tilde{\kappa}(S_i)$ is bounded between $0$ and $1$, it is clear that the estimator of $\Lambda^*$ following from this estimating equation is a non-decreasing step function where the jumps are at the ordered observed survival times: $0 = t_0 < t_1 < t_2 < \dots < t_K \leq \bar{\tau}$. However, using this estimating equation would involve a complex iterative optimization process that solves a $K$-dimensional system of equations, that is $$ \hat{\Lambda}_{\theta}(t_k \mid G = co)=\hat{\Lambda}_{\theta}(t_{k-1} \mid G = co)+\frac{\sum^n_{i=1} \left[ \tilde{\kappa}(S_i) \times \mathop{}\!d I_{i}(t_k) \right]}{\sum^n_{i=1} \left[\tilde{\kappa}(S_i) \times \tilde{I}_{i}(t_k)\exp(\psi_i(\theta,\hat{\Lambda}_{\theta}(t_k\mid G = co)))\right]}, $$ for $k=1,\dots,K$ and with $\hat{\Lambda}_{\theta}(t_0 \mid G = co)=0$. Following ZuckerDavidM2005APLM, we propose estimating $\Lambda^*$ as a step function with jumps at the ordered observed event times in the following way: $$ \Delta \hat{\Lambda}_{\theta}(t_k \mid G = co)=\frac{\sum^n_{i=1} \left[ \tilde{\kappa}(S_i) \times \mathop{}\!d I_{i}(t_k) \right]}{\sum^n_{i=1} \left[\tilde{\kappa}(S_i) \times \tilde{I}_{i}(t_k)\exp(\psi_i(\theta,\hat{\Lambda}_{\theta}(t_{k-1}\mid G = co)))\right]}, $$ with $\lim_{y \text{ } \uparrow \text{ } t_1}\hat{\Lambda}_{\theta}(y \mid G = co)=0$ and $\Delta \hat{\Lambda}_{\theta}(t_k \mid G = co) = \hat{\Lambda}_{\theta}(t_k \mid G = co) - \hat{\Lambda}_{\theta}(t_{k-1} \mid G= co)$. Therefore, the estimates $\hat{\Lambda}_{\theta}(t_1 \mid G = co),\dots,\hat{\Lambda}_{\theta}(t_K \mid G = co)$ can be obtained by a simple forward recursion. When the independence copula is specified, it is important to note that $\hat{\Lambda}_{\theta}$ cancels out from the denominator. This results in the estimator simplifying to a weighted version of the breslow1974covariance estimator. Moreover, the presence of $\hat{\Lambda}_{\theta}$ in the denominator when the specified copula deviates from the independence copula and the presence of $\tilde{\kappa}$ poses significant challenges in deriving the asymptotic properties of the estimator.
We set the expected value (conditional on $G=co$) of the first-order conditions (with respect to $\theta$) of the conditional likelihood evaluated at the true $\theta$ and $\Lambda$ to zero to arrive at the following equations:
where $$ U\big(S_i ,\theta^*,\Lambda^*\big) = \frac{\partial}{\partial \theta} \log f_{Y,\Delta_1,\Delta_{2} \mid Z,X,G,\theta^*,\Lambda^*}(Y_i,\Delta_{1,i},\Delta_{2,i} \mid Z_i,X_i,co), $$ Again, we can use Theorem 3.1 from ABADIE2003231 to rewrite equation this in the following way:
such that the score functions no longer depend on the unobserved variable $G$. Therefore, let $$M_n\big(\kappa,\theta,\Lambda\big) = n^{-1}\sum_{i=1}^n m\big(S_i,\kappa,\theta,\Lambda\big),$$ and $$m : \mathcal{S} \times\mathcal{K}\times \Theta \times \mathcal{L} \to \mathbb{R}^{\text{dim}(\Theta)} : \big(S,\kappa,\theta,\Lambda\big) \mapsto \kappa(S) \times U\big(S,\theta,\Lambda\big),$$ a measurable vector-valued function. Note that $\Theta$ is the parameter space for $\theta$ and $\mathcal{K}, \mathcal{L}$ are function spaces that will be defined in Section (ref). Further, let $M(\kappa, \theta,\Lambda) = \operatorname{\mathbb{E}}[m(S,\kappa, \theta,\Lambda)] $. Using what we have derived so far, we have the following estimating equations:
where $\hat{\theta}$ is defined as the solution to these weighted score equations. Note that this $Z$-estimator is only used for the theory, as in practice we will maximize the weighted logarithm of $L(\theta,\Lambda)$ with $\tilde{\kappa}(S_i)$ as the weights for each $i$. We do this by randomly generating $J$ starting values for $\theta$, denoted by $\tilde{\theta}_j$ with $j=1,\dots,J$, from a specified parameter space and calculating $\hat{\Lambda}_{\tilde{\theta}_j}$ for each $j$. The weighted log-likelihood is then maximized with respect to $\theta$ for each $j$, given $\hat{\Lambda}_{\tilde{\theta}_j}$, such that we have a set of estimates $\{\hat{\theta}_j\}_{j=1,\dots,J}$. Next, we continue with the $\hat{\theta}$ that has the highest weighted log-likelihood, and iterate calculating $\hat{\Lambda}_{\hat{\theta}}$ and maximizing the weighted log-likelihood until a convergence criterion is met or a maximum amount of iterations is reached.
In this section, it will be shown that the finite-dimensional parameter estimates $\hat{\theta}$ resulting from the $Z$-estimator defined in Section (ref) are consistent and asymptotically normal. The difficulty in establishing the asymptotic theory comes from our estimating equations involving both finite and infinite-dimensional parameters. This is complicated further by the fact that $\hat{\Lambda}_{\theta}$ depends on $\tilde{\kappa}$. Firstly, let
In addition to the previously made assumptions, we need the following regularity conditions to establish the asymptotic properties of the estimator:
Condition (ref) implies that the support of $X$, denoted by $\mathcal{X} \subset \mathbb{R}^m$ is bounded. Moreover, the densities of $Y$ or $X$ are also bounded and positive. Note that Assumption (ref) and Conditions (ref) and (ref) imply that $\exp\big(\psi(\theta,c)\big)$ is bounded on the interval $[\Psi_{\min},\Psi_{\max}]$ over $\theta \in \Theta, s \in \mathcal{S} $ and for all $ c \in \mathbb{R}_{\geq 0}$ and that Condition (ref) implies that $\kappa^*(S)$ and $\pi^*(X)$ are bounded away from 0 and 1 almost surely. Moreover, Condition (ref) implies that $h_{1,n} = o(n^{-1/2b} \wedge (\log n)^{\nicefrac{1}{(m_c+2d)}}n^{\nicefrac{-1}{(2m_c+4d)}})$ and $ h_{2,n} = o(n^{-1/2b} \wedge (\log n)^{\nicefrac{1}{(m_c+1+2d)}}n^{\nicefrac{-1}{(2m_c+2+4d)}})$. Further, we define $\mathcal{K}$ to be the class of all functions $\kappa(\cdot)$, defined over $\mathcal{S}$, which are bounded between 0 and 1 and that are at least $d$-th order continuously differentiable, where $d$ is defined by Condition (ref). Further, let $\mathcal{L}$ be the class of all functions $\Lambda(\cdot)$, defined over $[0,\bar{\tau}]$, for which $\Lambda(0)=0$, $\Lambda(\cdot)$ is non-decreasing, and $\Lambda(\bar{\tau}) < \Lambda_{\max}$, where $\Lambda_{\max}$ is defined by Lemma 3 in Section C of the Supplementary Material. We are now ready to state the following theorems, of which the proof can be found in Section D of the Supplementary Material.
To show the asymptotic properties of $\hat{\theta}$, we use results by chen2003estimation. However, to prove Theorem (ref), we need to extend their results to the case where the criterion function depends on two unknown functions, since $M_n$ depends on both $\tilde{\kappa}$ and $\hat{\Lambda}_{\theta}$. Checking the conditions of this extended version of the theorem is complicated further by $\hat{\Lambda}_{\theta}$ depending on $\Tilde{\kappa}$. Note that $\Omega$ has a very lengthy expression that, since we estimate $\Sigma$ through a naive bootstrap method, is of little use and therefore omitted.
In this section, various simulation studies are performed to investigate the finite sample performance of the proposed estimator. Firstly, we look at the performance of our estimator using multiple combinations of copulas and marginals. The proposed estimator is compared to two other estimators: one that does not account for unobserved heterogeneity (naive estimator) and one that uses the proposed method when $G$ is observed (oracle estimator). Secondly, we assess the performance of the proposed estimator under misspecification of the censoring distribution or the copula model. Lastly, we look at what happens to our estimates when we change either the proportion of compliers or the sample size.
The first steps of the data-generating process, which are the same for all designs, are as follows: $$ X = (X_1,X_2)^\top \text{ where } X_1 \sim \text{Bernoulli}(0.5) \text{ and } X_2 \sim U(0,1). $$ Further, we generate $G$ from a multinomial distribution with $\mathbb{P}(G=co)=2/3$ and $\mathbb{P}(G=at) = \mathbb{P}(G=nt)=1/6.$ Moreover, $W$ is generated from a Bernoulli$(\pi(X))$ distribution, where $$ \pi(X) = \frac{\exp\left(0.5X_1 + X_2 + 2X_1X_2 + \epsilon\right)}{1+\exp\left(0.5X_1 + X_2 + 2X_1X_2 + \epsilon\right)}, $$ with $\epsilon \sim \mathcal{N}(0,0.25^2).$ Note that by adding this $\epsilon$, a logistic regression model for $\pi(X)$ would be misspecified. It follows that we can determine $Z$ by using $G$ and $W$. The next step depends on the choice of the copula and censoring distribution. As an example, we consider the following design (Frank - Weibull): $$ \mathbb{P}(T \leq t, C\leq c\mid Z=z,X=x,G=g)=\mathcal{C}_{\xi_g}\left(F_{T \mid Z,X,G}(t\mid z,x,g),F_{C \mid Z,X,G}(c\mid z,x,g)\right),$$ with $\mathcal{C}_{\xi_g}$ a Frank copula, that is, $$ \mathcal{C}_{\xi_g}(u,v) = -\frac{1}{\xi_g}\log\left(1+\frac{\left(\exp(-{\xi_g} u)-1\right)\left(\exp(-{\xi_g} v)-1\right) }{\exp(-{\xi_g})-1}\right), \quad {\xi_g} \neq 0. $$ Furthermore, we specify the distribution of $T$ and $C$ by a proportional hazards and Weibull model respectively, that is, $$ F_{T \mid Z,X,G}(t\mid z,x,g)=1-\exp\{-\Lambda(t\mid G=g) \exp( z \alpha_g+x^\top \beta_g)\}, $$ $$ F_{C \mid Z,X,G}(c\mid z,x,g) = 1 - \exp\left(-\exp\left(\frac{\log(c) - \tilde{x}^{\top}\eta_g}{\nu_g}\right)\right), $$ where $\tilde{x} = (1,z,x^{\top})^{\top}$, $\beta_g^\top = (\beta_{1,g},\beta_{2,g})$ and $\eta_{g}^{\top} = (\eta_{0,g},\eta_{1,g},\eta_{2,g},\eta_{3,g},\nu_{g})$. Note that specifying the association parameter $(\xi_g)$ or Kendall's tau $(\tau_g)$, is equivalent because of the following relation: $$ \tau_g = 4\int_0^1\int_0^1\mathcal{C}_{\xi_g}(u,v)\textit{c}_{\xi_g}(u,v)dudv-1, $$ which is one-to-one for the (rotated) Archimedean and Gaussian copulas. This approach maintains the same model structure across all latent subgroups while allowing the parameter values to vary. Examples of this can be seen in Tables (ref) and (ref), which represent a low and high dependence scenario respectively.
To simulate the administrative censoring, let $A \sim U(0,15)$ for the low dependence scenario and let $A \sim U(0,50)$ for the high dependence scenario. Finally, we can generate the follow-up time $Y$ = min\{$T,C,A$\}, $\Delta_1 = \mathbbm{1}(Y=T)$ and $\Delta_2 = \mathbbm{1}(Y=C)$.
The data-generating process was repeated 500 times with a sample size of 1000 for each of the six simulation designs considered. For each of the designs, we compare our proposed estimator to two others. The naive estimator ignores the endogeneity issue such that everyone is considered to be a complier ($\tilde{\kappa}(S_i) = 1$ for each $i$). The oracle estimator uses the same estimation procedure as the proposed estimator but treats $G$ as if it were observed ($\tilde{\kappa}(S_i) = \mathbbm{1}(G_i=co)$ for each $i$). Therefore, our proposed method reduces to a one-step estimation procedure as we do not need to estimate $\kappa^*$. We report the bias, empirical standard deviation (ESD), root mean squared error (RMSE) and coverage rate (CR). Note that the CR indicates the percentage of simulations in which the true parameter value falls within the estimated 95% confidence interval. These coverage rates were computed using the warp-speed method described by warpspeed. This method calculates bootstrap confidence intervals by drawing only one bootstrap resample per Monte Carlo sample, significantly reducing computation time. However, this approach may result in less accurate coverage rates due to the reduced number of bootstrap resamples. To estimate the probability of being a complier, we modified the {\fontfamily{cmtt}\selectfont point.est.kernel} function from the {\fontfamily{cmtt}\selectfont CCQTEC} package wei2021estimation in {\fontfamily{cmtt}\selectfont R} to fit our setting. We use a multiplicative sixth-order Epanechnikov kernel and select the bandwidths from the range $\{0.01, 0.02,\cdots, 1\}$ using 10-fold cross-validation, with $a_{l,n} = 10n^{-1}$ and $a_{u,n} = 1-10n^{-1}$. In the second step, to estimate the parameters of interest, we maximize the weighted log-likelihood as explained in Section (ref) with $J=100$ and a maximum amount of iterations equal to 120.
The results for the scenario with low dependence can be found in Table (ref). Each of these three designs has around 30 to 40 percent dependent and 5 to 10 percent administrative censoring. Note that the design in the middle column is based on a data-generating process with negative dependence. The results show that our method has low bias across all designs, especially compared to the naive estimator where the bias is the highest for the main parameter of interest $\alpha$. Compared to the naive estimator, the proposed estimator exhibits a substantially lower RMSE for the parameter of interest $\alpha$. However, its RMSE is higher for the other parameters. Nonetheless, the coverage rates of the naive estimator frequently show substantial deviations from the nominal 95% level. In contrast, the coverage rates of the proposed estimator exhibit only minor deviations from the nominal level, likely due to the use of the warp-speed method. The results for the scenario with high dependence are similar and can be found in Section E of the Supplementary Material.
In this subsection, we evaluate the performance of the proposed method under misspecification of either the copula or the distribution of $C$. For each design, 500 datasets with a sample size of 1000 were generated. The parameter values used were those from the previously described low-dependence scenario. The Frank copula was used for all designs involving misspecification of the censoring distribution. When the copula was misspecified, a log-normal distribution was used for the censoring distribution. Detailed results on the impact of misspecification of the censoring distribution and copula are provided in Section E of the Supplementary Material. Overall, misspecification of the copula has a smaller effect on the bias and coverage rates compared to misspecification of the censoring distribution. However, this observation does not hold when the copula fails to model the correct direction of the dependence, as seen when data were generated using a Clayton(180) copula but a Clayton(90) copula was specified. The bias of the main parameter of interest, $\alpha$, remains relatively low across all designs, particularly compared to the naive estimator. For most designs, the bias and RMSE of $\alpha$ for the proposed estimator are comparable to those of the oracle estimator.
Finally, we examine the impact of the proportion of compliers and sample size on the bias and RMSE of the parameter of interest $\alpha$. To investigate the effect of the proportion of compliers, we generate 250 datasets for 5 different complier ratios $(1/10, 1/5, 1/3, 2/3, 1)$ with a sample size of 1000 using a Clayton(90) copula and a log-normal censoring distribution. For the effect of the sample size we also generated 250 datasets for each of the 5 sample sizes $(100,250,500,1000,1500)$ with a complier ratio of 2/3 and using the same copula and censoring distribution. The parameter values used to generate the data sets were those from the previously described low-dependence scenario. The results for the effect of the complier ratio and the sample size can be found in Figures (ref) and (ref) respectively. Our findings indicate that, regarding the bias, the proposed method is fairly robust to a low proportion of compliers. Even with only 10% compliers, the proposed estimator exhibits a bias that is comparable to the oracle estimator. Moreover, the proposed estimator substantially outperforms the naive estimator in terms of both bias and RMSE regardless of the complier ratio. As the proportion of compliers increases to 100% the bias and RMSE of all estimators converge, as would be expected. Regarding the effect of sample size, our proposed estimator maintains good performance with respect to bias, even when the sample size is as small as 100. In contrast, the bias of the naive estimator remains constant as the sample size increases. Furthermore, the RMSE of the naive estimator appears to plateau, showing little improvement with increasing sample size. Overall, these results suggest that our method performs effectively even in scenarios with a low proportion of compliers or small sample sizes, demonstrating its robustness under various conditions.
In this section, we present an empirical application of the proposed methodology to data from the Job Training Partnership Act (JTPA) study. A second data application uses data from the Health Insurance Plan of Greater New York (HIP) experiment and can be found in Section F of the Supplementary Material.
The data examined come from a large-scale randomized experiment known as the National Job Training Partnership Act (JTPA) Study and have been analyzed extensively by bloom1997benefits, abadie2002instrumental, frandsen2015treatment, wuthrich2020comparison, beyhum2024instrumental and crommen2024instrumental among others. The data that we investigate is the same as in abadie2002instrumental and wuthrich2020comparison, but the problem of interest differs. This is because we investigate the effect of the job training services on unemployment duration, while they look at the effect of JTPA services on the sum of earnings in the 30-month period after treatment assignment. The problem investigated is more similar to frandsen2015treatment and beyhum2024instrumental, but again differs as we allow for dependent censoring. Lastly, crommen2024instrumental analyzed the data under possible dependent censoring. However, their goal was to estimate a population treatment effect, while we focused on estimating a local treatment effect for the subgroup of compliers, which allows for more flexible modeling.
The JTPA study was conducted to assess the effectiveness of more than 600 federally funded programs established under the Job Training Partnership Act of 1982, aimed at improving the employability of eligible adults and out-of-school youths. These programs provided services such as classroom training, on-the-job training, and job search assistance. Funding for these initiatives began in October 1983 and continued into the late 1990s. Between 1987 and 1989, over 20000 adults and out-of-school youths who applied for JTPA were randomly assigned to either a treatment group, which was eligible for JTPA services, or a control group, which was not eligible for 18 months. However, due to some local program staff not strictly adhering to the randomization guidelines, about 3% of the control group members received JTPA services. It is important to note that we are comparing JTPA services to a mix of no services and other available services, as control group members could still access non-JTPA training. Participants were surveyed by data collection officers between 12 and 36 months after randomization, with an average survey period of 21 months. A second follow-up survey, involving a subset of 5468 participants, focused on the period between the two surveys and was conducted between 23 and 48 months after randomization. Figure 1 in Section G of the Supplementary Material plots a histogram of the observed follow-up time, where a higher censoring rate is indicated by darker shading.
We will focus our attention on the effect of JTPA job training programs on the sample of 3147 single mothers who reported being unemployed at the time of randomization and participated in the initial follow-up interview. The outcome of interest ($T$) is the time between randomization and employment. For participants who only participated in the initial interview, the outcome is fully observed if the individual is employed at the time of the survey $(\Delta_1=1)$. Otherwise, the outcome is censored at the date of the initial interview $(\Delta_2=1)$. For those who participated in the second follow-up interview, the outcome is fully observed if the individual is employed at the time of this second follow-up interview $(\Delta_1=1)$. If not, the outcome is considered to be independently censored at the second interview date $(\Delta_1=\Delta_2=0)$. For a graphical representation of the interview process, see Figure 2 in Section G of the Supplementary Material. Consequently, there may be a dependence between the time until employment and the censoring time if the decision to attend the second follow-up interview is influenced by the individual's employment status between the two interview dates. Note that this means that the individuals who were only invited to the first interview, and were still unemployed by this time, will be censored $(\Delta_1=1)$ instead of administratively censored $(\Delta_1=\Delta_2=0)$. This results from only observing who participated in the second interview, but not who was invited. Because we cannot know which type of censoring occurred, we opted to have them all be classified as censored $(\Delta_1=1)$.
The natural instrumental variable $W$ indicates whether an individual belongs to the control group ($W=0$) or the treatment group ($W=1$). We deem $W$ to be a valid instrument, as it is randomly assigned, moderately correlated with JTPA participation and influences the time to employment only through participation in a JTPA-funded program. The treatment variable $Z$, which is possibly endogenous, indicates whether the individual actually participated in a JTPA program ($Z=0$ for no participation and $Z=1$ for participation). This treatment variable can be confounded due to individuals moving themselves between the treatment and control groups in a non-random way. The covariates considered include the participant’s age (standardized), educational attainment (high school diploma or GED) and race (categorized as white or non-white). Approximately 31% of the total sample was assigned to the control group. The mean age of this subgroup of single mothers is approximately 28 years old, with 52% holding a GED or high school diploma and 47% identifying as white. Notably, 13% of the women in the control group managed to participate in JTPA services, in contrast to 3% of the entire control group. From the treatment group, around 70% participated in JTPA services. Additionally, the mean time to employment appears to be about 43 days shorter for individuals who participated in JTPA training. The censoring and independent censoring rates are comparable between both the control and treatment groups (20% censoring and 5% administrative censoring).
We start by splitting the data into two random samples of equal size. The first sample was used to select the specification for the copula and the distribution of $C$. Out of the possible 21 combinations, we select the model with the highest log-likelihood. We can simply use the log-likelihood since each combination of possible copulas and distributions has the same amount of parameters. The results of this selection procedure can be found in Section G of the Supplementary Material for both the proposed and naive estimator. It can be seen that the Frank copula with a log-logistic censoring distribution results in the highest log-likelihood for both of the estimators. The estimation results for both of these methods can be found in Tables (ref) and (ref). In these tables, we also include the estimation results of the model with the second-highest log-likelihood as a type of sensitivity analysis. Note that the $p$-value associated with the null hypothesis $\alpha = \alpha^*$, is calculated as $$ \frac{1}{B}\sum_{b=1}^B \mathbbm{1}\{\lvert \hat{\alpha}_{b}-\hat{\alpha} \rvert > \lvert \hat{\alpha} - \alpha^* \rvert \}, $$ where $B$ is the amount of bootstrap resamples and $\hat{\alpha}_{b}$ the bootstrap estimate based on the bootstrapped sample $b \in \{1,\dots,B\}$. The $p$-values for the other parameters can be computed similarly. In Tables (ref) and (ref), we test the null hypothesis that the true parameter equals zero, except for $\nu$ where we test the null hypothesis that $\nu = 1$.
The naive estimator appears to underestimate the effect of JTPA services on time until employment compared to the proposed estimator. At the 5% significance level, only the proposed estimator finds a significant effect of JTPA training in reducing unemployment duration, with the CCHR being 1.23. This finding suggests that the individuals participating in the treatment are likely to have a lower inherent ability to secure employment. For both estimators, variables such as age and having a high school diploma or GED do not significantly affect the duration of unemployment. Notably, being white is associated with a significant reduction in unemployment duration for the proposed estimator. Both the proposed and the naive estimator indicate that there is a strong positive dependence between $T$ and $C$, conditional on the treatment and the measured covariates.