Mortality-Air Pollution Associations in Low Exposure Environments (MAPLE): Phase 2

Abstract

Introduction: Mortality is associated with long-term exposure to fine particulate matter (particulate matter ≤2.5 μm in aerodynamic diameter; PM_2.5), although the magnitude and form of these associations remain poorly understood at lower concentrations. Knowledge gaps include the shape of concentration-response curves and the lowest levels of exposure at which increased risks are evident and the occurrence and extent of associations with specific causes of death. Here, we applied improved estimates of exposure to ambient PM_2.5 to national population-based cohorts in Canada, including a stacked cohort of 7.1 million people who responded to census year 1991, 1996, or 2001. The characterization of the shape of the concentration-response relationship for nonaccidental mortality and several specific causes of death at low levels of exposure was the focus of the Mortality-Air Pollution Associations in Low Exposure Environments (MAPLE) Phase 1 report. In the Phase 1 report we reported that associations between outdoor PM_2.5 concentrations and nonaccidental mortality were attenuated with the addition of ozone (O₃) or a measure of gaseous pollutant oxidant capacity (O_x), which was estimated from O₃ and nitrogen dioxide (NO₂) concentrations. This was motivated by our interests in understanding both the effects air pollutant mixtures may have on mortality and also the role of O₃ as a copollutant that shares common sources and precursor emissions with those of PM_2.5. In this Phase 2 report, we further explore the sensitivity of these associations with O₃ and O_x, evaluate sensitivity to other factors, such as regional variation, and present ambient PM_2.5 concentration-response relationships for specific causes of death.

Methods: PM_2.5 concentrations were estimated at 1 km² spatial resolution across North America using remote sensing of aerosol optical depth (AOD) combined with chemical transport model (GEOS-Chem) simulations of the AOD:surface PM_2.5 mass concentration relationship, land use information, and ground monitoring. These estimates were informed and further refined with collocated measurements of PM_2.5 and AOD, including targeted measurements in areas of low PM_2.5 concentrations collected at five locations across Canada. Ground measurements of PM_2.5 and total suspended particulate matter (TSP) mass concentrations from 1981 to 1999 were used to backcast remote-sensing-based estimates over that same time period, resulting in modeled annual surfaces from 1981 to 2016.

Annual exposures to PM_2.5 were then estimated for subjects in several national population-based Canadian cohorts using residential histories derived from annual postal code entries in income tax files. These cohorts included three census-based cohorts: the 1991 Canadian Census Health and Environment Cohort (CanCHEC; 2.5 million respondents), the 1996 CanCHEC (3 million respondents), the 2001 CanCHEC (3 million respondents), and a Stacked CanCHEC where duplicate records of respondents were excluded (Stacked CanCHEC; 7.1 million respondents). The Canadian Community Health Survey (CCHS) mortality cohort (mCCHS), derived from several pooled cycles of the CCHS (540,900 respondents), included additional individual information about health behaviors. Follow-up periods were completed to the end of 2016 for all cohorts. Cox proportional hazard ratios (HRs) were estimated for nonaccidental and other major causes of death using a 10-year moving average exposure and 1-year lag. All models were stratified by age, sex, immigrant status, and where appropriate, census year or survey cycle. Models were further adjusted for income adequacy quintile, visible minority status, Indigenous identity, educational attainment, labor-force status, marital status, occupation, and ecological covariates of community size, airshed, urban form, and four dimensions of the Canadian Marginalization Index (Can-Marg; instability, deprivation, dependency, and ethnic concentration). The mCCHS analyses were also adjusted for individual-level measures of smoking, alcohol consumption, fruit and vegetable consumption, body mass index (BMI), and exercise behavior.

In addition to linear models, the shape of the concentration-response function was investigated using restricted cubic splines (RCS). The number of knots were selected by minimizing the Bayesian Information Criterion (BIC). Two additional models were used to examine the association between nonaccidental mortality and PM_2.5. The first is the standard threshold model defined by a transformation of concentration equaling zero if the concentration was less than a specific threshold value and concentration minus the threshold value for concentrations above the threshold. The second additional model was an extension of the Shape Constrained Health Impact Function (SCHIF), the eSCHIF, which converts RCS predictions into functions potentially more suitable for use in health impact assessments. Given the RCS parameter estimates and their covariance matrix, 1,000 realizations of the RCS were simulated at concentrations from the minimum to the maximum concentration, by increments of 0.1 μg/m³. An eSCHIF was then fit to each of these RCS realizations. Thus, 1,000 eSCHIF predictions and uncertainty intervals were determined at each concentration within the total range.

