How low levels of fine air pollution relate to death rates in Canada over three census periods

We created three new, separate analytical cohorts from the 1991, 1996, and 2001 CanCHECs. Briefly, the CanCHECs are population-based, administrative data cohorts that link eligible census respondents (i.e., noninstitutional respondents to the mandatory Statistics Canada long-form census questionnaire that is distributed to 20% of all Canadian households) to their annual mailing address (1981–2016) and follow subjects for mortality. Information on a number of variables capturing the social and economic status of the subjects was available from the long-form census (Table 1).

The linkage was approved by Statistics Canada (linkage requests 037-2016 and 045-2015) and is governed by the Directive on Microdata Linkage (Statistics Canada 2017a). Eligible respondents were first linked probabilistically to tax records using sex, date of birth, postal code (PC), and spousal date of birth (if available). This initial linkage was necessary since linkage to the mortality database is based on the social insurance number (SIN), a unique personal identifier. The long form censuses did not capture the SIN, but they are available on tax records. The linkage rate to tax records near the time of cohort inception was approximately 80%, of which 99% were determined to be accurate matches (Christidis et al. 2018; Pinault et al. 2016a; Wilkins et al. 2008).

Mortality and PC history data were attached to the census-tax cohorts using Statistics Canada’s Social Data Linkage Environment (SDLE) Derived Record Depository (DRD) (Statistics Canada 2017b), a dynamic relational database. About 99.8% of all deaths that occurred in Canada between 1991 and 2016 were linked to the DRD before being linked to eligible census respondents. From this linkage, we obtained death date and underlying cause of death if it occurred between census day and 31 December 2016. Mortality data were coded by underlying cause of death according to the International Classification of Diseases, 9th Revision, prior to 2000 (ICD-9; WHO 1977), and 10th Revision post-2000 (ICD-10; WHO 2016).

We enhanced the cohort with a number of data elements characterizing the environment in which each subject lived, using PC histories from tax records, of which the primary source was income tax filings (1981 to 2016). We assigned a representative point (latitude and longitude) to each PC (Statistics Canada 2017c). In large cities, PCs often correspond to a single block face, though in rural areas, they can range over much larger areas. Similarly, the point estimates of PCs are accurate within 0.2km in urban centers and 5.6km in rural areas (Khan et al. 2018). These point estimates were used to derive estimates of air pollution and location-based contextual risk factors.

We note that these three linked cohorts are newly created using an enhanced linkage environment (SDLE) and thus are not identical to the CanCHEC cohorts used in previous publications (Crouse et al. 2015; Pinault et al. 2017).

We used annual ambient PM2.5 surfaces as our main exposure of interest at a 0.01˚×0.01˚ resolution (∼1km2) over North America for 1981–2016 (Meng et al. 2019; van Donkelaar et al. 2015). PM2.5 estimates for the years 1998–2012 were developed by relating satellite-based retrievals of total column aerosol optical depth to near-surface PM2.5 concentrations using the geophysical relationship simulated by a chemical transport model (CTM). These estimates were constrained using ground-based monitoring from the National Air Pollution Surveillance (NAPS) program stations, along with other North America–based measurements, land-use information, and simulated composition in a geographically weighted regression (V4.NA.01; van Donkelaar et al. 2015). For years outside this period, we used PM2.5 surfaces developed using a backcasting method (Meng et al. 2019) that applied observed annual trends in ground monitoring data for PM2.5 and coarser size fractions to adjust pregridded PM2.5 estimates backwards or forwards in time. We estimated a 3-y moving-average exposure window with 1-y lag for assigning PM2.5 exposures for consistency with previous studies, as ambient PM2.5 is regulated based on a 3-y time window in Canada (CCME 2012).

We assigned estimates of exposures to ambient ozone (O3; as a May–September daily maximum 8-h average) and nitrogen dioxide (NO2; annual) for inclusion in multipollutant models. Additionally, we estimated a measure of the combined oxidant capacity of O3 and NO2, expressed as Ox=2/3O3+1/3NO2 (Weichenthal et al. 2017). We estimated a 3-y average with 1-y lag for each of O3, NO2, and Ox for inclusion in the hazard models. Modeled O3 surfaces at 21-km spatial resolution were developed by Environment and Climate Change Canada (ECCC) for 2002–2015 using chemical transport modeling informed by surface observations (Robichaud and Ménard 2014; Robichaud et al. 2016). Estimates of ambient NO2 were based on a national land-use regression model (LUR) developed for 2006 (Hystad et al. 2011) with a spatial resolution of 100m. The LUR estimates were built using satellite-derived NO2 (with 10-km resolution), distances to highways and major roads, and roadway kernel density gradients as predictive variables.

