Instrumental variable estimation of dynamic treatment effects on a duration outcome

Instrumental variable estimation of dynamic treatment effects on a duration outcome

\if11 \fi

\if01 {

} \fi

{\it Keywords:} Dynamic treatment; Endogeneity; Instrumental variable; Censoring; Nonseparability. \\

\spacingset{1.4}

\setcounter{footnote}{0} \setcounter{equation}{0}

Introduction

Consider subjects in a certain state, such as unemployment or medical leave. The policymaker is sometimes interested in how the timing of some treatment affects the time the subjects remain in the state. An example comes from labor economics, where we may want to evaluate the effect of reducing unemployment benefits after some time on the unemployment duration. Also, in health economics, one may seek to find the optimal time to give therapy for workers on medical leave because of burnout. This setting corresponds to our empirical application.

{\color{black} The present paper considers the problem of estimating the causal effect of a treatment. The treatment is dynamic because it can be given at any point in time, and its effect can vary depending on its timing. Moreover,} the time $Z$ until the treatment is given is not randomly assigned. We solve the endogeneity issue thanks to an instrumental variable $W$ independent of the error term of the model and sufficiently related to $Z$. The outcome duration of interest $T$ is randomly right censored by a censoring time $C${\color{black}, assumed independent of the other variables of the model.} The timing of the treatment $Z$ is itself censored by $\min(T,C)$. {\color{black} The censoring of $Z$ by $T$ is endogenous, since the latter variables are dependent.} This context corresponds to studies where follow-up stops when the subjects leave the state of interest, or the treatment cannot be given to participants who leave the state of inflow. In the latter case, $Z$ corresponds to some latent duration to treatment, which is realized only when $Z\le T$. In the labor economics example mentioned above, this is justified because unemployment benefits cannot be reduced if the subject is no longer unemployed. In the illustration from health economics, there is no reason to treat cured workers.

The contributions of the paper are as follows. \textcolor{black}{We study a dynamic duration model where we are interested in the hazard rate of $T(z)$, where $T(z)$ is the potential outcome of the duration when the treatment happens at time $z$. We make a number of assumptions on this duration model. One assumption, called no anticipation, allows us to handle the issue that we never observe treatment times $Z$ which are larger than $T$ since $Z$ is right censored by $\min(T,C)$. A second assumption is a rank invariance condition common in the nonparametric instrumental variable (NPIV) literature. This assumption allows us to identify and estimate the hazard rate of $T(z)$ using an instrumental variable $W$ independent of the error term of the model. Specifically, we rewrite the duration model as a nonseparable NPIV model and adapt tools from the NPIV literature to study the identification of the regression function of the rewritten model in a fully nonparametric framework. The hazard rates of the potential outcomes are functionals of the regression function of the nonseparable model. Our identification results are not straightforward applications of the results from the nonseparable NPIV model literature since our nonseparable model is different from their standard quantile model because of the dynamic nature of the problem. For estimation purposes, we assume that the regression function follows a parametric model. We propose an estimation procedure and give conditions under which the semiparametric estimator is asymptotically normal. The finite sample properties of the estimator are assessed through simulations. We apply our methodology to evaluate the effect of the timing of therapy on the duration of medical leave for burnout.}

There exists an extensive literature on instrumental variable methods with randomly right-censored duration outcomes, where both the treatment and the instrument are time-independent. We only cite here some references to avoid lengthening the paper. Many works study semiparametric models, see, e.g., tchetgen2015instrumental for an additive hazard model, chernozhukov2015quantile for a quantile regression model, or martinussen2019instrumental for the Cox model. Other works (frandsen2015treatment, sant2016program, richardson2017nonparametric, blanco2019bounds, sant2021nonparametric) provide nonparametric estimation results for the average treatment effects on the compliers (see angrist1996identification) with binary treatment and instrument. Moreover, 2021nonparametric nonparametrically estimate the average quantile treatment effect over the whole population when both the treatment and the instrument are categorical. Also, centorrinoflorens2021 study nonparametric estimation when the treatment is continuous, and the model is additive. They allow the instrument to be categorical, even in this case. Our paper differs from this literature because it allows the treatment to be dynamic. Our dynamic setting is more complex than those previously analyzed since the treatment $Z$ is not always observed and endogenously censored by $\min(T,C)$. Among the papers above, 2021nonparametric is the most closely related to our paper since it relies on a nonseparable NPIV model as in CH. \textcolor{black}{There are two main differences with this paper. First, here $Z$ is dynamic, which creates endogenous censoring of $Z$ by $T$. Second, in the present context, $Z$ is a continuous variable. This makes the inverse problem ill-posed and requires a different identification analysis and a new estimation procedure.}