Sensitivity analyses were conducted to examine associations between PM_2.5 and mortality when in the presence of, or stratified by tertile of, O₃ or O_x. Additionally, associations between PM_2.5 and mortality were assessed for sensitivity to lower concentration thresholds, where person-years below a threshold value were assigned the mean exposure within that group. We also examined the sensitivity of the shape of the nonaccidental mortality-PM_2.5 association to removal of person-years at or above 12 μg/m³ (the current U.S. National Ambient Air Quality Standard) and 10 μg/m³ (the current Canadian and former [2005] World Health Organization [WHO] guideline, and current WHO Interim Target-4). Finally, differences in the shapes of PM_2.5-mortality associations were assessed across broad geographic regions (airsheds) within Canada.

Results: The refined PM_2.5 exposure estimates demonstrated improved performance relative to estimates applied previously and in the MAPLE Phase 1 report, with slightly reduced errors, including at lower ranges of concentrations (e.g., for PM_2.5 <10 μg/m³).

Positive associations between outdoor PM_2.5 concentrations and nonaccidental mortality were consistently observed in all cohorts. In the Stacked CanCHEC analyses (1.3 million deaths), each 10-μg/m³ increase in outdoor PM_2.5 concentration corresponded to an HR of 1.084 (95% confidence interval [CI]: 1.073 to 1.096) for nonaccidental mortality. For an interquartile range (IQR) increase in PM_2.5 mass concentration of 4.16 μg/m³ and for a mean annual nonaccidental death rate of 92.8 per 10,000 persons (over the 1991-2016 period for cohort participants ages 25-90), this HR corresponds to an additional 31.62 deaths per 100,000 people, which is equivalent to an additional 7,848 deaths per year in Canada, based on the 2016 population. In RCS models, mean HR predictions increased from the minimum concentration of 2.5 μg/m³ to 4.5 μg/m³, flattened from 4.5 μg/m³ to 8.0 μg/m³, then increased for concentrations above 8.0 μg/m³. The threshold model results reflected this pattern with -2 log-likelihood values being equal at 2.5 μg/m³ and 8.0 μg/m³. However, mean threshold model predictions monotonically increased over the concentration range with the lower 95% CI equal to one from 2.5 μg/m³ to 8.0 μg/m³. The RCS model was a superior predictor compared with any of the threshold models, including the linear model.

In the mCCHS cohort analyses inclusion of behavioral covariates did not substantially change the results for both linear and nonlinear models. We examined the sensitivity of the shape of the nonaccidental mortality-PM_2.5 association to removal of person-years at or above the current U.S. and Canadian standards of 12 μg/m³ and 10 μg/m³, respectively. In the full cohort and in both restricted cohorts, a steep increase was observed from the minimum concentration of 2.5 μg/m³ to 5 μg/m³. For the full cohort and the <12 μg/m³ cohort the relationship flattened over the 5 to 9 μg/m³ range and then increased above 9 μg/m³. A similar increase was observed for the <10 μg/m³ cohort followed by a clear decline in the magnitude of predictions over the 5 to 9 μg/m³ range and an increase above 9 μg/m³. Together these results suggest that a positive association exists for concentrations >9 μg/m³ with indications of adverse effects on mortality at concentrations as low as 2.5 μg/m³.

Among the other causes of death examined, PM_2.5 exposures were consistently associated with an increased hazard of mortality due to ischemic heart disease, respiratory disease, cardiovascular disease, and diabetes across all cohorts. Associations were observed in the Stacked CanCHEC but not in all other cohorts for cerebrovascular disease, pneumonia, and chronic obstructive pulmonary disease (COPD) mortality. No significant associations were observed between mortality and exposure to PM_2.5 for heart failure, lung cancer, and kidney failure.

