Evaluating Bias From Birth-cohort Effects in the Age-based Cox Proportional Hazards Model

Abstract

The nested case-control design is frequently used to evaluate exposures and health outcomes within the confines of a cohort study. When incidence-density sampling is used to identify controls, the resulting data can be analyzed using conditional logistic regression (equivalent to stratified Cox proportional hazards regression). In these studies, exposure lagging is often used to account for disease latency. In light of recent criticism of incidence-density sampling, we used simulated occupational cohorts to evaluate age-based incidence-density sampling for lagged exposures in the presence of birth-cohort effects and associations among time-related variables. Effect estimates were unbiased when adjusted for birth cohort; however, unadjusted effect estimates were biased, particularly when age at hire and year of hire were correlated. When the analysis included an adjustment for birth cohort, the inclusion of lagged-out cases and controls (assigned a lagged exposure of zero) did not introduce bias.

Since its use was first illustrated in 1977,¹ the nested case-control study design has been frequently used to evaluate exposures and health outcomes within occupational cohort studies (eg, lung cancer mortality among silica exposed workers²) and population-based cohort studies (eg, incident maternal breast cancer and insulin-like growth factor during early pregnancy³). Cases are cohort members who develop the outcome of interest, and incidence density or “risk-set” sampling is used to identify eligible controls. In addition to the case, each risk set includes all cohort members who are under observation and outcome-free at the time of the case's failure. The nested case-control study design has been shown to produce unbiased risk estimates under incidence-density sampling, in which controls are sampled at random (without replacement) from each risk set.^4–5 Variants of the nested case-control study design include randomly selecting m controls from the risk set for each case (simple 1: m design) or matching on additional factors such as sex and race (matched 1: m design).

Potential loss of information compared with the full cohort analysis is offset by gains in computational efficiency and reduced study costs. The nested case-control study design requires fewer evaluations of time-dependent covariates (eg, cumulative exposure or cumulative dose), which must be evaluated within risk sets at the case's failure time. Additional cost savings can result if further data collection is required, such as additional questionnaires or analysis of stored biologic samples.

When constructing risk sets, candidates for the time scale include calendar year, age, and time since cohort entry. Age is frequently used, because controlling for confounding from age is a primary concern.⁶ Analysis of the age-based risk sets via conditional logistic regression is asymptotically equivalent to a Cox proportional hazards regression analysis with age as the time scale.⁷ Use of age as the time scale precisely controls for age-related effects, but not necessarily for other time-related effects, which can be controlled in the analysis by including the appropriate terms in the regression model. This technique is often preferred over further matching of controls to cases.⁸

In studies of chronic diseases including cancer, exposure lagging (in which exposures occurring in the most recent L years are discounted) is often used to account for disease latency.⁹ For example, a nested case-control study of incident breast cancer among Finnish flight attendants (1:4 matched on year of birth [±1 year]) incorporated a 10-year exposure lag period for cumulative occupational radiation dose “to allow for the induction period of at least 10 years for radiation-induced solid tumors.”¹⁰ Another nested case-control study (1:5 matching on race and attained age) of male beryllium-exposed workers at a Reading, PA facility (hereafter, “the Reading cohort”) considered several metrics of beryllium exposure (duration, cumulative, average, and maximum) and 3 exposure lag periods (0, 10, and 20 years).¹¹ A subsequent reanalysis¹² of the nested case-control data identified confounding from birth cohort. Although this confounding was not considered in the original analysis, most positive associations between lagged beryllium exposure and lung cancer remained after adjusting for birth cohort.¹²

Critics of the study design have claimed that the results were biased, in part, due to the use of risk-set sampling in conjunction with exposure lagging.^13,14 Two simulation studies to evaluate the nested case-control study design in conjunction with exposure lagging did not observe bias.^15,16 Langholz and Richardson¹⁵ further described the theoretical basis for nested case-control studies both with and without exposure lagging. These studies have nonetheless been criticized because they did “not fully model important associations” in the Reading cohort.¹⁷

Here we report the results of a third simulation study that evaluates incidence-density sampling in conjunction with exposure lagging. The main study objective was to determine whether there is bias in cumulative exposure effect estimates from the Cox proportional hazards regression model in simulated cohorts in which risk was related to age, birth cohort, and lagged cumulative exposure. The simulated cohorts included associations among age at hire, year of hire, and exposure intensity. We also compared methods of adjusting for confounding by birth cohort and methods of treating lagged-out cases and controls, and we assessed the effect of incomplete follow-up.

