Occupational silica exposure and mortality from lung cancer and nonmalignant respiratory disease: G-estimation of structural nested accelerated failure time models.

BACKGROUND: Occupational exposure to crystalline silica is known to increase risks of both lung cancer and noninfectious nonmalignant respiratory diseases (NMRD). However, associations between silica exposure and survival times have not been described.
METHODS: In a longitudinal cohort of diatomaceous earth workers exposed to crystalline silica (primarily cristobalite) and followed from 1942 to 2011, we applied g-estimation of structural nested accelerated failure time models to adjust for time-varying confounding that could result in healthy worker survivor bias. A continuous measure of exposure was used in analyses estimating the hypothetical effect of banning exposure to silica on survival time. Since a ban is infeasible, sensitivity analyses examined the hypothetical effects of enforcing various Occupational Exposure Limits.
RESULTS: The estimated median number of years of life lost per worker (for all natural causes) due to silica exposure was 0.48 (95% confidence interval = 0.02, 1.01). For NMRD deaths, the corresponding estimate was 3.22 (0.82, 7.75) and for lung cancer deaths, 2.21 (0.97, 3.56). Cause-specific estimates were sensitive to the use of weights to adjust for competing events. Lung cancer mortality, which tended to occur at younger ages, was an important competing event for NMRD mortality. Sensitivity analyses supported the main results, but with larger estimates, and suggested that a strict limit would be nearly as effective as a complete ban on silica exposure.
CONCLUSIONS: Workplace exposure to crystalline silica in this industry appears to shorten survival times significantly, particularly for those who die of lung cancer or NMRD. More stringent exposure limits are probably warranted.

What this study adds

An alarming recent increase in severity of pneumoconioses in midlife Appalachian miners has been attributed partly to increased crystalline silica exposure. The Mine Safety and Health Administration acknowledges the need for an updated standard. Although silica is known to cause lung cancer and nonmalignant respiratory disease, its effect on survival time is unknown. Furthermore, most analyses lack rigorous adjustment for healthy worker survivor bias. This study applied g-estimation of accelerated failure time models to reduce healthy worker survivor bias and found that eliminating/reducing exposure would have lengthened life substantially for diatomaceous earth workers in California who died of respiratory causes.

Inhaled crystalline silica in the workplace is a recognized lung carcinogen[1-4] that causes many other adverse outcomes. The earliest identified was silicosis, a progressive fibrotic form of pneumoconiosis.[3,5,6] Occupational exposure to crystalline silica also reduces lung function[7-9] and increases the risk of other nonmalignant respiratory diseases (NMRD), including chronic obstructive pulmonary disease.[4,10,11] A recent alarming increase in pneumoconiosis prevalence and severity among Appalachian coal miners has been attributed at least partially to silica exposure.[12,13]

Quantifying the relation of silica exposure to outcomes other than silicosis is challenging. When occupational exposures adversely affect health, symptomatic workers may reduce their exposure by transferring to jobs with lower exposure or leaving employment entirely. The healthiest workers then accumulate the most exposure over time, giving rise to a paradoxical phenomenon called the healthy worker survivor effect (HWSE): higher cumulative exposure appears to be associated with longer survival and lower risk of the outcome.[14-16]

Traditional regressions yield biased results when prior exposure affects time-varying confounders, whether or not the confounders are included in the model.[17] The HWSE phenomenon is one example: health-related variables, even if unmeasured, may share causes with or directly affect the outcome and can also affect employment status, which in turn determines whether or not the worker continues to be exposed. Time-varying confounding by health and by employment status/duration is thus present, and it is hypothesized to be affected by prior exposure (Figure 1). Therefore, correct adjustment for it requires special methods,[17] e.g., g-estimation of a structural nested accelerated failure time model.[18,19] The recent occupational epidemiology literature includes several applications of this method.[20-26]

Figure 1.

Causal diagram showing exposure (crystalline silica) and outcome (death), as well as time-varying confounders: health (usually an unmeasured variable), employment status. To simplify the diagram, variables preceding the first exposure variables are not shown. Health after first exposure affects both the outcome and later exposures (by way of employment status). This means that health status is a confounder, and we need to adjust for it (or, since it is unmeasured, for its proxy employment status, which will block the backdoor path silica ← employment ← health → death). However, health and employment status are also affected by earlier exposures. If we analyze cumulative exposure to silica, adjusting for employment status or duration, we are adjusting for a variable that occurs after some of the exposure, which causes collider bias due to the pathway from earlier silica → employment ← health → death.

