Estimating Counterfactual Risk Under Hypothetical Interventions in the Presence of Competing Events: Crystalline Silica Exposure and Mortality From 2 Causes of Death.

Exposure to silica has been linked to excess risk of lung cancer and nonmalignant respiratory disease mortality. In this study we estimated risk for both these outcomes in relation to occupational silica exposure as well as the reduction in risk that would result from hypothetical interventions on exposure in a cohort of exposed workers. Analyses were carried out using data from an all-male study population consisting of 2,342 California diatomaceous earth workers regularly exposed to crystalline silica and followed between 1942 and 2011. We estimated subdistribution risk for each event under the natural course and interventions of interest using the parametric g-formula to adjust for healthy-worker survivor bias. The risk ratio for lung cancer mortality, comparing an intervention in which a theoretical maximum exposure limit was set at 0.05 mg/m3 (the current US regulatory limit) with the observed exposure concentrations, was 0.86 (95% confidence interval: 0.63, 1.22). The corresponding risk ratio for nonmalignant respiratory disease mortality was 0.69 (95% confidence interval: 0.52, 0.93). Our findings suggest that risks from both outcomes would have been considerably lower if historical silica exposures in this cohort had not exceeded current regulatory limits.

Exposure to crystalline silica has been linked to excess mortality risk, primarily from nonmalignant respiratory disease but also more recently from lung cancer (1–7). The International Agency for Research on Cancer classified silica as a confirmed (Category 1) human carcinogen in 1996 (8). Silica exposure occurs in multiple industries, including mining, shipbuilding, construction, and manufacturing, with millions of workers believed to be exposed in the United States and worldwide (9, 10). In the United States, the Occupational Health and Safety Administration was established in 1971, and the first permissible exposure limits (PELs) for silica exposures were adopted in the same year. A final rule was issued in 2016, designed to reduce lung cancer, silicosis, other nonmalignant respiratory disease and kidney disease in workers and setting a PEL for respirable crystalline silica of 0.05 mg/m3 averaged over an 8-hour shift (11).

The Mine Safety and Health Administration (MSHA) is the agency responsible for promulgating and enforcing safety and health standards for the mining industry in the United States. It requires metal and nonmetal mines to comply with exposure limits for airborne contaminants not otherwise regulated, based on the threshold limit values proposed by the American College of Government and Industrial Hygienists in 1973 (with a somewhat different policy for coal mines). The threshold limit values proposed by the American College of Government and Industrial Hygienists include a value for silica “designed to limit silica exposures to 0.10 mg/m3” and effectively half that value (0.05 mg/m3) for crystalline silica in the form of cristobalite and tridymite (12). MSHA has since issued a statement acknowledging that the standards are outdated and expressing the need to for a new rule to reduce miners’ exposure to respirable crystalline silica based on the Occupational Health and Safety Administration’s updated assessment (13). Although these limits are not enforceable as regulatory standards, the National Institute for Occupational Safety and Health first recommended an exposure limit of 0.05 mg/m3 averaged over an 8-hour shift in 1974 (14), and the American College of Government and Industrial Hygienists now recommends an even lower exposure limit of 0.025 mg/m3 (15).

Regulatory standards for workplace hazards that aim to reduce risk associated with exposures often depend on evidence from epidemiologic studies. Estimating absolute risk over time from observational studies can be challenging given that occupational studies may suffer from healthy-worker survivor bias taking the form of time-varying confounding affected by prior exposure (16), which cannot be addressed using traditional regression approaches and is generally expected to result in downward bias (17). This motivates the use of methods equipped to address this type of bias, collectively known as “g-methods” (18), such as the g-formula (19).

Another concern is the issue of competing risks, which can complicate the estimation of risk of mortality from one cause of interest when an exposure is associated with multiple causes of death, as is the case for crystalline silica. In the absence of independence of the competing events, failure to account for the change in risk of a competing event may lead to biased estimates of risk for the event of interest (20). The parametric g-formula can be used to generate effect estimates based on subdistribution functions of risk, which are often used to deal with competing events (20, 21).