METHODS

Cohorts were simulated based on methods developed by Richardson and Loomis,¹⁸ using SAS Software (version 9.2, SAS Institute Inc., Cary, NC). Each worker was randomly assigned a year of hire and age at hire (described below), a maximum exposure duration (exponential distribution [mean, 5 years]), and, to introduce censoring due to loss to follow-up, a maximum follow-up length (100 years minus an exponential random variable [mean, 15 years]). Three methods of assigning year of hire and age at hire were evaluated:

1. Year of hire was a uniform random variable (1935–1974) and age at hire was 18 years plus a gamma random variable (α = 1.6, β = 7.5).
2. Year of hire and age at hire as in (1) if hire year ≥1945; otherwise, age at hire was 18 years plus 8 times an exponential random variable (mean, 1 year) with probability 0.4 or uniform (18-65 years) with probability 0.6.
3. Year of hire and age at hire taken directly from the 3569 workers in the Reading cohort.¹¹

The second method was designed to produce an association between year of hire and age at hire that was similar to the Reading cohort.¹⁸ The third method was added because hire year and hire age distributions and associations would exactly mimic the Reading cohort. For each worker, birth year was calculated from the assigned hire age and hire year.

Exposure intensities were randomly assigned to each worker-year, using a lognormal distribution with a geometric mean of 2 units for workers hired in 1935, with a 2.5% decrease in each subsequent year, and a geometric standard deviation of 2.2. Supplemental analyses considered different methods of assigning the geometric mean: constant (one unit) for all calendar years and increasing (one-half unit in 1935 with a 2.5% increase in each subsequent year).

In the Cox proportional hazards model, the effect (φ) for a continuous exposure variable is the log hazard ratio for a 1-unit exposure increase. Four risk scenarios defined by the true exposure effect and lag period were considered:

1. φ = ln(1).
2. φ = ln(2) and lag = 0 years.
3. φ = ln(2) and lag = 10 years.
4. φ = ln(2) and lag = 20 years.

Each worker was followed from his hire year through his randomly assigned maximum follow-up. At each follow-up year, cumulative exposure under the assigned exposure lag period (CE_L, unit-years) was evaluated as:

where, intensity(t) was the randomly assigned exposure intensity for year t, and L was the exposure lag period. Each worker-year was assigned a probability of mortality from lung cancer h(t) and a probability of mortality from all other causes, c(t), that were based on published age- and year-specific death rates for US white males.¹⁹ Probability h(t) was a function of CE_L(t), φ, age(t), and birth year, conditional on survival to that age:

In equation (2), cumulative exposure was divided by 10 so that the exposure effect was the log hazard ratio for a 10 unit-year increase in cumulative exposure. Probability c(t) was a function of age(t) and birth year, conditional on survival to that age:

Models (2) (without the CE_L term) and (3) each accounted for 98.8% of the variability in the natural-log-transformed US white male death rates due to lung cancer and all other causes, respectively. These probabilities were used to assign Bernoulli random variables indicating the occurrence of death (lung cancer or all other causes combined) in a given follow-up year. Workers not assigned to die were censored at the end of their maximum follow-up. Workers assigned to die were cases if their earliest death was from lung cancer; otherwise, workers were censored at their earliest death.

Complete follow-up included all outcomes and person-time through death or loss to follow-up. Approximately 4% of the subjects were lost to follow-up under the null risk scenario. In addition, a study end year of 1990 was implemented to evaluate the effect of incomplete follow-up. In these analyses, workers who had not died or been lost to follow-up by 1990 were administratively censored at the study end date (31 December 1990). Approximately one-third of deaths were expected to occur through 1990 under the null risk scenario.

Five hundred cohorts were generated for each of the 36 scenarios (3 hire year/age methods, 3 exposure-intensity methods, and 4 risk scenarios). Each cohort, consisting of 4000 (hire year/age methods 1 and 2) or 3569 (method 3) workers, was analyzed using Cox proportional hazards regression (SAS PHREG procedure). Each analysis included all eligible controls (identified using incidence-density sampling by matching on attained age) rather than sample from the risk sets. Initial effect estimates were unadjusted for birth cohort. Three methods of adjusting for birth cohort were then compared:

1. Birth category (<1900 (reference), 1900–1919, 1920–1929, and 1930+), with cutpoints selected based on expected birth-cohort effects.
2. Birth year as a continuous variable.
3. Restricted cubic spline terms (ie, natural splines) for birth year (k = 4 knots at the 5th, 35th, 65th, and 95th percentiles of birth year for the cases).²⁰

