Chronic Exposure to Fine Particles and Mortality, 1974-2009
Chronic Exposure to Fine Particles and Mortality, 1974-2009
The Harvard Six Cities study population has been previously described (Dockery et al. 1993). Briefly, adults were randomly sampled from six cities in the eastern and midwestern United States between 1974 and 1977: in 1974, Watertown, Massachusetts; in 1975, Kingston and Harriman, Tennessee, and specific census tracts of St. Louis, Missouri; in 1976, Steubenville, Ohio, and Portage, Wyocena, and Pardeeville, Wisconsin; and in 1977, Topeka, Kansas. Information on age, sex, weight, height, educational level, smoking history, hypertension, and diabetes was collected by questionnaire at enrollment. All participants underwent spirometry tests at enrollment (Dockery et al. 1985) and chronic obstructive pulmonary disease (COPD) was defined as having
(FEV1 ÷ FVC) < 70%,
where FEV1 is forced expiratory volume in 1 sec, and FVC is forced vital capacity. This analysis, as in the previous analyses, was restricted to 8,096 white participants with acceptable pulmonary function measurements. The study was approved by the Harvard School of Public Health Human Subjects Committee and all participants signed an informed consent before participation.
Vital status and cause of death were determined by searching the National Death Index (NDI) for calendar years 1979–2009. Deaths before the NDI started in 1979 were identified by next of kin and Social Security records, and the cause of death was determined by a certified nosologist who reviewed death certificates (Dockery et al. 1993).
Survival times were calculated from enrollment until death or the end of follow-up (31 December 2009). For the 6 participants who were lost to follow-up before 1979, the censored survival times were calculated from enrollment to date of the last follow-up contact plus 6 months or the first day of the NDI (1 January 1979), whichever came first. For each cause of death category, participants who died from another cause were censored at time of death.
Annual PM2.5 concentration was assigned for each participant until death or censoring. PM2.5 concentration was measured in the participant's city by a centrally located monitor from 1979 to 1986–1988, depending on the city (Dockery et al. 1993). Therefore, the study has no spatial contrast on the within-city scale. PM2.5 concentrations for the years before monitoring started were assumed to be equal to the earliest monitored year. From the end of monitoring until 1998, PM2.5 concentration was estimated from PM10 (aerodynamic diameter < 10 μm) data from U.S. EPA monitors and visibility (extinction) data from the National Weather Service (Laden et al. 2006). From 1999 through 2009, direct measurements of PM2.5 were available from U.S. EPA monitors. For sensitivity analyses, we also predicted PM2.5 for 1999–2009 (correlation between predicted and measured was 0.97) using the formula applied to derive exposure estimates during the earlier period when PM2.5 was not measured.
We first replicated the original analysis separately for all-cause mortality, cardiovascular mortality as coded by the International Classification of Diseases, 9th Revision [ICD-9; World Health Organization (WHO) 1977] or the 10th Revision (ICD-10; WHO 1992), 400.0–440.9, I10.0–I70.9, respectively, lung-cancer mortality (ICD-9 162, ICD-10 C33.0–C34.9), and COPD mortality (ICD-9 490.0–496.0, ICD-10 J40.0–J47.0) for the 36-year follow-up from 1974 to 2009 using a Cox proportional hazards model with follow-up time as the time scale (Dockery et al. 1993; Laden et al. 2006). PM2.5 was included in each model as an annual time-dependent variable. The model was stratified by sex, age (1-year intervals) and time in the study (1-year intervals), so that each age/sex group had its own baseline hazard for each year of follow-up. The analysis was adjusted for potential confounders collected at baseline: smoking status (never, former, current), cumulative smoking (pack-years included separately for current and former smokers), educational level (< high-school, ≥ high school), and a linear and quadratic term for body mass index (BMI; kilograms per meter squared), using the Cox proportional hazards model formulated as follows:
hi s(t) = h0s(t) exp[β1Xi + β2Zi(t)], [1]
where hi is the instantaneous hazard probability of death for subject i in stratum s (defined by sex, age, and time in the study), h0s(t) is the baseline hazard function, Xi is the vector of time-independent variables, and Zi(t) is the vector of time-dependent variables. We evaluated models with 1-year (i.e., exposure during the year before death or censure) to 5-year lagged moving averages and chose the best fit model using Akaike's information criterion (AIC) (Akaike 1973). The best fit moving average was determined from participants who survived at least 5 years from enrollment, so that AIC criteria were evaluated among populations with comparable sizes. We then estimated mortality rate ratios (RR) associated with PM2.5 exposure during the best fit moving average on the whole sample size. Once the best exposure window was determined, we fit a penalized spline model using a cubic regression spline with 12 knots to estimate the shape of the concentration–response relation, and chose the optimal degree of freedom by minimizing AIC and evaluated nonlinearity with a Wald test. We investigated whether PM2.5 advanced date of death for participants with chronic conditions at enrollment. We also investigated the potential for effect modification of PM2.5 on mortality by smoking status at enrollment using interaction terms between such variables and PM2.5. Finally, we tested the hypothesis that the effect of PM2.5 changed over time by dividing the follow-up into four equally spaced time periods and testing interactions between period and PM2.5.
We performed sensitivity analyses using a second-degree polynomial distributed lag model to allow the effects of PM2.5 exposure to be distributed from 1 to 5 years before death or censor (Lepeule et al. 2006; Schwartz 2000); using predicted PM2.5 concentrations after 1999 instead of the measured PM2.5; considering only deaths from natural causes, with external causes of deaths (ICD-9 E800–E999, ICD-10 S00–T88 and V00–Y99) being censored at time of death; and considering only deaths that occurred in the state where the participants lived at enrollment. We next investigated the robustness of the results to alternative modeling assumptions by using a Poisson model with dummy variables for each year of follow-up, which is equivalent to a piecewise exponential proportionate hazard model with the baseline hazard changing each year (Laird and Oliver 1981):
log μit = log Eit + γtTt + β1Xi + β2Zi(t), [2]
where μit is the expected value of the death indicator for subject i at time t, Eit is the exposure duration of subject i at time t (log Eit being the offset), Tt is the vector of dummy variables for time by 1 year (piece-wise baseline hazard), Xi is the vector of the time-independent covariates, and Zi(t) is the vector of time-dependent variables. Using this Poisson survival analysis, we first compared the results to the Cox model and then relaxed the proportionate hazard assumption for sex, education, and cumulative smoking by including interaction terms of these variables with each year of follow-up. As an alternative to the previous analyses (Dockery et al. 1993; Laden et al. 2006), we used age in 5-year groups as the time scale, and adjusted the model for time trends (linear term). For specific causes of death, convergence issues led us to group age by 10 years. We then fit penalized spline models. Because RRs may vary over time and period-specific RRs may be biased, we used the Poisson model to calculate adjusted survival curves (Hernan 2010). We included product terms between PM2.5 and time in model 2 [Equation 2], thereby allowing the effect of PM2.5 to flexibly vary from year to year. We then predicted the survival probability for each year of follow-up for each participant under three scenarios using concentrations of PM2.5 throughout the entire follow-up period equal to 10, 15, or 20 μg/m.
p-Values < 0.05 were considered statistically significant. All analyses were repeated separately for all- and specific-causes of deaths. Analyses were conducted with SAS software, version 9.2 (SAS Institute Inc., Cary, NC) and R statistical software, version 2.12.2 (R Foundation for Statistical Computing, Vienna, Austria).
Methods
Study Population
The Harvard Six Cities study population has been previously described (Dockery et al. 1993). Briefly, adults were randomly sampled from six cities in the eastern and midwestern United States between 1974 and 1977: in 1974, Watertown, Massachusetts; in 1975, Kingston and Harriman, Tennessee, and specific census tracts of St. Louis, Missouri; in 1976, Steubenville, Ohio, and Portage, Wyocena, and Pardeeville, Wisconsin; and in 1977, Topeka, Kansas. Information on age, sex, weight, height, educational level, smoking history, hypertension, and diabetes was collected by questionnaire at enrollment. All participants underwent spirometry tests at enrollment (Dockery et al. 1985) and chronic obstructive pulmonary disease (COPD) was defined as having
(FEV1 ÷ FVC) < 70%,
where FEV1 is forced expiratory volume in 1 sec, and FVC is forced vital capacity. This analysis, as in the previous analyses, was restricted to 8,096 white participants with acceptable pulmonary function measurements. The study was approved by the Harvard School of Public Health Human Subjects Committee and all participants signed an informed consent before participation.
Mortality Follow-up
Vital status and cause of death were determined by searching the National Death Index (NDI) for calendar years 1979–2009. Deaths before the NDI started in 1979 were identified by next of kin and Social Security records, and the cause of death was determined by a certified nosologist who reviewed death certificates (Dockery et al. 1993).
Survival Time
Survival times were calculated from enrollment until death or the end of follow-up (31 December 2009). For the 6 participants who were lost to follow-up before 1979, the censored survival times were calculated from enrollment to date of the last follow-up contact plus 6 months or the first day of the NDI (1 January 1979), whichever came first. For each cause of death category, participants who died from another cause were censored at time of death.
Air Pollution Estimates
Annual PM2.5 concentration was assigned for each participant until death or censoring. PM2.5 concentration was measured in the participant's city by a centrally located monitor from 1979 to 1986–1988, depending on the city (Dockery et al. 1993). Therefore, the study has no spatial contrast on the within-city scale. PM2.5 concentrations for the years before monitoring started were assumed to be equal to the earliest monitored year. From the end of monitoring until 1998, PM2.5 concentration was estimated from PM10 (aerodynamic diameter < 10 μm) data from U.S. EPA monitors and visibility (extinction) data from the National Weather Service (Laden et al. 2006). From 1999 through 2009, direct measurements of PM2.5 were available from U.S. EPA monitors. For sensitivity analyses, we also predicted PM2.5 for 1999–2009 (correlation between predicted and measured was 0.97) using the formula applied to derive exposure estimates during the earlier period when PM2.5 was not measured.
Statistical Analysis
We first replicated the original analysis separately for all-cause mortality, cardiovascular mortality as coded by the International Classification of Diseases, 9th Revision [ICD-9; World Health Organization (WHO) 1977] or the 10th Revision (ICD-10; WHO 1992), 400.0–440.9, I10.0–I70.9, respectively, lung-cancer mortality (ICD-9 162, ICD-10 C33.0–C34.9), and COPD mortality (ICD-9 490.0–496.0, ICD-10 J40.0–J47.0) for the 36-year follow-up from 1974 to 2009 using a Cox proportional hazards model with follow-up time as the time scale (Dockery et al. 1993; Laden et al. 2006). PM2.5 was included in each model as an annual time-dependent variable. The model was stratified by sex, age (1-year intervals) and time in the study (1-year intervals), so that each age/sex group had its own baseline hazard for each year of follow-up. The analysis was adjusted for potential confounders collected at baseline: smoking status (never, former, current), cumulative smoking (pack-years included separately for current and former smokers), educational level (< high-school, ≥ high school), and a linear and quadratic term for body mass index (BMI; kilograms per meter squared), using the Cox proportional hazards model formulated as follows:
hi s(t) = h0s(t) exp[β1Xi + β2Zi(t)], [1]
where hi is the instantaneous hazard probability of death for subject i in stratum s (defined by sex, age, and time in the study), h0s(t) is the baseline hazard function, Xi is the vector of time-independent variables, and Zi(t) is the vector of time-dependent variables. We evaluated models with 1-year (i.e., exposure during the year before death or censure) to 5-year lagged moving averages and chose the best fit model using Akaike's information criterion (AIC) (Akaike 1973). The best fit moving average was determined from participants who survived at least 5 years from enrollment, so that AIC criteria were evaluated among populations with comparable sizes. We then estimated mortality rate ratios (RR) associated with PM2.5 exposure during the best fit moving average on the whole sample size. Once the best exposure window was determined, we fit a penalized spline model using a cubic regression spline with 12 knots to estimate the shape of the concentration–response relation, and chose the optimal degree of freedom by minimizing AIC and evaluated nonlinearity with a Wald test. We investigated whether PM2.5 advanced date of death for participants with chronic conditions at enrollment. We also investigated the potential for effect modification of PM2.5 on mortality by smoking status at enrollment using interaction terms between such variables and PM2.5. Finally, we tested the hypothesis that the effect of PM2.5 changed over time by dividing the follow-up into four equally spaced time periods and testing interactions between period and PM2.5.
Sensitivity Analyses
We performed sensitivity analyses using a second-degree polynomial distributed lag model to allow the effects of PM2.5 exposure to be distributed from 1 to 5 years before death or censor (Lepeule et al. 2006; Schwartz 2000); using predicted PM2.5 concentrations after 1999 instead of the measured PM2.5; considering only deaths from natural causes, with external causes of deaths (ICD-9 E800–E999, ICD-10 S00–T88 and V00–Y99) being censored at time of death; and considering only deaths that occurred in the state where the participants lived at enrollment. We next investigated the robustness of the results to alternative modeling assumptions by using a Poisson model with dummy variables for each year of follow-up, which is equivalent to a piecewise exponential proportionate hazard model with the baseline hazard changing each year (Laird and Oliver 1981):
log μit = log Eit + γtTt + β1Xi + β2Zi(t), [2]
where μit is the expected value of the death indicator for subject i at time t, Eit is the exposure duration of subject i at time t (log Eit being the offset), Tt is the vector of dummy variables for time by 1 year (piece-wise baseline hazard), Xi is the vector of the time-independent covariates, and Zi(t) is the vector of time-dependent variables. Using this Poisson survival analysis, we first compared the results to the Cox model and then relaxed the proportionate hazard assumption for sex, education, and cumulative smoking by including interaction terms of these variables with each year of follow-up. As an alternative to the previous analyses (Dockery et al. 1993; Laden et al. 2006), we used age in 5-year groups as the time scale, and adjusted the model for time trends (linear term). For specific causes of death, convergence issues led us to group age by 10 years. We then fit penalized spline models. Because RRs may vary over time and period-specific RRs may be biased, we used the Poisson model to calculate adjusted survival curves (Hernan 2010). We included product terms between PM2.5 and time in model 2 [Equation 2], thereby allowing the effect of PM2.5 to flexibly vary from year to year. We then predicted the survival probability for each year of follow-up for each participant under three scenarios using concentrations of PM2.5 throughout the entire follow-up period equal to 10, 15, or 20 μg/m.
p-Values < 0.05 were considered statistically significant. All analyses were repeated separately for all- and specific-causes of deaths. Analyses were conducted with SAS software, version 9.2 (SAS Institute Inc., Cary, NC) and R statistical software, version 2.12.2 (R Foundation for Statistical Computing, Vienna, Austria).
Source...