# A Comparison of Methods to Estimate the Hazard Ratio Under Conditions of Time-varying Confounding and Nonpositivity

In occupational epidemiologic studies, the healthy worker survivor effect refers to a process that leads to bias in the estimates of an association between cumulative exposure and a health outcome. In these settings, work status acts both as an intermediate and confounding variable and may violate the positivity assumption (the presence of exposed and unexposed observations in all strata of the confounder). Using Monte Carlo simulation, we assessed the degree to which crude, work-status adjusted, and weighted (marginal structural) Cox proportional hazards models are biased in the presence of time-varying confounding and nonpositivity. We simulated the data representing time-varying occupational exposure, work status, and mortality. Bias, coverage, and root mean squared error (MSE) were calculated relative to the true marginal exposure effect in a range of scenarios. For a base-case scenario, using crude, adjusted, and weighted Cox models, respectively, the hazard ratio was biased downward 19%, 9%, and 6%; 95% confidence interval coverage was 48%, 85%, and 91%; and root MSE was 0.20, 0.13, and 0.11. Although marginal structural models were less biased in most scenarios studied, neither standard nor marginal structural Cox proportional hazards models fully resolve the bias encountered under conditions of time-varying confounding and nonpositivity.

From the ^{a}Department of Epidemiology, Gillings School of Global Public Health, UNC-Chapel Hill, NC and ^{b}Department of Obstetrics and Gynecology and Duke Global Health Institute, Duke University

Submitted 21 October 2010; accepted 29 March 2011; posted 11 July 2011.

Supported in part by NIH grant R01CA117841 (D.B.R.), Fonds de Recherche en SantÃ© du QuÃ©bec through a Doctoral Research Award (A.I.N.), and NIH/NICHD grant K99-HD-06-3961 (D.J.W.).

Correspondence: Stephen R. Cole, Epidemiology, CB7415, University of North Carolina, Chapel Hill, NC 27599. E-mail: cole@unc.edu.

The healthy worker effect has long been recognized as a potential source of bias when estimating the association between an occupational exposure and health outcome such as mortality.^{1} Two aspects of the healthy worker effect can be distinguished from each other: the initial selection of healthy people into the work force, and the tendency for workers at increased risk of mortality to preferentially leave employment. The latter is known as the healthy worker survivor effect.^{2,3} Analytic methods to overcome the bias induced by the healthy-worker survivor effect began to appear in the early 1970s.^{2} In 1986, Robins^{4} identified this bias as one due to confounding of the association between cumulative exposure and mortality by time-varying work status, which is affected by prior exposure. Standard Cox proportional hazards regression models produce biased exposure-disease associations whether or not one adjusts for time-varying confounders affected by prior exposure.^{5â€“10} Marginal structural Cox proportional hazards regression models offer a viable alternative under such conditions.^{11}

In the context of occupational epidemiology, the healthy-worker survivor effect presents complications beyond what is encountered with standard time-varying confounding affected by prior exposure. Exposed persons who leave the workplace often have no chance of incurring work-based exposure at subsequent time points. If work status is a confounder, this situation results in a violation of the positivity assumption. The positivity assumption requires exposed and unexposed observations in all strata of the confounders at all time points.^{12â€“14} Violations of the positivity assumption can arise in diverse research settings.^{15,16} To make inferences that are not based on model interpolation or extrapolation, positivity is required.^{14}

Formally, the positivity assumption, which is also known as the experimental-treatment assignment assumption,^{17} is met when Pr(*X = x* | **L = l**) > 0 for all l, where Pr(**L = l**) â‰ 0, *X* is the exposure variable, and L is a vector of confounders.^{12} Violations of the positivity assumption (nonpositivity) are of 2 kinds: systematic and random. Systematic nonpositivity occurs when individuals cannot receive at least one level of the exposure within one or more of the confounder strata. Random nonpositivity occurs when no persons happen to be observed within one or more of the confounder strata.^{14} Although both types of nonpositivity can threaten the validity of inferences made with respect to exposure-outcome associations, as a structural feature of the scenario under study, systematic nonpositivity is of greater concern.

Marginal structural models are known to yield unbiased results under conditions of time-varying confounding when the positivity assumption is met.^{11} However, their performance under conditions of time-varying confounding with systematic nonpositivity (as occurs in the healthy-worker survivor effect) is unknown. Here, we use Monte Carlo simulation to compare the performance of hazard ratios derived from standard and marginal structural Cox models under conditions of time-varying confounding and nonpositivity in data that mimic sample sizes typically encountered in occupational epidemiology.

## METHODS

### Simulated Data

