The healthy worker survivor effect (HWSE) remains a potential source of bias even in occupational epidemiological studies with an internal comparison group. In this paper, we focus on how a study’s cohort definition can influence the HWSE and impact the exposure-response. We will do this with an example based on the Vermont granite workers.

The Vermont granite workers study has an important place in the literature involving silica and cancer.1^–4 Earlier analyses of the mortality experience of the Vermont granite industry relied on an external comparison group and job titles as surrogates for silica exposure.1 2 In 2004, Attfield and Costello incorporated quantitative exposure information into models for lung cancer mortality.4 A categorical analysis revealed an increasing trend in the dose-response, except for a decline in the highest exposure category (⩾6.0 mg-year/m³). The authors noted that 83% of those in the highest category of cumulative exposure had worked for at least 20 years by the start of follow-up, and that these self-selected survivors were likely to cause a healthy worker effect.

The HWSE results from less healthy workers transferring to jobs with lower exposure or leaving the workforce altogether. The healthier workers continue to be exposed (at higher levels), thus attaining higher cumulative exposures. As a result of this self-selection, the risk of disease at higher exposures appears to plateau or may even decrease. Sources of the HWSE bias and analytical approaches to reduce it in occupational epidemiology have been discussed previously.5^–7

Study design can impact the healthy worker effect. For example, cross-sectional studies of active workers over-sample longer term workers.8 Koskela et al9 used the term “cross-sectional follow-up” to describe follow-up studies of workers identified at one point in time. They found that these workers were healthier and had lower mortality rates than subjects captured as they enter the workforce (that is, open or inception cohort). Often, as is the case for the Vermont granite study, the cohort is a hybrid of these two approaches and includes both subjects already at work at the start of follow-up and others who entered the cohort when they were hired. We were interested in how the distinction between these workers in the study population exaggerated (or minimised) the HWSE in the Vermont granite industry study and impacted the association between silica exposure and lung cancer mortality.

We distinguish two types of subjects, referred to as prevalent hires and incident hires. Those employees already at work at the start of the start of follow-up (that is, part of a cross-sectional follow-up identified at one point in time) will be referred to as “prevalent hires”. “Incident hires” will describe those who first started working after the start of follow-up (that is, identified through the ongoing enrolment period of the inception cohort). The importance of this distinction is that to have made it into the Vermont granite workers study, prevalent hires must have endured the occupational exposure for some time and remained at work. Thus, prevalent hires may not be representative of all workers hired before the start of the study because those who were more susceptible to adverse health effects (either exposure related or unrelated) are more likely to have stopped working in the industry.10

The Vermont granite workers study provides an opportunity to examine the hypothesis that prevalent hires in the study population are a source of selection bias (HWSE) because more than half of the cohort was hired before the start of follow-up.

METHODS

Details on study design and exposure assessment have been described previously and are briefly reviewed here.1^–4

Cohort identification and follow-up

The most complete record of workers employed in the Vermont granite industry between 1950 and 1982 came from a medical surveillance programme. Granite workers who were x rayed at least once were eligible to be a part of this cohort study. A total of 5414 workers, all white males, were identified.

Information on vital status for 1950–82 was obtained from pension funds, the Vermont Division of Industrial Hygiene obituary notices of workers, and state vital records. For deceased workers, death certificates were obtained from the Quebec, New Hampshire, and Vermont Health Departments to determine the cause of death. When follow-up was extended to 1994, social security files on CD-rom were used to manually search for workers by last name to identify more deaths. A qualified nosologist coded the cause of death according to the Eighth Revision of the International Classification of Diseases. Study protocol was approved by the human subjects review committee at the Harvard School of Public Health.

Exposure assessment

Through their participation in the medical surveillance programme, granite workers provided work history information. Exposure assessment has been described in Attfield and Costello.4 Briefly, each worker’s quantitative silica exposure was estimated based on environmental surveys conducted between 1924 and 1977. Measured respirable silica dust concentrations for each job provided the basis for the exposure estimates. Cumulative exposure to respirable silica dust (mg-year/m³) was calculated for each granite worker by matching his work history with the estimated exposure for each job and year and summing the exposures over time. To be consistent with previous studies of lung cancer mortality and occupational exposures, we lagged exposure by 15 years.4