We temporally adjusted theandmodels to obtain exposure estimates over our study period (i.e., 1988–2015). Our adjustment was based on observed trends in ground monitoring data forandfrom the NAPS in Canada. For each of 24 census divisions (CDs) that had monitoring data available, we estimated yearly adjustment factors from the ratio of observed CD-average concentration in a specific year to the reference year(s) for which the original surfaces were estimated (i.e., 2006 forand 2002–2015 average for). We assigned adjustment factors for each PC from the closest CD. O 3 NO 2 NO 2 O 3 NO 2 O 3

We assigned contextual risk factors describing neighborhood-level characteristics and geographic identifiers using residential PC and data from the closest census (every 5 y from 1991 through 2016). We included in our analysis the Canadian Marginalization Index (CAN-Marg), population size of home community or city, an indicator of urban form, and regional airshed to capture risk factors beyond those captured at the subject level. We assigned these four categories of contextual covariates to residential PCs linked to census geography for each census year.

CAN-Marg is a publicly available index of neighborhood marginalization in Canada that was developed by Matheson et al. (2012) using an analysis of the 2001 and 2006 long-form census cycles. CAN-Marg consists of four dimensions that aim to capture different aspects of marginalization: material deprivation, residential instability, ethnic concentration, and dependency. Following the methodology of Matheson et al. (2012), we developed CAN-Marg using the 1991 and 1996 censuses. We assigned CAN-Marg to PC locations and then created quintiles (based on the cohort distribution) of the continuous values in the four Can-MARG dimensions in order to account for any potentially nonlinear associations with mortality.

We used a variable to describe the population size of a subject’s community (Pinault et al. 2017) (Table 2). We categorized geographic locations into the following: census metropolitan areas (CMAs) or census agglomerations (CAs; Statistics Canada 2003) with a population exceeding 1.5 million; 500,000–1.49 million; 100,000–499,999; 30,000–99,999; or 10,000–29,999, as well as non-CMAs/CAs. We note that although non-CMA/CAs are always rural areas, CMAs cover both the urban core of a city and the urban–rural fringe, such that some rural locations fall within a CMA/CA. As such, this variable does not perfectly delineate subjects living in rural vs. urban settings.

To further differentiate between the kinds of built environments and neighborhoods within communities, we created an urban form variable following the methodology developed by Gordon and Janzen (2013). This measure of urban form is informed by population density and the most frequently reported mode of transportation (active or transit) in each census tract as reported on each census cycle. The categories of this variable include an active urban core, transit-reliant suburb, car-reliant suburb, exurban, and non-CMA/CA. We note that mode of commute was not reported on the 1991 census cycle and was derived from the 1996 census.

We included airshed as a geographic covariate in our analysis (Crouse et al. 2016). Airsheds divide Canada into six regions (Western, Prairie, West Central, Southern Atlantic, East Central, and Northern) based on large-scale differences in air masses and meteorology. Airsheds can also be used to represent regional differences in mortality rates across Canada that remain uncaptured by other geographic covariates.

PC history was not available for each person in every year of follow-up, either because they did not file a tax return or from gaps in administrative data. Any gaps in PCs that had the same PC prior to and after the gap were assigned that PC for all years of the gap. After this imputation, 87.8% of person-years had an available PC. We imputed an additional 2.1% of person-years of missing PCs if the bounding PCs shared the first two characters (Finès et al. 2017; Pinault et al. 2017), totaling 89.9% of person-years with a PC.

After imputation, person-years were excluded if they did not have an assigned PC. Further exclusions of person-years occurred due to: immigrated to Canada less than 10 y before survey date (9,364,400 person-years), age during follow-up period exceeded 89 y (7,357,200), could not be linked to air pollution values (17,814,400), could not be linked to CAN-Marg values (25,973,900), could not be linked to CMA/CA size (25,613,100), could not be linked to airshed (25,545,500), 3-y moving average being informed by only 1 y of exposure (20,056,400), and year after subject death (17,936,100). The above are not mutually exclusive numbers of exclusions. The total available person-years for analyses were 150,996,500 after all exclusions (Figure S1).

