 Research article
 Open Access
 Published:
Time dependent hazard ratio estimation using instrumental variables without conditioning on an omitted covariate
BMC Medical Research Methodology volume 21, Article number: 56 (2021)
Abstract
Background
Estimation that employs instrumental variables (IV) can reduce or eliminate bias due to confounding. In observational studies, instruments result from natural experiments such as the effect of clinician preference or geographic distance on treatment selection. In randomized studies the randomization indicator is typically a valid instrument, especially if the study is blinded, e.g. no placebo effect. Estimation via instruments is a highly developed field for linear models but the use of instruments in timetoevent analysis is far from established. Various IVbased estimators of the hazard ratio (HR) from Cox’s regression models have been proposed.
Methods
We extend IV based estimation of Cox’s model beyond proportionality of hazards, and address estimation of a loglinear time dependent hazard ratio and a piecewise constant HR. We estimate the marginal timedependent hazard ratio unlike other approaches that estimate the hazard ratio conditional on the omitted covariates. We use estimating equations motivated by Martingale representations that resemble the partial likelihood score statistic. We conducted simulations that include the use of copulas to generate potential timestoevent that have a given marginal structural time dependent hazard ratio but are dependent on omitted covariates. We compare our approach to the partial likelihood estimator, and two other IV based approaches. We apply it to estimation of the time dependent hazard ratio for two vascular interventions.
Results
The method performs well in simulations of a stepwise timedependent hazard ratio, but illustrates some bias that increases as the hazard ratio moves away from unity (the value that typically underlies the null hypothesis). It compares well to other approaches when the hazard ratio is stepwise constant. It also performs well for estimation of a loglinear hazard ratio where no other instrumental variable approaches exist.
Conclusion
The estimating equations we propose for estimating a timedependent hazard ratio using an IV perform well in simulations. We encourage the use of our procedure for timedependent hazard ratio estimation when unmeasured confounding is a concern and a suitable instrumental variable exists.
Background
Confounding is a threat to all observational studies. The factors that influence the outcome of interest may also influence the selection of treatment. In randomized studies, intentiontotreat estimators are generally consistent estimators of the intentiontotreat effect but are not consistent estimators of the treatment effect. Instrumental variables (IV) influence treatment choice but otherwise have no effect on the outcome (exclusion restriction) and are independent of any other causes of the outcome (randomization assumption).
Estimation of treatment effects using instrumental variables for outcomes subject to right censoring has received the attention of many applied and several theoretical studies. Stukel et al. [1] proposed an adhoc estimator of the hazard ratio (HR) based on a linear model. MacKenzie et al. [2] proposed a hazard ratio estimator that assumes omitted covariates have an additive effect on the hazard. Tchetgen et al. [3] proposed an estimator of an additive hazards model based on two stage residual inclusion. They also argue that if the survival curve is close to unity (e.g. above 80%) for most of the followup that their additive hazards approach can be used as a good approximation of the multiplicative hazards approach. Li et al. [4] proposed a consistent estimator of a treatment satisfying an additive hazard model. Martinussen et al. [5] derived a consistent estimator for a structural Cox model. MartínezCamblor et al. [6] identified the role of frailties in estimation of the hazard ratio via the twostage residual inclusion algorithm if the treatment and omitted covariate jointly satisfy a Cox model. Wang et al. [7] derived an estimator of the marginal hazard ratio using a binary instrument.
The proportionality of hazards is a key assumption in estimation of the hazard ratio. Many test statistics have been proposed to test this assumption and estimate the hazard ratio as a function of time. The most simple approach to moving beyond proportionality of hazards is a piecewise constant function. The next most simple is a linear function of time. In this paper we extend the IV based estimator of the HR we proposed in 2014 MacKenzie et al. [2] beyond proportionality of hazards. We propose a method of estimating the HR as a timedependent function. We address two cases, (a) a piecewise constant HR and (b) a loglinear time dependent hazard ratio. We conduct Monte Carlo simulations that include the use of copulas to generate potential timestoevent. We demonstrate this method for timedependent hazard ratio estimation to compare the effect of two vascular interventions on survival.
Methods
Marginal structural Cox proportional hazards model
MacKenzie et al. [2] implemented a Cox proportional hazards model for the effect of treatment on the timetoevent that is marginal with respect to any omitted covariate. The term marginal is used to mean that if the treatment is applied independently of the omitted covariate, the distribution of the timetoevent conditioning only on the treatment satisfies a Cox PH model, which we now elaborate on.
Let {S_{x}} be the survival curve of the potential outcomes (timetoevent) if subjects receive level x of the exposure. In what follows we shall assume X is binary, e.g. 0 (no treatment) or 1 (treatment), although the methods described below are applicable to a continuously valued exposure. The structural Cox model is given by \(S_{x}(t)=S_{0}(t)^{\exp \{\beta _{X}\cdot x\}}\) and it is usually written in terms of the hazard function,
where λ_{0}(t) is the baseline hazard function. For any covariate, U, which also affects the timetoevent, let S_{x}(· ; u) be the survival curve for the potential timetoevent if a subject receives level x conditional on U=u. The marginal structural Cox proportional hazards model supposes that
where F_{U} stands for the distribution function of U on the studied population. The alternative to the marginal model is a conditional model that specifies a structural form for the conditional distribution given both the exposure X and omitted covariate U (e.g. multivariable Cox PH model). Due to noncollapsibility of Cox’s model [8], the marginal model that corresponds to this multivariable Cox model is not a Cox PH model.
Specifying that the marginal model is a Cox PH model has the following advantages. First, it makes no parametric assumption about how the omitted covariate U affects the distribution of the timetoevent: We allow for the possibility that the omitted covariate has a causal effect on the outcome without specifying a parametric model. Second, it does not require interpreting the treatment effect as conditional on a variable that is unknown (the omitted covariate). Third, this is the convention used in the reporting of randomized trials: It is unusual to analyze randomized trials using models that control for omitted covariates (such as Cox’s model with frailty), and the convention has been not to condition on measured characteristics either [9].
Time dependent hazard ratio
Proportionality of hazards is equivalent to a hazard ratio that does not vary over time. A constant hazard ratio can be a simple and good approximation, just as a linear model may be chosen when there is evidence of some nonlinearity. On the other hand there are times when the hazard ratio varies considerably and attention is focused on its time dependence. For instance, surgery may be associated with increased risks up front compared to endovascular procedures. There is a large literature on estimation of the hazard ratio as a function of time including smoothing splines [10], regression splines [11] or penalized regression splines [12] among other techniques.
Timetoevent, exposure and instrument
Suppose the observed data consists of independent and identically distributed observations, \(\left \{W_{i}, X_{i}, \Delta _{i}, T_{i}\right \}_{i=1}^{n}\) where W_{i} is the instrument (1≤i≤n) and X_{i} is the exposure level (either binary, ordinal or continuous). \(\Delta _{i}=1_{T^{0}_{i}<C_{i}}\) and \(T_{i}={min}(T^{0}_{i}, C_{i})\) are the event indicator and observed followup time respectively where \(T^{0}_{i}\) is the focus of interest, the timetoevent, and C_{i} is the censoring time. We assume that the potential censoring times {C_{i}(x)}_{x=0,1} are independent of the potential timetoevents \(\left \{T^{0}_{i}(x)\right \}_{\left \{x=0,1\right \}}\).
An instrument is defined by the following properties:

1.
It is associated with the exposure (W ⊥̸ ⊥X).

2.
It satisfies the following two restrictions:

(a)
There is no effect of the instrument on the outcome except through its effect on the exposure.

(b)
There are no confounders of the instrument and the outcome.

(a)
The latter two statements can be summarized as the potential outcomes are independent of the instrumental variable, {T^{0}(x)}_{{x=0,1}} ⊥ ⊥W.
To satisfy Assumption 1 the instrument may either be causal (i.e. it affects treatment choice) or noncausal (it is an effect of a variable that also affects exposure). Assumptions 2(a) and 2(b) are often combined by stating that conditional on X, the outcome and instrument are independent. The instrument, W, may be continuous, ordinal or binary.
For each subject there is one potential stochastic process {N_{i}(t;x),R_{i}(t;x)}_{t≥0,x} for each possible level of x where R_{i}(t;x) is the at risk indicator \(1_{T_{i}(x) \geq t}\), and N_{i}(t;x) is the counting process\(1_{T_{i} \leq t}\). Let R_{i}(t)=R_{i}(t;X_{i}) and N_{i}(t)=N_{i}(t;X_{i}).
In studies of instrumental variables another assumption that is typically introduced is homogeneity of the treatment effect and monotonicity of the effect of the instrument. In usual applications homogeneity of the treatment effect is the assumption that the effect is same for different levels of the omitted covariate, i.e. there is no interaction of the effect of the treatment and the omitted covariate. In the marginal model an estimate of the treatment effect will depend on the value of the omitted covariate. We do not specify the effect of the omitted covariate; we avoid this parameterization as the method we propose below relies on a first order linear approximation of the hazard in terms of the omitted covariate.
Another frequently addressed assumption in studies of instrumental variables is monotonicity of the effect of the instrument on the treatment. It assumes that in all subjects the instrument is positively associated with the treatment; for instance, the derivation of the complier average causal effect is dependent on it. In this study we make the monotonicity assumption, Pr [X≥xW=w_{1}]>Pr[X≥xW=w_{0}] if w_{1}>w_{0} for all x.
Estimating equations for a time dependent hazard ratio
Let λ_{x}(t) be the hazard function corresponding to survival in the population were all subjects in that population exposed to level x of the exposure, i.e., λ_{x}(t)=−dS_{x}(t)/S_{x}(t) which can also be written \(\mathbb E[dN_{i}(t;x)R_{i}(t;x)=1]\). We suppose that, at time t, the marginal (e.g. population) causal effect of having been exposed to level 1 relative to level 0 of the timevarying binary treatment X(t) is a change in the loghazard of β_{X}(t):
Further, suppose that the time dependent loghazard ratio is parameterized by β_{X}(t;θ). For instance, one simple parameterization is β_{X}(t;θ)=θ_{0}+θ_{1}·t.
If there is no selection bias, i.e. X is independent of any omitted covariate, U, the difference
which is known as a Martingale residual (Kalb eisch and Prentice [13]), has an expected value of zero and is independent of X\(\left (\Lambda _{0}(v)=\int _{0}^{v} \lambda _{0}(s)ds\right)\). This implies that
Moreover,
for any real function ϕ(·).
This equation suggests that an estimator of β_{X}(t;θ) is that \({\hat \beta _{X}(t;\theta)}\) for which
where τ is a suitably large point of time (e.g. max1≤i≤n{t_{i}}), K is the dimension of the parameter θ=(θ_{1},⋯,θ_{K}) and \(\left \{\phi _{k}(t)\right \}_{k=1}^{K}\) are suitably selected real functions. These functions should be chosen to minimize the variance of the resulting estimator which motivates the choices ϕ_{k}(t)=∂β_{X}(t;θ)/∂θ_{k}.
As Λ_{0}(t) is unknown we propose to solve for it using the estimating equations
for all t which is solved by the Breslow estimator (Lin [14])
Substitution of (4) into (2) yields the score equations for the partial likelihood (Cox [15]) for a timedependent hazard ratio (Abrahamowicz et al. [11]). That is, the estimating equations (2) are equivalent to the method of maximum partial likelihood estimation of a timedependent hazard ratio. If X is endogenous (confounded), then X_{i} is not necessarily independent of the Martingale residual
and therefore the estimating Eq. (2) are biased.
If W is an instrumental variable then it is approximately independent of the counting process conditional on X, or equivalently, the IV is independent of the Martingale residual. A heuristic justification starts by writing the causal hazard function conditional on an omitted covariate U that affects both the treatment and the timetoevent. We approximate the causal hazard function by a linear Taylor expansion, E [dN_{i}(t)R_{i}(t)=1]= exp{β_{X}(v;θ)·X_{i}}dΛ_{0}(v)+ϕ(t)[U−μ_{U}(t)] where μ_{U}(t) is the expected value of the omitted covariate among people at risk at time t. Then equating to zero the correlation of the instrument with the counting process leads to the following estimating equation θ,
If we substitute the Breslow estimator for Λ_{0}(v) into the latter equation it results in the estimating equation
If the timedependent function is constant, i.e. β_{X}(v;θ) equals a scalar θ this is the same estimating equation proposed in MacKenzie et al. [2].
IV estimation for piecewise proportional hazards
In applications, categorization of followup time is favored as an approach to assessing timedependence of the hazard ratio because of its simplicity. For instance, one can make statements such as, in the first month the hazard ratio was 1.5 but after one month was 0.5. The timeframe is categorized into intervals \(\left \{\tau _{i1},\tau _{i}\right \}_{i=1}^{K}\) for 0=τ_{0}<τ_{1}<…<τ_{K}=τ and within each window the hazard ratio is assumed to be constant, β_{X}(t;θ)=θ_{i} for τ_{i−1}≤t<τ_{i}. The partial likelihood estimator of a step function time dependent hazard ratio is equivalent to applying the partial likelihood estimator of the proportional hazards model within each time window. For instance, to estimate the hazard ratio in the interval [τ_{i−1},τ_{i}) exclude all subjects eliminated from risk before time τ_{i−1} and censor any events after time τ_{i} (1≤i≤K). For this piecewise constant hazard ratio the IV based estimator proposed in (6) is equivalent to applying the IV based method of MacKenzie et al. [2]. The R function supplied in that paper can be used accordingly to implement this approach for each time window (https://github.com/toddamackenzie/InstrumentalVariableHazardRatioEstimation).
IV estimation for linear time dependent log hazard ratio
Any parameterization of the hazard ratio could be implemented with the instrumental variable estimating equations we propose. We chose to illustrate the approach using a loglinear timedependent model of the hazard ratio, β_{X}(t;θ)=θ_{0}+θ_{1}·t. In this parameterization exp{θ_{0}} is the hazard ratio at inception and exp{θ_{1}} is the multiplicative change in the hazard ratio per unit time. The two parameters θ_{0} and θ_{1} can be estimated using the equations
and
respectively.
Monte Carlo simulations
We evaluated the behavior of the estimating equations we propose in (6) under two scenarios for the marginal timedependent hazard ratio; i) a three piece constant hazard ratio, and ii) a log linear timedependent hazard ratio. In the first scenario, the estimation method is equivalent to applying the approach of MacKenzie et al. [2] in early, middle and late followup settings, and that is how we implement and report the results. These simulations are novel in that the 2014 paper did not consider estimation of the hazard ratio later in followup, such as where the survival curve dips below 75%. Because a much greater proportion of the study will have either experienced the outcome event or been censored by that time, it is reasonable to expect that the performance of the method could differ over this timeperiod compared to the early followup period. For the loglinear function of time we consider a range of intercepts and a range of slopes. Specific details follow.
Comparators
In addition to evaluating the bias of our estimator for the stepwise constant and loglinear hazard ratio models, we evaluate the performance of the partial likelihood estimator. For the stepwise hazard model, because it is a proportional hazards model in each of the three periods, we also evaluate the performance of two other estimators designed for proportional hazards. First, we have evaluated the approach of Wang et al.[7] which estimates the hazard ratio for a marginal structural model like ours. Because Wang et al’s approach is designed for a binary instrument, we dichotomize the continuous instrument of our simulations at the median. Second, we have evaluated the performance of the estimator of MartínezCamblor et al. [6] which estimates the hazard ratios of a Cox model for the multivariable effect of the treatment and the omitted covariate; that is, it conditions on the omitted covariate.
Simulation methods
We have conducted extensive simulations to evaluate the bias of the estimator we propose. Our simulation uses a copula to generate timestoevent that depend on an omitted covariate and treatment in such a way that the marginal distribution of the potential outcomes, treated and untreated, satisfy a model with the specified timedependent hazard ratio. The settings and parameter values for the simulation are listed in Tables 1 and 2.
The steps of the simulation are:

1.
Randomly generate the bivariate random variable, (Y_{0},Y_{1}), from a bivariate copula., i.e. Y_{0} and Y_{1} are uniformly distributed but not independent.

2.
Let the omitted covariate be the standard normal deviate obtained as U=Φ^{−1}(Y_{0}).

3.
Randomly generate a continuous instrument, W, from the standard normal distribution.

4.
Randomly generate a binary exposure indicator X using a logistic link with intercept of zero (yields equal proportions of treated and untreated) and log odds ratios, α_{W} and α_{U}, for the association of W and U respectively.

5.
Obtain potential timetoevents T(0)=− log(Y_{0}) (unit exponential distribution) and \(T(1)=F_{1}^{1}(Y_{1})\) where F_{1} is the cumulative distribution function whose hazard equals the specified timedependent hazard ratio

6.
Right censor the timestoevent using a uniform distribution whose scale is set to achieve the specified censoring rate.
In our Monte Carlo simulation the confounding (endogeneity) is created by (i) setting the omitted covariate as a monotonic function of the potential outcome under no treatment, which is correlated with the potential outcome under treatment (as controlled by the bivariate copula using either Pearson or Kendal correlations of 0.50 or 0.99) and (ii) randomly generating the treatment indicator based on the omitted covariate (as well as the instrument). Therefore, the level of confounding by the omitted covariate is controlled by a single parameter, the odds ratio, exp (α_{UX}) relating U (the untreated timetoevent transformed to have a standard normal distribution) to X. For instance, if α_{UX}=0 there is no confounding. Before the generation of each dataset it was randomly drawn from a Uniform distribution on the interval log(1/10), log(10). Values above log(5) were considered positive confounding, and values below − log(5) were considered negative confounding.
The stepwise constant timedependent hazard ratio took values as specifed in Table 2. For each simulated dataset we randomly sampled with replacement from these grids for the early, mid and late followup hazard ratios. The early corresponded to approximately that time period up to which survival curve was above 90%, the middle 90% down to 75%, and late survival of 75% or less. The loglinear timedependent hazard ratio was generated based on the randomly drawn intercept and slope as specified in Table 2.
Complete R code is available at: https://github.com/toddamackenzie/InstrumentalVariableHazardRatioEstimation, for users interested in conducting the same simulations.
Simulation results for early, mid and late estimators of the HR
Figure 1 shows the results of the simulations when the instrument is strong, which we defined as an odds ratio between the treatment and the instrument that exceeds 5. Each red point represents the median of over 1000 maximum partial likelihood estimates of the stepwise hazard ratio. The partial likelihood estimators are clearly subject to bias from confounding, which decreases as followup increases. Each black point represents the median of over 1000 estimates using the instrument and estimating equations we propose. The green and blue points represent the median of the distribution of the estimators of Wang et al. [7] and MartínezCamblor et al. [6]. The instrumental variable based estimator we propose is unbiased as an estimator of the marginal hazard ratio in each of the three periods, with a slight tendency toward bias by confounding for large and small hazard ratios; that is, the magnitude of the bias increases as the HR moves away from null and in the direction opposite of the confounding. The estimators of Wang et al. [7] and MartínezCamblor et al. [6] demonstrate similarity and slightly better accuracy.
The bias was affected by the strength of the instrument; for each doubling of the odds ratio between the treatment and the instrument, the bias reduced by approximately 5%. Figure 2 shows the results restricted to instruments for which the odds ratio between treatment and instrument was less than 5. The results are similar to those reported for the scenario of a strong instrument. The bias of our estimator did not vary with respect to choice of copula, or the correlation between the treated and untreated potential outcomes (results not reported).
Simulation results for logLinear time dependent HR
Figure 3 shows the results of the simulations for estimation of the loglinear hazard ratio. The top row demonstrates that the intercept term of the loglinear hazard ratio is estimated with very little bias using the instrumental variable estimating equation we propose (blue points represent median of over 1000 estimates from that many simulated datasets), unlike the maximum partial likelihood estimators (points in red). Results from the estimation of slope parameter of the loglinear hazard ratio are found in the bottom row, which again show that the IV based estimating equation we propose yields unbiased estimators (blue points), unlike the partial likelihood estimator (in red) with the exception that the IV based estimator of the slope is biased when the true slope is zero. The estimators of Wang et al. [7] and MartínezCamblor et al.[6] are not designed for loglinear estimators or any timedependent hazard ratio beyond piecewise constant.
The simulation findings did not change with respect to the choice of the copula, the correlation between the treated and untreated potential outcomes, or the strength of the instrument.
Realworld application
We illustrate our method to address a comparative effectiveness question in patients with carotid stenosis. We use data from the nationwide Vascular Quality Initiative, VQI, (http://www.vascularqualityinitiative.org) to compare the composite outcome of death or stroke following intervention between those undergoing endarterectomy (CEA) and those undergoing carotid stenting (CAS). The data consists of 73,312 patients who received CEA and 12705 who received CA during the years 20032016. The number of events was 8,005 of which 730 occurred in the first 30 days, 2498 in months 1 through 12 and 4,777 after the first year. This example has been previously considered by Columbo et al. [16]. Figure 4 (top) shows the KaplanMeier estimates for the survival curves on both the CEA and the CAS groups.
The population of patients who undergo CEA may not be the same as those who undergo CAS, and therefore any estimator of comparative efficacy may be biased by confounding. Therefore we utilize an instrumental variable approach. The instrument we employ is the center level relative frequency of CEA versus CAS procedures over the 12 months prior to the current patient. A value of zero indicates that the patient received the surgery in a facility which during the prior 12 months, has not performed CEA. The validity of this instrument is argued in MartínezCamblor et al. [17]. Figure 4, bottom, depicts the instrument’s distribution (histograms at left and boxplot at right) in both groups.
Figure 5 shows estimates of the hazard ratio as a function of time. The left panel contains the partial likelihood estimators of the stepwise constant hazard ratio and the loglinear hazard ratio. After consulting with physicians and a visual inspection of the KaplanMeier survival curve, we chose the cutoffs of 1 month (30 days) and 6 months. In the first month the partial likelihood estimate of the hazard ratio comparing CEA to CAS is 0.46 (95%CI: 0.39 to 0.54), while it is 0.55 (95% CI: 0.48 to 0.62) in months 2 through 6, and considerably closer to unity thereafter at 0.75 (95%CI: 0.71 to 0.82). The corresponding estimates based on our instrumental variable approach are similar in the first interval, 0.44 (95%CI: 0.31 to 0.61), although the confidence interval is wider as expected for IV based estimators. The IV based estimator of the hazard ratio in the middle interval is closer to unity, 0.66 (95%CI: 0.50 to 0.88). The estimator in the final interval is 0.77 (95%CI: 0.65 to 0.91).
The model for the hazard ratio as a loglinear function of time indicates the same pattern of the hazard ratio moving toward unity as time progresses. Both the partial likelihood estimator and the IV based estimator indicate the hazard ratio begins at approximately 0.5 but the IV based estimator reaches unity earlier. The observed results at both early and late followup were quite close in the standard and proposed methodologies suggesting that, in average, the potential covariates affecting the survival have a minor impact on the effect of the treatment. Between the first and six months, the difference is around 20% indicating that in this period some omitted covariates spuriously enhance the observed effect of the treatment.
Discussion
We have proposed a method for using instruments in the estimation of timedependent hazard ratios. The framework is general enough to accommodate any parameterization for a timedependent hazard ratio. Like maximum partial likelihood estimation of Cox’s proportional hazards model it does not require a parameterization of the baseline hazard. We have illustrated our approach using two forms for a timedependent hazard ratio, (i) a piecewise constant hazard ratio and (ii) a log linear function of time.
Our approach focuses on timedependent hazard ratios that are marginal with respect to the omitted covariate. That is, just like estimators of treatment effects in randomized studies, the model we employ does not explicitly condition on the omitted covariate although it is motivated by a linear approximation of a hybrid hazard function in which the dependence on the omitted covariate is additive and the dependence on treatment and the observed covariates is multiplicative. Our approach is analogous to estimators of population averaged parameters such as generalized estimating equations. We encourage analysts to estimate both marginal models like ours and models that condition on the omitted covariate. A disadvantage of the conditional model is that it conditions on intangible characteristics, that is unmeasured characteristics, that somehow affected the treatment selection.
We conducted extensive simulations to evaluate our estimator of the piecewise hazard ratio and the loglinear hazard ratio. For the former, which reduces to a proportional hazards model for any of the time windows (in which the hazard ratio is constant) we compared our approach to the partial likelihood estimator, the approach of MartínezCamblor et al. [6] and the approach of Wang et al. [7]. For the loglinear hazard ratio model, we compared our approach to the partial likelihood estimator. The simulations indicated our estimator of a piecewise hazard ratio has some bias toward the null which is comparable to the bias of MartínezCamblor et al. [6] and Wang et al. [7]. Our simulations indicate our estimator of the loglinear hazard ratio has little bias, especially in comparison to maximum partial likelihood estimation when confounding by omitted covariates is present.
While we do not address it, this method generalizes readily to multivariable models without much additional complexity. That is, measured covariates can be incorporated into the model.
The method we propose has some notable limitations. We do not believe it is a consistent estimator because we derived our estimating equations by substituting for the baseline hazard the Breslow estimator without justification of why this is the ideal estimator using an instrumental variable. In addition, it is based on an approximation of the hazard as an additive function of the omitted covariate.
We believe our estimator can be improved by utilizing inverse probablity weighting in our estimating equations. In particular, like the approach of Wang et al. [7] and Huling et al. [18] one should weight the subjects at risk at any time t using the instrumental variable to make those subjects at risk generalizable to the entire population, i.e. those at risk at time t=0, at which time the instrument is orthogonal to the omitted covariates.
Conclusion
The estimating equations we propose for estimating a timedependent hazard ratio using an IV perform well in simulations. We encourage the use of our procedure for timedependent hazard ratio estimation when unmeasured confounding is a concern and a suitable instrumental variable exists.
Availability of data and materials
The code for this method and the simulations is available at: https://github.com/toddamackenzie/InstrumentalVariableHazardRatioEstimation. The data used in the example was accessed from Centers for Medicare and Medicad Servies under contract and is not available for sharing.
Abbreviations
 IV:

Instrumental Variable
 HR:

Hazard Ratio
References
 1
Stukel T, Fisher E, Wennberg D, Alter D, Gottlieb D, Vermeulen M. Analysis of observational studies in the presence of treatment selection bias: effects of invasive cardiac management on AMI survival using propensity score and instrumental variable methods. JAMA. 2007; 297(3):278–85.
 2
MacKenzie T, Tosteson T, Morden N, Stukel T, O’Malley A. Using instrumental variables to estimate a Cox’s proportional hazards regression subject to additive confounding. Health Serv Outcome Res Methodol. 2014; 14(12):54–68.
 3
Tchetgen E, Walter S, Vansteelandt S, Martinussen T, Glymour M. Instrumental variable estimation in a survival context. Epidemiol Camb Mass. 2015; 26(3):402.
 4
Li J, Fine J, Brookhart A. Instrumental variable additive hazards models. Biom. 2015; 71(1):122–30.
 5
Martinussen T, Vansteelandt S, Tchetgen Tchetgen E, Zucker D. Instrumental variables estimation of exposure effects on a timetoevent endpoint using structural cumulative survival models. Biom. 2017; 73(4):1140–9.
 6
MartínezCamblor P, MacKenzie T, Staiger D, Goodney P, James O’Malley A. An instrumental variable procedure for estimating Cox models with nonproportional hazards in the presence of unmeasured confounding. J R Stat Soc: Ser C: Appl Stat. 2019.
 7
Wang L, Tchetgen E, Martinussen T, Vansteelandt S. Learning Causal Hazard Ratio with Endogeneity. 2018. http://arxiv.org/abs/1807.053131807.05313.
 8
Martinussen T, Vansteelandt S. On collapsibility and confounding bias in Cox and Aalen regression models. Lifetime Data Anal. 2013; 19(3):279–96.
 9
Kahan B, Jairath V, Doré C, Morris T. The risks and rewards of covariate adjustment in randomized trials: an assessment of 12 outcomes from 8 studies. Trials. 2014; 15(1):139.
 10
Hastie T, Tibshirani R. Varyingcoefficient models. J R Stat Soc Ser B Methodol. 1993; 55(4):757–79.
 11
Abrahamowicz M, Mackenzie T, Esdaile J. Timedependent hazard ratio: modeling and hypothesis testing with application in lupus nephritis. J Am Stat Assoc. 1996; 91(436):1432–9.
 12
Gray R. Flexible methods for analyzing survival data using splines, with applications to breast cancer prognosis. J Am Stat Assoc. 1992; 87(420):942–51.
 13
Kalbfleisch J, Prentice R. The Statistical Analysis of Failure Time Data, 2nd: John Wiley & Sons; 2002. ISBN 9780471363576.
 14
Lin D. On the Breslow estimator. Lifetime Data Anal. 2007; 13(4):471–80.
 15
Cox D. Regression models and lifetables. J R Stat Soc Ser B Methodol. 1972; 34(2):187–220.
 16
Columbo J, MartínezCamblor P, MacKenzie T, Staiger D, Kang R, Goodney P, O’Malley A. Comparing longterm mortality after carotid endarterectomy vs carotid stenting using a novel instrumental variable method for risk adjustment in observational timetoevent data. JAMA Netw open. 2018; 1(5):181676.
 17
MartínezCamblor P, Mackenzie T, Staiger D, Goodney P, O’Malley A. Adjusting for bias introduced by instrumental variable estimation in the Cox proportional hazards model. Biostat. 2019; 20(1):80–96.
 18
Huling J, Yu M, O’Malley A. Instrumental variable based estimation under the semiparametric accelerated failure time model. Biom. 2019; 75(2):516–27.
Acknowledgements
The authors would like to acknowledge Phil Goodney and Jesse Columbo for facilitating access to the data.
Funding
This work was supported in part by a grant ME150328261 from the USA federally funded agency, PatientCentered Outcomes Research Institute, to develop statistical methodology for using instrumental variables in survival analysis. PCORI played no roles in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Author information
Affiliations
Contributions
TM concieved of the method, developed the math, performed the simulations and led the writing of the manuscript. PMC contributed to the theory, performed the statistical analysis for the example and contributed text and editing of the manuscript. JO contributed to the concept, simulations, example and writing of the manuscript. All three authors have approved the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
All data are collected under the auspices of an Agency for Healthcare Research and Quality designated Patient Safety Organization and were deidentified. Our study was approved by the Center for the Protection of Human Subjects at Dartmouth. All patient personal health information was protected, records, and outcomes were deidentified, and no testing or procedures were required for this study. Thus, the need for specific consent was waived.
Consent for publication
Not applicable.
Competing interests
There are no conflicts of interest to report.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visithttp://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
MacKenzie, T.A., MartinezCamblor, P. & O’Malley, A.J. Time dependent hazard ratio estimation using instrumental variables without conditioning on an omitted covariate. BMC Med Res Methodol 21, 56 (2021). https://doi.org/10.1186/s12874021012456
Received:
Accepted:
Published:
Keywords
 Causal inference
 Censoring
 Semiparametric model
 Marginal model