The link between cumulative exposure to crystalline silica (including both quartz and cristobalite) and lung cancer was confirmed in 2001 in a pooled epidemiologic analysis of data from 10 studies,[2] including a cohort of workers in the diatomaceous earth industry.[27-29] Diatomaceous earth is the fossilized remains of certain algae whose cell walls are mostly amorphous silica; it also contains a small proportion of quartz. After extraction from a quarry mine, processing includes calcining at high temperatures, which increases the proportion of cristobalite.[27] Processed diatomaceous earth is used as a filtration material for liquids and as an insecticide.

Workers in the diatomaceous earth cohort were exposed to crystalline silica, principally cristobalite, and followed for mortality from various causes. In traditional analyses, risks of both lung cancer and NMRDs increased with higher cumulative exposure.[4,27-29] However, the strength of these associations may have been underestimated due to bias from HWSE. If so, regulations based on those findings may not be stringent enough to protect worker health. The Occupational Safety and Health Administration (OSHA) recently promulgated a new standard for silica with a permissible exposure limit of 0.05 mg/m3,[30] but mines are regulated under the Mine Safety and Health Administration (MSHA) rather than OSHA. MSHA issued a statement in 2010 acknowledging that the industry standard, based on an assessment from 1973, is outdated.[31] Thus, more recent analyses may help in establishing a safer standard.

We applied g-estimation of structural nested accelerated failure time models to eliminate bias from HWSE in an examination of exposure to crystalline silica and survival time in this cohort of diatomaceous earth workers, focusing on mortality from natural causes, lung cancer, and NMRD.

Methods

Cohort description

The cohort is described in detail elsewhere.[27] Briefly, it includes 2342 white male workers from two diatomaceous earth mining and processing plants in Lompoc, California. Workers entered the cohort after employment for at least 1 year at either plant, including at least 1 day between January 1, 1942, and December 31, 1987. Work histories and silica exposure assessments were available from the beginning of plant operations (1902 and 1946 for the two plants) through 1994, with mortality follow-up from January 1, 1942, to December 31, 2011, based on National Death Index data, state driver’s license bureaus, and commercial credit bureaus. Workers lost to follow-up (N = 183) were censored the day after their last observed date of employment.

Industrial air monitoring measurements between 1962 and 1988 were used to estimate quantitative dust exposure, while data archived by the company provided information for the period 1948–1962.[32] Job-specific respirable dust exposure estimates were generated based on available measurements, and exposures before 1948 were extrapolated, accounting for changes over time.[29] Estimates for exposure to respirable crystalline silica were derived from the percent of silica contained in a given product and the exposure time to that product for each job.[29,32] Job-specific exposure intensities (mg/m3) were used to create the time-weighted average daily exposure intensity for each worker in each year, which was then lagged by 17 years to account for both the latency period for cancer and the unavailability of job history data after 1994. Prior work in this cohort suggests that a 10-year lag might be better for NMRD mortality,[33] but this would have reduced power, and we used the same lag to prioritize consistency between the two analyses.

Because two small operations in the plants involved chrysotile asbestos, asbestos exposures were derived from monitoring data and records of quantities of asbestos in mixed products from 1930 onwards. Exposures for earlier years were extrapolated.[29] Demographic information included hire year, duration of employment at study sites, dates of specific jobs held, and Latino ethnicity. Smoking status (ever/never, collected by the industry’s medical surveillance program starting in the 1960s[27]) was available for 50% of the cohort (N = 1171).

Separate analyses were conducted for mortality from all natural causes, lung cancer, and NMRD excluding pneumonia, influenza, and other infectious diseases.

Statistical methods

