A Theory Guide to Using Control Functions to Instrument Hazard Models

A Theory Guide to Using Control Functions to Instrument Hazard Models

\thispagestyle{empty}

\setcounter{page}{1}

\begingroup

Introduction

Hazard models (also known as survival models or duration models) specify the rate that events occur conditional on the process history (i.e., which events have occurred or not occurred and when). They are most commonly used in fields such as biostatistics, but have also found purchase in economics and finance, such as in household finance, labor economics, and health economics. Hazard models may also prove highly useful to the rapidly growing new field of genoeconomics -- the use of genes as instruments (also known as Mendelian randomization in biostatistics) to study economic contexts with instrumental variable methods. This modeling approach underpins causal inference for phenomena involving stochastic state transitions. Such phenomena include those for which the randomness is inherent -- e.g., the duration of unemployment spells is intrinsically random as a result of search-and-matching -- and those for which there are unobserved factors modeled as random shocks -- e.g., loan default can be modeled as a deterministic decision, but one that is made based on random factors observed by the agent but not the researcher. (See Kiefer, 1988, for some more examples of economic applications.)

Like many models however, hazard models can suffer from endogeneity, which can cause severe inconsistency and bias in parameter estimates. This is a common problem in practice -- endogeneity can arise due to measurement error in the regressors,\footnote{Strictly speaking, measurement error in the regressors is a form of endogeneity. This is because the measurement error, with which the mismeasured regressor is correlated, ends up in the composite error term -- introducing measurement error in a regressor therefore causes the mismeasured regressor to be correlated with the new, composite error term.} simultaneity, and omitted variables. Unfortunately, instrumental variable (IV) solutions to this endogeneity problem are non-trivial when the instrumented model is nonlinear. The naive generalized version of 2SLS, 2SPS -- two-stage predictor substitution, which consists of replacing endogenous regressors with their predicted values, is not generally valid in a nonlinear model. In this paper, I present a solution to this problem that uses the control function approach (also known in statistics as 2SRI -- two-stage residuals inclusion), which can be thought of as an alternative generalization of 2SLS.\footnote{The control function approach in a linear model is actually exactly equivalent to 2SLS when the same first-stage regressions are used.}

I focus on the use of the control function approach in a popular discrete-data analog of the continuous-time Cox proportional hazards model -- the Prentice and Gloeckler (1978) model -- though, the methods and discussion developed here can easily be applied to other hazard models (and binary choice models).\footnote{By “discrete data”, I mean both grouped-time data -- time-discretized data, such as monthly panel data -- as well as true discrete-time data -- e.g, data on the cycles of some process.} Proportional hazards models are the most commonly used hazard models because they can be distribution-free (i.e., where the baseline hazard is not assumed to follow any particular distribution) yet simple to estimate.\footnote{Other common types of models are additive, accelerated failure time, and proportional odds hazard models.} On the other hand, discrete-data hazard models are much less common than continuous-time models, but are important because they appropriately model the discrete nature of discrete data, which is by far the most common format of data in economics and many other areas. In contrast, continuous-time models that adjust for ties are just approximations to discrete-data hazard models; these inherently introduce misspecification bias. Obviating this misspecification bias motivates my focus on discrete-data hazard modeling.

The structure of the paper is as follows. In Section (ref), I contextualize the paper by explaining its place in and contribution to both the literature on hazard models with endogenous regressors as well as the literature on the control function approach. Section (ref) presents an introduction to the model, a discussion of the asymptotic theory (including the interpretation and selection of the control function as well as robustness/asymptotic inefficiency), and extensions of the core asymptotic theory. In Section (ref), I introduce and explain my Stata command, ivcloglog, which allows practitioners to easily instrument the Prentice and Gloeckler (1978) model via the control function approach. Finally, I conclude in Section (ref) by summarizing the the main findings of the paper. Appendix (ref) contains a theoretical proof of Theorem (ref), and Appendix (ref) contains a primer on grouped-time survival analysis.

Literature Review

In this section, I first present a survey of the endogenous regressor hazard model literature and then explain how my paper contributes to the econometrics and statistics literatures.

Hazard Models with Endogenous Regressors

Endogeneity is a familiar problem for statistical analyses, and the problems that it can create for hazard models are well-known (Prentice, 1982). Yet, the literature for non-additive hazard modeling often assumes exogenous regressors,\footnote{Additive hazard models, of course, are much easier to handle because they can be expressed as linear regression models, and therefore can be instrumented using 2SLS. For example, see J.\ Li et al.\ (2015), Martinussen, Vansteelandt, et al.\ (2017), and Tchetgen Tchetgen et al.\ (2015). (Also see Ying et al., 2019, for a paper that uses control functions instead.)} and works that do allow for endogenous regressors typically only allow for mismeasured regressors, e.g., Y.\ Li and Ryan (2004), Huang and C. Y. Wang (2000, 2006, 2018), and Song and C.\ Y.\ Wang (2014).\footnote{In such contexts, repeated observations of the mismeasured regressors are valid instruments. This is because, although they are subject to measurement error, the measurement error terms for each of the repeated observations are mutually independent (by assumption).} Works that allow for more general endogenous regressors -- “true” instrumental variable methods that go beyond just mismeasured regressors -- are much less common. Some of these are: Bijwaard (2008), Kianian et al.\ (2021), Loeys and Goetghebeur (2003), Loeys et al.\ (2005), Martinussen, S{\o}rensen, and Vansteelandt (2019), and L.\ Wang et al.\ (2022) -- see Choi and O'Malley (2017) for an overview of a bunch of IV methods. However, none of these provide a solution comparable in usefulness to that of 2SLS in linear models.\footnote{Another important (but less critical) point is that all of the mentioned works only consider continuous-time models and are therefore not suitable for discrete data.}