{\color{black}Some papers study the evaluation of the effect of the time to treatment on a survival outcome when the treatment is endogenous without relying on an instrumental variable. The first set of papers assumes conditional unconfoundedness of the treatment. See sianesi2004evaluation, AVdB2,lechner2009sequential, lechner2010identification, vikstrom2017dynamic and van2020policy, kastoryano2022dynamic for the literature in econometrics. hernan2010causal give an excellent overview of the methods available in the biostatistics literature to infer causal effects in this context. These studies do not rely on an instrumental variable to solve the endogeneity problem. Another approach is that of AVdB, where it is assumed that both the timing of the treatment and the duration outcome follow a mixed proportional hazards model.

Finally, like us, some papers use instrumental variables to estimate the effect of an endogenous time-varying treatment on a right-censored outcome. bijwaard2005correcting considers a setting where the instrument is binary, control group subjects are never treated and a mixed proportional hazard model is assumed. heckman2007 study a dynamic treatment effect model for outcomes that can be durations. However, they do not allow for censoring of $T$ by $C$ and of $Z$ by $T$, and they rely on a control function approach requiring additional structural assumptions, and no estimation theory is proposed. Another paper is VdBBM, which studies a one-sided noncompliance setting, where the time to treatment is either equal to the instrument or $\infty$ (untreated). {\color{black} They identify local average treatment effects on a subset of the population analogous to the compliers in static problems (see angrist1996identification). In contrast, the present work does not assume one-sided noncompliance, and our method allows us to estimate effects over the whole population. VdBBM also applies their estimation strategy to evaluate the causal effect of a reform of the unemployment insurance system in France. Our approach could also be applied to their dataset. Under our assumptions, we would identify the effect of the reform on the whole population rather than only part of it. Finally, note a recent line of work in biostatistics (tchetgen2018marginal, cui2017instrumental, michael2020instrumental) which studies the case where the instrument varies with time, ruling out the case with a static instrument studied in our paper. They have a local (time by time) estimation approach leveraging the time variation of the instrument. This estimation strategy can not be straightforwardly adapted to allow for a time-independent instrument.}}

{\color{black}Outline. }This paper is organized as follows. The model specification is given in Section 2. We study identification in Section 3. Section 4 is devoted to estimation and inference. Section 5 describes our simulations and the empirical application. Concluding remarks are given in Section 6. The proofs of the results and additional numerical experiments are in the online appendix.

The model

The duration model

\textcolor{black}{Let $T(z)$ be the potential outcome of the duration when the treatment time is set to $z\in\bar{{\mathbb{R}}}_+={\mathbb{R}}_+\cup\{\infty\}$. When $z=\infty$, $T(z)=T(\infty)$ corresponds to the duration that would have been realized if the subject of interest were never treated. The treatment time is a random variable $Z$, with support $\mathcal{Z}\subset{\mathbb{R}}_+$. We impose the consistency condition $T=T(Z)$. In our empirical application, studied in Section (ref), we want to evaluate the effect of a therapy for burnout on the duration of medical leave. We possess a dataset where each observation corresponds to a worker on medical leave for burnout. The variable $Z$ is the time until a therapy for burnout starts, $T(z)$ is the duration of medical leave that would have been realized if the therapy for burnout had been started at time $z$ and $T$ is the actual duration of medical leave. }