In sensitivity analyses, the addition of O₃ and O_x attenuated associations between PM_2.5 and mortality. When analyses were stratified by tertiles of copollutants, associations between PM_2.5 and mortality were only observed in the highest tertile of O₃ or O_x. Across broad regions of Canada, linear HR estimates and the shape of the eSCHIF varied substantially, possibly reflecting underlying differences in air pollutant mixtures not characterized by PM_2.5 mass concentrations or the included gaseous pollutants. Sensitivity analyses to assess regional variation in population characteristics and access to healthcare indicated that the observed regional differences in concentration-mortality relationships, specifically the flattening of the concentration-mortality relationship over the 5 to 9 μg/m³ range, was not likely related to variation in the makeup of the cohort or its access to healthcare, lending support to the potential role of spatially varying air pollutant mixtures not sufficiently characterized by PM_2.5 mass concentrations.

Conclusions: In several large, national Canadian cohorts, including a cohort of 7.1 million unique census respondents, associations were observed between exposure to PM_2.5 with nonaccidental mortality and several specific causes of death. Associations with nonaccidental mortality were observed using the eSCHIF methodology at concentrations as low as 2.5 μg/m³, and there was no clear evidence in the observed data of a lower threshold, below which PM_2.5 was not associated with nonaccidental mortality.

Statement Figure.

Shape of the association between outdoor PM_2.5 exposure and nonaccidental death. This plot shows how the risk of death changes over different PM_2.5 exposure concentrations. The relative risk of death compares the lowest observed PM_2.5 concentration (2.5 μg/m³) to all higher concentrations. (Adapted from Investigators’ Report Figure 29.)

Figure 1.

Schematic of the exposure development process for PM_2.5. GWR = geographically weighted regression.

Figure 2.

Location of collocated ground-based measurements of PM_2.5 and AOD. The background shows satellite-based estimates of PM_2.5 from (A) van Donkelaar et al. (2016) and (B) population density.

Figure 3.

Map of the 6 airsheds in Canada across the provinces and territories, with locations of large cities (black circles) and MAPLE PM_2.5 monitoring sites (red Xs).

Figure 4.

Average contribution of five major PM_2.5 chemical components to total PM_2.5 mass measured at MAPLE sampling sites in Downsview, Ontario (East Central Airshed); Halifax, Nova Scotia (Southern Atlantic Airshed); Kelowna, British Columbia (Western Airshed); Lethbridge, Alberta (Prairie Airshed); and Sherbrooke, Quebec (East Central Airshed).

Figure 5.

Comparison of mean PM_2.5 mass concentrations for 2000–2012 observed by in-situ ground-based monitors with (A) Phase 1 V4.NA.01 and (B) Phase 2 V4.NA.02-MAPLE satellite-derived PM_2.5 estimates for all North American monitor locations. Annotations include the root mean square difference (RMSD), line of best fit (y), coefficient of determination (R²) and the number of points (N). See Figure 6 for Canadian locations only.

Figure 6.

Comparison of mean PM_2.5 mass concentrations for 2000–2012 observed by in-situ ground-based monitors with (A) Phase 1 V4.NA.01 and (B) Phase 2 V4.NA.02-MAPLE satellite-derived PM_2.5 estimates for Canadian monitor locations and observed concentrations below 10 μg/m³. Annotations include the root mean square difference (RMSD), line of best fit (y), coefficient of determination (R²) and the number of points (N).

Figure 7.

Estimates of fine particulate matter (PM_2.5) annual means averaged over the entire study period (1981–2015). City-level estimates for the largest cities are shown in the circles.

Figure 8.

Estimates of fine particulate matter (PM_2.5) annual means averaged over decades within the study period (1981–1990, 1991–2000, 2001–2010, and 2011–2015). City-level estimates for the largest cities are shown in the circles.

Figure 9.

Estimates of O₃ annual means from 8-hour daily maxima, averaged over the entire study period (1981–2015). City-level estimates for the largest cities are shown in the insets.

Figure 10.

Estimates of O₃ annual means averaged over decades within the study period (1981–1990, 1991–2000, and 2001–2011).

Figure 11.

Estimates of O_x annual means calculated from annual O₃ and NO₂, averaged over the entire study period (1981–2015). City-level estimates for the largest cities are shown in the insets.

Figure 12.

Estimates of O_x annual means averaged over decades within the study period (1981–1990, 1991–2000, and 2001–2011).

Figure 13.