Many such papers make highly restrictive assumptions, severely reducing their applicability. This often entails requiring that the instruments and/or the treatments are binary variables. Indeed, Loeys and Goetghebeur (2003) and Loeys et al.\ (2005) require both to be binary, and L.\ Wang et al.\ (2022) and Kianian et al.\ (2021) require the instruments to be binary. The estimators in such works can also have some poor properties. For example, consider Kianian et al.\ (2021), who apply IPW to survival analysis.\footnote{IPW stands for inverse probability weighting, inverse propensity weighting, or inverse propensity-score weighting.} Developed by Horvitz and Thompson (1952), IPW is where the estimated probability of each treatment status (“propensity score”) is used to construct weights for observations, allowing the estimation of statistics corresponding to a (pseudo-)population with certain characteristics. For example, an average treatment effect can be calculated by comparing the outcomes of treated and untreated populations. The basic (i.e., singly robust) IPW methods involve using inverted propensity scores as weights. Trivially, this means that they are very sensitive to estimated propensity scores being close to or exactly zero -- this is the sample equivalent of the common support assumption being violated. Propensity scores close to zero result in extremely large weights, causing instability in the IPW estimates. This problem is pervasive enough to have motivated the development of stabilized IPW, where the weights are normalized or regularized (Robins, Hern\'an, Brumback, 2000), and doubly robust methods (i.e., AIPW -- augmented IPW), where outcome regression (i.e., “regular” regression methods) is combined with IPW. (See Kang & Schafer, 2007, for an overview of the latter.) Unfortunately, Kianian et al.\ (2021) use the Abadie (2003) IPW-type estimator, which is not stabilized, hampering their estimator's usefulness.

Much less restrictive proposals are given by Bijwaard (2008) and Martinussen, S{\o}rensen, and Vansteelandt (2019). Martinussen, S{\o}rensen, and Vansteelandt (2019) seem to have independently discovered essentially the same idea as Bijwaard (2008), so I will focus on discussing only the latter. Bijwaard (2008) is a seminal work that introduces his IVLR (“Instrumental Variable Linear Rank”) estimator; an excellent summary is Bijwaard (2009).\footnote{Specifically, the IVLR estimator covers the GAFT (generalized accelerated failure time) model, which nests as special cases the AFT (accelerated failure time) model as well as the Cox model and the MPH (mixed proportional hazard) model. (The MPH model is just a Cox model with “frailty” terms -- “random effects” in economics terminology.)} “Linear rank” refers to a rank test, an equivalent of the score test for coefficient significance in the Cox model. Ranks pop up here because they appear in the partial likelihood for a Cox model and in the likelihoods for other continuous-time models. The idea behind IVLR is that the instruments should not affect the hazard and that, based on this fact, a vector of rank test statistics will be centered at zero and asymptotically normal when evaluated at the true coefficient parameter. We can then get the coefficient estimates using what Bijwaard (2008) calls the “inverse rank estimation” approach of Tsiatis (1990): using the rank test statistics as moment functions by setting the sample versions of them to zero.\footnote{Bijwaard (2008) and Tsiatis (1990) refer to these moment functions as \textquotedbl estimating equations\textquotedbl .}

Bijwaard (2009) points out the IVLR estimator is similar to GMM and that there is a choice of weighting matrix. However, his IVLR estimator is not just similar to a GMM estimator -- it is a GMM estimator. I point this out because the ability to choose the weighting matrix means that IVLR is actually an IV-GMM estimator for hazard models. In particular, it implies that the standard properties for IV-GMM estimators carry over to his IVLR estimator. Couching Bijwaard's (2008) IVLR estimator in the IV-GMM framework allows us to understand it much more deeply. He specifically proposes that IVLR be implemented as efficient IVLR, where (an estimate of) the asymptotically efficient weighting matrix is used -- this makes efficient IVLR a feasible efficient IV-GMM estimator.

However, using feasible efficient IV-GMM entails robustness being sacrificed for the sake of inferential efficiency: it is well-known that there is an inherent tradeoff between efficiency and robustness for IV-GMM and that feasible efficient IV-GMM sacrifices the latter for the former. Therefore, I argue that an alternative instrumental variables technique for hazard models, more robust but less efficient than efficient IVLR, would be more practical, applicable, and attractive for empirical studies. After all, this tradeoff is especially well studied in linear models, where a strong preference for robustness over efficiency has developed, as evidenced by the fact that 2SLS is much more popular than feasible efficient IV-GMM. Thus, I propose using the control function approach instead, based on the fact that it satisfies these criteria -- standard results in the literature.

The Control Function Approach

This comprehensive theory guide is the first to provide a complete and rigorous theoretical foundation for the use of the control function approach to instrument hazard models, undergirding existing work. The first proposals on using the control functions in hazard models were made by Terza et al.\ (2008), Atiyat (2011), and Palmer (2013): Terza et al.\ (2008) and Atiyat (2011) proposed using the control function approach to instrument continuous-time hazard models, and Palmer (2013) proposed using the control function approach to instrument the Prentice and Gloeckler (1978) discrete-data model. Based on Palmer (2013), multiple applied works have used control functions to instrument the discrete-data Prentice and Gloeckler (1978) model (Coviello et al., 2015; Ganong & Noel, 2022; Liebersohn & Rothstein, 2023; Palmer, 2023). Similarly, inspired by Terza et al.\ (2008) and Atiyat (2011), there are also a number of applied works that have used control functions to instrument continuous-time models (e.g., Gore et al., 2010). Yet, a critical and pressing problem is that a complete theoretical justification is missing from the literature, as pointed out by Palmer (in personal communication) and by others in the literature (Choi & O'Malley, 2017).

Atiyat (2011) is a theory paper that explores the use of control functions in (continuous-time) hazard models, but ultimately does not provide a theoretical justification for their use. Specifically, he assumes that control functions can be used to instrument endogenous regressors but does not justify this or investigate when this is valid. Instead, Atiyat (2011) points towards Terza et al.\ (2008) for justification of this assumption

Unfortunately, Terza et al.\ (2008) provide little justification themselves. Firstly, they make the assumption that the first-stage residuals exactly equal any relevant unobserved confounders, which is an extremely strong assumption that cannot be justified in empirical applications! Secondly, their theoretical discussion of the control function approach is very brief and incomplete. In particular, they do not derive or present critical conditions that are required for the control function approach to be generally valid, such as what I call the instrument conditional independence assumption -- the control function equivalent of the 2SLS exclusion restriction!\footnote{The instrument conditional independence assumption is a generic requirement. An exception is that, when the second stage is probit, a weaker, conditional mean restriction can be used instead (Rivers & Vuong, 1988).} As another example, they do not prove how and under what conditions (e.g., the Lindeberg condition) their discussion would carry over to survival contexts.

Palmer (2013) proposed the use of the control function approach in discrete-data hazard models and utilized this to study loan default. However, the theoretical justification was heuristic and the task of developing rigorous theory on instrumenting hazard models was left as future work. This paper serves to fill in that gap in the literature. Moreover, this paper serves as a theory guide that complements and theoretically justifies Palmer (2023), which is an updated version of Palmer (2013) that also contains an empirical guide as well as a summary of the theory introduced in this paper.

This paper acts as a cornerstone for existing and future literature. It provides a rigorous theoretical basis for and advances the theory around the use of the control function approach to instrument hazard models, contributing to a nascent literature and validating multiple existing empirical studies. In addition, I also advance and broaden the theory around the use of control functions in discrete choice models -- after all, discrete-data hazard models can nominally be expressed as discrete choice models (Allison, 1982; Jenkins, 1995).\footnote{I focus on the Prentice and Gloeckler (1978) discrete-data proportional hazards model, but the methods and discussion developed here can be easily applied more generally to other hazard models and discrete choice models by modifying the primary model and associated control function specification assumption.} Existing theory papers on this mainly have a narrow focus on probit or Poisson regression due to convenient properties that are not found elsewhere (e.g., see Wooldridge, 2015); this paper relaxes the dependence on these properties. I also bring formal theoretical clarification to the Prentice and Gloeckler (1978) model by pointing out some of the unstated assumptions, such as the Lindeberg condition.

Theory

In this section, I introduce the version of the Prentice and Gloeckler (1978) discrete-data proportional hazards model instrumented via the control function approach. I discuss the Prentice and Gloeckler (1978) model specifically because it is the most popular discrete-data proportional hazards model and has certain advantages. Namely, it is simple, parsimonious, and easy to implement, and it is also flexibly parametric -- i.e., it is a distribution-free parametric model. Specifically, the Prentice and Gloeckler model estimates the baseline hazard with a stepped function, and this behavior is represented by the time fixed effects.\footnote{The model is not semi-parametric, but is similar in that it is distribution-free. On this basis, Han and Hausman (1990) also refer to flexibly parametric estimation as \textquotedblleft nonparametric\textquotedblright in a broad sense, but this usage has not been popularized.}

I then present sufficient conditions for consistent estimation and inference as well as a consistent estimator of the asymptotic variance. I also discuss practical considerations, such as the interpretation and selection of the control function and the robustness/asymptotic inefficiency of the estimator. After these is a discussion of extensions to multiple/transformed/discrete first-stage dependent variables.

The IV Prentice and Gloeckler (1978) Model with Control Functions

Model Overview

$l_{1}$ and $l_{2}$ denote the log-likelihoods corresponding to the first and second stages, respectively; their arguments are suppressed for brevity. $y_{it}$ denotes the outcome; here, it is an indicator for event occurrence -- a dummy variable that is 1 in the time period a failure occurs and 0 otherwise. $x_{it}$ denotes the continuous endogenous regressor. Let $z_{it}\coloneqq(z_{1it}',z_{2it}')$', where $z_{1it}\coloneqq(\varphi_{t}',\mathfrak{z}_{it}')$ is a transposed row vector of exogenous regressors and $z_{2it}$ is a transposed row vector of “excluded instruments” (i.e., “real” instruments that aren't the control variables $z_{1it}$). $\varphi_{t}$ represents time-period dummies, and $\mathfrak{z}_{it}$ represents other added control variables. As per Prentice and Gloeckler (1978), we implicitly assume that all regressors and instruments are constant within each time period, which is reflected in the notation.

The primary model is \[

\] where $u_{it}$ is a composite error term with some unknown distribution. $u_{it}$ represents a combination of unobserved heterogeneity that is potentially correlated with $x_{it}$ (which could induce omitted-variable bias) and the “true” error term, which has a standard Gumbel distribution. The primary model here is just the equivalent nominal cloglog (complementary log-log) representation of the Prentice and Gloeckler (1978) grouped-time proportional hazard model, albeit in a latent variable form.\footnote{Any Bernoulli model with $\Pr(y_{1}=1|x)=F_{-\upsilon}(x'\beta)$ for some function $F_{-\upsilon}$ has a latent index representation$y_{1}=\boldsymbol{1}(y_{1}^{*}>0)$, where $y_{1}^{*}\coloneqq x'\beta+\upsilon$ and $\upsilon$ has CDF $F_{\upsilon}(\cdotp)$ satisfying $F_{\upsilon}(a)\equiv1-F_{-\upsilon}(-a)$. In other words, we simply construct $\upsilon$ such that $F_{-\upsilon}(\cdotp)$ is the CDF of $-\upsilon$. This latent index representation is mechanically identical to the original Bernoulli model because $\Pr(y_{1}=1|x)=\Pr(y_{1}^{*}>0|x)=\Pr(-\upsilon<x'\beta|x)=F_{-\upsilon}(x'\beta).$ In the cloglog model, we have $F_{-\upsilon}(a)=1-e^{-e^{a}}$, which gives us $F_{\upsilon}(a)=e^{-e^{-a}}$, the CDF of the standard Gumbel distribution. ($\upsilon$ and $a$ are placeholder variables.)}\textsuperscript{,}\footnote{The Prentice and Gloeckler (1978) model has an equivalent nominal cloglog representation where we nominally treat each observation as IID. In general, discrete-time hazard models are identical to a discrete choice model with (1) an appropriate link function, (2) an indicator of the event being the outcome, (3) the regressors augmented by time fixed effects, and (4) with all observations for each entity only being kept until “failure” or censoring (Allison, 1982; Jenkins, 1995). However, the observations are not truly IID: observations for entity $i$ are cut off after “failure”. In this sense, the model is not a true cloglog model.} $\psi_{t}$ is a time-period fixed effect, and $(\beta_{1}',\beta_{2})'$ is the parameter vector of interest.

We can introduce the decomposition \[ u_{1it}\eqqcolon c(v_{it})+e_{it}, \] where we assume that the function $c(\cdot)$ is a function to be estimated, and that $v_{it}$ is the error of the auxiliary model \[ x_{it}=z_{it}'\pi+v_{it}\equiv z_{1it}'\pi_{1}+z_{2it}'\pi_{2}+v_{it}, \] where $\pi\coloneqq(\pi_{1}',\pi_{2}')'$. The decomposition just represents what happens with the new error term when the control function $c(v_{it})$ is introduced to the second stage model. For the purposes of this paper, we will assume that $c(v_{it})$ is a $Q$-th order polynomial with \[ c(v_{it})=p(v_{it})'\beta_{3}, \] where the vector of control function terms $p(v_{it})$ is defined by $p(v_{it})\coloneqq(v_{it},v_{it}^{2},\cdots,v_{it}^{Q})'$.\footnote{For expositional convenience, $p(v_{it})$ is defined as not containing a $1$ because $c(v_{it})$ would otherwise contain a constant, making it perfectly collinear with the time-period fixed effects (if all time-period fixed effects are included).} Using the decomposition, we have \[ y_{it}=\mathbf{1}(\psi_{t}+\mathfrak{\zeta}_{it}'\beta+e_{it}>0)\equiv\mathbf{1}(\psi_{t}+\mathfrak{z}_{it}'\beta_{1}+\beta_{2}x_{it}+p(v_{it})'\beta_{3}+e_{it}>0), \] where $\mathfrak{\zeta}_{it}\coloneqq(\mathfrak{z}_{it}',x_{it},p(v_{it})')'$ and $\beta\coloneqq(\beta_{1}',\beta_{2}',\beta_{3}')'$. Also, let $\gamma\coloneqq(\psi',\beta')'.$

Under certain conditions,\footnote{The inclusion of all non-collinear second-stage controls in the first stage is important not just for better efficiency but also robustness. Imagine adding some exogenous regressor $\tilde{z}_{1it}$ to the second stage but not the first stage. Since it is not in $z_{it}$, there is no guarantee that $e_{it}|x_{it},z_{1it},\tilde{z}_{1it},p(v_{it})\sim\mathrm{Gumbel}(0,1)$ -- see Appendix (ref). This is just like in 2SLS, where including the exogenous regressor $\tilde{z}_{1it}$ in the first stage guarantees that it will be orthogonal to $v_{it}$. (In 2SLS, $v_{it}$ shows up in the new, composite second-stage error term.) Obviously, we do not need to add $\tilde{z}_{1it}$ to the first stage (which would entail dropping some instruments collinear with $\tilde{z}_{1it}$) if $\tilde{z}_{1it}$ is perfectly collinear with the first-stage variables $z_{1it}$: after all, $v_{it}$ would be exactly the same either way. More formally, using the fact that $(x_{it},z_{it})\raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}} e_{it}|p(v_{it})$ and the perfect collinearity of $\tilde{z}_{1it}$ with $z_{it}$, we have that $e_{it}|x_{it},z_{1it},p(v_{it})=e_{it}|x_{it},z_{it},p(v_{it})=e_{it}|x_{it},z_{it},\tilde{z}_{1it},p(v_{it})=e_{it}|x_{it},z_{1it},\tilde{z}_{1it},p(v_{it})$. Therefore, when $\tilde{z}_{1it}$ is perfectly collinear with $z_{1it}$, we still have that $e_{it}|x_{it},z_{1it},\tilde{z}_{1it},p(v_{it})\sim\mathrm{Gumbel}(0,1)$ even if $\tilde{z}_{1it}$ is omitted from the second stage.} \[ e_{it}|x_{it},z_{1it},p(v_{it})\sim\mathrm{Gumbel}(0,1). \] This just means that the primary model is a cloglog model when the control function terms $p(v_{it})$ are included alongside the original second-stage regressors $x_{it}$ and $z_{1it}$. Therefore, this implies that including the terms of $p(v_{it})$ as controls allows us to recover an exogenous Prentice and Gloeckler (1978) model when the base assumptions for it are satisfied. Consequently, if $p(v_{it})$ were known, we could just do that to get consistent estimates of $\beta_{1}$ and $\beta_{2}$, but we cannot because $v_{it}$ is unobserved.\footnote{Strictly speaking, just like other common binary outcome models, we are actually estimating the scaled parameters (the ratio of the unscaled parameters to the standard deviation of $e_{it}$) -- the unscaled are not identified (Wooldridge, 2010). After all, $y_{it}^{*}$ has a purely ordinal meaning, so the scale of $\beta$ and $e_{it}$ in $y_{it}^{*}=\psi_{t}+\mathfrak{\zeta}_{it}'\beta+e_{it}$ is entirely arbitrary. The convention for probit and logit is to set the variance of $e_{it}$ to 1, causing $e_{it}$ have a standard normal/logistic distribution and the scaled parameters to equal the unscaled. However, the convention for cloglog is instead to set the variance of $e_{it}$ to $\sfrac{\pi^{2}}{6}$, which is the variance of the standard Gumbel distribution. Technically, this means that cloglog point estimates do not truly estimate $\gamma$, but this does not matter (except in simulations).} Instead, if we have an estimate $\hat{\pi}$ of the auxiliary model parameters, we can consider substituting $v_{it}$ with $\hat{v}_{it}\coloneqq x_{it}-z_{it}'\hat{\pi}$ when estimating.\footnote{Strictly speaking, $\hat{v}_{it}$ is not an estimate of $v_{it}$ because $v_{it}$ is not a parameter.} To see why, first substitute the unobserved variable $v_{it}$ for the observed variables in $z_{it}$ and unknown parameter vector $\pi$ using the auxiliary model. This gives us \[ y_{it}=\mathbf{1}(\psi_{t}+\mathfrak{z}_{it}'\beta_{1}+\beta_{2}x_{it}+p(x_{it}-z_{it}'\pi)'\beta_{3}+e_{it}>0). \] This is a classic two-step Z-estimation setup: were it not for $\pi$ being unknown, we would be able to directly estimate $\psi_{t}$ and $\beta$ using MLE.\footnote{M-estimation theory refers to extremum estimation theory in general, whereas Z-estimation theory (coined by Van der Vaart, 1998) is a special case where first-order condition optimization is possible. Classical GMM, such as Hansen (1982) or Newey and McFadden (1994), covers ergodic or IID data, and its theorems are special cases of finite-dimensional Z-theorems. Non-classical GMM simply refers to finite-dimensional Z-theorems that use different CLTs to achieve analogous results. One famous example of a finite-dimensional Z-theorem is Theorem 3.3 in Pakes and Pollard (1989), which is a slightly more general version of Theorem 3 in Huber (1967) that allows the sample moment functions to be non-smooth and even discontinuous. It forms the theoretical basis of the Simulated Method of Moments (SMM), also known as the Method of Simulated Moments (MSM), where moment functions can be simulated to avoid the intractability or non-existence of an analytical form.} Under appropriate assumptions, we can show that we can replace $\pi$ with $\hat{\pi}$ and instead do \emph{QMLE} -- quasi-maximum likelihood estimation.\footnote{Technically, this is specifically \emph{2SCMLE} -- two-stage conditional maximum likelihood estimation (Rivers & Vuong, 1988).}

Notice that, just like in sequential classical GMM (Newey, 1984; Newey & McFadden, 1994), the sequential nature of this estimator means that we can “stack” the moment functions of the first and second stages (namely, the first-stage score and second-stage quasi-score). This is because the just-identification of both stages implies that the vector of moment conditions in each can be solved for exactly, implying that combining them into one vector of moment conditions yields the exact same estimator.\footnote{“Just-identification” is used here in the general sense (one moment condition for each parameter) rather than the more common meaning for IV (the numbers of “excluded instruments” and endogenous regressors being equal).} In other words, we can treat the estimator like a one-step Z-estimator, allowing us to apply a Z-theorem to find sufficient conditions and get results.

Moment Functions

Let $\theta\coloneqq(\pi',\psi',\beta')'$. The moment functions are \[ \hat{m}(\theta)\coloneqq\mathopen\mathclose\bgroup\originalleft(

\aftergroup\egroup\originalright)\equiv\mathopen\mathclose\bgroup\originalleft(

\aftergroup\egroup\originalright)\equiv\frac{1}{n}\stackrel[i=1]{n}{\sum}g_{i}(\theta), \] where \[ \hat{m}_{1}(\theta)\equiv\hat{m}_{1}(\pi)\coloneqq\frac{1}{n}\frac{\partial l_{1}}{\partial\pi}\equiv\frac{1}{n}\stackrel[i=1]{n}{\sum}\stackrel[t=1]{T}{\sum}z_{it}'v_{it}\equiv\frac{1}{n}\stackrel[i=1]{n}{\sum}\underset{g_{1i}(\theta)}{\underbrace{\stackrel[t=1]{T}{\sum}z_{it}'(x_{it}-z_{it}'\pi)}} \] and\footnote{Strictly speaking, Prentice and Gloeckler (1978) actually use the unnormalized score (i.e., $n\times\hat{m}(\theta)$) in their paper. I normalize by $n$ here to emphasize the stochastic convergence of the score.} \[ \hat{m}_{2}(\theta)\equiv\hat{m}_{2}(\psi,\beta,\pi)\coloneqq\frac{1}{n}\mathopen\mathclose\bgroup\originalleft(

\aftergroup\egroup\originalright)\equiv\frac{1}{n}\mathopen\mathclose\bgroup\originalleft(

\aftergroup\egroup\originalright)\equiv\frac{1}{n}\stackrel[i=1]{n}{\sum}\underset{\mathclap{\coloneqq g_{2i}(\theta)\equiv\frac{1}{T}\stackrel[t=1]{T}{\sum}g_{2it}(\theta)}}{\underbrace{\mathopen\mathclose\bgroup\originalleft(

\aftergroup\egroup\originalright)}}, \] where $n$ is the number of entities (indexed by $i$), $T$ is the number of time periods, and $K$ is the number of regressors. The elements of $g_{2i}(\theta)$ are defined by \[ \frac{\partial l_{2i}}{\partial\psi_{t}}=

;\qquad\frac{\partial l_{2i}}{\partial\beta_{k}}=\mathfrak{\hat{\zeta}}_{is,k}\mathopen\mathclose\bgroup\originalleft(\delta_{i}b_{is}-\stackrel[j=1]{s-1}{\sum}h_{ij}\aftergroup\egroup\originalright), \] where $\delta_{i}$ is the censoring indicator, $t$ represents a time period, $s$ is the time period of survival (with subscript $i$ suppressed for clarity), $\mathfrak{\hat{\zeta}}_{it,k}$ is the $k$-th element of $\mathfrak{\hat{\zeta}}_{it}$ (i.e., the observation for the $k$-th regressor, excluding the time fixed effects, for entity $i$ and time $t$$)$, and where \[ b_{is}\coloneqq\frac{h_{it}e^{-h_{it}}}{1-e^{-h_{it}}};\qquad h_{it}\coloneqq e^{\psi_{t}+\mathfrak{\zeta}_{it}'\beta}. \] $\hat{m}_{1}(\theta)$ is a score function, but $\hat{m}_{2}(\theta)$ is technically a quasi-score function because $p(\hat{v}_{it})$ is a vector of generated regressors. $g_{1i}(\theta)$ is the entity-specific first-stage score contribution,$g_{2i}(\theta)$ is the entity-specific second-stage quasi-score contribution, and $g_{2it}(\theta)$ is the observation-specific second-stage quasi-score contribution; $g_{i}(\theta)\coloneqq(g_{1i}(\theta)',g_{2i}(\theta)')'$. $\hat{\theta}\coloneqq(\hat{\pi}',\hat{\psi}',\hat{\beta}')'$ is the solution of \[ \hat{m}(\hat{\theta})=0. \]

Root-$n$ Asymptotic Theory

Here, I present assumptions that are sufficient to infer that standard GMM asymptotic properties hold, such as consistency and asymptotic normality for $\hat{\theta}$.

$G$ has the structure

where $\xi_{it}\coloneqq(\varphi_{t}',\mathfrak{\zeta}_{it}')'\equiv(\varphi_{t}',\mathfrak{z}_{it}',x_{it},p(\hat{v}_{it})')'$ is a transposed row vector of second-stage regressors, $z_{it}$ is a transposed row vector of auxiliary model instruments (i.e., the same first-stage regressors as before), $\gamma$ is the coefficient vector on $\xi_{it}$, and $\mkern3mu\mathchar'26\mkern-12mu d_{it}$ and $d_{it}$ are defined by \[ \mkern3mu\mathchar'26\mkern-12mu d_{it}=-d_{it}\stackrel[q=1]{Q}{\sum}\beta_{3q}v_{it}^{q-1};\qquad d_{it}=

, \] with $Q$ being the order of the control function polynomial and $\beta_{3q}$ being the coefficient on $v^{q}$.

$\widehat{G}$ has the structure

$\tilde{n}$ is the total number of observations -- i.e., the number of pairs $(i,t)$. $\widehat{\Xi}$ and $Z$ are the matrix versions of $\hat{\xi}_{it}\coloneqq(\varphi_{t}',\mathfrak{\hat{\zeta}}_{it}')'$ and $z_{it}$ -- i.e., $\hat{\xi}_{it}'$ and $z_{it}'$ stacked; $\hat{\xi}_{it}$ only differs from $\xi_{it}$ in that it contains $p(\hat{v}_{it})$ rather than $p(v_{it})$. $\mathrlap{\hskip1.1ex\widehat{\vphantom{D}}}\mkern2mu\rule[0.75ex]{0.75ex}{0.06ex}\mkern-8mu D$ and $\widehat{D}$ are defined by \[ \mathrlap{\hskip1.1ex\widehat{\vphantom{D}}}\mkern2mu\rule[0.75ex]{0.75ex}{0.06ex}\mkern-8mu D\coloneqq\mathrm{diag}\mathopen{}\mathclose\bgroup\originalleft(\mathrlap{\hskip1ex\hat{\vphantom{d}}}\mkern3mu\mathchar'26\mkern-12mu d_{it}\aftergroup\egroup\originalright);\qquad\widehat{D}\coloneqq\mathrm{diag}\mathopen{}\mathclose\bgroup\originalleft(\hat{d}_{it}\aftergroup\egroup\originalright), \] with \[ \mathrlap{\hskip1ex\hat{\vphantom{d}}}\mkern3mu\mathchar'26\mkern-12mu d_{it}=-\hat{d}_{it}\stackrel[q=1]{Q}{\sum}\hat{\beta}_{3q}\hat{v}_{it}^{q-1};\qquad\hat{d}_{it}=

. \]

The invertibility of $\widehat{G}$ is implied by the invertibility of $Z'Z$ and $\widehat{\Xi}'\widehat{\Xi}$.

The Interpretation and Selection of the Control Function

Let $\hat{c}(\hat{v}_{it})$ denote the estimated control function, which is composed of the first-stage residuals $\hat{v}_{it}$ fed into some estimate $\hat{c}(\cdot)$ of $c(\cdot)$. In this paper, I have focused on $\hat{c}(\hat{v}_{it})\coloneqq\hat{p}(v_{it})'\hat{\beta}_{3}$, but the discussion here still applies to non-polynomial control functions. Note that the unobserved control function $c(v_{it})$ represents the unobserved heterogeneity in the “true model”. We thus are assuming that the unobserved heterogeneity takes the functional form of $c(\cdot)$ and that it is only a one-dimensional function of $v_{it}$.

The general interpretation of the estimated control function $\hat{c}(\hat{v}_{it})$ is that it is a good proxy for the unobserved heterogeneity $c(v_{it})$. In models where the second stage is linear or probit, a conditional mean restriction is required rather than a conditional independence assumption: consequently, $c(v_{it})$ just needs to flexible enough to have the same mean as the unobserved heterogeneity, conditional on $v_{it}$ (e.g., see Rivers and Vuong, 1988). Then, the estimated control function $\hat{c}(\hat{v}_{it})$ is estimating the mean of the unobserved heterogeneity.

However, that specific interpretation is not true for non-probit (nominally) discrete choice models because a conditional independence assumption is required. Instead, $c(v_{it})$ must have the same distribution as the unobserved heterogeneity; hence, a correct interpretation of $\hat{c}(\hat{v}_{it})$ is that it is an optimal approximation in distribution to the unobserved heterogeneity.\footnote{The control function approach can also be understood as the flexibly parametric counterpart of the Blundell and Powell (2003, 2004) semi-parametric approach, which does not parametrically restrict the relationship between $u_{it}$ and $e_{it}$. (See Imbens & Newey, 2009, for the generalized framework and Wooldridge, 2015, for a summary.)} After all, when we carry out MLE-type estimation, what we are really doing is fitting the (GLM) pseudo-residuals $\hat{e}_{it}$ such that their empirical distribution is an optimal fit to the assumed distribution of $e_{it}$, with the norm being the Kullback-Leibler divergence.\footnote{Equivalently, with the estimated control function included in the fitted model, we attempt to make the hazard function as close to proportional as possible in some metric. See (ref).} By the additivity property of the Kullback-Leibler divergence and the fact that the distribution of $e_{it}$ does not depend on $c(v_{it})$, we see that $\hat{c}(\cdot)$ is fitted such that $\hat{c}(\hat{v}_{it})$ is an optimal approximation in distribution to $c(v_{it})$ in the Kullback-Leibler sense, subject to the functional form restrictions of $c(\cdot)$.

A key takeaway is that more flexibility in the functional form of $c(\cdot)$ can increase improve accuracy by allowing $\hat{c}(\hat{v}_{it})$ to better approximate the distribution of the unobserved heterogeneity. Though, too flexible a functional form can result in concerns of overfitting and the curse of dimensionality; thus, the practitioner is looking for a sweet spot in the middle. This is essentially a machine learning-style dictionary selection problem, and there are many techniques from the machine learning literature that can be used to help with this selection process, e.g., cross-validation.

Finally, the other key takeaway is that this implies that an appropriate robustness check is to examine how $\hat{\theta}$ changes in response to different functional forms being specified for $c(\cdot)$, e.g., by varying the degree of a polynomial control function. After the selection process, if $\hat{\theta}$ is not sensitive to slight changes to the functional form, then this suggests that the functional form of $c(\cdot)$ is near-optimal.

Asymptotic Inefficiency

Control function methods are asymptotically inefficient (this includes 2SLS) -- compared with feasible efficient IV-GMM methods, they trade efficiency for robustness. The parallels of QMLE to true MLE can be leveraged to give an intuitive and in-depth explanation of the cause of this inefficiency.

Rather than a joint likelihood, we use a conditional likelihood. In other words, we see that the conditional information matrix does not reflect all of the information available: we pass only the first-stage estimates into the second stage, and any other information about the excluded instruments is lost from the perspective of the second stage. This intuitively explains why the second-stage asymptotic variance matrix would be inflated.

I give a more in-depth explanation here. We effectively conduct constrained optimization of the true likelihood when we do QMLE. Note that adding more information tends to sharpen (i.e., increase the magnitude of the curvature) of the objective function; thus, we are constrained in not using that discarded information to sharpen the objective function. Therefore, the magnitude of the curvature of the objective function is limited. Now, note that $G_{22}$ is the expected gradient of the second-stage quasi-score -- the negative of the expected Hessian of the second-stage conditional log-likelihood (i.e., second-stage quasi-log-likelihood), and thus embodies the negative of the curvature. The limiting of the curvature therefore limits the size of the determinant of $G_{22}$ since the determinant of a Hessian, evaluated at a critical point of a function, equals the Gaussian curvature at that point. Because the conditional information matrix equality implies that the second-stage asymptotic variance equals $G_{22}^{-1}$, using the fact that $det(A^{-1})=\sfrac{1}{det(A)}$ for an arbitrary matrix $A$, we can see that this ultimately results in the (root-$n$-scaled) asymptotic variance matrix being inflated in the sense that its determinant is limited in how close to zero it can get.

Extensions

Multiple Endogenous Regressors

Under very similar assumptions, it is trivial to extend Theorem (ref) to allow for multiple endogenous regressors (strictly speaking, multiple first-stage dependent variables), i.e., when $x_{it}$ is a vector and there are multiple first-stage regressions -- one for each element of $x_{it}$. The only noteworthy assumption here is Assumption (ref). The overall, multivariate control function $c(v_{it})$ is allowed to simply be composed of the separate, univariate control functions for each first-stage regression, i.e., just having separate control functions for each each endogenous regressor. However, note that $c(v_{it})$ is also allowed to be a non-separable function of all first-stage residual terms -- e.g., a multivariate polynomial with all residual terms as arguments. Thus, when there are multiple endogenous regressors, there is additional flexibility in the choice of functional form for the control function $c(v_{it})$.

The next thing to note is what the asymptotic variance $V\equiv G^{-1}\Omega(G^{-1})'$ and its estimator $\widehat{V}\equiv\widehat{G}^{-1}\widehat{\Omega}(\widehat{G}^{-1})'$ become. $\Omega$ and $\widehat{\Omega}$ have the same structure as before, but $G$ and $\widehat{G}$ do not. Denoting the number of first-stage dependent variables by $\kappa$, $G$ is given by

On the other hand, $\widehat{G}$ is given by

where $Z^{(1)}$ is the matrix of instruments for the first auxiliary equation (i.e., for the first endogenous regressor) and $Z^{(\kappa)}$ is the matrix of instruments for the last auxiliary equation (i.e., for the last endogenous regressor).

When identical instruments are used for all first-stage regressions, the expression simplifies to \[ \mathopen\mathclose\bgroup\originalleft[

\aftergroup\egroup\originalright], \] where $I_{\kappa}$ is an identity matrix, $\boldsymbol{1}_{\kappa}'$ is a row vector of ones, and $\otimes$ represents the Kronecker product.

Discrete Endogenous Regressors That Are Functions of Continuous Variables

In the primary model, the endogenous regressors $x_{it}$ can actually be replaced by $f(x_{it})$ for any (measurable, potentially vector-valued) function $f$ without needing to change the auxiliary model or the control function: Theorem (ref) still carries over.\footnote{This function of $x_{it}$ being measurable just means that it is a random variable itself.} This property is characteristic of the control function approach and applies more generally to other models (e.g., linear or probit). Ultimately, this is because $(x_{it},z_{it})\raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}} e_{it}|p(v_{it})\implies(f(x_{it}),z_{it})\raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}} e_{it}|p(v_{it})$: it guarantees that $f(x_{it})$ will also be conditionally independent of the new error term when the control function is added as a regressor, just like $x_{it}$. In other words, when your endogenous regressors are $f(x_{it})$, you only need to run one first-stage regression for each element of $x_{it}$.

As an example, suppose that your model has five endogenous regressors of the form $x_{1it}$, $x_{1it}^{2}$, $x_{2it}$, $x_{2it}^{2}$, and $x_{1it}x_{2it}$, for some observed scalar, continuous variables $x_{1it}$ and $x_{2it}.$ You only need to run two first-stage regressions -- one for $x_{1it}$ and one for $x_{2it}$ -- and include the resulting two control functions (or one non-separable, multivariate control function).

Note that the untransformed $x_{it}$ is allowed to be absent from the primary model. This can particularly useful -- for example, suppose that your model has one endogenous regressor, $f(x_{it})=\boldsymbol{1}(x_{it}>0)$, for some observed scalar, continuous variable $x_{it}$. If we want to directly instrument $f(x_{it})$, since $\boldsymbol{1}(x_{it}>0)$ is a binary variable, the first stage is not allowed to be linear (see below). We could estimate a nonlinear model for the first stage and use the residuals from that to construct a control function, but this requires stronger assumptions (Wooldridge, 2015). Instead, because $x_{it}$ is observed rather than latent, you can just use $x_{it}$ as a first-stage dependent variable rather than $\boldsymbol{1}(x_{it}>0)$. Then, construct a control function as normal. For an example of this, see Palmer (2013, 2023), who does this to instrument an indicator for a loan being underwater. A loan is underwater if its loan-to-value ratio exceeds 100%, so the loan-is-underwater indicator can be instrumented by using loan-to-value ratio as a first-stage dependent variable.

Linear versus Nonlinear Auxiliary Models

It is a well-known fact that we can always run a linear first stage if the second stage is linear (Angrist & Pischke, 2009; Kelejian, 1971). Even if the true auxiliary model (“first stage”) is nonlinear, assuming that the true primary model (“second stage”) is linear, under typical 2SLS assumptions, the resulting estimates will still be consistent. This is why running a nonlinear first stage in such a scenario is considered a “forbidden regression” (Angrist & Pischke, 2009; Hausman, 1975) -- if you did so, the resulting estimates would be much less robust to misspecification since their consistency is only guaranteed under strong assumptions. Unfortunately, when the primary model is nonlinear, this result no longer applies. Given that the second stage is cloglog, this means we must carefully consider the form of the first stage.

Using a linear first stage is not necessary if the first-stage dependent variables are continuous. So, why is this assumed in the Section (ref) model? Well, as with all two-stage IV methods, we want an estimator of the conditional expectation function of the first stage in order to exploit some sort of orthogonality condition (here, conditional independence) between the instruments and the second-stage error term. OLS is an estimator with the most parsimonious functional form: OLS provides the best linear approximation to the conditional expectation function. So, even if the true first stage is not truly linear, OLS will often provide a good approximation to reality when the first-stage dependent variables are continuous.

On the other hand, if the first-stage dependent variables are discrete, then we know that a linear (i.e., a linear probability model) is a poor approximation to reality if the predicted probabilities are substantially different from 0.5. This motivates the use of a nonlinear first stage. Theorem (ref) can easily be extended to allow for this: because the form of the first stage is not central to the proof (see Appendix (ref)), a nonlinear first stage is fine in principle -- we just need to make additional assumptions. For more detail, see Wooldridge's (2015) work on the control function approach for discrete first-stage dependent variables.

Stata Implementation

An Automated Stata Command

The GMM estimation procedure described in Section (ref) can be easily and conveniently implemented using the Stata command ivcloglog, which is available from Boston College's SSC (Statistical Software Components archive). As with all commands available through SSC, Stata users can install ivcloglog by simply running ssc install ivcloglog in Stata. (\texttt{\textbf{ssc}} is a built-in Stata command.)

Before covering how to estimate the instrumented version of the Prentice and Gloeckler (1978) model, I will first give an overview of how to estimate the original, basic version (i.e., with all regressors exogenous). Since the Prentice and Gloeckler (1978) model can be nominally expressed as a cloglog model with time fixed effects (Allison, 1982; Jenkins, 1995), the setup required for estimating it in Stata is very simple:

(For a user-created command, pgmhaz8, that automates these steps and estimates a model with exogenous regressors, see Jenkins, 2004.) The process for estimating a model with endogenous regressors is identical, except with the use of the command ivcloglog instead of cloglog.

ivcloglog checks for collinear variables (in the instruments $Z$ and augmented second-stage regressors $\widehat{\Xi}$) or perfect predictors (in the augmented second-stage regressors $\widehat{\Xi}$) and excludes them from the estimation. Using the remaining variables, ivcloglog obtains the point estimates using regress (for the first-stage) and \texttt{\textbf{cloglog}} (for the second-stage) and then calculates an appropriate variance-covariance matrix estimate (VCE) that covers both the first and second stages using the formulae detailed in Section (ref). (The OLS first stage means that all endogenous regressors are assumed to be either continuous or a known, possibly discontinuous function of the first-stage dependent variables.) By default, the VCE is simply obtained through the \texttt{\textbf{gmm}} command.

However, Stata may fail to do so through the gmm command due to numerical difficulties. In particular, although the earlier collinearity-checking guarantees that $\widehat{G}$ is exactly invertible (see Part 2B), if $\widehat{G}$ is highly near-singular, Stata may nonetheless declare that $\widehat{G}$ is numerically singular and refuse to obtain the VCE $\widehat{V}=\widehat{G}^{-1}\widehat{\Omega}(\widehat{G}^{-1})'$.\footnote{Technically, it is more efficient to find $\widehat{V}$ by using LU decomposition to solve $\widehat{G}\widehat{V}\widehat{G}'=\widehat{\Omega}$ for $\widehat{V}$ rather than directly calculating $\widehat{V}=\widehat{G}^{-1}\widehat{\Omega}(\widehat{G}^{-1})'$, and this is presumably what Stata actually does. In any case, Stata will not obtain $\widehat{V}$.} For example, this can occur when the dataset has many observations and is very sparse. The presence of sparse columns in $Z$ and $\widehat{\Xi}$ means that changing the corresponding fixed effect parameters affects the moment functions very little, resulting in near-zero columns in the (sample) matrix of moment function derivatives $\widehat{G}$. Stata can treat these as actually being columns of zeros.

This sparsity issue can be particularly common when estimating any version of the Prentice and Gloeckler (1978) model due to the need to include time fixed effects, and can be compounded if other fixed effects are included as controls. This is because each level of a fixed effect represents a dummy variable (i.e., a binary indicator variable), and we must directly include these dummy variables as regressors -- since the second stage is nonlinear rather than linear, you cannot simply instead “absorb” the fixed effects by applying a within-transformation to the whole dataset.

The option difficult_vce, which requests that alternative code be used to calculate the VCE, resolves this problem. With this option, the code will construct $\widehat{G}$ and $\widehat{\Omega}$ and then solve for $\widehat{V}$ with a tolerance threshold of 0 (i.e., floating-point numbers are treated as 0 only if they are exactly 0). Since $\widehat{G}$ is exactly invertible, this prevents Stata from declaring that $\widehat{G}$ is numerically singular.

An alternative, fast procedure for obtaining the second-stage VCE for generic control function models is described by Terza (2016).

Conclusion

This paper introduces a theory-complete IV methodology for the estimation of hazard models and discrete choice models) that have endogenous regressors, including when there is measurement error in the regressors. Although there are existing methods for dealing with mismeasured regressors in hazard models, these are not “true” IV methods: they do not allow for other types of endogenous regressors, and are therefore limited in their usefulness. On the other hand, there do exist proposals on IV methodologies to tackle the problem of endogeneity in hazard models, but none give an adequate treatment of the subject: these are either highly heuristic and not theoretically grounded or dangerously unrobust. Similarly, existing papers on the use of control functions to instrument discrete choice models are narrowly restricted to probit and Poisson regression, and a solution for generic discrete choice models is missing.

To fill this gap, I rigorously demonstrate how the IV relevance condition and instrument conditional independence assumption (the control function analog of the 2SLS exclusion restriction), together with minor conditions and the assumption of a sufficiently flexible functional form for the control function, allow causal effects to be recovered from hazard models and discrete choice models with generic endogenous regressors.

This paper contributes to a nascent literature, and theoretically justifies and undergirds already existing empirical work that uses the proposed estimator. It supplies a theory guide, thus serving as a crutch for existing and future empirical work, and this is compounded by the guidance, presented in this paper, on my Stata command, ivcloglog, which allows practitioners to easily implement the proposed estimator. Furthermore, this paper serves as a starting point for future theoretical work to further improve the robustness of IV methodologies in hazard models. All in all, this paper opens up promising new theoretical and empirical avenues for practitioners to investigate.