Cumulative exposure, truncated at the case's age at death, was lagged by zero, 10, or 20 years for cohorts generated with exposure lag periods of zero, 10, or 20 years, respectively. In lagged analyses, 2 methods of treating lagged-out workers were compared:

1. Excluding cases who died in the first 10 (or 20) years of follow-up and, when risk sets were determined, excluding controls with fewer than 10 (or 20) years of follow-up at the case's age at death.
2. Including lagged-out cases and controls (assigned a cumulative exposure of zero unit-years).

Model fit was measured by Akaike's Information Criterion and Schwarz's Bayesian Criterion.²⁰

Exposure effect estimates were summarized for each simulation scenario. Summary statistics included the mean of the estimated hazard ratios and the empirical 5% type I error rates (ie, the fraction of cohorts for which the 95% confidence interval for the log hazard ratio did not contain the true value). We used a 2-sided z-test for proportions to determine if the empirical 5% type I error rate was significantly different from the nominal 0.05.

RESULTS

As expected, year of hire and age at hire were not correlated in cohorts simulated under hire year/age method 1 (mean correlation: 0.00), whereas hire year and hire age were negatively correlated in cohorts simulated under hire year/age method 2 (mean correlation, −0.18) and method 3 (correlation, −0.13) (Table 1). Also as expected, birth year and unlagged cumulative exposure were not correlated in cohorts simulated with constant exposure intensity (mean correlations, 0.01 to 0.03). In contrast, birth year and unlagged cumulative exposure were negatively correlated in cohorts simulated with decreasing exposure intensities (mean correlations, −0.23 to −0.10) and positively correlated in cohorts with increasing exposure intensities (mean correlations, 0.14 to 0.20). For scenarios with a positive lag, the correlation between birth year and lagged cumulative exposure increased with the lag period. For example, for simulations based on the Reading cohort and a decreasing exposure intensity, mean correlations between birth year and cumulative exposure were −0.10, −0.02, and 0.09 for lag periods of 0, 10, and 20 years, respectively.

TABLE 1: Mean Correlation Coefficients for Hire Year and Hire Age SYMBOL ₁ and for Birth Year and Cumulative Exposure SYMBOL ₂, by Hire Year/Age Method, Exposure Intensity, Hazard Ratio, and Lag Period

Figure

Cox proportional hazards regression analyses of cohort data simulated under the decreasing exposure intensity method are summarized in Table 2 for complete follow-up. For hire year/age method 1, there was no confounding by birth cohort and the empirical 5% type I error rates were not inflated. For hire year/age method 2, some confounding was observed for unadjusted analyses, but only for the positive, unlagged scenario (mean hazard ratio 1.95 compared with 2.0; 10.8% empirical type I error rate). For birth-cohort- adjusted analyses, the empirical 5% type I error rates were preserved for category of birth and restricted cubic spline adjustments, but with some minor residual confounding for continuous year-of-birth adjustment. For hire year/age method 3 (Reading cohort), results were similar to method 2, but minor bias was also observed for unadjusted analyses in the positive, 20-year lagged scenario when lagged-out workers were included in the analysis (mean hazard ratio 2.06 compared with 2.0) with an empirical 5% type I error rate of 10.2%. No bias was observed in any analysis that included an adjustment for birth-year category or included restricted cubic spline terms for birth year, regardless of lagging or the treatment of lagged-out workers; some bias was observed for continuous year-of-birth adjustment under a 20-year lag period.

TABLE 2: Mean Hazard Ratios and Empirical Type I Error Rates (α = 0.05) for Cohorts Simulated Under a Decreasing Exposure Intensity for Complete Follow-up (Through Death or Loss to Follow-up), by Hire Year/Age Method, Hazard Ratio, and Lag Period

Similar results were observed when outcomes and person-time were administratively censored in 1990 (Table 3). However, greater bias was observed in unadjusted analyses in the positive, 20-year lagged scenario for analyses that included lagged-out workers (mean hazard ratio 2.13 compared with 2.0; 21.4% empirical 5% type I error rate) or excluded lagged-out workers (mean hazard ratio 2.07 compared with 2.0; 9.6% empirical 5% type I error rate). For these analyses, including restricted cubic spline terms for birth year produced unbiased estimates.