Predicted relative risk of nonaccidental mortality by PM_2.5 concentration (μg/m³) for the 1991, 1996, and 2001 CanCHEC cohorts. Predictions are (left column) unadjusted for either O₃ or O_x, (center column) adjusted for O₃, and (right column) adjusted for O_x. Mean predictions are displayed as solid blue lines with 95% CIs as grey-shaded areas. Green x-axis tick marks show RCS knot locations.

Figure 14.

Predicted relative risk of nonaccidental mortality by PM_2.5 concentration (μg/m³) for the Stacked CanCHEC cohort by number of knots used in the RCS fit [3–18 knots; Rows 1–4]. Mean predictions are displayed as solid blue lines with 95% CIs as grey-shaded areas. Green x-axis tick marks show RCS knot locations. The pink curve at 9 knots has the lowest BIC value. From left to right the last row of panels shows -2 log-likelihood values, AIC, BIC, and the highest concentration below which the RCS lower confidence limit is less than one.

Figure 14.

Predicted relative risk of nonaccidental mortality by PM_2.5 concentration (μg/m³) for the Stacked CanCHEC cohort by number of knots used in the RCS fit [3–18 knots; Rows 1–4]. Mean predictions are displayed as solid blue lines with 95% CIs as grey-shaded areas. Green x-axis tick marks show RCS knot locations. The pink curve at 9 knots has the lowest BIC value. From left to right the last row of panels shows -2 log-likelihood values, AIC, BIC, and the highest concentration below which the RCS lower confidence limit is less than one.

Figure 15.

RCS-based predicted relative risk by PM_2.5 concentration (μg/m³) for the mCCHS cohort by: (left column) unadjusted and (right column) adjusted for behavioral covariates. Corresponding predictions adjusted for either O₃ (middle row) or O_x (bottom row) are also presented. Mean predictions are displayed as solid blue lines with 95% CIs as grey-shaded areas. Green x-axis tick marks show RCS knot locations.

Figure 16.

RCS-based predicted relative risk by PM_2.5 concentration (μg/m³) for the Stacked CanCHEC cohort by cause of death. Mean predictions are displayed as solid blue lines with 95% CIs as grey-shaded areas. Green x-axis tick marks show RCS knot locations.

Figure 17.

RCS-based predicted relative risk by PM_2.5 concentration (μg/m³) for the Stacked CanCHEC cohort by aggregated causes of death: (top row) nonaccidental; (middle row) cardiovascular; (bottom row) respiratory. From left to right the columns show: PM_2.5 concentration unadjusted for O₃, adjusted for linear O₃, adjusted for RCS of O₃; and O₃ RCS adjusted for RCS PM_2.5. Mean predictions are displayed as solid blue lines with 95% CIs as grey-shaded areas. Green x-axis tick marks show RCS knot locations.

Figure 18.

RCS-based predicted relative risk by PM_2.5 concentration (μg/m³) for the Stacked CanCHEC cohort by aggregated causes of death: (top row) nonaccidental; (middle row) cardiovascular; (bottom row) respiratory. The left column shows PM_2.5 concentration unadjusted for O_x; the right column shows PM_2.5 concentration adjusted for linear O_x. Mean predictions are displayed as solid blue lines with 95% CIs as grey- shaded areas. Green x-axis tick marks show RCS knot locations.

Figure 19.

RCS-based predicted relative risk of nonaccidental mortality by PM_2.5 concentration (μg/m³) for the Stacked CanCHEC cohort by (top row) O₃ or (bottom row) O_x. From left to right the columns show: the first, second, and third tertiles. Mean predictions are displayed as solid blue lines with 95% CIs as grey-shaded areas. Green x-axis tick marks show RCS knot locations.

Figure 20.

Threshold model comparisons with linear and RCS models. (A) -2 log-likelihood values (orange dots) by threshold concentration; threshold concentrations for models representing statistically significantly improved fit (P< 0.05; orange dots lying below purple dashed line) over the linear model (red dot) using the likelihood-ratio test; and RCS -2 log-likelihood value (black dot). (B) Distribution of estimated threshold values under the ensemble model. (C) Predicted relative risk by PM_2.5 concentration (μg/m³) for nonaccidental mortality in the Stacked CanCHEC cohort; mean predictions (solid blue line with 95% CIs as grey-shaded areas); RCS mean predictions (solid red line) and 95% CIs (dashed red lines).

Figure 21.