Our primary statistical model relating exposure to mortality was the Cox proportional hazards model (Cox 1972). Participants were at least 25 y of age at the beginning of each cohort, and the time axis was the year of follow-up until 2016. Person-years before a census year and after a subject’s death year were excluded from analysis. Events were determined by year of death for nonaccidental causes. The Cox model baseline hazard function was stratified by age (5-y groups), sex, and immigrant status (yes or no). This latter strata variable was included since immigrants to Canada live longer, on average, than do Canadian-born citizens (Ng 2011). We excluded immigrants living in Canada for less than 10 y at cohort commencement due to the healthy immigrant effect (Ng 2011) and lack of knowledge of their historical air pollution exposures. Each subject was censored at 89 y of age, either at the start of each cohort or during follow-up, due to evidence suggesting an increased mismatch between home address and the tax return mailing address with increasing age (Bérard-Chagnon 2017). We postulate that relatives of elderly people were completing their tax returns. Each of the three CanCHEC cohorts (1991, 1996, and 2001) was examined separately. Estimates of the cohort-specific hazard ratios were then pooled to form a single summary hazard ratio. We also conducted a test for differences in the hazard ratios between cohorts (Cochran 1950).

We fit two covariate adjustment models for each cohort. The first was based on a directed acyclic graph (DAG; Figure S2) and consisted of all the geographically based predictors: CAN-Marg (four dimensions), airshed, urban form, and CMA/CA size. The second model, denoted as “Full,” additionally included the subject-level predictors (income, education, occupational class, Indigenous status, visible minority status, employment status, and marital status), which are not a priori causes of outdoor PM2.5 concentrations, but which may contribute to confounding owing to a chance imbalance across the PM2.5 distribution.

We also conducted analysis by categories of: immigrant status (yes or no), sex (male or female), and age during follow-up (, 65–74, or) for each cohort separately, again pooling the cohort-specific hazard ratio estimates among the three cohorts. In addition, we examined theassociation, adjusting for,, orby cohort. < 65 ≥ 75 y PM 2.5 O 3 NO 2 O x

The main purpose of the current study was to describe the association between PM2.5 and mortality in a manner that can be used for risk and benefits assessment. The standard approach is the log-linear (LL) model that relates the logarithm of the hazard ratio to exposure in a linear manner: logHR(PM2.5)=βPM2.5. Here, β represents a change in relative risk per unit change in concentration estimated using the Cox model. Nasari et al. (2016) developed the SCHIF in order to extend the LL model to nonlinear transformations of exposure, T(PM2.5), with the form: SCHIF(PM2.5)=θT(PM2.5). Nasari et al. (2016) proposed a specific family of transformations based on a sigmodal function that could accommodate a variety of shapes they suggested would be suitable for risk and benefits assessment. Here, θ represents a change in risk per unit change in T(PM2.5). The SCHIF can then be used in benefits assessment in a manner similar to the LL model after a suitable transformation of concentration.

The SCHIF approach has two major limitations: The first is in defining an appropriate number of transformations of a sigmodal function that can capture all shapes of interest; the second is that the method requires considerable computational capacity if the selected family is very large. This can be a serious limitation when cohort sizes are very large, such as with the CanCHECs.

Spline methods have also been proposed to characterize the shape. RCS with a few knots have been used (Beelen et al. 2014) in addition to smoothing splines (Di et al. 2017). However, the manner in which splines are presented by graphic representation of the mean predictions and uncertainty bounds over the concentration range limits their use in risk and benefits assessment, as these assessments typically require a differentiable algebraic function in addition to a quantitative estimate of uncertainty by concentration.

We developed and applied a new method that combines the flexibility of splines and the ease of use of the SCHIF in benefits assessment. Our method involves three steps. The first step is a data reduction step in which we fit a RCS with a large number of knots in order to characterize the shape of the concentration–response relationship in sufficient detail. RCS can easily be fit to large cohorts, as they only involve a series of transformations of concentration. Here, we have converted millions of person-years of data into a few hundred observations of RCS predictions over the observed concentration range. In step 2, we smooth the potential erratic predictions due to the large number of knots using a MISS, and in step 3, we fit our SCHIF function to the MISS predictions. In addition, we model the uncertainty in the spline fit as a cubic polynomial in concentration in a manner that assigns all uncertainty to theparameter in the SCHIF model, but unlike the LL model, uncertainty can vary with concentration. We now have a differentiable algebraic function of both relative risk and its uncertainty by concentration. This approach also allows for visualization of the SCHIF as well as its representation of the underlying data (as summarized by the RCS). θ