In the present analysis, we used this method to assess the cumulative risk of lung cancer and nonmalignant respiratory disease mortality, both strongly associated with silica exposures, under hypothetical interventions on crystalline silica exposures, including interventions in compliance with MSHA regulations in a cohort study of California diatomaceous earth workers (3). We assessed counterfactual risk under hypothetical interventions while addressing time-varying confounding affected by prior exposure, and we estimated subdistribution functions of risk accounting for changes in counterfactual risk of competing causes of death.

METHODS

Study population

Analyses were carried out using data from a cohort of diatomaceous earth workers in California; the cohort is described in greater detail elsewhere (3). Briefly, the study population consisted of 2,342 male workers from 2 diatomaceous earth plants in Lompoc, California. Inclusion criteria were cumulative employment for at least 1 year at either plant and having worked for 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 for one plant and 1946 for the other) through 1994. Mortality follow-up was based on data from the National Death Index, state driver’s license bureaus, and commercial credit bureaus (4) and was available from January 1, 1942, to December 31, 2011, for a maximum follow-up time of 70 years. Complete mortality follow-up was not available for 183 participants. These subjects were considered alive until their last observed date of employment and censored afterward. Demographic information on the cohort included work history (hire year, duration of employment at study sites, and dates of specific jobs held) and ethnicity. Information on smoking status (ever/never) was also available for 50% of the cohort (n = 1,171).

Exposure assessment

Quantitative dust exposure estimates were determined primarily from industrial air monitoring measurements made between 1962 and 1988, with company-archived data providing some additional information for the period of 1948–1962 (22). Job-specific respirable crystalline silica (mostly in the form of cristobalite in this cohort) and respirable dust exposure estimates were generated based on available measurements. Exposures before 1948 were based on extrapolated job-specific exposures that accounted for interventions to reduce dust exposures and other changes over time (1). The estimates for crystalline silica were derived from the percent of silica contained in a given diatomaceous earth product and the exposure time to that product for a given job. Detailed work history was available through 1994, by which time 88% of the cohort had terminated employment at the participating plants. Job-specific exposure concentrations were combined with job duration and summed for cumulative exposures to silica and dust (milligrams per cubic meter–years) (1, 22).

Asbestos exposures were also derived, because 2 small operations involving chrysotile asbestos occurred over the study period. Estimates were based on monitoring data and records of quantities of asbestos in mixed products from 1930 onward, while data was extrapolated to determine exposures for earlier years (1, 23).

Statistical analyses