The Figure presents a causal directed acyclic graph^{18} representing the healthy worker survivor effect^{19} on which our simulated data are based. Baseline exposure status *X*(0), and an unmeasured confounder *U*, were each generated as independent Bernoulli random variables with a probability of 0.5. To mimic a study of an occupational cohort in which people enter follow-up at start of employment, work status at baseline *W*(0) was always set to 1; as a constant, this variable may be ignored without consequence. Follow-up work status was denoted by *W*(1) and was defined as a Bernoulli random variable with a probability of being at work of Pr[*W*(1) = 1] = expit[Î¸ + Î³*U* âˆ’ Î²*X*(0)], where expit(Â·) = exp(Â·) / [1+ exp(Â·)], where Î¸ represents the intercept (chosen such that the marginal probability of *W*(1) is approximately 0.5), and where Î³ and Î² represent the association between the unmeasured confounder *U* and *W*(1) and between baseline exposure *X*(0) and *W*(1), respectively. To induce nonpositivity, follow-up exposure status *X*(1) was defined as zero if *W*(1) = 0, and a Bernoulli random variable with a probability of 0.5 if *W*(1) = 1. Survival time *T* followed an exponential distribution with a survival function exp(âˆ’Î»*t*) and a rate parameter of Î» = exp[âˆ’Î±Ï‡ + Î”*U*], where Î± is the log hazard ratio of the association between the cumulative exposure Ï‡ = *X*(0) + *X*(1) and *T*, and Î”, the log-hazard ratio of the association between *U* and *T*. Finally, the data were administratively censored to induce a 0.25 event proportion.

### Statistical Analysis

Performance was assessed using 3 regression models. Models 1 and 2 were a standard unadjusted Cox model, and a standard Cox model adjusted for work status, respectively. Model 3 was a marginal structural Cox model weighted by the inverse probability of exposure history, which in our scenario (Figure) amounts to weighting for the inverse probability that *X*(1) = *x* given work status *W*(1).^{13} The hazard ratios derived from each of these models were compared with the true hazard ratio. An estimate of the true hazard ratio was obtained using a marginal structural Cox model weighted by the inverse probability that exposure status *X*(1) = *x* given the unmeasured confounder *U*, generated with a sample size of 20 million. According to the causal structure specified in the Figure, weighting for *U*, rather than *W*(1) as in model 3, allows us to ascertain the unbiased (or true) marginal association between cumulative exposure Ï‡ and outcome *T*. Of course, in real research settings, one cannot weight for *U* because it is unobserved. Parameters for model 3 and the true hazard ratio were estimated using stabilized inverse probability of exposure weights, calculated as

where *Z* = *W*(1) when weighting for *W*(1), and *Z* = *U* when weighting for *U.* SAS version 9.2 was used for all analyses (SAS Institute, Inc, Cary, NC).

### Simulation Summary

We investigated the behavior of the 3 Cox models (ie, crude, adjusted, and weighted) under 13 different specifications of Î±, Î², Î”, and Î³ representing the odds ratios or hazard ratios of 1, 2, 5, and 9, respectively. We defined the base-case scenario as: (i) a hazard ratio of 9 for Î”, the association between *U* and *T*; (ii) an odds ratio of 9 for Î³, the association between *U* and *W*(1); (iii) an odds ratio of 9 for Î² the association between *X*(0) and *W*(1); and (iv) a hazard ratio of 2 for Î±, the association between cumulative exposure Ï‡ and *T*. Using Monte Carlo simulation, we assessed the performance of hazard ratios derived from the 3 modeling strategies based on 20,000 trials, each with a sample size of 1000 to mimic modest-sized occupational cancer cohort studies.

To compare the 3 models across different specifications of Î±, Î², Î”, and Î³, we calculated 3 summaries: (i) bias, computed as the mean of the estimator across 20,000 trials minus the true value; (ii) 95% confidence interval (CI) coverage, computed as the proportion of times the 95% CI from the 3 models contained the true value; and (iii) root mean squared error (root MSE), computed as the square root of the sum of the squared bias and the simulated variance for the estimator. All summaries were calculated for Î±, the association between cumulative exposure Ï‡ and the outcome *T*. Confidence intervals for the standard Cox models (crude and work-status-adjusted) were estimated using model-based standard errors. Confidence intervals for the marginal structural Cox model were estimated using robust standard errors. The simulated variance was computed as the square of the simulated standard deviation of the estimator across all 20,000 trials.

## RESULTS

Each regression method resulted in biased estimates under some conditions. In general, the weighted association was least biased, the crude association most biased, and the adjusted association intermediate. Improved performance of the marginal structural Cox model relative to standard Cox regression may be expected, as the weighted association accounts for time-varying confounding affected by prior exposure, while neither the crude nor adjusted associations do so. The adjusted Cox model tended to be less biased than the crude Cox model, suggesting that, as expected,^{20} the collider stratification bias induced in the adjusted model is smaller than the confounding in the crude model.