Predicted relative risk by PM_2.5 concentration (μg/m³) for the Stacked CanCHEC cohort by airshed. Mean predictions are displayed as solid blue lines with 95% CIs as grey-shaded areas. Green x-axis tick marks show RCS knot locations.

Figure 22.

RCS predictions over PM_2.5 concentration range based on 9 knots with and without adjustment for proximity to healthcare resources. Mean and uncertainty bounds with adjustment (blue line and blue-shaded area) and without adjustment (black line and grey-shaded area).

Figure 23.

RCS sensitivity analyses for nonaccidental mortality in the Stacked CanCHEC, using BIC on the fully adjusted model for the best-fitting spline for two airsheds that show a large dip in the concentration–response shape and removal of population groups for which linkage to mortality may have been less complete. Black line: full cohort; red line: immigrants removed from cohort; blue line: Indigenous respondents removed from cohort.

Figure 24.

RCS predictions over PM_2.5 concentration range based on 9 knots specifying age as 1-year or 5-year strata. Mean and uncertainty bounds with 1- and 5-year strata (blue line and blue-shaded area) (black line and grey-shaded area).

Figure 25.

RCS predictions for nonaccidental mortality over the PM_2.5 concentration range for the full cohort, person-years with PM_2.5 <12 μg/m³, and person-years with PM_2.5 <10 μg/m³.

Figure 26.

Comparison of standard and simulation-based methods to determine mean and 95% confidence intervals. (A) Relative risk estimates and (B) PAF Estimates (%) per 0.1-μg/m³ are also presented. Mean and 95% CI: Standard = dashed red line and pink-shaded area; Simulated = dashed blue line and blue dotted lines

Figure 27.

2001 cohort nonaccidental death RCS mean relative risk predictions over the PM_2.5 concentration range, with eSCHIF relative risk predictions. (A) Relative risk estimates; (B) PAF estimates. Mean and 95% CI: RCS = solid red line and pink-shaded area; eSCHIF = solid blue line and dashed blue lines. Green x-axis tick marks show RCS knot locations.

Figure 28.

2001 cohort nonaccidental death RCS mean relative risk predictions adjusted for O₃ over the PM_2.5 concentration range, with eSCHIF relative risk predictions. (A) Relative risk estimates; (B) PAF estimates. Mean and 95% CI: RCS = solid red line and pink-shaded area; eSCHIF = solid blue line and dashed blue lines. Green x-axis tick marks show RCS knot locations.

Figure 29.

Stacked cohort nonaccidental death RCS (9 knots) mean relative risk predictions over the PM_2.5 concentration range, with eSCHIF relative risk predictions. (A) Relative risk estimates; (B) PAF estimates. Mean and 95% CI: RCS = solid red line and pink-shaded area; eSCHIF = solid blue line and dashed blue lines. Green x-axis tick marks show RCS knot locations.

Figure 30.

Stacked cohort nonaccidental death ensemble RCS mean relative risk predictions over the PM_2.5 concentration range, with eSCHIF relative risk predictions. (A) Relative risk estimates; (B) PAF estimates. Mean and 95% CI: RCS = solid red line and pink-shaded area; eSCHIF = solid blue line and dashed blue lines.

Commentary Figure 1.

MAPLE study cohorts. The five cohorts included CanCHEC 1991, CanCHEC 1996, CanCHEC 2001, Stacked CanCHEC, and CCHS. The Stacked CanCHEC and CCHS cohorts included groups of respondents that entered the study in different census or survey years (dashed lines). Participants were followed until death or the end of the study.

Commentary Figure 2.

Airsheds of Canada. The associations between PM_2.5 and mortality were also examined by airshed because regional geographical features and weather conditions influence ambient air quality.

Commentary Figure 3.

Ambient PM_2.5 exposure and nonaccidental mortality in the Stacked CanCHEC cohort. Ambient PM_2.5 exposure with a 10-year moving average and 1-year lag was associated with higher total nonaccidental and select cause-specific mortality rates using a linear model and controlling for individual- and community-level sociodemographic variables.

Commentary Figure 4.

Concentration–response curves for ambient PM_2.5 exposure and relative risk of nonaccidental mortality in the Stacked CanCHEC cohort. The RCS and extended SCHIF 95% CIs are wider at low and high PM_2.5 levels to reflect greater uncertainty in the hazard ratio at these levels of exposure relative to the mean concentration. (Adapted from Investigators’ Report Figures 20 and 29.)