Classification of subjects by date of hire

The Vermont granite cohort is comprised of both incident hire and prevalent hire subjects. A schematic for potential subjects in the cohort is presented in figure 1. Incident hires—that is, those hired during the follow-up period (starting in 1950 and later)—were captured as they entered the workforce. Prevalent hires worked before the start of follow-up and continued working until the start of follow-up in 1950. Those who started and stopped working before 1950 will be missing. Thus, the probability of being missing is related to length of employment as well as to health status.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1 Schematic of Vermont granite workers cohort definition. The circles with lines extending to the right represent possible person-time experiences of Vermont granite workers, which ended either as a lung cancer death (×) or being censored (|). Study follow-up began in 1950 and ended in 1994. The original design captured subjects who had worked anytime in between, regardless of date of hire.

Data analysis

Cox proportional hazards regression was used to estimate the mortality rate ratio (MRR) and 95% confidence interval (CI) for the relationship between silica and lung cancer mortality (SAS version 9.1). Risk sets were generated using age as the timeline. We controlled for calendar year achieved risk age because lung cancer mortality was increasing over the time period of this study and silica dust exposures were decreasing. Calendar year was included as an ordinal variable representing increasing decade at risk age. We also adjusted for age at hire with a continuous variable. Smoking information was not available.

Cumulative exposures typically have non-linear relationships with occupational disease—for example, a plateau or decrease in risk among those with the highest exposures. Such patterns are generally attributed to the healthy worker effect.6 7 To avoid a linearity assumption, we initially analysed cumulative silica exposure (15-year lag) as a categorical variable. Quartiles based on the exposure distribution of the cases in the original cohort were used to define the categories. To compare the silica dose-response by hire status, the same categories for silica were also analysed after stratifying by incident or prevalent hire status. We tested for statistical interaction by comparing the log likelihood for the model with silica exposure, hire status, and their cross products to the log likelihood from a model without the cross products. However, there were differences in the distribution of exposure for incident and prevalent hires. Therefore, we also determined categories separately by hire status and analysed these for an association with lung cancer mortality. In categorical analyses stratified by hire status, cut-points were determined by the distribution of the cases in each stratum. A trend test was conducted for categorical analyses by including an ordinal term representing the categories with increasing cumulative silica exposure. All statistical tests were two-sided.

To more fully characterise the dose-response relationship for lung cancer mortality and silica, semi-parametric models based on penalised splines were fit using the approach suggested by Thurston, et al.17 We used penalised splines to allow for non-linearity because they are insensitive to the location or number of knots and have been applied in Cox models of occupational cohorts.16^–18 Spline models were generated using R 1.9.0. Currently, R is unable to generate a penalised spline for an interaction term in a Cox model. Therefore, penalised spline models were fit to prevalent and incident hires in a stratified analysis. Akaike’s Information Criterion (AIC) was used to guide the selection of the optimal number of degrees of freedom.19 A sensitivity analysis examined the impact of moving the date of hire relative to the start of follow-up. Exposure-response splines were generated for three cohort definitions: workers hired in 1940 or later (10 years before start of follow-up), hired in 1930 or later (20 years before start of follow-up), and workers hired in 1920 or later (30 years before start of follow-up).

RESULTS

A total of 5414 subjects were initially identified as having worked anytime between 1950 and 1982 and having at least one x ray. This surveillance programme is reported to have included 98% of the eligible workforce.20 Thus, 98% of workers hired from 1950 and later (the incident hires) were captured. However, there is no complete record of all workers who started and stopped work before 1950; therefore, we cannot say what proportion of people who worked before 1950 was represented by the prevalent hires in this study.