There were exceptions to this general pattern, however. First, in Table 1, rows 8 and 11, there was no time-varying confounding because either *e* ^{Î³} = 1 or eÎ” = 1, respectively. In these 2 scenarios, all associations were unbiased. Second, in Table 1, rows 5 and 6, the time-varying confounder was not affected by prior exposure (ie, *e* ^{Î²} = 1), or was only modestly affected by prior exposure (ie, *e* ^{Î²} = 2). In these 2 scenarios, the adjusted association was less biased than the weighted association.

The same pattern can be observed with the root mean squared error. The weighted associations were more accurate (ie, smaller root MSE) than the adjusted, which in turn were more accurate than the crude (Table 2). There were again 2 exceptions. First, when either *e* ^{Î³} = 1 or *e* ^{Î”} = 1, root MSE for the weighted association was slightly higher than the crude or adjusted associations. Because there was no bias in these settings, this higher root MSE reflects the effect of lower precision in the weighted model. This lower precision is the expected result of unnecessary protection against time-varying confounding, which is not present in these 2 settings. Second, when either *e* ^{Î²} = 1 or *e* ^{Î²} = 2, the adjusted association was most accurate. These findings coincide with the patterns of bias observed in Table 1.

Finally, as demonstrated in Table 3, the 95% CI coverage for the weighted and adjusted associations proved on average to be roughly equivalent. For instance, across the 13 scenarios explored, the mean coverage for the crude, adjusted, and weighted associations were 64.8%, 89.2%, and 88.8%, respectively. However, the 95% CI coverage for the weighted association was less than 90% only for the 3 scenarios where the association between baseline exposure and subsequent work status was null to moderate. These are exactly the scenarios where the time-varying confounder is affected by prior exposure only moderately to not at all.

## DISCUSSION

Our findings indicate that neither standard nor marginal structural Cox models adequately deal with the bias induced by time-varying confounding with nonpositivity in the context of the healthy-worker survivor effect. As the results in Table 1 suggest, in the absence of time-varying confounding all methods explored recover the correct association, irrespective of nonpositivity. Furthermore, when time-varying confounding is not affected by prior exposure (ie, when the association between baseline exposure and subsequent work status was null or small), standard adjustment appears to provide a reasonable estimate of the true association. Consequently, when the healthy-worker survivor effect is thought to be operating, researchers could assess the association between exposure and subsequent work status. Stronger associations would suggest use of nonstandard methods, whereas null or modest associations might warrant the use of standard Cox models adjusted for work status. Finally, although the marginal structural Cox model proved more robust to increases in the strength of the association between the unmeasured confounder and either work status or mortality, the degree of this improvement over the standard work status adjusted model was modest.

Previous simulation research has considered the healthy-worker survivor effect under settings of no true exposure effect, and evaluated the ability to control for it using standard methods.^{21} As in that work, our simulation scenarios illustrate negative bias in estimates under a null exposure effect due to the healthy-worker survivor effect (Table 1, row 2). Further, we show that although weighted methods perform better under a null exposure effect, both standard and weighted methods give biased results. Previous simulation studies have not considered the performance of methods in settings with an actual exposure effect. With the few exceptions mentioned previously, marginal structural models consistently perform better in these settings.

The healthy worker survivor effect has been recognized by epidemiologists for more than 40 years,^{22} yet most of the proposed solutions inadequately address the biases that may arise under such conditions.^{2} For example, although restricting the analysis to survivors of a given period is likely to decrease confounding bias in certain settings,^{23} it may induce a selection bias such that the net bias is increased after restriction.^{19} Furthermore, proposed solutions such as exposure lagging^{24} and work status adjustment^{22} necessitate the unrealistic assumption that work status is not a risk factor for mortality.^{2} However, as a time-varying confounder affected by prior exposure, adjusting for work status may induce collider-stratification bias.^{19,20,25} Moreover, a key assumption underlying exposure lagging is that time off work is equivalent to time on work at zero exposure.^{26} Yet, if workers leave their jobs for health-related reasons (whether induced by exposure or some other cause) that also predispose them to higher rates of the event under study, then the risk of the event in persons no longer at work would be higher than the risk of the event in those at work but unexposed. Consequently, time off work is not likely to be equivalent to time on work under no exposure.^{26} Despite these limitations, however, exposure lagging may successfully be used to specify the appropriate exposure given an assumed empirical induction period.

Though marginal structural models are able to control the bias induced by time-varying confounders that are affected by prior exposure,^{11} as we described here, they do not account for the bias induced by nonpositivity. Because persons who have left the workplace cannot receive work-based exposures at subsequent time-points, the model used to estimate the inverse probability weights encounters a zero cell in the stratum where *W*(1) = 0 and *X*(1) = 1. Thus, the inverse of the probability of being exposed at time *t* = 1 when one has left the workplace is undefined. As a result, using the inverse probability weights fails to fully eliminate the bias induced by nonpositivity.