We applied the parametric g-formula to assess cumulative risk of lung cancer and nonmalignant respiratory disease mortality under hypothetical interventions on crystalline silica exposures. The process for the parametric g-formula and all models considered is described in greater detail in Web Appendix 1 (available at https://academic.oup.com/aje), along with a directed acyclic graph depicting the hypothesized relationships in this study (Web Figure 1). Briefly, we fitted parametric models for each of the following: the outcome (lung cancer and nonmalignant respiratory disease mortality in separate analyses), competing events (mortality from all causes other than the event of interest in each of the analyses), loss to follow-up, exposure (silica) and other time-varying covariates (total respirable dust, asbestos (all lagged by 15 years)), and employment status. Each model was fitted conditional on prior exposure and covariate histories including baseline covariates: age fitted as a cubic spline, calendar year as a multicategory variable with 1 level per decade, an indicator variable for Hispanic ethnicity, smoking status (ever/never/missing), and, as continuous terms, cumulative exposures accrued prior to beginning of follow-up. Silica and respirable dust exposures were predicted in a 2-step process: first, logistic models were used to predict a binary indicator for whether exposures were greater than 0; for the instances where exposures were predicted to be greater than 0, actual exposure values were predicted based on linear models for the log-transformed exposures. Approximately 15% of actively employed person-time was exposed to asbestos, with a wide range of exposures among the exposed, resulting in unstable prediction models for quantitative exposures to asbestos, particularly in some of the bootstrap samples. Thus, our estimation process relied on binary indicators of any asbestos exposure in each year, predicted using logistic regression models. In models for the outcome and censoring events, exposures to silica and respirable dust were fitted as separate cubic spline terms for exposure in the prior year and cumulative exposure up to the year before.

We then generated a large Monte Carlo pseudosample (n = 50,000), based on the observed distributions of the baseline covariates. In this pseudosample we simulated exposure and covariate values at each age (up to age 90 years) using the parameters of the models for the exposure and covariates. Under no intervention, the simulation uses values for the exposure and time-varying covariates that are predicted from the model and then predicts the risk under the observed (natural) course of events. For all other interventions, the exposure values were changed from the predicted values according to the specific intervention. We considered 2 interventions, one setting a hypothetical maximum exposure limit on average daily crystalline silica exposures equivalent to the current MSHA effective PEL for cristobalite of 0.05 mg/m3 at each time for the duration of follow-up, and one setting crystalline silica exposures to zero. (The more conservative PEL corresponding to cristobalite was chosen, because the majority of crystalline silica exposure in this cohort was in the form of cristobalite; we did not have information on proportion of silica form (i.e., cristobalite vs quartz) on the level of individually assigned exposures.) For the intervention setting a hypothetical exposure limit of 0.05 mg/m3, all predicted silica exposure values above 0.05 mg/m3 were replaced with 0.05 mg/m3 and otherwise remained unchanged. We compared both interventions with no intervention (i.e., the observed natural course (what actually happened)).

Risk, with age as the time scale (defined as the probability of death by the event of interest by age t (24)), was calculated for each intervention using an estimator for the subdistribution of the event of interest, in the presence of competing risks (20, 25). Subdistribution cumulative incidence functions can be seen as parts of a composite outcome that is the complement of survival. In the absence of competing events (as in the case of all cause-mortality), risk can be seen as a complement function of survival , so that and reaches 1 as follow-up time approaches infinity . In the presence of competing risks (as with multiple causes of death) we assume that the complement of survival is now a composite outcome (e.g., and and ); here we estimate for 2 mortality causes of interest, in separate analyses, while accounting for mortality from all causes other than each cause of interest . Thus, using these functions we can estimate absolute population-average risk of a particular cause of death as the analogue of the proportion of cases with the cause of interest in a hypothetical closed cohort with no loss to follow-up, where death from multiple causes is possible. We repeated the entire above process in 200 bootstrap samples. The standard deviation of the estimates from the bootstrap samples was used as an estimate of the standard error for the estimates of risk differences and the log of the risk ratios (26), which were then used to generate 95% confidence intervals. For additional details on the general application of the parametric g-formula, refer to Taubman et al. (27) or Cole et al. (28).

We also performed a sensitivity analysis with multiple imputation for missing smoking status using the “MI” procedure in SAS (SAS Institute, Inc., Cary, North Carolina) and application of the parametric g-formula in each of the imputed data sets. Values for smoking were imputed in 50 data sets using logistic regression, with the aforementioned baseline variables, cumulative exposure to silica, respirable dust and asbestos, death by the event of interest or competing event, and loss to follow-up, as covariates included in the imputation. Other sensitivity analyses included exclusion of deaths from infectious lung diseases (pneumonia, influenza, and acute respiratory infections) from the definition of nonmalignant respiratory disease as an outcome of interest. All analyses were carried out in SAS, version 9.4 (SAS Institute, Inc.).

RESULTS

Demographic and baseline characteristics of the study population are presented in Table 1. The observed numbers of death from the causes of interest were 113 for lung cancer and 165 for nonmalignant respiratory disease (127 excluding all infectious lung diseases). Mean age at lung cancer death was 66.6 (standard deviation, 10.2) years and at nonmalignant respiratory death, 70.7 (11.8) years. Mean annual average daily exposure to crystalline silica during actively employed person-time in this cohort was 0.16 (standard deviation, 0.21) mg/m3 while the mean cumulative exposure at end of follow-up was 2.16 (standard deviation, 3.51) mg/m3-years. Figure 1 summarizes the distribution of annual average daily exposures for this study, with indicators for the median and the current PEL, which corresponded to the 38th percentile of the observed average daily exposures during actively employed person-time.

Table 1.

Characteristics of a Cohort of 2,342 California Diatomaceous Earth Workers Exposed to Crystalline Silica and Followed for Mortality Between 1942 and 2011

Characteristic	No. of Participants	%	Median (Range)	Mean (SD)
Hispanic	546	23.3
Ever-smokersa	861	73.5
Age at beginning of follow-up, years			27 (17–61)
Year of hire			1952 (1908–1986)
Year of birth			1927 (1881–1966)
Duration of employment, years			5 (1–50)
Duration of follow-up, years			39 (1–70)
Total deaths	1,219	52.0
Lung cancer deaths	113	4.8
Nonmalignant respiratory deaths	165	7.0
Cumulative silica exposure, mg/m3-years				2.16 (3.51)
Cumulative asbestos exposure, fibers/mL-years				1.44 (4.44)

Abbreviation: SD, standard deviation.

a Smoking data were available for 1,171 participants. Number and percentage of ever-smokers is based on this subset of participants.

Figure 1.

Distribution of annual average daily crystalline silica exposure (in mg/m3) during active employment in a cohort of 2,342 California diatomaceous earth workers followed between 1942 and 2011. The vertical lines represent the median value of 0.12 mg/m3 (gray line) and the value corresponding to the current Mine Safety and Health Administration permissible exposure limit of 0.05 mg/m3 (black line).

Observed risk curves for both events of interest, contrasted with the simulated natural course (no intervention), suggested that our models predict the observed data well (Web Figure 2). The median predicted annual average daily silica exposure during active person-time under the simulated natural course was 0.09 mg/m3, with an interquartile range, of 0.02–0.20, compared with the a median of 0.12 mg/m3 (interquartile range, 0.01–0.20) in the observed data. Table 2 presents estimates of risk of lung cancer and nonmalignant respiratory disease mortality under hypothetical interventions on crystalline silica exposures, along with risk ratio and risk difference compared with the natural course. The risk ratio for lung cancer mortality risk under an intervention setting a hypothetical exposure limit of 0.05 mg/m3 compared with the observed risk was 0.86 (95% confidence interval: 0.63, 1.22). The corresponding risk ratio for nonmalignant respiratory disease under the same intervention was 0.69 (95% confidence interval: 0.52, 0.93). Interventions setting exposure to zero resulted in risk ratios of 0.82 (95% confidence interval: 0.53, 1.26) for lung cancer and 0.63 (95% confidence interval: 0.43, 0.91) for nonmalignant respiratory disease mortality. Based on these estimates, the attributable fraction (estimated as the ratio of the risk difference between the risk under an intervention of always unexposed and the risk in the exposure (here the observed risk) over the risk in the exposed) of mortality due to the observed silica exposures in this population was estimated at 18% for lung cancer mortality and 37% for nonmalignant respiratory disease mortality. (The attributable fraction reported can also be estimated as (risk ratio −1)/risk ratio (29), where the risk ratio in the formula is the inverse of the risk ratio comparing the intervention of always unexposed with the natural course, for each event of interest.) Risk curves (over age) for both causes of death are depicted in Figure 2. Risk curves over age for the competing events for each cause of death for the natural course and under no exposure are depicted in Web Figure 3.

Table 2.

Lung Cancer and Nonmalignant Respiratory Disease Mortality at Age 90 Years Under the Natural Course and Hypothetical Interventions on Crystalline Silica in a Simulation Using Data From 2,342 California Diatomaceous Earth Workers Exposed to Crystalline Silica and Followed for Mortality Between 1942 and 2011

Intervention	Cumulative Risk, %	Rangea	RR	95% CI	RD	95% CI
Lung Cancer Mortality
Simulated natural course	7.2	5.1–15.4	1.00	Referent	0.0	Referent
Annual average silica ≤ 0.05 mg/m3	6.2	3.8–10.9	0.86	0.63, 1.22	−1.0	−3.4, 1.4
Annual average silica = 0 mg/m3	5.9	3.6–11.8	0.82	0.53, 1.26	−1.3	−4.0, 1.4
Nonmalignant Respiratory Disease Mortality
Simulated natural course	12.0	6.9–15.1	1.00	Referent	0.0	Referent
Annual average silica ≤ 0.05 mg/m3	8.3	5.0–13.1	0.69	0.52, 0.93	−3.7	−7.2, −0.2
Annual average silica = 0 mg/m3	7.5	4.9–13.8	0.63	0.43, 0.91	−4.5	−8.3, −0.7

Abbreviations: CI, confidence interval; RD, risk difference; RR, risk ratio.

a Range of risk estimates from 200 bootstrap samples.

Figure 2.

Risk of lung cancer (A) and nonmalignant respiratory disease (NMRD) mortality (B) in a cohort of 2,342 California diatomaceous earth workers followed between 1942 and 2011, under the natural course (solid lines) and under hypothetical interventions for annual average daily crystalline silica exposures of ≤0.05 mg/m3 or 0 mg/m3 (unexposed).

Results using a definition of nonmalignant respiratory disease mortality that did not include infectious disease were generally very similar to the wider definition of the outcome (Table 3). Table 4 summarizes estimates of risk, risk ratio, and risk difference for both outcomes of interest under the same interventions on exposure based on analyses from 50 multiply imputed data sets for missing smoking data as well as the range of estimates from all 50 data sets. The mean point estimates were similar to analysis using an ever/never/missing classification for smoking status and a relatively narrow range of estimates across all data sets.

Table 3.

Nonmalignant Respiratory Disease Mortality (Excluding Infectious Disease) at Age 90 Years Under the Natural Course and Hypothetical Interventions on Crystalline Silica in a Simulation Using Data From 2,342 California Diatomaceous Earth Workers Exposed to Crystalline Silica and Followed for Mortality Between 1942 and 2011

Intervention	Cumulative Risk	Rangea	RR	95% CI	RD	95% CI
Simulated natural course	9.4	5.1–16.8	1.00	Referent	0.0	Referent
Annual average silica ≤ 0.05 mg/m3	6.6	3.2–11.1	0.70	0.46, 1.02	−2.8	−6.4, 0.7
Annual average silica = 0 mg/m3	5.9	3.0–11.3	0.63	0.42, 0.96	−3.5	−6.5, −0.3

Abbreviations: CI, confidence interval; RD, risk difference; RR, risk ratio.

a Range of risk estimates from 200 bootstrap samples.

Table 4.

Averages and Ranges of Estimates From 50 Multiply Imputed Data Sets for Missing Smoking Status in a Simulation Using Data From 2,342 California Diatomaceous Earth Workers Exposed to Crystalline Silica and Followed for Mortality Between 1942 and 2011

Intervention	Risk	Minimum, Maximuma	RR	Minimum, Maximuma	RD	Minimum, Maximuma
Lung Cancer Mortality
Simulated natural course	7.5	7.0, 8.1	1.00	Referent	0.0	Referent
Annual average silica ≤ 0.05 mg/m3	6.3	5.9, 6.5	0.84	0.76, 0.87	−1.2	−2.0, −0.9
Annual average silica = 0 mg/m3	6.0	5.3, 6.3	0.79	0.71, 0.81	−1.5	−2.5, −1.3
Nonmalignant Respiratory Disease Mortality
Simulated natural course	11.8	11.5, 12.4	1.00	Referent	0.0	Referent
Annual average silica ≤ 0.05 mg/m3	8.2	8.0, 8.5	0.69	0.66, 0.73	−3.6	−4.1, −3.2
Annual average silica = 0 mg/m3	7.4	7.2, 7.8	0.63	0.59, 0.67	−4.4	−5.1, −3.8

Abbreviations: RD, risk difference; RR, risk ratio.

a Range of estimates from all 50 multiple imputation data sets.

DISCUSSION

We applied the parametric g-formula to estimate the counterfactual cumulative risk due to 2 different causes of death in a cohort of workers in the diatomaceous earth industry under hypothetical interventions on exposure to crystalline silica. We observed considerable risk reduction in mortality from both causes under a hypothetical intervention setting historical exposures to levels in compliance with PELs for crystalline silica currently applicable to nonmetal mines, although some excess risk appears to persist at these levels. Compliance with the current PELs would have resulted in 62% of the actively employed person-time in this cohort being intervened on to reduce exposure. The majority of excess risk due to the 2 causes of interest would have been prevented with this intervention. The remaining excess risk would still exceed the threshold of 1/1,000 deaths, which has been identified as “significant risk” by the US Supreme Court (30). In addition, the attributable fractions based on our estimates are considerable for both causes but especially for nonmalignant respiratory disease mortality.

Under assumptions of conditional exchangeability (no unmeasured confounding), counterfactual consistency (i.e., every individual’s counterfactual outcome under their observed exposure history is equal to their observed outcome), correct model specification, and no information bias, the g-formula can generate estimates of health outcomes that we would expect under interventions on exposure, such as introduction of regulatory exposure limits (19). The g-formula is suited to situations where time-varying confounders affected by exposure are present, as is the case with healthy-worker survivor bias in occupational data (31).

Our results are based on estimates of risk had historical exposures been compliant with theoretical exposure limits based on current PELs. In addition, they are based on an observed population and do not rely on assumptions of transportability of any exposure effects to external populations (as would be the case in life-table calculations of risk, where effect estimates based on observational studies are applied in hypothetical populations with baseline risk comparable to that observed in the general population). Interventions considered in the current study were of the nature of a theoretical maximum and did not allow exposures above this maximum but did allow for a range of exposures in the population under the maximum. In addition, employment status was predicted as a function of the observed data, counterfactual exposure, and covariate histories under each intervention, rather than implicitly being intervened on, as in “worst-case scenario” risk assessment estimates with fixed exposures and durations of employment for everyone in the target population (5, 6). Our results are therefore more representative of real-world settings where length of employment tenure varies, and compliance with a certain exposure limit still results in an exposure distribution below the limit rather than a fixed exposure value. They are not, however, directly comparable to previously reported risk estimates (5, 6).

In the assessment of interventions in occupational studies, counterfactual consistency may depend on how well the actual interventions necessary to reduce exposure levels are represented (with respect to the outcome) by the mechanisms that result in low exposures in the observed data. The assumption of positivity also stipulates that the probability of any counterfactual exposure history of interest is nonzero within levels of covariates. Positivity violations were not a major issue, as even in the case of the hypothetical intervention where exposures are completely removed in our study, 18% of the observed, actively employed person-time was actually unexposed. From the regulatory standpoint, agencies such as the Occupational Safety and Health Administration and MSHA are required to promulgate standards that reduce significant risk to the extent that it is feasible, with both technological and economic feasibility considerations considered (32), so potential limits are bounded by additional factors beyond health concerns. However, even if we assume that cost was not an issue, an intervention setting crystalline silica exposures to zero is not considered feasible for this industry given the ubiquitous nature of the exposure. Results under this intervention in the current study, however, are useful as they allowed us to estimate potential attributable fractions in the observed population.

In the present study, we estimated absolute risk with respect to the interventions of interest. Estimates of risk are of direct relevance to risk quantification and management, and measures of effect based on risk, such as the risk ratio and risk difference, are collapsible. (Collapsibility here refers in general terms to the property of a measure of association of an exposure and outcome, which is constant across levels of a covariate, where the marginal measure computed ignoring the covariate will equal the constant stratum-specific estimates across the levels of the covariate (33).) In cases of noncollapsibility, conditional estimates of outcome parameters may differ from marginal ones, even in the absence of confounding (34–37). The property of collapsibility makes it preferable to use risk differences and risk ratio as effect measures, rather than contrasts based on the rate, hazard, or odds, which are not collapsible parameters.

We also addressed the issue of competing risks by estimating subdistribution cumulative incidence of each of the events of interest (20, 21). Under the assumptions highlighted above, we estimated risk of a particular event of interest (e.g., lung cancer mortality) under an intervention, while accounting for changes in the composite or net cumulative incidence function of all competing events over time. In studies of exposures such as silica, which may exert effects on more than one cause of death, estimation of effects of specific causes of death becomes complicated. While it is common practice to treat competing events in survival analysis as right-censored observations, interpretation of effect estimates under this approach can be limited (25), and competing events in survival data can lead to biased effect estimates because of collider bias (38, 39).

The subdistribution estimates also have advantages over estimates that address competing events with use of inverse probability of censoring weights, which essentially assess the effect of a joint intervention composed of the intervention with respect to the exposure of interest together with an intervention eliminating the competing event (40). While it may sometimes be useful to envision a population where no one experiences a competing event, in many situations, including the current study, this is not a well-defined population (38, 41). The identifiability of target parameters where no one experienced competing events also requires strong assumptions, such as independence of competing events given measured covariates (42, 43). Estimation of subdistribution risk enabled us to quantify counterfactual risk in a population where competing events were present and also affected by exposure, which is more realistic and more informative for a public health framework. We did, however, assume that no participants were lost to follow-up in the simulations for the parametric g-formula. This assumption does not suffer from the same limitations of interpretability as no one dying of competing events, because it deals with actual right-censored observations and simply envisions a study without loss to follow-up (40).

Results for nonmalignant respiratory disease were comparable on the relative scale, and greater risk reductions were observed on the absolute scale, when using a wide definition of the outcome compared with a more restrictive definition excluding infectious respiratory disease cases. Our finding suggests that silica exposure increases risk of infectious, as well as noninfectious, respiratory disease. Evidence from both in vitro and in vivo studies suggests that silica exposures result in inflammatory responses and immune dysfunction, including toxicity to macrophages and T-cell dysregulation (44–48). This may in turn compromise lung defenses against pneumonia and other infectious diseases, and mortality from these diseases could therefore be a source of excess mortality risk associated with silica exposure.

A limitation of the study was the lack of refined adjustment for smoking. Although smoking information was available for 50% of the workers in the present study, information on smoking intensity or duration were not available; thus we relied on ever/never-smoker classification for confounding control. Smoking should certainly be considered a common cause of both events of interest (as well as all-cause mortality) and also appears to be a positive confounder in this cohort (4). Inability to fully adjust for smoking could result in potential residual confounding. The use of qualitative rather than quantitative information regarding asbestos exposures was a similar limitation. Work history files were available through 1994, with approximately 12% of the cohort still employed at that time. Potential exposures accrued after 1994 are not quantified for this small portion of the population, although given the 15-year lag applied in the analyses the potential window for this misclassification is small. A general limitation of the parametric g-formula is the increased sensitivity to model misspecification due to the increased number of parametric models. As a means of validation, we compared our simulated results under no intervention with the observed data and showed a close match (Web Figure 2), providing assurance that our choices of models were reasonable. Furthermore, our results appeared robust to modeling assumptions regarding exposure models.

In summary, we estimated cumulative risk of lung cancer and nonmalignant respiratory disease mortality, 2 causes of death linked to occupational silica exposure, under hypothetical interventions while addressing issues of both time-varying confounding affected by prior exposure and competing events. Our findings suggest that mortality risk from both outcomes would have been considerably lower if historical occupational silica exposures in this cohort had not exceeded current regulatory limits.

Click here for additional data file.