Six subjects were excluded due to incomplete demographic data (for example, date of birth). There were 213 lung cancer deaths among the remaining 5408 Vermont granite workers in the original cohort. Upon determination of risk sets for cases, 227 subjects did not contribute to any risk set either because they were censored before the youngest case, or they were older than the oldest risk age at the start of the follow-up period. A total of 5181 subjects were included in the risk sets, and 52% of these were prevalent hires.

Table 1 compares descriptive characteristics for prevalent and incident hires, which together comprise the original data set. Overall, subjects were approximately evenly divided between prevalent hires (n = 2603) and incident hires (n = 2578) while 74% of the lung cancer deaths were prevalent hires (prevalent hire cases = 157, incident hire cases = 56). On average, there was a 30-year difference in year of hire for prevalent and incident hires (1932 compared to 1961, respectively). Incident hires started working in the industry at a slightly older age than the prevalent hires and worked about 20 years less by the end of follow-up than prevalent hires. The duration of follow-up was not very different between the two groups, which shared the earliest possible date for the start of follow-up of 1 January 1950.

The most notable difference between prevalent and incident hires was the accumulation of exposure. Overall, prevalent hires had a higher mean cumulative silica exposure than the incident hires. This was due in part to their longer work tenure, but it was also influenced by the implementation of industry-wide control technologies in the 1940s, which reduced silica dust concentrations. Among incident hires, there was no difference in mean silica exposure between cases and non-cases (both 0.3 mg-year/m³). Among prevalent hires, however, the mean silica exposure was higher among non-cases (3.7 mg-year/m³) than cases (2.9 mg-year/m³).

The categorical models for cumulative exposure are presented in table 2. In the model for the original cohort, there was an elevated mortality rate with higher cumulative exposure, achieving more than a doubling in the mortality rate (MRR 2.54, 95% CI 1.58 to 4.06, p for trend  = 0.001) before the dose-response began to decline. Stratifying by hire status, the mortality rate for prevalent hires increased to approximately a doubling relative to the lowest exposed category (MRR 2.14, 95% CI 1.18 to 3.90) and there was a statistically significant trend (p = 0.02). However, for the incident hires, 68% of cases fell within the first quartile of exposure of the original cohort, only two (4%) in the third quartile, and none in the highest. Although a test for interaction can be conducted (results were not significant), it was not informative of a difference in the dose-response because there was so little overlap in the exposure distributions of the prevalent and incident hires. When including prevalent and incident hires in the same categorical model, the incident hires therefore served as the referent group for the prevalent hires (71.7% of person-time in the referent group was from the incident hires).

Table 2 also presents stratified categorical analyses with the cutpoints for the quartiles determined separately for incident and prevalent hires. For prevalent hires, the mortality rate ratio increased to 1.90 (95% CI 1.18 to 3.08) and then declined in the highest exposure category. Because the incident hires had lower power due to the number of cases and a narrower exposure range, we analysed tertiles of exposure. For incident hires, we observed a dose response which reached a mortality rate ratio of 2.51 (95% CI 1.14 to 5.55) with a statistically significant trend (p = 0.02) over the low exposure range (0.4 to 2.0 mg-year/m³).