We have compared the performance of crude, adjusted, and marginal structural Cox models under conditions of systematic nonpositivity, using the healthy-worker survivor effect as an example. We would expect our results to generalize to any research context with a causal structure similar to the one outlined in the Figure. Whether these results generalize to settings where nonpositivity is random,^{14} remains uncertain. However, recent evidence suggests that, for marginal structural models, increased random nonpositivity results in less efficient estimators.^{27} A lack of positivity may also arise in other research settings, such as when assessing the independent effects of neighborhood deprivation on preterm birth^{15} or the verbal ability of children.^{16} It is uncertain how marginal structural models would deal with nonpositivity bias encountered in such settings.

The present findings should be interpreted with limitations in mind. First, our scenarios were characterized by data with only 2 exposure time points. The bias may be amplified or attenuated under conditions where more than 2 time points are present. Second, in our simulations, when bias is present, we cannot partition the bias into components due to time-varying confounding and due to nonpositivity. However, the results presented in Table 1, rows 5â€“7, do support the speculation that the inverse relationship between the bias for the weighted association and the magnitude of the effect of baseline exposure on subsequent work status is due to the relationship between nonpositivity and time-varying confounding. Specifically, when the effect of baseline exposure on subsequent work status is small (or null), the time-varying confounder is affected only weakly (or not at all) by prior exposure, and the constant nonpositivity bias plays a proportionately larger role in biasing the estimates. Finally, the choice of parameterizations used to specify the magnitude of effects was only a small portion of all possible scenarios. As in all Monte Carlo research, results apply only to the scenarios studied.^{28} Future research might assess these methods using a broader set of scenarios based on established cohorts, and using methodological approaches better suited to handle nonpositivity, such as history-restricted marginal structural models,^{29} g-estimation of a structural nested accelerated failure-time model,^{6} or the parametric g-formula.^{30,31}

Despite these limitations, our approach does offer strengths. To our knowledge, this is the first assessment of the performance of marginal structural models in the context of the healthy-worker survivor effect. Our causal structure and simulated data are characterized by only 2 biases, time-varying confounding and nonpositivity, allowing for a more clear understanding of the role that these biases play in occupational epidemiology. Finally, we compared the performance of a commonly used resolution to the healthy-worker survivor effect (adjustment for work status) to a less common, but perhaps more theoretically justified, approach (marginal structural models).

In conclusion, neither standard nor marginal structural Cox proportional hazards models can fully resolve the bias encountered under conditions of time-varying confounding with nonpositivity, such as in the healthy-worker survivor effect.

## ACKNOWLEDGMENT

We thank Chanelle Howe for helpful comments on earlier versions of the manuscript.

## REFERENCES

*Epidemiology*. 1994;5:189â€“196.

*Encyclopedia of Environmetrics*. New York: John-Wiley & Sons; 2002.

*Math Model*. 1986;7:1393â€“1512.

*The Encyclopedia of Biostatistics*. Chichester, United Kingdom: John Wiley and Sons; 1998.

*Biometrika*. 1992;79:321â€“334.

*Stat Med*. 1993;12:1605â€“1628.

*Pneumocystis carinii*pneumonia on the survival of AIDS patients [erratum in

*Epidemiology.*1993;3:189].

*Epidemiology*. 1992;3:319â€“336.

*AmJ Epidemiol*. 1998;148:390â€“401.

*Epidemiology*. 2000;11:550â€“560.

*J Epidemiol Commun Health*. 2006;60:578â€“586.

*Am J Epidemiol*. 2008;168:656â€“664.

*Am J Epidemiol*. 2010;171:674â€“677.

*Am J Epidemiol*. 2010;171:664â€“673.

*Proc Natl Acad Sci USA*. 2008;105:845â€“852.

*Am J Epidemiol*. 2005;162:382â€“388.

*Biometrika*. 1995;82:669â€“688.

*Epidemiology*. 2004;15:615â€“125.

*Epidemiology*. 2003;14:300â€“306.

*Am J Epidemiol*. 1996;143:202â€“210.

*Am J Epidemiol*. 1982;116:177â€“188.

*Br J Prev Soc Med*. 1976;30:225â€“230.

*Radiat Res*. 1979;79:122â€“148.

*Int J Epidemiol*. 2010;39:417â€“420.

*Occup Environ Med*. 1996;53:455â€“462.

*Int J Biostat*. 2010;6.

*Epidemiology*. 1997;8:453â€“456.

*Am J Epidemiol*. 2010;171:1233â€“1243.

*J Chronic Dis*. 1987;40(suppl 2):139Sâ€“161S.

*Int J Epidemiol*. 2009;38:1599â€“1611.