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.
72,146 characters · 20 sections · 38 citation commands
Instrumental variable estimation of the proportional hazards model by presmoothing
\newif\ifband \newif\ifcorr \bandfalse \corrfalse \spacingset{1.2}
\spacingset{1.4}
The proportional hazards (henceforth, PH) model, first introduced in Cox1972's seminal paper, has become a popular tool for analyzing time-to-event data, such as the unemployment duration of job seekers or the duration until the bankruptcy of some firms. It is a semiparametric model of the hazard rate, which has the convenient property that the nonparametric component of the model - called the baseline hazard - can be profiled out in the likelihood. The resulting partial likelihood allows for a simple and stable estimation procedure of the parametric part of the model.
In practice, observational data may be affected by unmeasured confounding that influences both the treatment and the outcome, creating an endogeneity issue. In econometrics, a common solution to the endogeneity problem is the use of an instrumental variable (henceforth, IV) approach. There exists several different IV models, which differ in terms of the unobserved heterogeneity that they allow.
In this paper, we address the problem of endogeneity in the PH model by embedding the problem in the instrumental variable quantile regression (henceforth, IVQR) model introduced in chernozhukov2005iv. Indeed, by making a rank invariance assumption, we reformulate the proportional hazards model into a semiparametric version of the instrumental variable quantile regression model of chernozhukov2005iv. This allows us to use some of the tools from the literature on the IVQR model. The IVQR model and its associated estimation procedures have become popular in econometrics research, see Chapter 9 in koenker2017handbook for a recent review. By constraining the treatment effects to only vary by quantile, this approach allows estimating quantile treatment effects over the full population. This is in contrast to the local average treatment effects framework of angrist1996identification, which does not restrict the heterogeneity of treatment effects but requires a monotonicity assumption to hold and only allows to estimate treatment effects on subpopulations. See wuthrich2020comparison for a detailed comparison of the two approaches.
One of the main difficulties generated by the IVQR model is estimation. The causal regression function of interest can be shown to be the solution of a system of nonlinear integral equations. A natural approach to estimation is to solve an empirical version of this system. However, as noted in chernozhukov2006instrumental, when the causal regression function follows a (semi)parametric model this leads to a severely nonsmooth and nonconvex objective function. For the linear-in-parameters semiparametric quantile regression model, some ingenuous solutions to this problem have been developed in chernozhukov2006instrumental, kaplan2017smoothed and kaido2021decentralization among others. However, to the best of our knowledge, there does not exist a general computationally advantageous approach to the estimation of semiparametric IVQR models.
In this paper, we focus on the case where the endogenous variable and the instrument are both discrete and there are (possibly continuous) exogenous covariables. To simplify the exposition, we assume that the endogenous variable and the instrument have the same number of modalities. This setting is very relevant in practice since it includes randomized experiments with noncompliance (where both the instrument and the endogenous variable are binary), which is the main focus of chernozhukov2005iv, or, more recently, wuthrich2019closed. In the remainder of the paper, we often call “treatment” the discrete endogenous variable because of this link. As common with time-to-event data, we assume that the outcome variable is randomly right-censored, where the censoring time is conditionally independent of the outcome.
We propose a new solution to the problem of estimating IVQR models with a semiparametric quantile regression function. We apply this new methodology to the PH model but stress that it could be useful in other semiparametric models, such as the proportional odds model or distribution regression. This approach is based on a presmoothing strategy. Our procedure can be summarized in three steps.
Step 1 is simplified by the fact that the treatment and instrument are discrete, making the inverse problem well-posed. Moreover, as we argue in the paper, in the practically relevant case of randomized experiments with one-sided noncompliance, where both the treatment and the instrument are binary and there is full compliance in the control group, our nonparametric estimator possesses a closed-form expression. For Step 1, we use beran1981nonparametric's estimator of the conditional survival function to construct an empirical version of the system of equations generated by the model and then solve the estimated system to obtain the causal regression function of interest. This allows us to deal with random right-censoring, a common feature of duration data. The fact that we can use the partial likelihood estimator in Step 3 is very convenient since the latter has a convex objective function and is implemented in traditional statistical software. We stress again that the proposed procedure could be applied to other models (such as the proportional odds model or distribution regression) by replacing Step 3 with standard estimators of these models (under exogeneity). The main advantage of our procedure is to bring us back (in Step 2) to the case of exogenous data, which is well-studied and for which most semiparametric models were originally motivated.
We show that the resulting estimator is asymptotically normal. We also develop a theory for the nonparametric estimator of Step 1 which is new to the literature. Simulations demonstrate that our estimator possesses good finite sample properties. Finally, the procedure is applied to evaluate the causal effect of reemployment bonuses using data from the Illinois Unemployment Incentive Experiment.
We call our methodology presmoothing because it is related to literature in (semi)parametric statistics advocating for estimation procedures that first use nonparametric estimators before projecting the nonparametric estimators on (semi)parametric spaces in a second step, see cristobal1987class,akritas1996use or, more recently, musta2022presmoothing. The main takeaway from this statistical literature is that, asymptotically, the two-step estimator recovers the property of the “oracle” second-step estimator (which assumes that the nonparametric object estimated in the first step is known). Presmoothing is particularly useful in cases where estimating directly the semiparametric model is complicated while nonparametric estimation is simpler. One of the contributions of this article is to bring this presmoothing tool to econometrics.
Finally, before outlining the different sections of the paper, we note that there already exist quite a few papers studying instrumental variable methods for duration models. To the best of our knowledge, no other paper considers the PH model embedded in the IVQR model. Please note, that for reasons of space, we do not cite all relevant works. bijwaard2005correcting consider a mixed proportional hazards model in which the unobserved heterogeneity enters multiplicatively, that is differently from the IVQR model. wang2022instrumental considers a PH model with endogeneity but they make an additive (rather than quantile) homogeneity of treatment effects assumption. As argued in chernozhukov2005iv, the restrictions on unobserved heterogeneity in the IVQR model can be seen as weaker. beyhum2022nonparametric studies the nonparametric IVQR model with random right censoring with discrete treatment and instrument but no covariates. The present paper generalizes beyhum2022nonparametric by allowing for continuous covariables in nonparametric estimation and studying the three-step estimation of the PH model. hong2003inference,chen2018sequential,wang2021moment,beyhum2023instrumental study the linear-in-parameters IVQR model under random right-censoring but do not cover the PH model. The way we reformulate the PH duration model into an IVQR model is related to the reformulation of the dynamic duration model into a dynamic IVQR model in beyhum2023variable.
The rest of the paper is organized as follows. In Section (ref), we specify the model. Identification results are derived in Section (ref). Section (ref) is devoted to the estimation theory. In Section (ref) we discuss the asymptotic theory of the proposed estimator. In Section (ref), the finite sample performance of the proposed method is investigated through simulations. In Section (ref), we illustrate the method by means of an empirical application of the method. All the technical details are deferred to the Appendix.
We are interested in the effect of an endogenous treatment $Z$ on a non-negative outcome $T$. We consider a treatment $Z=(Z_1,\dots,Z_{d_Z})^\top$ expressed as a vector of $d_Z$ binary variables and taking $L$ different values $\{z_1,\dots,z_L\} = \mathcal{Z}$. In our data example e.g., we have a categorical treatment with three levels that can be defined by means of two binary (dummy) variables, so $d_Z=2$. Since there are three levels, $L=3$ in this case and $z_1=(0,0)^\top, z_2=(1,0)^\top$and $z_3=(0,1)^\top$. There are some exogenous covariables $X$ with compact support $\mathcal{X}\subset\mathbb{R}^{d_X}$. Let $T(z,x)$ be the potential outcome duration under treatment equal to $z\in\mathcal{Z}$ and covariables equal to $x\in\mathcal{X}$. Let us assume that $T(z,x)$ is a continuous random variable and also define by $$\lambda(z,x,t):= \lim_{dt \to 0}\frac{{\mathbb{P}}(T(z,x)\in[t,t+dt]|T(z,x)\ge t)}{dt} $$ the structural hazard of $T(z,x)$ at time $t\in{\mathbb{R}}_+$. We impose the consistency condition $T=T(Z,X)$. We also assume that $T(z,x)$ follows a proportional hazards model, that is there exists a baseline hazard $\lambda_0:{\mathbb{R}}_+\xrightarrow{}{\mathbb{R}}_+$ and a vector of regression coefficients $(\beta_{0,z},\beta_{0,x})\in\mathbb{R}^{d_Z+d_X}$ such that
Remark that this proportional hazards assumption implies that the support of $T(z,x)$ does not depend on $z,x$ and is equal to that of $\lambda_0$, or equivalently of $T$. We denote by $\mathcal{T}$ this support of $\lambda_0$. Let us now define the structural cumulative hazard $$\Lambda(z,x,t):=\int_0^t \lambda(z,x,s)ds=\Lambda_0(t)\exp((z,x)^\top\beta_0),$$ where $\Lambda_0(t):=\int_0^t \lambda_0(s)ds$ is the cumulative baseline hazard. Let us also introduce $Q(z,x):= \Lambda(z,x,T(z,x))$, which is the structural hazard of $T(z,x)$ evaluated at $T(z,x)$. The random variables $\{Q(z,x)\}_{(z,x)\in\mathcal{Z}\times \mathcal{X}}$ can be thought of as the unobserved heterogeneity of the model. We make the following Assumption.
By Assumption (ref) $(ii)$, the structural cumulative hazard $\Lambda(z,x,\cdot)$ is strictly increasing on $\mathcal{T}$. Hence, we can define the inverse of $\Lambda(z,x,\cdot)$ on $\mathcal{T}$, which we denote by $\Lambda(z,x,\cdot)^{-1}$. This means that $\Lambda(z,x,\cdot)^{-1}(s)$ is the unique element $t$ of $\mathcal{T}$ such that $\Lambda(z,x,t)=s$ (by definition the image of the cumulative hazard on $\mathcal{T}$ is ${\mathbb{R}}_+$ so that $\Lambda(z,x,\cdot)^{-1}(s)$ is defined on all ${\mathbb{R}}_+$). Because of this, we have $T(z,x)= \Lambda(z,x,\cdot)^{-1}(Q(z,x))$. Assumption (ref) $(i)$ then implies that for two subjects $i$ and $j$ , $T_i(z,x)>T_j(z,x)$ implies $Q_i(z,x) = Q_i>Q_j = Q_j(z,x)$, which leads to $T_i(z',x')>T_j(z',x')$, for all $(z,x),(z',x')\in\mathcal{Z}\times \mathcal{X}$. Hence, under Assumption (ref) $(i)$ the rank in the outcome of any two subjects is the same across all potential outcomes. Because of this, Assumption (ref) $(i)$ is a rank invariance assumption as in chernozhukov2005iv. 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, note that the rank invariance assumption does not restrict the possible values of the structural hazard $\lambda(z,x,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,x)$. It only imposes limits on the joint distribution of potential outcomes, that is the distribution of $(T(z,x))_{(z,x)\in \mathcal{Z}\times \mathcal{X}}$. This rank invariance assumption could be relaxed into a rank similarity assumption as in chernozhukov2005iv while keeping all results valid. Remark that Assumption (ref) allows us to write
We assume to observe a categorical instrumental variable $W$. To simplify the theoretical arguments, we also suppose that the instrument has the same number of modalities as the treatment. Therefore, $W$ has support $\mathcal{W} = \{w_1,\dots,w_L\}\subset \mathbb{R}^{d_W}$. Differently from the variable $Z$, the instrument $W$ does not need to consist of binary variables. Importantly, assuming that the cardinalities of the support of $Z$ and $W$ are equal does not limit our discussion to a specific case. First, it would be possible, but burdensome, to extend both identification and estimation results to the case where the number of modalities of $W$ is greater or equal to the number of modalities of $Z$. Second, it is always possible to obtain the same number of modalities between $Z$ and $W$ by aggregation of modalities of $W$. We impose the following condition.
where `$\protect\mathpalette{\protect\independenT}{\perp}$' stands for statistical independence. This assumption formalizes the exogeneity of the instrument $W$ and covariables $X$. The duration $T$ is randomly right censored by a random variable $C$ with support in ${\mathbb{R}}_+$ so that we do not observe $T$ but $Y=\min(T,C)$. The observables are $(Y, \delta, Z, X, W)$, where $\delta =I(T\le C)$. We impose the following standard independent censoring assumption:
In this section, thanks to Assumption (ref), we reformulate the duration model into a semiparametric IVQR model. Let us define, for $u\in[0,1]$,
and $U=1-\exp(-Q)$. By standard arguments in duration analysis, $Q$ follows a unit exponential distribution, so that $U$ follows a standard uniform distribution. Then, inverting equation (ref), we obtain that
where $\mathcal{U}[0,1]$ stands for the uniform distribution on $[0,1]$. Remark that by Assumption (ref) $(ii)$ and the inverse function theorem, $\varphi(z,x,\cdot)$ is strictly increasing and differentiable with continuous derivative on $[0,1]$, for all $(z,x)\in\mathcal{Z}\times \mathcal{X}$. As claimed in the introduction, this is an IVQR model as in chernozhukov2005iv, where the causal quantile regression function $\varphi$ follows the semiparametric model given in (ref), which is generated by the proportional hazards assumption on the structural hazard. From now on, this paper will analyze directly this semiparametric IVQR model. \\
Notation. Finally, we introduce some additional notation. Define $V = (Z,X) \in \mathcal{V} = \mathcal{Z}\times \mathcal{X}\subset\mathbb{R}^{d_V}$, where $d_V=d_Z+d_X$. For $v = (z,x)$, we write $\varphi(v,\cdot)$ for $\varphi(z,x,\cdot)$. For $l = 1,\dots,L$, we also write $\varphi_l^x(\cdot)$ for $\varphi(z_l,x,\cdot)$ and $\varphi^x(\cdot)$ for $(\varphi_l^x(\cdot))_{l=1}^L $.
In this section, we provide an identification result of the parameter of interest $\beta_0$. The identification can be obtained in two steps. First, we analyze the identification of $\varphi$. Second, we study the identification of $\beta_0$ given the first step. Note that we keep the analysis of identification brief since it does not constitute the main contribution of the paper (which is the estimation procedure).
Let us first recall the standard characterization of $\varphi$ as the solution to a system of nonlinear equations, that is, for every $x\in\mathcal{X}$, $\varphi^x(u)$ is a solution of the following system of equations in $(\theta_l)_{l=1}^L\in\mathbb{R}^{L}_{+}$:
where $ F(t,z|x,w) = P(T\le t, Z =z|X = x,W=w)$. Indeed, by the fact that $\varphi(z,x,\cdot)$ is strictly increasing and the independence between $U$ and $(W,X)$, for all $k\in\{1,\dots,L\}$, it follows that:
chernozhukov2005iv give conditions under which there is a unique solution to (ref). In the context of the PH model, there is an additional layer of complexity due to the fact that $F$ is not identified everywhere because of censoring. For $z\in\mathcal{Z}, x\in\mathcal{X}$ and $w\in\mathcal{W}$, let us denote by $c_{z,x,w}$ the upper bound of the support of the distribution of $C$ given $Z=z,X=x, W=w$. Also, given two real numbers $c,d$, we denote by $c\wedge d$ the value $\min(c,d)$. Moreover, let $\bar t$ be the upper bound of $\mathcal{T}$. Thanks to Assumption (ref) and by standard “Kaplan-Meier type” arguments from the duration analysis literature, $F(\cdot, z|x,w)$ is identified on $[0,\bar t\wedge c_{z,x,w}]$, and $\bar t\wedge c_{z,x,w}$ is itself identified. Let us now introduce $$\bar u^x= \inf_{z\in\mathcal{Z}, w\in\mathcal{W}} \varphi(z,x,\cdot)^{-1}(\bar t\wedge c_{z,x,w}),$$ where $\varphi(z,x,\cdot)^{-1}(t)=\inf\{u\in\mathbb{R}_+:\varphi(z,x,u)\ge t\} $ is the pseudo-inverse of $\varphi(z,x,u)$. For $u\le \bar u^x$, the left-hand-side of (ref) will be identified at $(\theta_\ell)_{\ell=1}^L =\varphi^x(u)$ for all $x\in\mathcal{X}$, so that identification can follow along the lines of chernozhukov2005iv. We provide below an identification result assuming that the solutions to the system (ref) are unique, a property for which sufficient conditions can be found in chernozhukov2005iv.
The proof of identification of $\beta_0$ is constructive, meaning that we can use its rationale to estimate $\beta_0$ in the subsequent section. The proof uses the knowledge of $\varphi$ to solve the problem of identifying $\beta_0$ under exogeneity.
Formally, let $\bar u:=\inf_{x\in\mathcal{X}} \bar u^x$. From Theorem (ref), $\varphi$ is identified for $(z,x,u)\in \mathcal{Z}\times\mathcal{X}\times[0,\bar u]$. As $\bar u$ is not known precisely in practice, we choose a point $\bar U$ such that $0<\bar U<\bar u$. We then define two random variables $U_g$ and $U_g^c$ such that $(U^g,U_g^c)\protect\mathpalette{\protect\independenT}{\perp}(Z,X)$. The subscript `g' highlights that these random variables are generated. The variable $U_g$ is drawn from a uniform distribution $\mathcal{U}[0,1]$, is independent of $U^c_g$, and serves as the generative process. $U_g^c$ acts as a censoring variable for $U_g$ with support in $[0,\bar U]$, ensuring $\min(U_g,U_g^c)$ does not exceed $\bar U$, thus respecting the identification constraints. Hence, any continuous distribution on $[0,\bar U]$ for $U_g^c$ is acceptable.
We then define the random variables $T_g=\varphi(Z,X,U_g)$, $C_g =\varphi(Z,X,U_g^c)$, $Y_g= T_g\wedge C_g$ and $\Delta=I(T_g\le C_g)$. By Theorem (ref), the distribution of $(Y_g,\Delta, Z,X)$ is identified. As shown in Section (ref), the random variable $T_g$ follows a PH model with parameter $\beta_0$, where exogeneity holds as $U_g$ is independent of $(Z,X)$. Furthermore, the censoring mechanism is conditionally independent as $T_g\protect\mathpalette{\protect\independenT}{\perp} C_g|Z,X$ by construction. Hence, $\beta_0$ is identifiable from the distribution of $(Y_g,\Delta, Z,X)$ under standard assumptions for the identification of the PH model (e.g., see tsiatis1981large).
This argument leads us to the following theorem.
For later theoretical discussions, it will be beneficial to set $\tilde U = \min(U_g,U_g^c)$. We will also employ the equalities $Y_g = \varphi(Z,X,\tilde U)$ and $\Delta = I(U_g\le U_g^c)$.
In this section, we discuss the estimation of the parameter vector $\beta_0$, given an independent and identically distributed (iid) sample $\{(Y_i,\delta_i,Z_i,X_i,W_i)\}_{i=1}^n$ of observables. We first outline a naïve, but ultimately unsuitable approach to estimation before presenting the three steps of the estimation procedure mentioned in the introduction. For sake of simplicity, we restrict ourselves to the case $d_X=1$, so that the already involved conditions we state next do not depend on the parameter $d_X$.
A naïve approach to estimation would proceed as follows. By substituting equation (ref) into (ref), we obtain $$ \sum_{l=1}^L F\left(\left.\Lambda_0^{-1}\left(\frac{-\log(1-u)}{\exp((z_l^\top,x^\top)\beta_0)}\right),z_l\right|x,w_k\right)= u \quad\text{for}\, k=1,\dots,L,\,u\in[0,\bar u).$$ This leads to
One could estimate $F$ by a conditional beran1981nonparametric's estimator (see Section (ref) for its definition) and then obtain an empirical analog to equation (ref). This empirical criterion could be minimized over $\Lambda_0$ and $\beta_0$. As noted in the introduction, this approach would result in a highly nonsmooth and nonconvex problem.
We first build a nonparametric estimator of $\varphi$. The estimator is based on (ref), where $F$ is first estimated via a conditional Beran's estimator \\
Step 1.1: Nonparametric estimator of $F$ under random right censoring
Remark that $F(t,z|x,w) = F(t|z,x,w)p_{z,x,w}$, where $F(t|z,x,w) = P(T\le t|Z=z,X=x,W=w)$ and $p_{z,x,w}= P(Z=z|X=x,W=w)$. We propose to estimate $F(t|z,x,w)$ non-parametrically using the beran1981nonparametric estimator, given by
where $\eta_j(t) = I(Y_j\le t, \delta_j = 1)$ and $B_{h}^{z,w}(x-X_k,Z_k,W_k)$ is a sequence of non-negative weights adding up to 1. In our case, we adopt the Nadaraya-Watson type weights, which are specified as follows:
where $K(\cdot)$ is a univariate kernel function, and $h=h_n$ is a bandwidth depending on $n$ and converging to zero as $n\xrightarrow{}\infty$. To ensure that the final criterion function is smooth, we smooth $\tilde F$ with respect to time. Consider a further kernel $\tilde K$ and define $H(t) = \int_{-\infty}^t \tilde K(u)\text{d}u$, with a bandwidth $\epsilon = \epsilon_n$ depending on $n$ and converging to zero as $n\xrightarrow{} \infty$. We obtain
Consider now the following estimator for the quantity $p_{z,x,w}$:
The final estimator for $F(t,z|x,w) $ is given by
Step 1.2: Final nonparametric estimator of $\varphi$\\ Let $\bar T<\infty$ be an upper bound on $\max_{v\in\mathcal{V}}\varphi(v,\bar U)$, let $\mathcal{F}_{Z}^{\bar U, \bar T} = \{f:\mathcal{Z}\times [0,\bar U]\xrightarrow{} [0,\bar T]\}$ and let $\mathcal{F}_{\uparrow}^{Z,W}$ be the set of maps from $\mathbb{R}_+\times \mathcal{Z}\times \mathcal{W}$ to $\mathbb{R}_+$ that are continuous and increasing in the first argument. Let $x\in\mathcal{X}$ be fixed. We denote by $F^x(t,z|w)$ the function $F(t,z|x,w)$. We have shown that for each $x\in\mathcal{X}$, $\varphi^x$ belongs to the set of solutions to the equation:
where $A$ is an operator that associates to an element $(\tilde\varphi^x,\tilde F^x)\in \mathcal{F}_{Z}^{\bar U, \bar T}\times \mathcal{F}_{\uparrow}^{Z,W}$, a map from $[0,\bar U]$ to $\mathbb{R}^L$, denoted by $A(\tilde\varphi^x,\tilde F^x)$ and defined as $$ A(\tilde\varphi^x,\tilde F^x)(u) = \Big(\sum_{l=1}^L \tilde F^x(\tilde \varphi^x_{l}(u), z_l|w_k) - u\Big)_{k=1}^{L}, $$ where we wrote $\tilde\varphi^x_l(u)$ for $\tilde\varphi^x(z_l,u)$. Thus, we define the estimator of $\varphi^x$ by
where $\|\cdot\|$ denotes the Euclidean norm.
We now show how to use the nonparametric estimator obtained in Step 1 to generate “proxy” observations for which the treatment is exogenous. For each $i=1,\dots,n$, we generate a random observation of $(U_g,U_g^c)$, namely $(U_{g,i},U_{g,i}^c)$, from which we obtain $\tilde U_i = \min(U_{g,i},U_{g,i}^c)$ and $\Delta_i = I(U_{g,i}\le U_{g,i}^c)$. Then, the proxy observations correspond to $\{(\hat Y_{g,i}, \Delta_i,V_i)\}_{i=1}^n$, where $\hat Y_{g,i}=\hat \varphi(V_i,\tilde U_{i}) = \hat\varphi^{X_i}(Z_i,\tilde U_i)$.
The last step involves obtaining an estimator $\hat\beta$ for $\beta_0$ using the standard partial likelihood estimator of Cox1972 based on the “proxy” observations obtained in Step 2. Denote by $\mathcal{B}$ a compact set that contains $\beta_0$ as an internal point. Thus, $\hat \beta$ is obtained as the minimizer of the standard score function for a PH model, and so
We now discuss the optimization program (ref), noticing that the objective function $\|A(\theta,\hat F^x)\|$ might have multiple local minima. Therefore, the computation of the estimator $\hat \varphi^x$ can require starting the optimization algorithm at different initialization points. However, when the system of equations provided in (ref) reduces to a triangular system, the estimator is easier to compute. This means that, after possibly relabelling the points in $\mathcal{Z}$ and $\mathcal{W}$, for a fixed value $u\in[0,\bar U]$ the system takes the form:
This case occurs, for instance, in the presence of one-sided noncompliance, that is, $Z$ and $W$ are binary with support $\mathcal{Z}=\mathcal{W}=\{0,1\}$, and $P(Z=1|W=0,X) = 0$. In fact, the corresponding system takes the form:
Another example, in which $Z$ and $W$ are not binary, corresponds to the empirical application discussed in Section (ref).
When the system is of the form provided in (ref), the estimator has a closed-form expression based on the following inductive algorithm:
where $\hat F(\cdot,z|x,w)^{-1}(u) = \inf\{t\in\mathbb{R}_+:\hat F(t,z|x,w)\ge u\} $ denotes the pseudo-inverse of the function $\hat F(\cdot,z|x,w)$. Note that, based on the assumptions on the model, $F(\cdot,z|x,w)$ is strictly monotone increasing on $[0,\bar T]$, for each $z\in\mathcal{Z},x\in\mathcal{X},w\in\mathcal{W}$. In addition, by its definition, $\hat F(t,z|x,w)$ is monotone increasing in $t\in[0,\bar T]$, and this implies the inversion procedure in (ref) is well-posed.
In conclusion, the computation of $\hat \varphi^x$ only requires the inversion of some functions, which is easily available in any programming software. For instance, in R, the inversion of a function can be performed using uniroot.
We summarize the algorithm for the estimation of $\hat \beta$ as follows:
In this section, we will discuss the asymptotic results of the estimator $\hat \beta$. First, we will consider its consistency. Then, we will discuss its asymptotic normality.
Define the following quantities:
Define also
and
Consider the following regularity assumption, where $\alpha\in(0,1)$ corresponds to the order of Lipschitz continuity needed for later theoretical results. A further specification of the value of $\alpha$ is given in Appendix (ref).
A valid example of the values for $(u,w,\alpha,\nu,\pi)$ is, for instance, $(0.143,0.144,0.001,4,3)$. With this choice, we can also construct a valid example for $K(x)$ and $\tilde K(x)$. Consider the function $\bar K(x)= I(|x|\le 1)\frac{693}{512}(1-x^2)^5$. Then, it is easy to check that $K(x) = \tilde K(x) = \frac{3}{2}\bar K(x) + \frac{1}{2}x\bar K'(x)$ satisfies condition $(ii)$.
We obtain the following result, which concerns the consistency and the rate of convergence of $\hat F(t,z|x,w)$ to $F(t,z|x,w)$. In addition, it states the consistency of $\hat f(t,z|x,w)$ in $t$ and $x$ up to order $2+\alpha$. These results will be used to ensure the regularity and the rate of convergence of the estimator $\hat \varphi$, which is crucial for the rate of convergence of $\hat\beta$. Lastly, we show that the estimators $\hat F(t,z|x,w)$ and $\hat f(t,z|x,w)$ admit an iid representation that will be employed to obtain an analogous representation for $\hat \varphi$, which will lead to the asymptotic normality of $\hat{\beta}$.
To proceed, denote by $\|\cdot\|_\infty$ the infinity norm. The argument over which the supremum is taken may vary throughout the paper, and it is specified whenever any doubt might arise. Introduce the following assumption:
Assumption (ref) is related to the shape of the objective function $A$ and ensures that there is a unique solution to the system of equations (ref) in $[0,\bar U]$ and, hence, that there is a unique minimum to the program (ref) when $\hat F$ estimates $F$ well enough.
Denote the Fréchet derivative of $A$ in its first argument at the point $(\tilde\varphi^x,\tilde F^x)\in \mathcal{F}_{Z}^{\bar U, \bar T}\times \mathcal{F}_{\uparrow}^{Z,W}$ by
where the subscript $l,k$ indicates the $(l,k)$th entry of the matrix for $l,k=1,\dots,L$, and $\tilde f^x$ is the derivative in the first argument of $\tilde F^x$. Consider the following assumption.
Note that Theorem (ref) does not specify the asymptotic covariance matrix of $\sqrt{n}(\hat\beta-\beta_0)$ since the proof is based on showing that the estimator can be written as an empirical process over a Donsker class of functions. This result ensures asymptotic normality, but to obtain an explicit expression of the covariance matrix, one would need a complete representation of the empirical process. This is possible, but, as can be easily seen from the proof, the expression would be so involved that it would lack meaningful insights into the distribution and be difficult to use for practical inference purposes.
In this section, we evaluate the finite sample performance of our proposed method using Monte Carlo simulation. Our main objective is to compare the performance of the partial likelihood estimator proposed by Cox, an estimator that does not take endogeneity into account, with our new estimator. We will consider several designs, of which the specifications are given in Table (ref). Since our method allows for the inclusion of both discrete and continuous exogenous variables $X$, we will consider several univariate distributions for $X$: the continuous Beta design, the continuous uniform design, and the discrete Bernoulli design.
For all designs, we set the structural baseline cumulative hazard function $\Lambda_0$ to that of an exponential variable with mean 1, i.e. $\Lambda_0(s) = s$, and fix the value of $\beta_0 = (\beta_{0,z},\beta_{0,x})^\top$. This results in the potential outcome $T(z,x)$ following an exponential distribution with a rate (reciprocal of the mean) equal to $\exp(x\beta_{0,x})$ if $z = 0$ and $\exp(\beta_{0,z}+x\beta_{0,x})$ if $z=1$. The specifications of the design are given in Table (ref). It is worth noting that the variable $Z$ is endogenous as it depends on the unobserved heterogeneity term $U$, and that $W$ is a proper instrumental variable.
In each design, the probability of $Z = 1$ given a certain value of the variable $W$ is as follows: $P(Z = 1|W = 0) = 0$, and $P(Z = 1|W = 1)$ is around 0.54, mimicking the empirical application considered in the next section. Moreover, in each design, $X$ and $W$ are independent.
We set the censoring variable, $C$, to be distributed as an exponential variable with rate (reciprocal of the mean) $\lambda$. We then consider different levels of censoring, specifically 20% and 40%, by setting $\lambda = 0.43$ and $\lambda= 1.15$, respectively, for the discrete Bernoulli design, $\lambda = 0.30$ and $\lambda = 0.82$, respectively, for the continuous uniform design, and $\lambda = 0.33$ and $\lambda = 0.87$, respectively, for the continuous Beta design. Finally, with $Y=\min(T,C)$ and $\delta = I(T\le C)$ we generate an iid sample of size $n = 500,1000$ observations having the same distribution as $(Y,\delta,Z,X,W)$.
We employ the Epanechnikov kernel to smooth the variable $X$ when estimating $F^x$ and, for each $X=x$, we select the bandwidth with a data-driven direct plug-in estimate of the Integrated Mean Squared Error-optimal bandwidth. This corresponds to the default bandwidth selection proposed in the R package nprobust, for which we refer to calonico2019nprobust. Our simulation studies showed that smoothing in the argument $t$ during the estimation procedure for $F^x$ leads to a little deterioration of the result and significantly increases the estimation time. Therefore, we consider smoothing in the argument $t$ to be a theoretical requirement that permits to show the asymptotic normality of the estimator, but that can be safely avoided during the method's application. Similarly, we chose to use the Epanechnikov kernel for our simulations, which is a function of order $\nu =2$, even though Assumption (ref) specifies an order of $\nu \ge 4$ for the kernel, due to its widespread acceptance in the field.
The upper bound of $C$ is infinite, allowing for any value of $\bar U\in(0,1)$ to be theoretically chosen, and so we set $\bar U = 0.9$. As we explain in Appendix (ref), we can choose $U_g^c$ having a degenerate distribution with a unique point mass at $\bar U$.
The results given in Table (ref) are based on $N=500$ replications. We present the average bias, standard deviation, and Mean Squared Error (MSE) for each component of $\beta = (\beta_z,\beta_x)^\top$. Additionally, we report the coverage of $95\%$ bootstrap confidence intervals. The intervals are constructed by estimating the standard deviation of the proposed estimator based on a naive bootstrap resampling procedure with replacement, and then plugging in the standard deviation in a normal approximation of the bounds of confidence intervals. To accelerate the simulations, we employed the method proposed in giacomini2013warp. To evaluate the estimator on a component-wise basis, we also include the Root Mean Squared Error (RMSE) as a measure to assess the overall power of the estimator. The value is computed following the formula $RMSE = \sqrt{N^{-1}\sum_{j=1}^N \|\hat\beta^{(j)}-\beta_0\|^2}$, where $\hat\beta^{(j)}$ is the estimation of simulation $j$. For comparison purposes, the same quantities are also given for the partial likelihood estimator.
The results show that the proposed estimator has low bias. Additionally, as the sample size increases, the bias, and coverage of the estimator approach their theoretical values of 0 and 0.95 respectively. The performance of the estimator improves when the proportion of censored observations is lower. It can be argued that the estimator performs better under the discrete Bernoulli design, which may be attributed to the simpler structure of this design, as the exogenous component is discrete rather than continuous. As expected, the standard partial likelihood estimator is biased.
The Illinois Unemployment Incentive Experiment was a controlled social experiment conducted by the Illinois Department of Employment Security in 1984 to evaluate if cash bonuses reduce the duration of unemployment. New claimants for Unemployment Insurance (UI) were randomly assigned to one of three groups: the Job Search Incentive Experiment group (JSIE), the Hiring Incentive Experiment group (HIE), or the control group. The JSIE group was eligible for a \$500 bonus if they found a job of at least 30 hours per week within 11 weeks of the start of their unemployment period and held the job for 4 months. The Hiring Incentive Experiment group had the same eligibility requirements, but the \$500 bonus was given to the hiring company instead. A detailed description of the experiment can be found in woodbury1987bonuses.
It is worth noting that to be part of one of the three aforementioned groups, claimants must be between 20 and 55 years old and have a valid unemployment insurance claim. The duration of unemployment $T$ was recorded as the number of weeks in which participants received unemployment benefits, resulting in a discrete data set. However, the true unemployment duration is continuous but only observed by intervals. This creates a partial identification issue, not within the scope of this paper. The data is treated as continuous. It is important to note that unemployment insurance is granted for 26 weeks, so observations can only be made up until the end of the UI claim and are therefore right censored at 26 weeks. Approximately 40% of the observations are censored.
To clarify the connection between our model and the experiment, we examine the impact of the cash bonus in the JSIE (or HIE) experiment. Although group assignment is randomly determined, participant agreement to participate in the experiment is not independent. This is because, to participate in the program, job seekers must also read the experiment description and sign an agreement form at the start of their follow-up, and this can be influenced by the participant’s motivation to find work or by other personal attributes, such as specific skills. A framework with two treatments, as indicated by the following variable $Z=(Z_1,Z_2)^\top$, will be considered:
Therefore, $Z=(0,0)^\top$ corresponds to the control group, and $Z=(1,0)^\top$, $Z=(0,1)^\top$, correspond to the JSIE and HIE group, respectively. To address selection bias we utilize the group assignment variable $W$ as instrumental variable:
We consider for this analysis a subset of 1,543 non-white women out of the 12,101 available individuals and we report in Table (ref) the sample sizes for each $(W, Z)$ combination. JSIE and HIE experiments have refusal rates of 20% and 35%, indicating significant selection biases.
We use the proposed method to estimate the vector of proportional hazards coefficients. Similar to Section (ref), we estimate the conditional distribution function using the Epanechnikov kernel, and the same selection method for the bandwidth using data-driven direct plug-in estimates of the Integrated Mean Squared Error-optimal bandwidth. To avoid numerical problems, we set $\bar U = 0.5$ and $\bar T = 26$. Table (ref) displays the estimation results of the proposed estimator compared to the partial likelihood estimator, where the standard deviation is computed using 500 bootstrap resamples with replications, and the confidence intervals are constructed using the estimated standard deviation and a normal approximation.
The results in Table (ref) show a remarkable difference between the proposed estimator and the standard partial likelihood estimator. Specifically, the estimated values for $\beta_{z,1}$ and $\beta_{z,2}$, using the proposed estimator, are both positive and significant at the 95% confidence level. This provides empirical evidence that both treatments increase the hazard rate of finding a job, that is they reduce unemployment duration. Conversely, the values for these coefficients estimated using the standard method are close to zero and not significant at the same confidence level, leading to conclude that there is no evidence of any effects of the treatments on the duration.
Additionally, the estimated value for $\beta_{x}$ using the proposed procedure is positive but not significantly different from zero at the 95% confidence level. This could be interpreted as there being no evidence of the effect of the age of non-white women on the duration of unemployment. In contrast, the estimated value for that coefficient is negative and significantly different from zero when the standard partial likelihood estimation procedure is used. This would mean that the older the woman, the more difficult it is to find a job.
In conclusion, the proposed estimator and the partial likelihood estimator yield different results leading to notably different cause-effect interpretations. It could be argued that our estimator provides a more reliable estimation of the regression coefficients compared to the partial likelihood estimator, especially when assessing the influence of the treatments on the duration outcome variable.
The authors thank Gerda Claeskens, Elia Lapenta and Juan Carlos Pardo-Fern´andez for comments that improved the paper. Jad Beyhum undertook most of this work while employed by CREST, ENSAI. Ingrid Van Keilegom acknowledges support from the FWO and F.R.S.-FNRS under the Excellence of Science (EOS) programme, project ASTeRISK (grant No. 40007517).