\textcolor{black}{We assume that $T(z)$ is a continuous random variable and let $\lambda(z,\cdot):{\mathbb{R}}_+\mapsto{\mathbb{R}}_+$ be its hazard rate, that is $$\lambda(z,t)= \lim_{dt \to 0}\frac{{\mathbb{P}}(T(z)\in[t,t+dt]|T(z)\ge t)}{dt}.$$ We call $\lambda$ the “structural hazard". Remark that $\lambda(z,\cdot)$ differs from the hazard rate of $T$ conditional on $Z=z$ in general (this is the endogeneity issue). The main goal of the paper is to identify and estimate this structural hazard. }

\textcolor{black}{As mentioned in the introduction, one of the problems considered in this paper is that $Z$ is censored by $T$. When $Z>T$, we never know to which potential outcome $T$ corresponds. To circumvent this issue, we impose the so-called no anticipation assumption (see AVdB, VdBBM, van2020policy), which is standard in the literature on dynamic treatment effects. The no anticipation assumption can be formally stated as follows.

This assumption means that the hazard rate $\lambda(z,t)$ of $T(z)$ does not depend on $z$ when $z>t$. Under Assumption (ref), all observations for which $Z>T$ correspond to the same structural hazard (equal to $\lambda(\infty,\cdot)$), which allows us to solve the aforementioned issue that we do not observe $Z$ when $Z>T$. Denote by $\Lambda(z,t)=\int_0^t \lambda(z,s)ds$ the structural cumulative hazard of $T(z)$. By integration and derivation, Assumption (ref) is also equivalent to $\Lambda(z,t)=\Lambda(z',t)$ when $t<\min(z,z')$, which is the way the no anticipation assumption is stated in AVdB (except that the strict inequality $t<\min(z,z')$ is replaced by the weak inequality $t\le\min(z,z')$ in AVdB). We illustrate the no anticipation assumption in Figure (ref), where we draw the structural hazard rates under treatment levels $z=5$ and $z=7$ (in the caption, we use $I(\cdot)$ to denote the indicator function). In this illustration, the treatment increases the hazard rate, but our model also allows the treatment to reduce the hazard rate. Before the treatment, the hazard rate does not depend on the treatment time, while after the treatment it does. Finally, note that the no anticipation assumption does not rule out that subjects actually anticipate receiving future treatment. It rather means that the counterfactuals that we are interested in correspond to a setting where the intervention does not change the anticipations but only the actual value of the treatment (see the discussion in AVdB for more details).}

\textcolor{black}{To be able to identify treatment effects over the full population under endogeneity, we impose constraints on the unobserved heterogeneity of the model. For $z\in\bar{{\mathbb{R}}}_+$ and $t\in{\mathbb{R}}_{+}$, recall that $\Lambda(z,t)=\int_0^t \lambda(z,s)ds$ is the structural cumulative hazard under treatment $z$ and let $U(z)=\Lambda(z,T(z))$ be the hazard of $T(z)$ evaluated at $T(z)$. The $\{U(z)\}_{z\in\bar{{\mathbb{R}}}_+}$ can be thought of as the unobserved heterogeneity of the model. In our empirical application, it could correspond to the underlying mental health of the subject. We impose the following conditions on them.

Condition (ii) allows the hazard rate $\lambda(z,\cdot)$ to be discontinuous at the time of treatment $z$. This permits behaviors similar to that of Figure (ref), where treatment makes the hazard rate “jump". We show in the online appendix (Lemma (ref)) that Assumption (ref) (ii) implies that the cumulative hazard $\Lambda(z,\cdot)$ maps the support of $T(z)$ to ${\mathbb{R}}_+$ and is strictly increasing on the support of $T(z)$. Therefore we can define $\Lambda(z,\cdot)^{-1}$, the inverse of the mapping $\Lambda(z,\cdot)$ restricted to the support of $T(z)$ (that is $\Lambda(z,\cdot)^{-1}(u)$ is the unique element $t$ of the support of $T(z)$ such that $\Lambda(z,t)=u$). This implies that $T(z)=\Lambda(z,\cdot)^{-1}(U(z))$ is strictly increasing in $U(z)$. Then, Assumption (ref) (i) yields that, for two subjects $i$ and $j$, $T_i(z)>T_j(z)$ implies $U_i(z) = U_i>U_j = U_j(z)$, which leads to $T_i(z')>T_j(z')$, for all $z,z'\in\bar{{\mathbb{R}}}_+$. In other words, Assumption (ref) (i) implies that the rank in the outcome of any two subjects is the same across all potential outcomes. Assumption (ref) (i) is therefore a rank invariance assumption as in CH. This assumption restricts the heterogeneity of the treatment effects on the duration: the treatment can change the quantiles of the distribution of the potential outcomes but it cannot change the rank that a subject has in this distribution. Moreover, remark that the rank invariance assumption does not restrict the possible values of the structural hazard $\lambda(z,t)$ beyond the continuity condition in Assumption (ref) (ii) and therefore, in this sense, rank invariance is not a constraint on the marginal distribution of $T(z)$. It only imposes limits on the joint distribution of potential outcomes, that is the distribution of $(T(z))_{z\in \mathcal{Z}}$. Note also that we could relax the rank invariance assumption into a rank similarity assumption as in CH while keeping all results valid.}

\textcolor{black}{Another important fact is that Assumption (ref) (ii) implies that for all $z\in\bar{{\mathbb{R}}}_+$, $U(z) \sim\text{Exp}(1)$, where $\text{Exp}(1)$ is the unit exponential distribution (see Lemma (ref)). In virtue of Assumption (ref) (i), we can therefore write $\Lambda(Z,T)=U,$ with $U\sim \text{Exp}(1)$. }

\textcolor{black}{The variables $U$ and $Z$ may be dependent, which creates an endogeneity issue, that is $\lambda(z,\cdot)$ may differ from the hazard rate of $T$ given $Z=z$. In the context of the empirical application, the endogeneity may, for instance, be due to the fact that subjects with worse burnout are more likely to be treated early. There exists an instrument $W$ allowing to solve this issue. In the real data, the instrument is related to the medical center to which workers on medical leave are assigned (see Section (ref) for more details). The support of $W$ is denoted $\mathcal{W}$. For simplicity, we limit ourselves to the case where $W$ is scalar, that is $\mathcal{W}\subset {\mathbb{R}}$. \textcolor{black}{We impose the following assumption:

} The duration $T$ is randomly right censored by a random variable $C$ with support in $\bar{{\mathbb{R}}}_+$, so that we do not observe $T$ but $Y=\min(T,C)$. In the application, $C$ corresponds to the duration during which subjects on medical leave are followed in the data (the follow-up stops after two years of medical leave or at the end of 2020). The observables are $(Y, \delta, \tilde Z, \tilde D, W)$, where $\delta =I(T\le C)$, $\tilde Z=\min(Z,Y)$, is a censored version of $Z$ and $\tilde D= I(Z\le Y)$ is a censored version of $D=I(Z\le T)$, the treatment indicator. Note that $D=1$ for treated observations only (in the sense that they receive treatment before the end of their spell). In the burnout data, we have $\delta =1$ if the duration of medical leave is observed and $0$ otherwise, $\tilde Z$ is the minimum between the treatment time, the censoring time and the duration of medical leave, $D$ is equal to $1$ if the subject is treated before the end of its medical leave and $0$ otherwise and $\widetilde{D}=1$ if we observe the treatment time in the dataset and $0$ otherwise.}

Reformulation as a nonseparable NPIV model

\textcolor{black}{As mentioned in the introduction, this paper makes use of tools from the NPIV model to solve the dynamic problem. Recall that for all $z\in\bar{{\mathbb{R}}}_+$, $\Lambda(z,\cdot)^{-1}(u)$ is the unique element $t$ of the support of $T(z)$ such that $\Lambda(z,t)=u$ and that $I(\cdot)$ is the indicator function. The reformulation as a nonseparable NPIV model is based on the following lemma. }

Equation (ref) means that there exists mappings $\varphi_0$ and $\varphi_1$ such that $T$ is equal to $\varphi_0(U)$ when $Z<\varphi_0(U)$ and equal to $\varphi_1(Z,U)$ otherwise. Intuitively, this signifies that there are two regressions functions: one of the “not yet treated" corresponding to $\varphi_0$, and one for the “already treated" which is $\varphi_1$. Writing the NPIV model as functions of these $\varphi_0$ and $\varphi_1$ allows us to define a parameter space which is a vector space (see Section (ref)). Note that, by definition, $\varphi_0$ and $\varphi_1(z,\cdot)$ are strictly increasing. Moreover, the fact that $\varphi_1(z,\varphi_0^{-1}(z))=z$ for all $z\in{\mathbb{R}}_+$ (see the end of Lemma (ref)) implies that $\varphi$ is strictly increasing too.

\textcolor{black}{Equation (ref) defines a nonseparable NPIV model similar to that of CH. Hence, we can use tools from the literature on nonseparable NPIV models to obtain identification results on $\varphi$. Throughout the paper, we will therefore identify and estimate $\varphi$ rather than $\lambda$ because this approach simplifies the mathematical analysis. Using Lemma (ref) and Assumption (ref), one can show that $\varphi(z,\cdot)=\Lambda(z,\cdot)^{-1}$ (see the proof of Lemma (ref) in the online appendix for more details). The structural hazard $\lambda(z,t)$ is therefore the derivative of the inverse of $\varphi(z,\cdot)$ at $t$. As a result, $\varphi$ is a one-to-one transformation of $\lambda$ and identification of $\varphi$ implies identification of $\lambda$ and therefore of many quantities of interest in duration models (such as the survival function or the cumulative hazard). Note that, the present model involves two additional complexities with respect to the original model of CH. First, $T$ is right censored by $C$. Second, $Z$ is right censored by $\min(T,C)$ and to (partially) solve this problem we have imposed a no anticipation assumption which implies that $\varphi$ follows the specific functional form of the right-hand side of (ref), which is different from the one used in the standard quantile model of CH. As we will see later, this implies that we have to adapt results of the nonseparable NPIV model to our nonstandard model.}

Identification

We study identification in a fully nonparametric setting. First, we show that our model generates a system of integral equations. Then we derive identification results based on this system of equations.

Identification equation

In order to formulate the identification equation, we introduce the following reduced form quantities. For $t\in{\mathbb{R}}_+ $, $w\in\mathcal{W}$, let

Moreover, for a mapping $\phi:\mathcal{Z}\mapsto {\mathbb{R}}_+$ and $w\in\mathcal{W}$, we define $$F_1(\phi,w)= {\mathbb{P}}(T\le \phi(Z), D=1,W\le w).$$ The following theorem states the system of equations that we use to obtain identification results.

\textcolor{black}{In Section (ref) of the online Appendix, we provide an alternative characterization of the model in terms of reduced-form (conditional) hazard rates and survival functions of the observed duration $T$, which are natural quantities in the context of duration models. }

Identification without censoring

In this section, we discuss the simple case where there is no censoring, that is $C=\infty$ a.s. Then, $Y=T$ and $\tilde D= D$, which implies that $F_0(t,w)$ is identified for all $t\in{\mathbb{R}}_+$ and $w\in\mathcal{W}$. Moreover, when $D=1$, we have $\tilde Z=Z$, and, hence, $F_1(\phi,w)$ is identified for all $\phi:\mathcal{Z}\mapsto {\mathbb{R}}_+$ and $w\in\mathcal{W}$. \textcolor{black}{Therefore, $F_0,F_1,F_W$ in (ref) are all identified in the absence of censoring.} Hence, uniqueness of the solutions to (ref) implies identification. This allows us to derive identification results in the next two subsubsections.

\textcolor{black}{We focus on identification of $\varphi(\cdot,u)$ for a given $u\in{\mathbb{R}}_+$. At this point, it is useful to define the set to which $(\varphi_0(u),\varphi_1(\cdot,u))$ belongs. Let the parameter space $\mathcal{P}$ be the set of $(\psi_0,\psi_1)$ such that $\psi_0\in {\mathbb{R}}$, and $\psi_1$ is a bounded mapping from $\mathcal{Z}$ to ${\mathbb{R}}$. This set $\mathcal{P}$ is a vector space and we endow it with the norm $\|\cdot\|_{\mathcal{P}}$, where $\left\|(\psi_0,\psi_1)\right\|_{\mathcal{P}}^2= E[(\psi(Z))^2|U=u],$ with $\psi(z)=\psi_0I(z> \psi_0)+\psi_1(z)I(z\le\psi_0).$ The function $\psi$ is the mapping "induced" by $(\psi_0,\psi_1)$. It is similar to the object that we want to identify ($\varphi(\cdot,u)$). The norm $\left\|(\psi_0,\psi_1)\right\|_{\mathcal{P}}$ is finite because $\psi_1$ is bounded. We assume that $\psi_1$ is bounded, because, in practice, to identify $\varphi_1(\cdot,u)$ in the presence of censoring with finite support, we will have to assume that $\varphi_1(\cdot,u)$ is bounded by the upper bound of the support of the censoring variable (see the discussion in Section (ref) and Assumption (ref)(ii) later in the paper). The fact that we restrict the parameter space to bounded functions $\psi_1$ also allows us to weaken the completeness conditions for identification.}

\textcolor{black}{Remark also that since (i) $\mathcal{P}$ is not a standard $L^2$-space with respect to a continuous distribution and (ii) what we wish to identify is not $(\varphi_0(u),\varphi_1(\cdot,u))$ but $\varphi(\cdot,u)$, we cannot rely on the high-level identification theory for nonseparable NPIV models as in CH,chen2014local. Therefore to obtain our identification results, we adapt the proofs of these papers to our specific model.}

Local identification

\textcolor{black}{We start by local identification.

} \textcolor{black}{Notice that we only seek to identify $\varphi(z,u)$ for all $z\in\mathcal{Z}$ and not $(\varphi_0(u),\varphi_1(z,u)), z\in \mathcal{Z}$. This is because it is not possible (and not interesting) to identify $\varphi_1(z,u)$ for $z>\varphi_0(u)$ because $\varphi_1(z,u)$ never generates the data for such $z$.} The main assumption for local identification is the following bounded completeness condition.

Intuitively, this condition means that $Z$ and $W$ are sufficiently dependent given $U=u$. Assumption (ref) is implied by the bounded completeness of $Z$ given $W,U=u$, that is for all {\color{black} bounded} functions

Such a condition is imposed in cazals2016nonparametric. In separable NPIV models, the related bounded completeness condition

for $m$ belonging to some class of functions, is often imposed newey2003,Darolles. Condition (ref) is a “conditional on $U=u$" version of condition (ref). Several authors have provided sufficient conditions for bounded completeness as stated in (ref) newey2003,xavier2011, hu2018nonparametric, andrews2011. These sufficient conditions are restrictions on the family of distributions $\{Z|W=w\}_{W=w}$ where “$Z|W=w$" stands for the distribution of $Z$ given $W=w$. To obtain sufficient conditions for (ref), it therefore suffices to impose the sufficient conditions from newey2003,xavier2011, hu2018nonparametric, andrews2011 on the family of distributions $\{Z|W=w,U=u\}_{W=w}$ rather than $\{Z|W=w\}_{W=w}$. As a last remark, notice that restricting conditions (ref) to bounded functions makes it more likely to hold (see the aforementioned papers for further details).

In addition to Assumption (ref), we also impose some regularity conditions. Since these conditions are technically involved, they are stated in the online appendix (see Assumption (ref) in Section (ref)). Remark that these regularity conditions require Gâteaux differentiability of some operator but no Fréchet differentiability is needed. \textcolor{black}{We have the following local identification result.

We have shown here local identification on a segment $\mathcal{N}_{\epsilon}$. As noted in chen2014local, in nonparametric nonlinear structural models (as the one of the present paper), it is often not possible to derive local identification results when $\mathcal{N}$ is an open ball (in the topology defined by $\|\cdot\|_{\mathcal{P}}$). We therefore focused on a smaller set, which is not an open ball.}

Global identification

\textcolor{black}{Next, we discuss global identification, which we define as follows.

We adapt the theory of CH to the case of the present model. We introduce $\epsilon=T-\varphi_0(u)(1-D)- \varphi_1(Z,u)D.$ Let $f_{\epsilon|D,W}(\cdot|0,w)$ be the density of $\epsilon$ given $D=0$ and $W=w$, $f_{\epsilon|D,Z,W}(\cdot|1,z,w)$ be the density of $\epsilon$ given $D=1$, $Z=z$ and $W=w$ and $f_{Z|D,W}(\cdot|d)$ be the density of $Z$ given $D=d$ (their existence is guaranteed under our assumptions). Let us make the next Assumption:

For $(\Delta_0,\Delta_1)\in\mathcal{P}$, we define

We make the following hypothesis:

This is a type of bounded strong completeness condition. It is the counterpart of Assumption L1$^*$ of CH in our model. This assumption holds when $Z$ and $W$ are sufficiently dependent given $U$. To the best of our knowledge, the only sufficient conditions known for this type of assumption correspond to condition L2$^*$ in CH. In Section (ref) of the online appendix, we give sufficient conditions for Assumption (ref) which are in the spirit of condition L2$^*$ in CH. These conditions include assuming that a family (in $w$) of distributions is boundedly complete and therefore relate Assumptions (ref) to standard bounded completeness conditions (as in (ref)). We have the following global identification result.

}

Identification with censoring

Let us now consider the case where $T$ is right censored. In this case, $F_0$ and $F_1$ may not be identified everywhere because of censoring. We make the following assumption on the censoring:

This assumption could be relaxed. For instance, the fact that $C$ and $U$ are independent given $Z,W$ would suffice for identification. However, we impose the stronger Assumption (ref) to simplify the exposition and the estimation.

Let $c_0$ be the upper bound of the support of $C$. Identification of $F_0(t,w)$ and $F_1(\psi,w)$ for $\psi:\ \mathcal{Z}\mapsto [0,t]$ is only possible for $t\in[0,c_0]$. Hence, we can only check that $\varphi(\cdot,u)$ is the solution of the identification equation for $u\in[0,u_0]$, where $u_0 =\inf\{u\in{\mathbb{R}}_+ : \varphi(z,u)< c_0 \: \forall z\in\mathcal{Z}\}.$ Define $G(t)= {\mathbb{P}}(C\ge t)$, the survival function of $C$. Since $F_0(t,w)= E[I(T\le t,D=0,W\le w)]$, \textcolor{black}{using the law of iterated expectations and Assumption (ref), we can show that}

for all $t\in[0,c_0)$. Similarly, for all $\phi:\mathcal{Z}\mapsto [0,c_0)$, it holds that

The proof of (ref) and (ref) is given in Section (ref) of the online appendix. By standard arguments from the survival analysis literature, $G(t)$ is identified for all $t\in[0,\sup\{t\in{\mathbb{R}}_+:\ {\mathbb{P}}(T\ge t)>0\}]$ \textcolor{black}{(on this interval $G(t)$ is equal to the population analog of the Kaplan-Meier estimator of the survival function of $C$ which identifies it)}. Hence, the quantities on the right-hand side of (ref) and (ref) are identified and so are $F_0(t,w)$ for all $t\in[0,c_0)$ and $ F_1(\phi, w)$ for all $\phi:\mathcal{Z}\mapsto [0,c_0)$. Therefore, identification results on $\varphi(\cdot,u)$ for $u\in[0,u_0]$ can be obtained as in the case without censoring. \textcolor{black}{Note however that the fact that we can identify $\varphi$ only up to $u_0$ has several implications. First, it means that average treatment effects ($E[T(z)-T(z')]$) are not identified. Second, only some quantile treatment effects ($\varphi(z,u)-\varphi(z',u)$) are identified. Third, the structural hazard $\lambda(z,t)=(\varphi(z,\cdot)^{-1})'(t)$ is identified only for $t\le \varphi(z,u_0)$.}

Estimation

Parametric regression function

Our strategy is to first estimate $F_0,F_1$ and $F_W$, and then to solve an estimate of equation (ref) for $\varphi_0$ and $\varphi_1$, which is obtained by plug-in of the estimates of $F_0,F_1$ and $F_W$. Since (ref) is a {complicated} integral equation, this method is unlikely to deliver precise nonparametric estimates of $\varphi$ on datasets of reasonable size. As a result, we decide to assume that $(\varphi_0,\varphi_1)$ follows a parametric model $\{\varphi_{\theta0},\varphi_{\theta1}\}_{\theta\in\Theta}$, that is $\varphi_0=\varphi_{\theta_*0}$ and $\varphi_1=\varphi_{\theta_*1}$ for some $\theta_*\in \Theta$, where $\Theta\subset{\mathbb{R}}^K$ is the parameter set. Here, for all $\theta\in\Theta$, $\varphi_{\theta0}$ is a mapping from ${\mathbb{R}}_+$ to ${\mathbb{R}}_+$ and $\varphi_{\theta1}$ is a mapping from $\mathcal{Z}\times {\mathbb{R}}_+$ to ${\mathbb{R}}_+$ For all $\theta\in\Theta$, we can also define $\varphi_{\theta}:\mathcal{Z}\times{\mathbb{R}}_+\mapsto{\mathbb{R}}_+$ such that $\varphi_{\theta}(z,u)=\varphi_{\theta0}(u)I(z>\varphi_{\theta0}(u))+\varphi_{\theta1}(z,u)I(z\le \varphi_{\theta0}(u))$ for all $z\in\mathcal{Z},u\in{\mathbb{R}}_+$. This approach has three additional advantages. First, it avoids the need for regularization, since the parameter set has finite dimension. Second, parametric shapes enable to summarize simply the properties of $\varphi$. Finally, they allow to estimate $\varphi(\cdot,u)$ even for $u> u_{0}$, which is useful when the study has insufficient follow-up. Note that, the model remains semiparametric since $F_0$ and $F_1$ are not parametrically constrained. Let us give some examples of parametric models for $\varphi$. \\

Example 1: Weibull model. A first example comes from the Weibull distribution. Recall that $T(z)$ is the potential outcome of $T$ when the treatment time is set to $z$. We assume that before $z$, the hazard rate of $T(z)$ corresponds to that of a Weibull distribution with parameters $\theta_{00},\theta_{01}$. After $z$, the hazard rate $T(z)$ is that of a Weibull distribution with parameters $\theta_{10},\theta_{11}$. The structural hazard of $T(z)$ at time $t$ is therefore given by

By inverting the cumulative hazard, it can be shown that

Example 2: Log-normal model. The second example comes from the log-normal distribution. Before $z$, the hazard rate of $T(z)$ is assumed to be equal to that of a log-normal distribution with mean $\theta_{00}$ and variance $\theta_{01}$. After $z$, the hazard rate of $T(z)$ corresponds to that of a log-normal distribution with mean $\theta_{10}$ and variance $\theta_{11}$. The structural hazard rate of $T(z)$ at time $t$ is then given by

where $\phi$ and $\Phi$ are respectively the density and the cumulative distribution function of a standard normal distribution. Inverting the cumulative hazard, we obtain

where $R_z(\theta) = \left[1 - \Phi \left(\frac{\log(z) - \theta_{10}}{\theta_{11}} \right)\right]\left[1 - \Phi \left(\frac{\log(z) - \theta_{00}}{\theta_{01}} \right)\right]^{-1}.$

Estimation of the integral equation

For $\theta\in \Theta$, let $$M_\theta(u,w) = F_0(\varphi_{\theta 0}(u),w) +F_1(\varphi_{\theta 1}(\cdot,u),w) - (1-e^{-u})F_W(w)$$ be the value of the identifying equation (ref) in $(\varphi_{\theta0},\varphi_{\theta1})$. In order to estimate $\theta_*$ using (ref), it is necessary to estimate the unknown operator $M$. Assume that we possess an i.i.d. sample $\{Y_i, \delta_i,\tilde Z_i, \tilde D_i, W_i\}_{i=1}^n$. We estimate $F_0$ and $F_1$ using (ref) and (ref). Let

The Kaplan-Meier estimator of $G(t)$ is given by $\widehat{G}(t) = \prod_{s< t}\left(1-\frac{dN(s)}{Y(s)}\right),$ where $dN(s)=N(s)-\lim\limits_{s'\to s,s'<s}N(s')$. In turn, $F_0$ and $F_1$ are estimated by

Then, $F_W$ is estimated by $\widehat{F}_W(w)= n^{-1}\sum_{i=1}^n I(W_i \le w)$. Finally, the estimator of $M$ is $$\widehat{M}_\theta(u,w)=\widehat{F}_0(\varphi_{\theta0}(u),w) + \widehat{F}_1(\varphi_{\theta1}(\cdot,u), w)- (1-e^{-u})\widehat{F}_W(w). $$ Remark that, although $Z$ and $W$ are continuous random variables, we avoid smoothing because we use an unconditional identification equation.

Estimator of $\theta_*$

It is computationally impossible to solve the estimated identifying equation for every $u$ and $w$. Hence, we solve it on a grid $0\le u_1<\dots<u_m<u_0$ of values of $u$, where $m\in {\mathbb{N}}_*$ is fixed, and at $w=W_1,\dots, W_n$. The estimator of $\theta_*$ is

where $\widehat{L}(\theta)= (nm)^{-1}\sum_{i=1}^n\sum_{j=1}^mp(u_j) \widehat{M}_\theta(u_j,W_i)^2$ and $p(\cdot)$ is some weighting function, typically $p(u)=e^{-u}$. {\color{black}This type of estimators is akin to minimum distance from independence estimators as in brown2002weighted. We can not rely directly on the theory of brown2002weighted. Indeed, in our case, because of censoring, $G$ has to be estimated in a first step. Moreover, weights on $W$ are chosen according to the empirical measure, while they are set according to a fixed measure selected by the researcher in brown2002weighted. This avoids the need to choose weights on $W$ and might lead to greater efficiency. }

Asymptotic normality

Now, we state the conditions that we impose to show asymptotic normality. The first assumption ensures identification, that is

where $L(\theta)=m^{-1}\sum_{j=1}^mp(u_j)E[M_\theta(u_j,W)^2].$

Condition (i) was studied in the previous section, whereas condition (ii) restricts the choices of $\Theta$ and $u_m$, and \textcolor{black}{(iii) is a constraint on the parametric family that guarantees that $\theta_*$ is identified when $\varphi(\cdot,u_j)$ is known for all $j\in\{1,\dots, m\}$. Together (i) and (iii) ensure that $\theta$ is identified from (ref), which yields (ref). Condition (ii) implies that $\varphi_{\theta_*1}(\cdot,u)=\varphi_1(\cdot,u)$ is bounded when $c_0<\infty$ (that is censoring has finite support).} Let $||\cdot||$ denote the Euclidean norm in ${\mathbb{R}}^K$. We also impose some regularity conditions:

These standard and mild conditions depend simultaneously on the regularity of the mapping $\theta\mapsto \varphi_\theta$ and on the distribution of $(U,Z,W)$.

Bootstrap

Since the asymptotic variance matrix $\Sigma$ has a complicated expression (see the proof of Theorem (ref)), we rely on the nonparametric bootstrap for inference. Let $\{Y_{bi}, \delta_{bi},\tilde Z_{bi}, \tilde D_{bi}, W_{bi}\}_{i=1}^n$ be the bootstrap sample drawn with replacement from the original sample $\{X_i=(Y_i, \delta_i,\tilde Z_i, \tilde D_i, W_i)\}_{i=1}^n$. Let also $\widehat{\theta}_b$ be the value of the estimator computed on the bootstrap sample $b$. The following result allows to build confidence intervals with the na\"ive bootstrap using the empirical distribution of the estimates in the bootstrap samples.

Numerical experiments

Simulations

For a sample of size $n$, and for $i = 1,\dots,n$, we generate the instrument $W_i$ and $U_i$ from two independent exponential distributions with parameter equal to one. We also generate an additional error term $r_i \sim Exp(1)$. The treatment time $Z_i$ is taken equal to $Z_i = \sqrt{2r_i U^{\alpha}_i W^{\beta}_i},$ where the parameter $\alpha$ controls the level of endogeneity, and the parameter $\beta$ controls the strength of the instrument. When $\alpha = 0$, the treatment time, $Z_i$, is independent of $U_i$, and when $\beta=0$, the instrument cannot explain any of the variation in the treatment time. Finally, $$ T_i = \varphi(Z_i,U_i) = \varphi_0(U_i) I(\varphi_0(U_i) < Z_i) + \varphi_1(Z_i,U_i) I(\varphi_0(U_i) \ge Z_i).$$

For $\varphi$, we consider the two parametric models of Section (ref). For the Weibull model (example 1), we take the true parameter vector $\theta = (\theta_{00},\theta_{10},\theta_{01},\theta_{11})^\top = (1,2,1.5,2)^\top$. For the log-normal model (example 2), the true parameter vector is $\theta = (\theta_{00},\theta_{10},\theta_{01},\theta_{11})^\top = (0,1,1,1)^\top$. For both designs, we consider $\alpha \in \{ 0.25,0.75 \}$ to vary the level of endogeneity, and $\beta \in \{ 0.5,1 \}$ to vary the strength of the instrument. We further take censoring into account as follows:

The sample size is fixed to $n \in \{ 500,1000,3000 \}$. We thus have a total of $3 \times 2^4 = 48$ simulation schemes, and we run $R = 1000$ replications for each scheme. The optimization algorithm is started at 100 random values. Each of these starting values yields a local minimum. The local minimum which leads to the lowest value of the objective function corresponds to the estimate. We use a grid $u_1,\dots,u_{100}$ of $100$ values of $U$, where $u_1$ (respectively $u_{100}$) is the 0.025 (respectively, 0.975) quantile of the unit exponential distribution and the points are equally spaced.

We report the bias and standard error of the estimator and the coverage of the bootstrap percentile confidence intervals at level $90\%,\ 95\%$ and $99\%$. As the computation time increases with the sample size, the coverage of the confidence intervals are evaluated using the Warp-Speed method of giacomini2013, where only one bootstrap resampling is used for each simulated sample. We also provide the average number of treated units observed, $\bar{D}$, and the average number of uncensored observations, $\bar{\delta}$.

We summarize the simulations results for the Weibull design with censoring in Table (ref). Between $33$ and $43\%$ of observations are treated. With a strong instrument ($\beta=1$), the bias is low, even when $n=500$. The coverage of the confidence intervals improves when $n$ grows and are relatively close to nominal for $n=3000$. When the instrument is weaker ($\beta=0.5$), the estimator exhibits good performance in terms of bias when $n=1000$ or $n=3000$. \textcolor{black}{The results for the Weibull design without censoring and the log-normal model (with and without censoring) are reported in Section (ref) of the online appendix.}

\afterpage{ \thispagestyle{empty}

}

Empirical application

We use data from a large Belgian insurance company. This insurer offers a product to other companies, which consists in paying for the salaries of their client' workers who are on medical leave because of burnout. This insurance product also contains the possibility of following a free therapy for these workers suffering from burnout. The goal of the company is to reduce the duration of medical leave (our duration variable in this application) and, hence, it would like to know the effect of the start of the therapy (our treatment) on the duration of medical leave.

We take advantage of the treatment assignment mechanism to evaluate the causal effect of the treatment. The employees on medical leave are assigned to one of several partner institutions for medical care, which are responsible for carrying out the therapy. These partners evaluate the patient and decide to propose or not the therapy on the basis of the expected medical benefits of the therapy. The employee may then accept or decline to be treated, which suggests that $Z$ is endogenous. The exact date of the start of the therapy varies depending on the availabilities of the partner and the patient, which makes the treatment time-varying.

The assignment by the insurance company of the partner was only based on geographical distance, which should make it exogenous. Moreover, partners are more or less likely to offer a therapy to assigned patients, hence the propensity of the partner to offer treatment has an impact on the treatment time $Z$. We use the proportion of patients who are given the possibility to be treated by the partner in a given year as the instrumental variable. This choice of instrument is common in medical studies (see brookhart2006evaluating,chen2011use for reviews).

The sample consists of 838 individuals who entered medical leave for burnout between 2017 and 2020. Their age at the start of the medical leave was between 30 and 39 years old. They are observed until one of the following events happens: their medical leave ends, their medical leave exceeds 2 years, or the study is ended (at the end of 2020). In the latter two cases the individual is censored. Since the duration of follow-up depends on external factors, we expect censoring to be uninformative. Around 41% of the observations are censored and 48% are treated before censoring. The average (respectively, median) duration of medical leave for uncensored observations is 189 days (respectively, 159 days). For the observations for which the treatment time is uncensored, the average is 112 days and the median is equal to 90 days.

For both the Weibull model and the log-normal model, we estimate the parameters using 100 random starting values around the values corresponding to the na\"ive fit of the Weibull distribution or the log-normal distribution to the data of durations and censoring indicators. We use a grid $u_1,\dots,u_{100}$ of $100$ values of $U$, where $u_1$ is the 0.025 quantile of the unit exponential distribution and the points are equally spaced. The upper bound $u_{100}$ is chosen such that

where $c_0$, the maximum follow-up time, is equal to 2 years. Condition (ref) ensures that Assumption (I) (ii) is satisfied for some $\Theta$ in a neighborhood of $\widehat{\theta}$. Following this approach, we chose $u_{100}$ equal to the 0.9 (respectively, 0.8) quantile of a unit exponential for the Weibull (respectively, log-normal) model. The curves of the hazard rates for the individuals who are never treated ($z=\infty$) and those who received treatment at time 0 ($z=0$) corresponding to the estimated parameters are plotted in Figure (ref) for the Weibull model and in Figure (ref) for the log-normal model.

These hazard rates are computed by plug-in of the estimates in (ref) and (ref). By definition of our Weibull and log-normal models, for an arbitrary value of $z$ the structural hazard rate at $t$ under treatment at time $z$ is equal to that of the never treated for $t \le z$ and is equal to that of the treated at time $0$ for $t>z$. Hence, in Figures (ref) and (ref), the estimated structural hazard under treatment $z$ corresponds to the red curve before $z$ and then “jumps” to the blue dashed curve.

We see also that the therapy appears to increase the hazard rate for all possible treatment timings. As a result, the treatment should be administered at time $0$ in order to minimize the duration of medical leave.

Bootstrap confidence intervals for the hazard rates are given in Section (ref) of the online appendix. Note that the estimated hazards exhibit different shapes under the two models. This is to be expected since the models assume different parametric forms. However, the treatment significantly increases the hazard rate in both models (see the bootstrap confidence intervals in the online appendix). The robustness of this conclusion to the choice of the model constitutes statistical evidence supporting the efficacy of the therapy.

Concluding remarks

This paper develops an instrumental variable approach to estimate the causal effect of the time until a treatment is started on a possibly right-censored duration outcome. Therefore, the treatment $Z$ corresponds to the jump of a counting process with single jump. As an extension, it would be of interest to consider procedures where $Z$ is a more general process. Another possible research direction could be the development of a control function approach in the context of the present paper. \textcolor{black}{This might allow us to use a discrete instrument, as in xavier2015,torgovitsky2015.}