Risk of Long-term Exposure to Fine Particulate Matter
The study cohort is a subset of the 1991–2001 Canadian census mortality follow-up study (Wilkins et al. 2008). Persons were eligible for the census mortality cohort if they were ≥ 25 years of age; were a usual resident of Canada on the census reference day (4 June 1991); were not a long-term resident of an institution such as a prison, hospital, or nursing home; and had been among the 20% of Canadian households (~ 3.6 million respondents) selected randomly for enumeration with the mandatory long-form questionnaire. Subjects in the census cohort were linked to the Canadian Mortality Database (Statistics Canada 2005b) from 4 June 1991 to 31 December 2001 using deterministic and probabilistic linkage methods (Wilkins et al. 2008). A random selection of 125,100 linked records were excluded from the cohort so that the final sample represented no more than 15% of the Canadian population, as stipulated in the record linkage protocol, which left ~ 2.7 million subjects in the cohort (i.e., ~ 76% of the 3.6 million respondents). The 1991–2001 Canadian census mortality follow-up study received approval by the Statistics Canada Policy Committee (reference no. 012-2001) after consultation with Statistics Canada Confidentiality and Legislation Committee, Data Access and Control Services Division, and the Federal Privacy Commissioner. This approval is equivalent to that of standard research ethics boards.
In the present study, we included only subjects who were nonimmigrants (i.e., only those granted Canadian citizenship by birth, which left ~ 2.1 million subjects for the analysis) because we were interested in capturing the exposure experience of longer-term residents of Canada. Immigrants to Canada have unknown prior exposures and are more likely to live in areas that are characterized by higher ground-level concentrations of PM2.5 (e.g., large cities such as Toronto, Montreal, and Vancouver) than those born in Canada (Villeneuve et al. 2011).
Immigrants, and especially recent immigrants, tend to have better health and health behaviors than the Canadian-born population (Ali et al. 2004; McDonald and Kennedy 2004) and to live longer than the nonimmigrant population (Wilkins et al. 2008). Because of more limited information on long-term exposure and because of immigrant-related traits that complicate the analysis, the immigrant subset of the population was excluded from our analysis and will be the subject of a separate, future analysis.
Mortality data included underlying cause of death [coded to the International Classification of Diseases, 9th Revision (ICD-9; WHO 1977), for deaths before 2000 and to ICD-10 (WHO 1992) for those deaths registered from 2000 onward] and date of death. Additionally, location of residence of each subject at baseline was aggregated to 1991 enumeration areas. Enumeration areas range in size from approximately 650 dwellings in urban areas to < 100 dwellings in rural areas. In 1991 there were 45,710 enumeration areas across the country.
Observations of concentrations of PM2.5 from ground-based stations were available for the follow-up period of the cohort in only 11 cities. We calculated the mean annual concentration in these cities averaged over the 1987–2001 period (i.e., 5 years before baseline and the full 10 years of follow-up) and assigned exposure to each cohort member living in the corresponding 11 census divisions [for details regarding the compilation of the historical ground-based observations, see Supplemental Material (http://dx.doi.org/10.1289/ehp.1104049)]. This subsample represented 43% of the full cohort.
A second set of exposure estimates was created for the full cohort using estimates of PM2.5 derived from satellite remote sensing observations during the period 2001–2006 [van Donkelaar et al. 2010; for additional descriptions of the satellite data retrievals, see Supplemental Material (http://dx.doi.org/10.1289/ehp.1104049)]. These estimates were available on a grid with a spatial resolution of 10 km × 10 km. Previous analyses showed that the satellite-derived estimates of PM2.5 were in close agreement (Pearson correlation coefficient r = 0.77, slope = 1.07, n = 1,057) with ground-based measurements both in Canada and in the United States (van Donkelaar et al. 2010). We overlaid the surface on a map layer representing the 1991 enumeration area boundaries and computed the mean concentration of PM2.5 within the boundaries of each enumeration area across the country. We thus assigned exposure to the satellite-derived estimates of PM2.5 to all cohort members by linking the exposure surface to their enumeration area of residence in 1991. The sparsely inhabited northern territories of Canada were excluded from our study because of the absence of estimates of PM2.5.
Geographical context can influence risk of mortality and may confound the association between mortality and air pollution because areas characterized by some population groups (e.g., those of lower income, the unemployed) and by some environmental characteristics (urban vs. rural) tend also to be characterized by higher ambient concentrations of air pollution (Crouse et al. 2009; Jerrett et al. 2001). To derive ecological covariates, we aggregated census data from 1991 describing socioeconomic and demographic characteristics to the smaller "neighborhood" scale of census tracts as well as to the larger "community" scale of census divisions, which may correspond to (or be of comparable size to) counties. To adjust for regional variations in these variables across Canada, we subtracted the census division mean from the values of each census tract. We thus compiled ecological variables describing the proportion of unemployed adults (≥ 15 years of age), the proportion of adults who had not completed high school, and the proportion of individuals in the lowest income quintile (as represented by the ratio of family income divided by the low-income cutoff) at both scales for each cohort member. The low-income cutoff, which varies by community and family size, is defined by Statistics Canada to identify those who need to spend a greater proportion of their income on basic necessities than does an average family of similar size and thus provides an indicator of deprivation that is adjusted for regional variation in the cost of living (Statistics Canada 2001). We also created a categorical variable that indicated the population size of the subject's home community given that those who live in rural areas tend to have poorer health than those who live in more urban areas (Eberhardt and Pamuk 2004).
We estimated hazard ratios (HRs) using a standard Cox proportional hazards model as well as with a nested, spatial random-effects Cox model (Krewski et al. 2009; Ma et al. 2003), in which random effects were represented by clusters, as defined below. The baseline hazard function for both models was stratified by single-year age groups and sex. Follow-up time was measured in days, calculated from 4 June 1991 to 31 December 2001.
The standard Cox survival model assumes that the survival time of each subject is statistically independent from that of other subjects after controlling for the mortality risk factors included in the model. These risk factors, however, may not completely explain differences in mortality between subjects and may introduce dependencies among subjects. One approach to explain some of these dependencies is to include information on the location of each subject, such as their community and neighborhood. Subjects living in the same community or neighborhood may be more likely to share mortality risk factors than subjects residing in different locations.
In our nested, spatial random-effects Cox model, therefore, we defined two levels of spatial clusters: a first cluster level defined by census divisions, and a second cluster level defined by census tracts within census divisions. We assumed that if two census divisions were adjacent, then they were correlated, and if not adjacent, uncorrelated. A similar assumption was made for the census tracts within each census division. Census tracts in different census divisions were assumed to be uncorrelated.
We examined the sensitivity of the estimates of the effects of air pollution on mortality due to the proportional hazard assumption using methods presented by Sheppard et al. (2011). We developed survival models for each year of follow-up and summarized the HRs for PM2.5 over the multiple years, and we found that the HR and its associated standard error was similar to those achieved by a single, multiyear model, which suggested that our estimate of the PM2.5 HR and the corresponding standard error was not sensitive to the proportional hazards assumption.
We also examined the sensitivity of the PM2.5 association with mortality to inclusion of selected sets of covariates: those measured at the subject level (i.e., those covariates listed in Table 1), an indictor of urban size, and the contextual covariates measured at the census division and census tract levels.
We examined the shape of the relationship between PM2.5 and mortality using natural cubic spine functions with one, two, three, or four degrees of freedom (df) for both the standard Cox survival model and the spatial random-effects model, and for several causes of death. We then examined plots of the concentration–response curves and used the Bayesian information criterion (BIC) to assess the relative goodness of fit for these models.
We developed models for mortality from ischemic heart disease (ICD-9: codes 410–414; ICD-10: I20–I25), cerebrovascular disease (ICD-9: 430–434, 436–438; ICD-10: I60–I69), cardiovascular disease (ICD-9: 410–417, 420–438, 440–449; ICD-10: I20–I28, I30–I52, I60–I79), circulatory disease (ICD-9: 390–459; ICD-10: I00–I99), and all nonaccidental causes (ICD-9 codes < 800; ICD-10 codes starting with letters A through R). HRs and 95% confidence intervals (CIs) were calculated for an increment of 10 μg/m in estimated concentrations of PM2.5.
Methods
The Study Cohort
The study cohort is a subset of the 1991–2001 Canadian census mortality follow-up study (Wilkins et al. 2008). Persons were eligible for the census mortality cohort if they were ≥ 25 years of age; were a usual resident of Canada on the census reference day (4 June 1991); were not a long-term resident of an institution such as a prison, hospital, or nursing home; and had been among the 20% of Canadian households (~ 3.6 million respondents) selected randomly for enumeration with the mandatory long-form questionnaire. Subjects in the census cohort were linked to the Canadian Mortality Database (Statistics Canada 2005b) from 4 June 1991 to 31 December 2001 using deterministic and probabilistic linkage methods (Wilkins et al. 2008). A random selection of 125,100 linked records were excluded from the cohort so that the final sample represented no more than 15% of the Canadian population, as stipulated in the record linkage protocol, which left ~ 2.7 million subjects in the cohort (i.e., ~ 76% of the 3.6 million respondents). The 1991–2001 Canadian census mortality follow-up study received approval by the Statistics Canada Policy Committee (reference no. 012-2001) after consultation with Statistics Canada Confidentiality and Legislation Committee, Data Access and Control Services Division, and the Federal Privacy Commissioner. This approval is equivalent to that of standard research ethics boards.
In the present study, we included only subjects who were nonimmigrants (i.e., only those granted Canadian citizenship by birth, which left ~ 2.1 million subjects for the analysis) because we were interested in capturing the exposure experience of longer-term residents of Canada. Immigrants to Canada have unknown prior exposures and are more likely to live in areas that are characterized by higher ground-level concentrations of PM2.5 (e.g., large cities such as Toronto, Montreal, and Vancouver) than those born in Canada (Villeneuve et al. 2011).
Immigrants, and especially recent immigrants, tend to have better health and health behaviors than the Canadian-born population (Ali et al. 2004; McDonald and Kennedy 2004) and to live longer than the nonimmigrant population (Wilkins et al. 2008). Because of more limited information on long-term exposure and because of immigrant-related traits that complicate the analysis, the immigrant subset of the population was excluded from our analysis and will be the subject of a separate, future analysis.
Mortality data included underlying cause of death [coded to the International Classification of Diseases, 9th Revision (ICD-9; WHO 1977), for deaths before 2000 and to ICD-10 (WHO 1992) for those deaths registered from 2000 onward] and date of death. Additionally, location of residence of each subject at baseline was aggregated to 1991 enumeration areas. Enumeration areas range in size from approximately 650 dwellings in urban areas to < 100 dwellings in rural areas. In 1991 there were 45,710 enumeration areas across the country.
Assignment of Concentrations of PM2.5
Observations of concentrations of PM2.5 from ground-based stations were available for the follow-up period of the cohort in only 11 cities. We calculated the mean annual concentration in these cities averaged over the 1987–2001 period (i.e., 5 years before baseline and the full 10 years of follow-up) and assigned exposure to each cohort member living in the corresponding 11 census divisions [for details regarding the compilation of the historical ground-based observations, see Supplemental Material (http://dx.doi.org/10.1289/ehp.1104049)]. This subsample represented 43% of the full cohort.
A second set of exposure estimates was created for the full cohort using estimates of PM2.5 derived from satellite remote sensing observations during the period 2001–2006 [van Donkelaar et al. 2010; for additional descriptions of the satellite data retrievals, see Supplemental Material (http://dx.doi.org/10.1289/ehp.1104049)]. These estimates were available on a grid with a spatial resolution of 10 km × 10 km. Previous analyses showed that the satellite-derived estimates of PM2.5 were in close agreement (Pearson correlation coefficient r = 0.77, slope = 1.07, n = 1,057) with ground-based measurements both in Canada and in the United States (van Donkelaar et al. 2010). We overlaid the surface on a map layer representing the 1991 enumeration area boundaries and computed the mean concentration of PM2.5 within the boundaries of each enumeration area across the country. We thus assigned exposure to the satellite-derived estimates of PM2.5 to all cohort members by linking the exposure surface to their enumeration area of residence in 1991. The sparsely inhabited northern territories of Canada were excluded from our study because of the absence of estimates of PM2.5.
Contextual Variables
Geographical context can influence risk of mortality and may confound the association between mortality and air pollution because areas characterized by some population groups (e.g., those of lower income, the unemployed) and by some environmental characteristics (urban vs. rural) tend also to be characterized by higher ambient concentrations of air pollution (Crouse et al. 2009; Jerrett et al. 2001). To derive ecological covariates, we aggregated census data from 1991 describing socioeconomic and demographic characteristics to the smaller "neighborhood" scale of census tracts as well as to the larger "community" scale of census divisions, which may correspond to (or be of comparable size to) counties. To adjust for regional variations in these variables across Canada, we subtracted the census division mean from the values of each census tract. We thus compiled ecological variables describing the proportion of unemployed adults (≥ 15 years of age), the proportion of adults who had not completed high school, and the proportion of individuals in the lowest income quintile (as represented by the ratio of family income divided by the low-income cutoff) at both scales for each cohort member. The low-income cutoff, which varies by community and family size, is defined by Statistics Canada to identify those who need to spend a greater proportion of their income on basic necessities than does an average family of similar size and thus provides an indicator of deprivation that is adjusted for regional variation in the cost of living (Statistics Canada 2001). We also created a categorical variable that indicated the population size of the subject's home community given that those who live in rural areas tend to have poorer health than those who live in more urban areas (Eberhardt and Pamuk 2004).
Statistical Methods
We estimated hazard ratios (HRs) using a standard Cox proportional hazards model as well as with a nested, spatial random-effects Cox model (Krewski et al. 2009; Ma et al. 2003), in which random effects were represented by clusters, as defined below. The baseline hazard function for both models was stratified by single-year age groups and sex. Follow-up time was measured in days, calculated from 4 June 1991 to 31 December 2001.
The standard Cox survival model assumes that the survival time of each subject is statistically independent from that of other subjects after controlling for the mortality risk factors included in the model. These risk factors, however, may not completely explain differences in mortality between subjects and may introduce dependencies among subjects. One approach to explain some of these dependencies is to include information on the location of each subject, such as their community and neighborhood. Subjects living in the same community or neighborhood may be more likely to share mortality risk factors than subjects residing in different locations.
In our nested, spatial random-effects Cox model, therefore, we defined two levels of spatial clusters: a first cluster level defined by census divisions, and a second cluster level defined by census tracts within census divisions. We assumed that if two census divisions were adjacent, then they were correlated, and if not adjacent, uncorrelated. A similar assumption was made for the census tracts within each census division. Census tracts in different census divisions were assumed to be uncorrelated.
We examined the sensitivity of the estimates of the effects of air pollution on mortality due to the proportional hazard assumption using methods presented by Sheppard et al. (2011). We developed survival models for each year of follow-up and summarized the HRs for PM2.5 over the multiple years, and we found that the HR and its associated standard error was similar to those achieved by a single, multiyear model, which suggested that our estimate of the PM2.5 HR and the corresponding standard error was not sensitive to the proportional hazards assumption.
We also examined the sensitivity of the PM2.5 association with mortality to inclusion of selected sets of covariates: those measured at the subject level (i.e., those covariates listed in Table 1), an indictor of urban size, and the contextual covariates measured at the census division and census tract levels.
We examined the shape of the relationship between PM2.5 and mortality using natural cubic spine functions with one, two, three, or four degrees of freedom (df) for both the standard Cox survival model and the spatial random-effects model, and for several causes of death. We then examined plots of the concentration–response curves and used the Bayesian information criterion (BIC) to assess the relative goodness of fit for these models.
We developed models for mortality from ischemic heart disease (ICD-9: codes 410–414; ICD-10: I20–I25), cerebrovascular disease (ICD-9: 430–434, 436–438; ICD-10: I60–I69), cardiovascular disease (ICD-9: 410–417, 420–438, 440–449; ICD-10: I20–I28, I30–I52, I60–I79), circulatory disease (ICD-9: 390–459; ICD-10: I00–I99), and all nonaccidental causes (ICD-9 codes < 800; ICD-10 codes starting with letters A through R). HRs and 95% confidence intervals (CIs) were calculated for an increment of 10 μg/m in estimated concentrations of PM2.5.
SHARE