Penalised splines were fit and predicted MRR plotted against the range of exposure (figs 2 through 4). For the original cohort (fig 2 and table 3), the MRR increased linearly up to a cumulative exposure of 4.0 mg-year/m³, where the MRR reached a maximum of 2.50 (95% CI 1.81 to 3.47) and then the MRR declined below the null. For prevalent hires (fig 3 and table 3), the shape of the dose-response curve was similar to the original cohort, although the maximum MRR was lower (1.97, 95% CI 1.52 to 2.56). Restricting to incident hires (fig 4 and table 3), the MRR increased up to 2.19 (95% CI 1.27 to 3.76) at a much lower cumulative exposure of 0.6 mg-year/m³. This was followed by a drop in the MRR, which increased dramatically, albeit with very wide confidence intervals, at the highest exposure of 2.0 mg-year/m³.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2 Mortality rate ratio for cumulative silica exposure and lung cancer using the original cohort definition for the Vermont granite workers study. Mortality rate ratio estimated by a Cox proportional hazards model with penalised splines. Model used age as the timeline and controlled for calendar year and age at hire. The original cohort was comprised of subjects who worked in the Vermont granite industry anytime between the start of follow-up on 1 January 1950 and 31 December 1982, regardless of date of hire (for example, includes prevalent and incident hires).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3 Mortality rate ratio for cumulative silica exposure and lung cancer among only the prevalent hires in the Vermont granite workers study. Mortality rate ratio was estimated by a Cox proportional hazards model with penalised splines. Model used age as the timeline and controlled for calendar year and age at hire. Prevalent hires were already at work prior to the start of follow-up on 1 January 1950 and were still at work on that date.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4 Mortality rate ratio for cumulative silica exposure and lung cancer among only the incident hires in the Vermont granite workers study. Mortality rate ratio was estimated by a Cox proportional hazards model with penalised splines. Model used age as the timeline and controlled for calendar year and age at hire. Incident hires started working in the Vermont granite industry between 1 January 1950 (the start of follow-up) and 31 December 1982.

To further examine the impact of including prevalent hires, we conducted a sensitivity analysis. Figure 5 and table 4 present what occurred when an increasing proportion of prevalent hires was included in the analysis. Penalised splines containing the incident hires with prevalent hires hired within 10 years before the start of follow-up (hired 1940 or later), 20 years before the start of follow-up (hired 1930 or later), and hired as much as 30 years before the start of follow-up (hired 1920 or later). The steepest slope was among those hired 1940 and later. As the date of hire for cohort eligibility was moved further away from the start of follow-up (for example, earlier calendar date) and allowed for a higher proportion of prevalent hires, the exposure range increased and became more skewed (fig 5).

• Download figure
• Open in new tab
• Download powerpoint

Figure 5 Comparison of dose-response curves for silica exposure and lung cancer mortality from the Vermont granite industry defined by a different date of hire eligibility: 1940 or later, 1930 or later, and 1920 or later. Follow-up for all subjects began in 1950 regardless of date of hire. Mortality rate ratios estimated with Cox proportional hazards models using penalised splines. Models used age as the timeline and controlled for calendar year and age at hire.

In table 4, as the percentage of prevalent hires decreased, the maximum MRR occurred at a lower cumulative exposure (from 4.55 mg-year/m³ for those hired 1920 or later (45.7% prevalent hires) to 2.25 mg-year/m³ among those hired 1940 or later (27.2% prevalent hires)). Furthermore, the maximum MRR was somewhat higher with fewer prevalent hires, although the lower limits of the confidence intervals were similar.

DISCUSSION

Prevalent hires were older at time of lung cancer death, younger at hire, and worked longer which allowed them to accumulate higher silica exposures than the incident hires. With penalised splines, we observed that the exposure-response increased and then declined in a non-linear fashion, whether incident hires and prevalent hires were combined or analysed separately. The main limitations to restricting the study to incident hires were the loss in power as the number of cases decreased by 74% and the dramatically reduced exposure range. A sensitivity analysis indicated that as the date of hire criterion moved closer to the start of follow-up, the exposure-response slope became steeper, allowing for a doubling in the mortality rate ratio to be reached at a lower cumulative exposure.

Two recent analyses of the Vermont granite industry cohort that included the extended follow-up period reached opposing conclusions about whether there was an association between lung cancer mortality and silica dust exposure.3 4 21 We demonstrated that HWSE bias may be exaggerated in the original cohort, which was used in both of the previous analyses, with the presence of prevalent hires potentially weakening the association between silica and lung cancer. When we accounted for this selection bias, the lung cancer mortality rate increased with increasing cumulative silica up to 2.25 mg-year/m³, even when including subjects who started working up to 10 years before start of follow-up. Indeed, if we used the strictest cohort definition limited to incident hires only, all of whom started working after dust control measures had been implemented, there was evidence of an exposure-response relationship, exhibited with the steepest slope and statistically significant association between silica and lung cancer mortality when cumulative exposure reached 0.6 mg-year/m³.