TABLE 3: Mean Hazard Ratios and Empirical Type I Error Rates (α = 0.05) for Cohorts Simulated Under a Decreasing Exposure Intensity for Incomplete Follow-up (Through the Earlier of Death, Loss to Follow-up, or Administrative Censoring in 1990) by Hire Year/Age Method, Hazard Ratio, and Lag Period

Results for constant and increasing exposure intensity scenarios (eAppendix, eTables 1-4 [https://links.lww.com/EDE/A452]) were similar, with inflated empirical type I error rates and some bias for unadjusted analyses. The maximum observed bias for unadjusted analyses was approximately 15% (mean hazard ratio 2.31 compared with 2.0) for the scenario based on the Reading cohort, increasing exposure intensity, 20-year lag period (including lagged-out workers), with follow-up through 1990.

Based on Akaike's Information Criterion, analyses that included a birth-cohort adjustment (any of the methods evaluated) generally fit better than unadjusted analyses (results not shown). This was particularly true for cohorts simulated under hire year/age methods 2 and 3 and for complete follow-up. The best-fitting method varied among the adjustment methods evaluated, but the restricted cubic spline method was generally superior to the others. According to Schwarz's Bayesian Criterion, which incurs a greater penalty for increasing the number of model parameters, unadjusted analyses were generally better fitting for cohorts simulated under hire year/age method 1. For cohorts simulated under hire year/age methods 2 and 3, the best fitting method varied, but continuous metrics (birth year, restricted cubic splines) were generally better fitting than the categorical metrics (category of birth year).

DISCUSSION

Risk-set sampling in conjunction with exposure lagging was evaluated by Langholz and Richardson¹⁵ in a simulation study based on the Colorado Plateau uranium miners cohort. In this study, no bias was observed when a simple random sample of controls was selected from the risk set for each case identified using attained age alone.¹⁵ Hein et al¹⁶ simulated occupational cohorts with time-dependent exposures to evaluate incidence-density sampling in the presence of exposure lagging. In these simulations, exposure effect estimates from analyses that matched on attained age alone were unbiased, and, in lagged analyses, including lagged-out cases and controls (assigned zero lagged exposure) did not introduce bias.¹⁶ Yet, in a rejoinder to these 2 studies, Deubner and Roth¹⁷ suggested that the Reading cohort data “display a complex web of associations” and expressed a reservation “that these simulations do not fully model important associations” in the data. Although the Hein et al¹⁶ simulated cohorts were fairly realistic, with loss-to-follow-up censoring and with mechanisms for assigning dates and causes of death based on age and cumulative exposure, they did not include associations among hire year, hire age, and exposure intensity that have been observed in the Reading cohort.¹¹

In our simulations, birth-cohort effects were integrated into the risk model and associations among other time-related variables were introduced by the distributions of hire year, hire age, and exposure intensity. That is, there was potential for bias from confounding factors; however, the presence of birth-cohort effects did not necessarily require adjustment in the analysis to obtain unbiased exposure effect estimates. Indeed, for simulated cohorts with age at hire independent of year of hire, unadjusted analyses fit the best and resulted in little bias, regardless of the method of assigning exposure intensities. For simulated cohorts with hire age negatively correlated with hire year, adjustment for birth cohort was sometimes necessary to obtain unbiased estimates. In these simulations, however, failure to control for birth-cohort effects did not result in substantial bias. The maximum bias observed for any scenario was only 15% in the hazard ratio. Furthermore, the need for adjustment was not limited to scenarios for which the exposure-response relation was lagged. Incorporating an exposure lag period requires consideration of the treatment of lagged-out cases and controls. In these simulations, the inclusion of lagged-out cases and controls did not introduce bias as long as the analysis included an adjustment for birth cohort. Excluding lagged-out cases and controls reduced, but did not necessarily eliminate, bias associated with unadjusted analyses.

In these simulations, controls were not “matched” to cases on birth year; rather, birth-cohort confounding was controlled by including terms in the regression model because adjusting for differences in potential confounders is generally recommended over matching.⁸ Options for adjusting for continuous confounders include categorization, the inclusion of higher-order terms, and the use of regression splines. In these simulations, the type I error rate was always preserved and results were unbiased when restricted cubic splines were used to adjust for birth cohort. Generally, including terms for category of birth year was unbiased, but including birth year as a linear term in the model occasionally resulted in bias (likely from residual confounding). This was not unexpected because birth-cohort effects are highly nonlinear. Here, the risk model (based on actual death rates for US white men) included first- and second-order terms for birth year. Restricted cubic splines have been promoted as one method to minimize residual confounding when adjusting for a continuous covariate.^21–23

Compared with Hein et al,¹⁶ the cohorts simulated here were more realistic because the risk model incorporated calendar-year effects and because exposure intensity was varied over time and within-worker. In addition, relations among year of hire and age at hire approximated (method 2) or represented exactly (method 3) relations in the Reading cohort. Greater degrees of bias were observed for simulated cohorts based on the Reading cohort compared with method 2, likely due to differences in the distribution of hire year (uniform in the simulated cohorts, but highly right skewed in the Reading cohort).

Simulation studies have been used to evaluate the nested case-control study design in a variety of settings.^24–26 In fact, Deubner et al²⁷ conducted their own “empirical evaluation” of incidence-density sampling in conjunction with exposure lagging using the same Reading cohort data. In their evaluation, “case” status was randomly assigned to 142 cohort subjects, all subjects were assigned an average exposure of 20 units, and incidence-density sampling was used to identify 5 controls from the risk set for each case. Because they observed that cumulative exposure was higher among “cases” compared with “controls” when a 20-year lag period was implemented, they concluded that “the study design produced a biased case-control lagged exposure difference under the null hypothesis and could not distinguish qualitatively between null and alternate [sic] hypotheses.”²⁷

This “evaluation,” however, was fundamentally flawed.^28,29 Langholz and Richardson²⁸ noted that the empirical evaluation detected a “real association between length of employment and age at exit from the study” because case status was randomly assigned to workers but age at exit and duration of employment were not. Because their explanation did not consider birth cohort, we provide here an alternative. Because Deubner et al²⁷ randomly assigned case status to the entire cohort, workers who had not died by the study end date of 31 December 1992 were allowed to be selected as cases. Each person in the risk set for a case who was alive at the end of follow-up must have been born earlier than the case. For example, if a case was 50 years old and alive at the study end date, then every control in the risk set must have been either (i) alive as of 31 December 1992 and at least 50 years old on that date or (ii) deceased prior to this date and at least 50 years old at death. In both situations, the control would have had to have been born before 31 December 1942 (the case's birthdate). Thus, as long as a significant proportion of the cohort had not died by study end, controls will have an average date of birth earlier than cases using the methods of Deubner et al.²⁷ In the Reading cohort, approximately 55% of the cohort was still alive as of the study end date; consequently, the empirical evaluation of Deubner et al²⁷ actually introduced a birth-cohort effect, that was not considered in their analysis. We replicated their results (data not shown) and then considered an alternative analysis that randomly assigned lung cancer case status to “deceased,” rather than all, cohort members. In the alternative analysis, there was no association between lagged exposure and case status. This flaw in the design of the empirical evaluation²⁷ negates the authors' conclusions regarding the nested case-control study design with lagged exposures.

Furthermore, the “P values” reported in Table 2 in Deubner et al²⁷ are invalid because they were obtained by performing t tests on the combined data set from their 10 trials. Combining data from the 10 trials artificially inflates the sample size by a factor of 10, and hence decreases the P values, thus making the results appear to be more significant than they might otherwise be. It would be more correct to perform and summarize 10 t tests, or, even better, 10 Cox regressions. Additionally, Monte Carlo simulation studies are generally conducted with more than 10 simulations.

The nested case-control study design has been recommended as an efficient study design.^30–32 Theoretical arguments have shown that the nested case-control study design is unbiased under simple random sampling of controls (without replacement) from the risk set for each case, as long as incidence-density sampling is used.^7,33 In addition, the nested case-control study design has been shown to be unbiased for other sampling strategies that incorporate matching on potential confounders or counter-matching on an exposure surrogate (as long as the appropriate sampling weights are included).^34,35 Although these simulations analyzed the full cohort rather than simple nested case-control samples, the results lend support to the unbiased properties of exposure effect estimates obtained from incidence-density sampling.

In summary, for simulated occupational cohorts with birth-cohort effects and time-related associations among year of hire, age at hire, and exposure intensity, incidence-density sampling with matching on attained age was shown to produce unbiased exposure effect estimates as long as the analysis included an adjustment for birth cohort. The results of this study lend support to the use of incidence-density sampling and the nested case-control study design. We observed biased effect estimates when there was no adjustment for birth cohort, particularly when year of hire and age at hire were correlated. Effect estimates were unbiased when adjusted for birth-cohort effects, and, as before, including lagged-out cases and controls did not introduce bias.