In 2001, Sanderson et al at the National Institute for Occupational Safety and Health (NIOSH) published results from a nested case-control study of beryllium exposure and lung cancer mortality.^{1} Increased lung cancer mortality was observed among workers with higher estimated beryllium exposures when cumulative exposure was log transformed and lagged by 10 or 20 years. Five controls were selected for each case using incidence-density sampling by matching on race and attained age. The study design was immediately criticized by Deubner et al in a Letter to the Editor, which stated “ … the method of selecting control subjects introduced a serious bias” and “Matching on attained age biases the sample with respect to age entered into the cohort (ie, age of hire), since the older the age of initial eligibility, the more likely is survival to any subsequent age.”^{2} Deubner et al expressed concern that controls selected by matching on attained age would be more likely to have no exposure under a lag period because they had an older mean hire age.^{2} The mean hire age for controls was 37 years compared with 33 years for cases.^{1}

In 2007, Levy et al published a reanalysis of the NIOSH beryllium case-control data in which controls were matched to cases using a nonstandard criterion.^{3} Specifically, this nonstandard criterion required eligible controls to match cases on race, attained age, and age at death or censor (within 3 years of the case's age at death).^{3} Levy et al advised “Our discussion regarding the matching issue applies only to the use of risk set sampling in combination with exposure lagging … we would recommend closer matching of cases to controls in studies such as the one considered here, which does consider latency … ”^{3} The extra matching criterion reduced the exposure odds ratio (OR) from 1.06 (95% confidence interval [CI] = 1.00–1.12) to 0.99 (0.91–1.08).^{3} NIOSH investigators expressed concern about this extra matching criterion in a Letter to the Editor, stating “Information about the control that occurs after the time (age at death) of the case may not be used to select controls without introducing a potential bias.”^{4} In response, Levy et al stated, “There is, to our knowledge, nothing stated in the literature that one can not perform analysis on subgroups of the risk set on which the controls are selected so long as the findings are extrapolated to the appropriate target population.”^{5} Furthermore, in response to a 2008 NIOSH reanalysis of the beryllium case-control data,^{6} Deubner et al continued to advocate “closer matching of cases and controls on age at censor” to “reverse” the effect of “study design artifacts” on the exposure odds ratio.^{7}

To address this criticism of the nested case-control study design, we performed a simulation study to determine whether matching on the case's age at death in addition to attained age introduces bias in the risk estimates. We invoked statistical simulation since it can be difficult to ascertain methodological artifacts or bias in existing cohorts. Simulation studies have the advantage that the true exposure-response relation is known, and various statistical methods can be evaluated side-by-side.

## BACKGROUND

Occupational cohort mortality studies assemble cohorts of workers to estimate the effect of exposure on mortality.^{8} If additional exposure or covariate information is required from the workers, it may become necessary for both feasibility and financial reasons to study a subset of the cohort, known as a nested case-control sample, especially when the cohort is large or the outcome is rare. Workers who die of the outcome of interest are considered cases. Not all workers will die of the outcome of interest during the study period; workers may die of other causes, become lost to follow-up, or be alive at the study end date. For each case, the risk set includes all workers “at risk” at the time of the case's death (ie, workers who entered the cohort prior to the case's death and remain alive at the time of the case's death). Incidence-density sampling involves randomly sampling controls from these risk sets. This approach is frequently used for nested case-control studies because it has been shown to produce unbiased risk estimates.^{9–10} Some consequences of incidence-density sampling are that controls in the risk set for a particular case may become future cases and workers can appear in multiple risk sets.

The selection of an appropriate time scale (eg, follow-up time, calendar time or age) for defining the risk sets depends on the nature of the study. For occupational cohort studies with cancer outcomes, controls are usually matched to cases based on attained age because age is one of the most important risk factors for cancer mortality and is related to exposure.^{11} Matching on attained age requires that eligible controls be at risk for the outcome of interest (ie, under observation and alive) at the case's age at death. When controls are matched to cases based on attained age and the entire risk set is used for each case, nested case-control studies analyzed using conditional logistic regression yield results that are equivalent to Cox proportional hazards regression with age as the time scale.^{12} When age is used as the time scale, however, it may be necessary to include terms in the model for additional time-related variables, such as calendar year.^{9,12}

Additional matching criteria such as sex, race, or birth year may be used to reduce confounding, or these variables may be controlled in the analysis by including the appropriate regression variables. Covariate information may be time-varying (eg, cumulative exposure, tobacco use) or independent of time (eg, sex, race); however, within risk sets, covariate information must be evaluated at the time (ie, age) of selection in order for results to be unbiased.^{13} The additional matching criterion proposed by Levy et al requires that eligible controls' age at death or censor match the case's age at death, requiring the investigator to know something about the control that happens after the case's age at death. Consequently, matching on age at death may produce biased results.^{13}

The primary objective of this paper was to use simulated occupational cohorts to determine whether matching the control's age at death or censor on the case's age at death introduces bias in nested case-control studies. We considered risk scenarios for unlagged cumulative exposure and cumulative exposure under a 10-year lag period (ie, discounting exposures occurring in the most recent 10 years). When a lag period is incorporated into the analysis, workers who die or are otherwise censored before the end of the lag period are considered “lagged-out”; therefore, a secondary objective was to determine whether the risk estimates depend on the treatment of lagged-out workers. A tertiary objective was to examine risk estimates when the lag period was not correctly specified in the analysis.

## METHODS

Simulations were conducted using SAS Software (version 9.1.3, SAS Institute Inc., Cary, NC).

### Cohort Simulation

The algorithm for simulating the cohorts was based on methods developed by Richardson and Loomis who simulated occupational cohorts with time-dependent exposures to assess the impact of exposure categorization for grouped analyses of cohort data.^{14} We considered 10 simulation scenarios defined by the distribution of exposure intensity, the exposure effect, and the lag period. We evaluated 2 exposure intensities (constant [1] and exponential [mean, 1 unit]) and 3 exposure effects (null, positive, and negative). For non-null exposure effects, we evaluated lag periods of 0 and 10 years. We generated 500 cohorts for each simulation scenario and 4000 workers for each cohort. Each worker was randomly assigned values for age at first exposure (18 years plus an exponential random variable [mean, 10 years]), maximum exposure duration (exponential [mean, 25 years]), and maximum follow-up time (40 years minus an exponential random variable [mean, 5 years]); these distributions were used by Richardson and Loomis to simulate cohorts similar to an existing nuclear industry worker cohort.^{14}

Each worker was followed for their randomly assigned maximum follow-up time. At each follow-up year, we evaluated the worker's current age (years) and cumulative exposure (intensity multiplied by exposure duration) under the lag period (CE_{L}, unit-years). Each worker-year was assigned a probability of mortality from the outcome of interest as a function of CE_{L}, the exposure effect (φ), and age, conditional on survival to that age:

The exposure effect (φ) was *ln*(1), *ln*(1.5), and *ln*(0.67) for the null, positive, and negative models, respectively. Each worker-year was also assigned a probability of mortality from all other outcomes as a function of age alone, conditional on survival to that age:

Specific parameters for these conditional probabilities (hazard rates) were used by Richardson and Loomis.^{14}

These conditional probabilities were used to assign Bernoulli random variables for the outcome of interest and all other outcomes combined. A Bernoulli value of one indicated that a “death” occurred in that year. Workers not randomly assigned to “die” were censored at the end of their randomly assigned maximum follow-up time. Workers randomly assigned to “die” were cases if their earliest “death” was from the outcome of interest, otherwise they were censored at their earliest death. The final cohort data included beginning and ending ages for exposure, exposure intensity, beginning and ending ages of follow-up, and case status at the end of follow-up.

### Data Analysis

Exposure effects were estimated using conditional logistic regression (SAS PHREG procedure). Controls were selected using incidence-density sampling.^{15} Two methods for constructing risk sets were compared: matching controls to cases based on attained age alone and matching controls to cases based on both attained age and age at death (ie, eligible controls had to additionally die or be otherwise censored within 3 years of the case's age at death). All eligible controls were included; however, in supplemental analyses, 5 controls were randomly selected without replacement from the risk set for each case (5:1 matching). Model fit was measured by the −2 log likelihood statistic. Initially, cumulative exposure was lagged by 0 and 10 years for cohorts generated with exposure lags of 0 and 10 years, respectively (ie, the lag period was correctly specified). Next, a 10-year lag period was applied to cohorts generated with a lag of 0 years and vice versa to explore misspecification of the lag period. Two methods of treating lagged-out workers were compared. In the first, cases that died in the first 10 years of follow-up were excluded and, when risk sets were determined, workers with less than 10 years of follow-up at the case's age at death were considered ineligible. In the second, lagged-out cases and controls were included and assigned a cumulative exposure of zero. Cumulative exposure was truncated at the case's age at death.

### Expression of Results

We summarized characteristics (age at first exposure, truncated cumulative exposure, and age at death or censor) of the 5:1 sampled risk sets for matching based on attained age alone and matching based on attained age plus age at death or censor. We summarized estimated exposure effects for each simulation scenario. For analyses that included all eligible controls, the estimated exposure effect can be interpreted as the log hazard ratio (HR) for a 10 unit-year increase in cumulative exposure. Summary statistics included the mean and standard deviation of the estimated HRs, bias (the mean percent difference between the empirical and theoretical HRs), and the 5% and 1% rejection frequencies (for evaluating the type I error rate). For supplemental analyses that employed a 5:1 matching criterion, the estimated exposure effect can be interpreted as the log OR for a 10 unit-year increase in cumulative exposure.

### Additional Analyses

In these simulations, workers were followed for no more than 40 years. The effect of varying the maximum length of follow-up was evaluated in supplemental analyses by specifying maximum follow-up times of 20, 30, 50, and 60 years for the simulation scenario of greatest interest (exponential intensity, positive effect, and 10-year lag). Additional supplemental analyses also considered varying the distributions of age at first exposure and maximum exposure duration for this scenario.

## RESULTS

Characteristics of the sampled risk sets are summarized in Table 1. In the absence of an exposure effect, we expected cases and controls to have similar mean exposures. For constant intensity, the mean cumulative exposure was similar for cases and controls selected by matching on attained age alone (mean difference = −0.022 unit-years; 95% confidence interval [CI] = −0.098 to 0.053); however, controls selected by matching on attained age plus age at death were younger at first exposure (−5.26 years; −5.37 to −5.14) with higher cumulative exposures (mean difference = 1.86 unit-years, 95% CI = 1.79 to 1.94) compared with cases. Similar results were observed under the null model for exponential intensity.

For positive and negative exposure effects, we expected cases to have mean exposures compared with controls that were higher and lower, respectively. For the positive scenarios, the mean cumulative exposure for cases was substantially higher compared with controls selected by matching on attained age alone (eg, for constant intensity, the mean difference was 4.52 unit-years [4.46 to 4.58] and 2.61 unit-years [2.56 to 2.66] for 0- and 10-year lag periods, respectively). Likewise, for the negative scenarios, the mean cumulative exposure for cases was lower compared with controls selected by matching on attained age alone (eg, for constant intensity, the mean difference was −4.12 unit-years [−4.21 to −4.04] and −2.41 unit-years [−2.47 to −2.34] for 0- and 10-year lag periods, respectively).

For all scenarios, controls selected by matching on attained age plus age at death were markedly different from controls selected by matching on attained age alone. Under the extra criterion, controls were, on average, approximately 5 years younger at first exposure with higher, on average, truncated cumulative exposures compared with controls selected by matching on attained age alone. For example, under the constant intensity, null model, the mean age at first exposure difference was −4.95 years (−5.06 to −4.85) and the mean truncated cumulative exposure difference was 1.84 unit-years (1.79 to 1.89).

Risk set construction is illustrated in Figure 1 for a case that died at age 64 years from a cohort simulated under the null model with exponential intensity. Age at risk end, necessarily greater than age at first exposure, is plotted against age at first exposure so points are bounded below by the identity line (y = x). In these simulations, no worker was followed for more than 40 years, so points are bounded above by the line y = x + 40. Eligible controls for matching on attained age lie in the upper left quadrant defined by the case's age at death. Under the additional criterion that eligible controls have an age at death or censor within 3 years of the case's age at death, the pool of eligible controls is reduced (Fig. 2).

Simulation results evaluating the 2 methods for constructing risk sets are presented in Table 2. All eligible controls were included and the lag period was specified correctly. For unlagged analyses, selecting controls by matching on attained age resulted in low bias (mean percent difference = −1.5% to 0.1%) and generally preserved the 5% type I error rate (3.2% to 5.6%). For lagged analyses, selecting controls by matching on attained age and including lagged-out workers resulted in low bias (mean percent difference = −0.5% to 0.9%) with slightly higher empirical 5% type I error rates (3.8% to 6.4%); excluding lagged-out workers resulted in slightly higher biases (mean percent difference = 0.2% to 4.1%). For every scenario, selecting controls by matching on attained age plus age at death or censor resulted in downwardly biased estimates and grossly inflated empirical type I error rates. The bias was greater for constant exposure intensities (mean percent difference = −32% to –8.4%) compared with exponential exposure intensities (−13.7% to −3.8%). For the scenario of greatest interest (positive model under a 10-year lag), the empirical 5% type I error rates were 100% under the extra matching criterion.

Simulation results evaluating misspecification of the lag period are presented in Table 3. All eligible controls were included, and controls were selected by matching on attained age alone. Incorporating a 10-year lag period in the analysis when the true lag was zero resulted in positively biased HRs for the positive exposure effect and negatively biased HRs for the negative exposure effect. The lag period that produced the highest estimate did not necessarily produce the best fitting model. For a positive exposure effect with a true lag of 0 years, the 10-year lagged analysis produced higher exposure estimates in nearly every cohort; however, the unlagged analysis usually provided a better fit to the data based on the −2 log likelihood. Similar results were observed for the negative models. On the other hand, when the true lag was 10 years, the unlagged analysis resulted in lower exposure estimates for the positive model and higher exposure estimates for the negative model.

Similar results (not shown) were observed for all scenarios in supplemental analyses that employed a 5:1 matching criterion rather than include all eligible controls. Additional supplemental analyses for the scenario of greatest interest (exponential intensity, positive exposure effect, and a 10-year lag period) evaluated the amount of bias under different distributions for maximum follow-up duration, age at first exposure, and maximum exposure duration. Similar results (not shown) were observed for matching on attained age alone when the maximum follow-up was increased from 40 up to 60 years; however, biases associated with matching on attained age plus age at death increased from approximately −10% (Table 2) to approximately −17%. Similar results (not shown) were observed when the distributions of age at first exposure and maximum exposure duration were varied.

## DISCUSSION

Lubin and Gail described biases in relative risk estimation associated with various control selection methods and concluded that “case-control analyses of cohort data lead to biased estimates of the relative risk … unless the controls for a case at time *t* are randomly selected from all those at risk at *t*.”^{9} Our simulations confirm this result. Incidence-density sampling based on attained age as the time scale did not introduce a systematic bias, whereas adding age at death or censor as a matching criterion resulted in downwardly biased estimates. The amount of bias was greater in simulations where exposure intensity was constant compared with exponentially distributed. In these simulations, risk is modeled as a function of exposure duration when intensity is constant and of cumulative exposure when intensity is exponentially distributed. The latter scenario is of greater interest, as it represents a more realistic exposure distribution. Levy et al reported an approximate 7% reduction in the lung cancer OR in the NIOSH beryllium cohort when controls were selected by additionally matching on age at death.^{3} This “reduction” is in line with the amount of bias that we observed under the extra matching criterion in cohorts simulated under an exponential intensity.

Matching on age at death, in addition to attained age, resulted in downwardly biased exposure estimates for every scenario (null, positive, and negative). However, our findings are limited to the actual simulation scenarios that were evaluated. In addition, we cannot determine, based on these simulations, whether matching on attained age plus age at death will always result in bias, or whether conditions exist for which the extra matching criterion might produce unbiased results. However, we varied the distributions of the randomly assigned age at first exposure, exposure intensity, maximum exposure duration, and maximum follow-up duration in supplemental analyses and observed similar results. In fact, we never encountered a scenario for which the extra matching criterion was unbiased.

The exact cause of the bias observed with the extra matching criterion cannot be identified using a simulation study; however, a consequence of matching on attained age plus age at death appears to be that sampled controls have a younger age at first exposure and consequently a higher truncated cumulative exposure compared with controls selected by matching on attained age alone. For incidence-density sampling based on attained age, eligible controls lie in a “triangle” defined by age at risk end, age at exposure begin, and the case's age at death (Fig. 1). Under the extra criterion, eligible controls must be selected from the “bottom band” of this triangle and are more likely to be younger at first exposure (Fig. 2). Within risk sets, cumulative exposure is truncated to the case's age at death; therefore, selecting controls with a younger mean age at first exposure results in the potential for a higher mean truncated cumulative exposure, which explains the downward bias observed with the additional criterion. Indeed, when we compared age at first exposure and cumulative exposure for controls selected by matching on attained age alone and attained age plus age at death, mean age at first exposure was substantially lower and mean truncated cumulative exposure was substantially higher in the latter compared with the former in every scenario evaluated.

Deubner et al^{2} and Levy et al^{3} expressed concern that controls selected using incidence-density sampling with 5:1 matching on race and attained age were older at first exposure than cases by about 4 years in the 2001 analysis of the NIOSH beryllium cohort.^{1} In these simulations, matching on attained age alone selected controls that were, on average, similar to cases when the underlying exposure-response relation was null, but older than cases at first exposure when the underlying exposure-response relation was positive (and younger when the association was negative). Overall, results obtained by matching on attained age alone were generally unbiased, so it appears that this concern is unwarranted.

Matching on attained age can, however, result in bias when the lag period is not specified correctly. When risk was assigned based on a true lag of 0 years, unnecessarily lagging in the analysis biased effect estimates away from the null. Not lagging in the analysis when risk was assigned based on a true lag of 10 years, however, biased effect estimates toward the null. Biases resulting from incorrect assumptions about the true lag period are introduced by exposure misclassification.^{16} The selection of an appropriate exposure lag can be difficult. Researchers may report results from several lag periods, results for a single lag period selected a priori based on theoretical considerations, or results from the lag period that produces the best fitting model.^{17–19} Salvan et al caution against using the highest-estimate criterion in which the lag that produced the highest estimate of risk is selected.^{16} Our results lend support to this caution. In addition, selection of the best fitting model may protect against bias from lagging cumulative exposure when the true lag period is zero. For the positive-unlagged scenario, the 10-year lagged analysis produced higher effect estimates but the unlagged analysis provided a better fit to the data in all (exponential intensity) or nearly all (constant intensity) of the simulated cohorts.

Additionally, we evaluated 2 methods of treating lagged-out workers. When risk was assigned based on a lag of 10 years, effect estimates from analyses that included lagged-out workers were generally similarly or less biased than those that excluded lagged-out workers. The treatment of lagged-out workers is an important design consideration, especially when the exposure-response relation is nonlinear and a log transformation is required.^{6} When risk was associated with a lagged cumulative exposure, we did not consider the effect of a log transformation nor did we consider other misspecifications of the lag. Additional work is needed to evaluate the effects of both lagging and logging exposure on risk estimates.

Although time-related (eg, calendar year, birth cohort) effects can be important in cohort studies, calendar year did not influence risk in these simulations, nor did it influence exposure intensity. Thus, it was not considered in any analyses. The original 2001 analysis of the NIOSH beryllium cohort did not consider additional time-related effects^{1}; however, the 2008 reanalysis by Schubauer-Berigan et al adjusted for birth cohort to account for known differences in smoking rates by birth year.^{8} Additional work is needed to evaluate the effects of time-related effects on the risk estimates when selecting controls by incidence-density sampling with matching on attained age.

Finally, we did not consider varying the length of the additional match criterion. All analyses used 3 years as suggested by Levy et al.^{3} We expect that a closer matching criterion (<3 years) would produce even greater bias and a relaxed criterion would produce less bias.

In summary, incidence-density sampling with matching based on attained age plus age at death was found to introduce a downward bias. This approach is not recommended for nested case-control studies. Although the results from the simulated cohorts cannot resolve the debate on the form and magnitude of the relation between beryllium exposure and lung cancer mortality in the NIOSH beryllium cohort, these results indicate that the correct relation cannot be ascertained using the extra matching criterion suggested by Levy et al.^{3} Incidence-density sampling with matching based on attained age alone was generally unbiased and including lagged-out workers did not introduce bias. When the true lag was zero, lagging cumulative exposure by 10 years biased exposure estimates away from the null. Not lagging cumulative exposure when the true lag was 10 years biased exposure estimates toward the null. Consequently, the decision to lag cumulative exposure in an analysis should be made based on a priori information and model fit.

## REFERENCES

*Am J Ind Med*. 2001;39:133–144.

*Am J Ind Med*. 2001;40:284–285.

*J Occup Environ Med*. 2007;49:96–101.

*J Occup Environ Med*. 2007;49:708–709.

*J Occup Environ Med*. 2007;49:709–711.

*Occup Environ Med*. 2008;65:379–383.

*Occup Environ Med.*Available at: [http://oem.bmj.com/cgi/eletters/oem.2007.033654v1]. Accessed November 16, 2007.

*J R Stat Soc A*. 1977;140:469–491.

*Biometrics*. 1984;40:63–75.

*Biometrics*. 1986;42:293–299.

*J Am Stat Assoc*. 1983;78:1–12.

*J Chronic Dis*. 1986;39:379–388.

*Occup Environ Med*. 2004;61:930–935.

*Occup Environ Med*. 2004;61:e59.

*Epidemiol*. 1995;6:387–390.

*J Natl Cancer Inst*. 2007;99:357–364.

*Radiat Res*. 2007;167:222–232.

*Occup Environ Med*. 2004;61:2–7.