These findings are consistent with other occupational studies of lung cancer where a subset of incident hires was examined.22 23 For example, in a study of a metal production industry in Colorado, the strongest standardised mortality ratio for lung cancer mortality was achieved at a lower cumulative cadmium exposure for incident hires than for the prevalent hires.22 In another study of Tacoma smelter workers, only incident hires exhibited a linear dose-response with increasing arsenic exposure, resulting in a stronger dose-response at lower levels of exposure.23

Although HWSE will be minimised in a cohort composed of incident hires only, it is important to recognise that this bias will continue to be present because those who are less healthy will still be more likely to remove themselves from higher exposure earlier than those who are healthier. The use of lagging exposure has been suggested as a means of reducing HWSE.24 In addition to the 15-year lag, we also examined models without a lag and with shorter lag periods (5 and 10 years), and the same pattern of stronger dose-response for incident hires persisted. Adjusting for time since hire or employment status has also been promoted as a means of controlling for this ongoing HWSE. However, there was no change in the estimates when time since hire or employment status was included in the model. Further, controlling for these time-varying factors, which in addition to their confounding role may also be intermediate variables, would induce bias.25^–27 As a result, structural nested models have been proposed to address HWSE.28 29 A greater discussion of methods that address the ongoing HWSE in occupational cohorts is beyond the scope of this analysis.

An alternative explanation for the observed differences in the exposure-response relationships between prevalent and incident hires is measurement error. The cumulative exposure estimates were derived from six environmental surveys conducted between 1924 and 1977. It could be argued that exposure measurements improved with time, thereby reducing exposure misclassification and its associated bias towards the null. As a result, the stronger association with fewer prevalent hires could be an artifact of less measurement error. However, even if measurement error decreased with time, the bias induced by including prevalent hires would still be present. Furthermore, to the extent that measurement error decreased over time, focusing on incident hires will reduce error.

Incident hires were older, on average, than prevalent hires at the time of first employment. This difference was larger for cases than non-cases. Consider two subjects, both hired at age 25. The prevalent hire started work in 1940 and the incident hire started in 1960. The prevalent hire will be 20 years older than the incident hire in any given year of follow-up. If both subjects eventually die of lung cancer at age 65, only the prevalent hire will be identified as a case during follow-up (1950–94). As age is a risk factor for lung cancer, incident hires who start work older and die from lung cancer are more likely to be identified as cases than incident hires who start work younger.30 To address potential confounding by age at hire, we adjusted for it in the models.

In this analysis, we could not control for smoking. The absence of smoking information is a common challenge for occupational studies based on historical data. It has been demonstrated that the magnitude of uncontrolled confounding due to smoking is modest in occupational studies of lung cancer.31^–34 Most of this attention has focused on smoking as an unmeasured positive confounder, including discussion in the original analysis of this cohort.4 If, however, smoking cases are more likely to leave work than non-smoking cases, smoking in fact could be operating as a negative confounder in this study.

The impact of healthy worker bias in studies of respiratory cancer has been questioned as symptoms do not typically appear until shortly before diagnosis.35 Some have argued that the bias does exist and provided empirical evidence.6 24 36 When the exposure suspected of causing lung cancer is known to cause other adverse respiratory effects, which are, in turn, associated with lung cancer, then HWSE becomes even more plausible. In the case of silica, the association with lung cancer has been found to be stronger among silicotics, suggesting that silicosis may be on the causal pathway to lung cancer.37 In addition, silica causes deficits in pulmonary function and x ray changes consistent with silicosis, which can lead to difficulty breathing, persistent cough, and bronchitis.38 If these earlier health effects are related to silicosis, then they can bias associations with lung cancer because leaving work is related to both lower exposure (truncated exposures) and outcome.

We conclude that HWSE bias may be reduced in occupational cohort studies by restricting date of hire to be close to the start of follow-up. Further, we urge that authors of occupational epidemiological studies clearly state what the dates of hire and follow-up were for study subjects to ensure that this bias is limited in any analysis.