Specifically, we selected 15 knots defined at the 2nd, 4th, 10th, 14th, 18th, 22nd, 26th, 50th, 74th, 78th, 82nd, 86th, 90th, 94th, and 98th percentiles of theperson-year distribution. We selected a large number of knots covering both the lower and upper quartiles in order to capture a variety of desired shapes. From this first step, we obtain estimates of the logarithm of the RCS hazard ratio (logRCS) and the associated standard error at several hundred concentrations between the minimum and the 99th percentile of the exposure distribution. We do not include predictions above the 99th percentile, since RCS are linear beyond the highest knot concentration. This linear form can have some influence on the shape of the SCHIF throughout the concentration range, and especially over the higher concentrations, since the SCHIF is a single algebraic function. We also fixed the logRCS to zero at the minimum concentration, and its associated standard error was also set to zero. PM 2.5

In step 2, we smooth the potentially erratic logRCS predictions with a MISS in order to obtain predictions suitable to model with the SCHIF algebraic function, which itself is monotonically increasing (Pya and Wood 2015). The SCHIF hazard ratio function has the form: logSCHIF(z)=θf(z)l(z)with l(z)=1[1+exp (−z−μτr)] a logistic function in concentration. Here, θ,μ, and τ are unknown parameters to be estimated from the data, r is the range in the translated exposure, and z=PM2.5−min(PM2.5) such that logSCHIF(0)=0. The function f(z) can take two forms: f(z)=z (linear) and f(z)=log(z+1) (log). We have constructed the SCHIF to be similar to the LL model, logLL(z)=βz, by writing: logSCHIF(z)=θT(z), where T(z)=f(z)l(z) is a specific transformation of concentration.

The linear formcan model both linear and sublinear associations, while the log formcan model supralinear associations with mortality. Both forms can accommodate S-shaped functions through. Sets of valuesare selected that define the shape of. Larger values ofresult in larger ranges of concentration for which a sublinear association is modeled at lower concentrations due to the property of the logistic function. Larger values ofgenerate shapes forwith less curvature. By limiting the ranges for, we can limit the amount of curvature in the SCHIF. f z ( ) z = f l o g ( ) z ( ) z + 1 = l ( ) z ( ) μ τ , l ( ) z μ τ l ( ) z ( ) μ τ ,

A linear regression model was constructed using each transformation as the single predictor and the MISS prediction as the response. Using the MISS predictions, we were then able to select a wide range of values of the parameters to examine a wide variety of shapes that is not possible by modeling the subject-level cohort data. We selected values of μ ranging from 0 to r by integers, and τ ranging from 0.1 to 1 by 0.1 increments. For each set of parameters and the two forms of f(z), we obtained an estimate of θ and its standard error. We then created a single SCHIF curve by a weighted average of all the SCHIF curves examined, with weights determined by the fit of each curve on the MISS values. However, as the model averaged predictions at each concentration are themselves a potentially complicated function, these predictions can be summarized as a single algebraic function. Specifically, we fit a generalization of the SCHIF model logSCHIF(z)=θlog(zα+1)1+exp[−(z−μ)/v]to the mean SCHIF predicted curve over the concentration range. We added an additional parameter α to model the combination of the linear and log forms of f(z) used in the fitting step. The function log(zα+1) is nearly linear in z for large values of α. We collapse the product τr into a single parameter v to simplify the reporting of the parameter estimates.

In the LL model, all uncertainty in the hazard ratio is assigned to the single unknown parameter,. We aim to make a similar characterization of uncertainty in the SCHIF predictions, where all the uncertainty is ascribed to the parameter. We do this by considering a model of the standard error in the RCS predictions. However, unlike the LL model, RCS standard errors can vary in a nonlinear manner with concentration. We therefore consider a model for the standard error as a function of concentration of the form:, withdenoting our standard error model of, dependent on concentration. We select a general model that can accommodate a variety of shapes such as a cubic polynomial with the form:. β θ s s T e R C S e θ ( ) z ( ) z ( ) z = × s e θ ( ) z θ s z e θ σ 0 σ 1 σ 2 σ 3 ( ) z = + + + z 2 z 3

Finally, we construct pooled SCHIF models among the three cohorts in the following manner: Let vc(z) be the variance of the logarithm of the SCHIF prediction, logSCHIFc(z), at concentration z for cohort c=1,2,3. We construct a meta-analytic summary of the SCHIF predictions among the three cohorts as: logSCHIFPooled(z)=∑c=13wclogSCHIFc(z),where wc(z)=[1/vc(z)]/∑c=131/vc(z). For the variance of logSCHIFPooled(z), we include the variation in predictions among the cohorts in addition to the sampling uncertainty for each cohort as: ∑c=13wc2(z){vc(z)+[logSCHIFc(z)−logSCHIFPooled(z)]2}.