Since both lung cancer and NMRD mortality analyses can be subject to HWSE,[20] we conducted path analysis to determine whether HWSE could cause bias in this cohort: was employment status associated with future exposure and outcomes, and also affected by prior exposure? Since workers could not be exposed to silica in the Lompoc plants after leaving employment, employment status affected future exposure. We therefore only checked whether employment was associated with the outcomes, independent of cumulative exposure, and whether exposure to silica affected time to employment termination. For the former, we used proportional hazards models, adjusting for covariates and total cumulative exposure, to assess whether employment duration was independently associated with mortality. For the latter, we used g-estimation of an accelerated failure time model (see below) to evaluate whether previous exposure might have caused workers to terminate employment earlier than they otherwise would have. After confirming that employment was a time-varying confounder affected by prior exposure, we proceeded to assess the relation between exposure and the outcomes of interest using g-estimation to adjust correctly for this confounding.

For each outcome (mortality from each of the causes of interest, and employment termination for the path analysis described above), survival time is defined as the time from cohort entry to the date of that outcome. We modeled counterfactual survival times that would have been observed under no exposure, starting from entry into the cohort, as a function of observed exposures and observed survival times, assuming a log-linear exposure-response in a structural nested accelerated failure time model:

In this equation, is observed survival time, denotes obser ved exposure (mg/m3) in year , and represents the unknown constant parameter of the model; when we use its true value, represents the counterfactual survival time under no exposure (see eAppendix; http://links.lww.com/EE/A22 for more intuition). We estimated using g-estimation, which was developed by Robins[18] and is explained in pedagogical detail by Hernán et al.[19] This parameter can be interpreted as the negative log of the ratio of median survival times, comparing what would have happened if everyone had been exposed at an intensity of one unit every year of follow-up (regardless of employment status) to what would have happened if everyone had always remained unexposed. We subtracted this ratio from 1 to obtain the relative difference in median survival times, but these measures of etiologic effect are quantified with reference to an implausible scenario, as workers are not exposed to silica after leaving employment. We therefore used our estimate to calculate the median number of years of life that would have been saved per worker (compared to the observed scenario) under a ban on exposure to crystalline silica in the Lompoc plants, which has a more practical interpretation.[26]

Combining exposure over time into cumulative exposure would conflate the possibly different effects of intensity and duration of exposure.[34] This model takes into account a worker’s entire quantitative history of exposure to silica (i.e., at each time ) without summarizing it in a single number (cumulative exposure).

G-estimation leverages the assumption that there are no unmeasured confounders of the exposure–outcome relation or, equivalently, that counterfactual outcomes are statistically independent of observed exposures, conditional on measured covariates. This assumption implies that if we knew the counterfactual survival times, we could include them (or any function of them) along with all measured covariates in a model predicting observed exposures, and the coefficient of the counterfactual survival times would be 0. Our g-estimation procedure entailed choosing “candidate” values (within a reasonable search interval) for the unknown parameter in the structural model and using them to calculate corresponding values for the counterfactual unexposed survival times from the structural model, observed exposures, and observed outcome times. Then, for each value of , we tested whether the candidate counterfactual survival times were conditionally independent of observed exposures. Our estimate is the value of that makes have coefficient equal to zero in a traditional multivariate regression predicting exposure . We predicted the level of exposure in categories defined by quartiles of the observed distribution, using ordinal logistic regression adjusted for the measured confounders: Latino ethnicity, smoking (ever/never/missing), age, calendar year, time taken off work in the previous year, exposure to silica in the previous year, previous cumulative exposures to dust and asbestos, and employment duration prior to start of follow-up.[19]

Not all workers died during follow-up; the analysis adjusted for administrative censoring by end of follow-up (see eAppendix; http://links.lww.com/EE/A22).[18] We considered two additional types of censoring. All analyses were adjusted for censoring by loss to follow-up by applying inverse probability of censoring weights, in order to assess what would have been observed in the absence of loss to follow-up. Analyses were conducted with and without further weights to adjust for censoring by competing risks. For all natural mortality, the only potential competing risks are external causes of death; for lung cancer mortality, we considered death from NMRD a competing risk, and vice versa, because these two outcomes share many causes and are known to be related to silica exposure. The weights were equal to the inverse of the probability of remaining uncensored (i.e., not being lost to follow-up, and not dying from specific competing risks: the product of the estimated probabilities) from that time forward.[19] These probabilities were predicted on the basis of Latino ethnicity, age, calendar year, smoking, cumulative prior exposures to silica, dust, and asbestos, and employment duration prior to cohort entry.

As sensitivity analyses, we repeated each analysis, replacing the quantitative exposure metric with each of a series of binary exposure metrics defined by silica concentrations above versus below cutoffs[25] corresponding to the median observed annual average daily exposure (0.117 mg/m3), the OSHA standard (0.05 mg/m3),[30] the American Conference of Governmental and Industrial Hygienists’ recommended limit (0.025 mg/m3),[35] and a ban (0 mg/m3).

Analyses were carried out in SAS 9.4 (SAS Institute, Inc., Cary, NC). The Committee for the Protection of Human Subjects at the University of California, Berkeley, approved the study.

RESULTS

Table 1 contains demographic information about the all-male, white cohort. Smoking data were available for half the cohort, but this half contained a lower proportion of Latinos. Among workers for whom smoking was recorded, Latinos were less likely to be ever-smokers than non-Latinos were. About 10% of the cohort was hired before follow-up began; 40% of workers survived to the administrative end of follow-up.

Table 1

Demographic and workplace characteristics of California diatomaceous earth industry cohort of white men (N = 2342)

Table 2 summarizes the exposures assessed during active employment. About 73% of the person-years were after employment termination. Cumulative exposure to silica (lagged by 17 years) was under 5 mg/m3 years for 90% of the study population.

Table 2

Exposures during active employment in the diatomaceous earth industry in Lompoc, California (person-years, except for cumulative exposures, which are the final cumulative exposures for each worker), with 17-year lag

Path analysis found significant associations between employment duration and natural-cause mortality, lung cancer, and NMRD, confirming that it is a confounder. Furthermore, exposure shortened time to leaving work: workers would have terminated employment a median of 0.47 (95% confidence interval = 0.11, 0.93) years later if they had never been exposed. Thus, time-varying confounding affected by prior exposure was present, necessitating the use of g-methods to avoid bias from HWSE.

For mortality from all natural causes, if everyone had been exposed to crystalline silica at an average daily intensity of 0.1 mg/m3 (approximately the 32nd percentile of nonzero exposures), every year from entering the cohort to the end of follow-up, survival time from cohort entry would have been at least 3.0% (95% confidence interval = 0.1, 5.7) shorter for half the workers than it would have been under no exposure. The corresponding estimate for lung cancer from analysis without weights for competing risks (9.8% [5.3, 14.1]) was stronger than that for all natural causes. However, no estimate was found when adjusting for competing risks from NMRD deaths (the estimating function did not cross 0). Estimates for NMRD were sensitive to the use of weights to adjust for censoring by competing deaths from lung cancer. With adjustment, NMRD mortality was also strongly associated with silica exposure: half the workers would have died an estimated 10% (3.2, 16.3) sooner if always exposed at 0.1 mg/m3 than they would have under no exposure.

Figure 2 presents estimates of the median number of years of life that would have been saved per worker who died of various causes if they had never been exposed to silica, compared to their observed survival times under the exposures they actually experienced.

Figure 2.

Median number of years of life that could have been saved per worker if exposure to crystalline silica had been eliminated starting in 1925, among workers who died of various causes during follow-up. All analyses adjusted for the following confounders: Latino ethnicity, age, calendar year, smoking (ever/never/missing), time taken off work in the previous year, exposure to silica in the previous year, previous cumulative exposures to dust and asbestos, and employment duration prior to start of follow-up. Diamonds represent estimates from analyses that were adjusted for censoring by loss to follow-up. The hollow circle represents the estimate from an analysis additionally adjusted for censoring by deaths from external causes. The solid circle represents the estimate from an analysis additionally adjusted for censoring by deaths from lung cancer. No estimate was found in the lung cancer analysis additionally adjusted for censoring by deaths from nonmalignant respiratory disease (NMRD).

Results from sensitivity analyses with binary exposure metrics agreed qualitatively with the main analyses using a quantitative exposure variable: lowering exposures would have increased survival times among workers who died of the causes of interest. Estimates were stronger than those seen in the main analyses (Table 3).

Table 3

Results from sensitivity analyses with binary exposure metrics: median number of years of life that would have been saved per worker who died of the specified causes, under various exposure limits

DISCUSSION

Previous analyses applying traditional methods (including both regressions and standardized mortality ratios comparing rates in workers to those in the general population) have found reasonably consistent associations between exposure to silica and lung function,[7-9] as well as rates or risks of both lung cancer[2-4] and NMRD.[4-6,10,11] With survival time as the metric, our analysis does not allow direct quantitative inference about the effect of silica exposure on cause-specific mortality risks because the model assumes that those who died of NMRD or lung cancer would have died of those causes even if they had not been exposed to silica. (This assumption is false for silicosis as a single outcome; our composite NMRD outcome includes COPD and pneumoconiosis mortality.) Among workers for whom smoking was measured, only two lung cancer deaths and five NMRD deaths occurred in nonsmokers. Most of the cases would have therefore had elevated risks for these outcomes even without exposure to silica, making the assumption reasonably plausible. Controlling bias from the HWSE and assuming these workers would have died of the same causes even if they had not been exposed, we found that silica exposure was associated with earlier death—a new finding distinct from previously reported increased risks. However, these data cannot distinguish between earlier onset and more aggressive forms of the disease. Exposure to silica had a stronger, though less precise, relationship with NMRD than with lung cancer. Sensitivity analyses using binary exposure metrics suggested that setting an exposure limit at 0.025 mg/m3 would have saved nearly as many years of life overall as a complete ban in this cohort, though this was less true for lung cancer.

We presented estimates of relative differences of median survival times comparing two hypothetical interventions (exposing everyone to crystalline silica at an average daily intensity of 0.1 mg/m3 in every year versus never exposing anyone). Since workers are never exposed after they leave work, this estimate is not what would be observed under any feasible intervention. Nevertheless, it is informative in its own right as an estimate of the etiologic effect of a static intervention,[36,37] and it resembles the type of measurement used by OSHA in setting standards (i.e., estimated mortality following continuous exposure at a given level for 40 years). We also presented estimates of the median number of years of life that were lost due to exposure to silica, i.e., estimates of the effect of an intervention that, while perhaps still not feasible, is meaningful from the regulatory point of view. These estimates compare what would have happened if silica exposure had been banned to what actually happened. This effect measure depends not only on the magnitude of the etiologic effect but also on the distribution of exposure among those who died of the cause of interest. Thus, even if the etiologic effect of silica exposure (the true value of ) were the same in another population of similar workers, the number of years of life lost due to silica exposure would be different if the population had a different exposure distribution.

Estimates from the analysis of mortality from natural causes were not sensitive to adjustment for censoring by deaths from external causes. This is consistent with the plausible assumption that mortality from external causes either does not share risk factors with mortality from natural causes or is not associated with exposure. An advantage of analyzing mortality from all natural causes is that, regardless of exposure, everyone will eventually experience the outcome (unless they first die from external causes), so the assumption that the accelerated failure time model requires is met. We found that workers would have lived about a half-year longer (median) if they had never been exposed to silica. This finding is presumably driven mostly by diseases already linked to silica exposure; indeed, the estimates for lung cancer and NMRD were stronger.

Lung cancer and NMRD share their most important risk factors, so we expected each to act as an informative censoring event in the analysis of the other. Each of the two cause-specific analyses was sensitive to adjustment for competing deaths from the other cause using inverse probability of censoring weights. The unweighted estimate for NMRD was essentially null, even though the weighted estimate was strong. Combined, these results are consistent with lung cancer death acting as an informative censoring event in the analysis of NMRD: without adjustment for competing risks, the association between exposure to silica and NMRD mortality is difficult to detect.

By contrast, the analysis of lung cancer did not yield an estimate when weighting to adjust for censoring by NMRD death, though the unweighted estimate was convincing. The lung cancer deaths generally occurred at younger ages than deaths from NMRD, perhaps indicating that the people at risk of lung cancer did not live long enough for NMRD to progress to death. In that case, NMRD mortality did not act as a competing event for lung cancer mortality, whereas lung cancer deaths could censor deaths from NMRD, which were more likely to occur later in life.

Estimates from analyses without inverse probability of censoring weights for competing risks are easier to interpret; weighting the analysis is equivalent to estimating the effect in a pseudo-population in which death from the competing event cannot occur. Although preventing all deaths from lung cancer is unrealistic, the NMRD mortality analysis in a pseudo-population in which lung cancer mortality does not occur allows us to estimate the direct effect of exposure when the pathway of disease prevention via premature death is not included. Our results suggest that the impact of censoring by lung cancer death should be considered in analyses of NMRD mortality.

Other ways to deal with competing risks and HWSE without envisioning interventions to remove competing risks include analyzing the data by applying either the parametric g-formula[38,39] or targeted maximum likelihood estimation[40,41] to adjust for HWSE, while estimating contrasts based on subdistribution functions of risk.[42] Although these methods can compare dynamic interventions on exposures occurring while at work, they require additional assumptions and interventions on employment status to estimate an etiologic effect. Furthermore, unlike the model presented here, these methods do not model survival time directly. An application of the parametric g-formula to estimate subdistribution functions of risks in this cohort found qualitatively similar results, with silica linked more strongly to risk of NMRD mortality than lung cancer.[43] Our substantive conclusions (linking silica exposure to shorter survival times from both causes, with the same relative strengths observed for risks in Neophytou et al.[43]) are thus robust to different modeling assumptions and effect measures, lending further confidence to both results.

Results can be interpreted causally only if certain assumptions hold: consistency (see eAppendix; http://links.lww.com/EE/A22), correct model specification, and conditional exchangeability.

Our results were somewhat sensitive to the exposure metric used in the specification of the structural model. Sensitivity analyses with binary exposure metrics yielded larger estimates; the reasons for this are unclear. One possible explanation is that the main results are underestimated due to violations of parametric assumptions (if, e.g., the relation with survival time is not log-linear). Another possibility is that the relevant exposure is misclassified by the binary cutoffs. See eAppendix; http://links.lww.com/EE/A22 for more details on parametric modeling assumptions and the sensitivity analyses.

Conditional exchangeability (see eAppendix; http://links.lww.com/EE/A22) is the explicit basis for g-estimation.[19] If an important confounder was unmeasured or excluded in the exposure model, then results would be biased. In this study, only crude smoking information was available and only for half the cohort, so there may be residual confounding by smoking or bias due to using a “missing” category in the analysis.[44] Exposures tended to be higher among those for whom the information was missing, but this was mainly because these were employees who left employment earlier during follow-up, when exposure levels were higher. If workers diagnosed with silicosis were more likely to quit smoking, then smoking status or intensity at time might be affected by prior exposure. Detailed smoking data were unavailable,[27] but ever-smoking is not a strong predictor of exposure. Among workers whose smoking data were available, exposure levels while employed were similar for ever-smokers and nonsmokers. Results from other analyses in this cohort did not change much when multiple imputation was used to fill in missing smoking data, compared to using ever/never/missing categories.[43] Thus, as in many occupational studies, even though smoking is the most important risk factor for the outcomes under study here, its association with exposure is probably too weak to make it a strong confounder.[45-47]

One limitation of this analysis is that because of the 17-year lag, we consider an intervention on exposure beginning in 1925, but follow-up begins in 1942; workers are only included if some of their employment occurred after the beginning of 1942. We cannot know if some excluded workers, who left employment before the start of follow-up, would have remained employed long enough under the intervention considered to be included in the cohort. Our results only apply to the workers included, who may be less susceptible than the excluded workers to the health effects of exposure to silica. However, in that case our results would be biased downward; the true effect might be slightly stronger in a cohort that includes all the workers. Only about 10% of our cohort had been employed for more than a year prior to the start of follow-up, though these workers accounted for 16% of natural cause mortality and lung cancer mortality, and 24% of NMRD mortality.

While it has been known for some time that workplace silica exposure increases the risks of diseases that shorten life, our results quantify the number of years of life lost due to exposure to crystalline silica. Applying a method that controls bias due to the HWSE, we estimated that survival times would have been significantly longer under a hypothetical intervention banning exposure and that limiting exposures to 0.025 mg/m3 might be nearly as effective for these workers. Given that MSHA has recognized the need to re-evaluate industry standards for crystalline silica,[31] this analysis and additional research should help guide policy.

Conflict of interest statement

The authors declare that they have no conflicts of interest with regard to the content of this report.