In order to obtain an algebraic function for the pooled SCHIF, we used nonlinear regression to estimate the SCHIF parameters, withdefining the data for the regression. We also modeled the standard error of the pooled SCHIF in a manner similar to that for each cohort separately. The variance of the pooled SCHIFs is a function of both the variance of each cohort-specific SCHIF prediction and the squared difference between the cohort-specific SCHIF predictions and the pooled SCHIF prediction. This latter term captures the uncertainty in both the shape and magnitude of the hazard ratio predictions among the three cohorts. l o g S C H I F P o o l e d ( ) z

Our study has a number of strengths and limitations. These analyses include a large number of deaths (nearly 1.5 million); population representativeness; a number of subject-level mortality predictors, such as occupation, income, and education, that are often not available in such large administrative-based cohorts (Di et al. 2017); and annual location information through linkage to income tax filings back to 1981. Our temporally resolved air pollution estimates could then be assigned to each year of follow-up. We also developed a method to combine the flexibility of splines to characterize the shape of the relationship between PM2.5 and mortality with the usefulness of the SCHIF in risk and benefits assessments. This new approach is feasible with respect to computing resources for very large cohort studies.

Despite these important strengths, the study includes a number of limitations. Specifically, the cohorts lack individual information on behavioral mortality risk factors such as smoking, diet, and obesity. We have previously examined the influence of further adjustment for these missing risk factors using statistical methods of indirect adjustment (Crouse et al. 2015; Erickson et al. 2019; Shin et al. 2014) and direct adjustment in a cohort based on annual health surveys in Canada where such information was available (Pinault et al. 2016b). For both of these approaches, we found that the PM2.5–mortality association remained positive. These observations indicate that it is unlikely that the positive associations we observed are entirely due to lack of inclusion of missing risk factors.

We captured information on social and economic status of each subject only at cohort inception. Thus, information on marital status, income, education, occupation, and employment may have changed over time, and we were not able to model any potential influence of such changes on the–mortality association. The 1991 cohort may have been most influenced by this lack of temporal adjustment given its 25 y of follow-up. PM 2.5

As we examined a 3-y moving-averageexposure window lagged by 1 y in all analyses, for the 2001 cohort, the initial exposure window assigned to the 2001 follow-up year was a 1998–2000 average and thus did not include backcasted exposures. The 1991 cohort used more years of backcasted exposures (1988–1997) than the 1996 cohort (1993–1997). Backcasted predictions may have induced greater exposure error in the two earlier cohorts. PM 2.5

We employed three very different exposure models for,, and. Theexposure model included remote sensing information coupled with CTM predictions. Land-use information informed any bias in the CTM predictions for the chemical components of total mass, but were not included as direct predictors of ground measurement data. Themodel used both land-use and remote sensing information, while themodel fused ground data with CTM predictions. The spatial resolution of the models was also different, withhaving the finest resolution ofandat, while themodel was at a resolution of. It was therefore difficult to directly compare the associations with mortality between pollutants and even more difficult with multiple-pollutant models subject to potential differential exposure error. However, the observed stronger associations withmay be due to a lower exposure error in both the original surfaces and the temporal adjustments, sinceis known to be more spatially homogenous than, whileis known to have the largest spatial variation. Additionally, although the quantitative estimates presented here reflect the large and diverse geographic area of Canada, they may not apply to places around the world with notably different sources and compositions of, for example, in areas where a much higher proportion of theis from dust and sand. PM 2.5 NO 2 O 3 PM 2.5 NO 2 O 3 NO 2 100 m PM 2.5 1 km O 3 21 km O 3 O 3 PM 2.5 NO 2 PM 2.5 PM 2.5

In summary, we found positive associations betweenand mortality in all three cohorts, and a similar supralinear concentration–response relationship. Lower uncertainty bounds were elevated above unity for all concentrations in the 1991 cohort, abovein the 1996 cohort, and abovein the 2001 cohort, suggesting our confidence in identifying concentrations where there exists a positive response is declining as concentrations decline over time. Therefore, interest exists in examining subsequent CanCHEC cohorts as they become available with the administration of more recent long form censuses. PM 2.5 2 μ g m / 3 5 μ g m / 3