Association of long-term PM2.5 exposure with mortality using different air pollution exposure models: impacts in rural and urban California.

Most PM2.5-associated mortality studies are not conducted in rural areas where mortality rates may differ when population characteristics, health care access, and PM2.5 composition differ. PM2.5-associated mortality was investigated in the elderly residing in rural-urban zip codes. Exposure (2000-2006) was estimated using different models and Poisson regression was performed using 2006 mortality data. PM2.5 models estimated comparable exposures, although subtle differences were observed in rate ratios (RR) within areas by health outcomes. Cardiovascular disease (CVD), ischemic heart disease (IHD), and cardiopulmonary disease (CPD), mortality was significantly associated with rural, urban, and statewide chronic PM2.5 exposures. We observed larger effect sizes in RRs for CVD, CPD, and all-cause (AC) with similar sizes for IHD mortality in rural areas compared to urban areas. PM2.5 was significantly associated with AC mortality in rural areas and statewide; however, in urban areas, only the most restrictive exposure model showed an association. Given the results seen, future mortality studies should consider adjusting for differences with rural-urban variables.

Introduction

Exposure to particulate matter (PM) air pollution is a leading risk factor for premature mortality globally (Lim et al. 2013). An extensive body of epidemiological research has established a strong association between chronic exposures to fine PM less than 2.5 μm in diameter (PM2.5) and cardiovascular and all-cause non-accidental mortality (Dockery et al. 1993; Pope et al. 2002; Pope 2004; Laden et al. 2006; Miller et al. 2007; Chen et al. 2008; Puett et al. 2009; Ostro et al. 2010; Crouse et al. 2012; Lepeule et al. 2012; Cesaroni et al. 2013; Hoek et al. 2013). In fact, the U.S. EPA (2009) Integrated Science Assessment (ISA) recently concluded that a causal relationship exists between long-term PM2.5 exposure and cardiovascular effects and mortality (U.S. EPA 2009). These epidemiological studies have been conducted mainly in urban areas with very few rural areas included since most regulatory air pollution monitors are sited in populated areas. However, the responses to PM2.5 exposure in rural areas may be different from those in urban areas. For example, there are differences in population characteristics (Hart et al. 2005; Johnson et al. 2005) and access to health care (Laditka et al. 2009; Probst et al. 2011) between rural and urban areas, which may result in increased susceptibility to PM2.5 impacts. These factors and others may drive a phenomenon described as the “non-metropolitan mortality penalty,” a recently observed pattern, in which mortality rates are higher in rural compared to urban areas (Cosby et al. 2008; Cossman et al. 2010).

In addition, differences exist between particle chemical and physical composition in rural vs. urban areas, which may be due to variations in sources and seasonality differences (Chow et al. 1994; Motallebi et al. 2003; Eiguren-Fernandez et al. 2004; Rinehart et al. 2006; Chen et al. 2007; Hu et al. 2014). Different methodologies used to estimate long-term PM2.5 exposure may also influence the statistical estimation of the health impacts because each model results in a unique pattern of non-differential misclassification bias inherent to its structure, especially in rural areas where air pollution monitors are sparser. For these reasons, a comparison is warranted of the long-term health impacts of ambient PM2.5 exposure in rural vs. urban areas.

The elderly is the target population in this study because the health outcome of interest is cardiovascular disease related, and advanced cardiovascular disease (CVD) is rare in populations younger than 65 years of age. Air pollution-related health studies that have the youngest populations tend to show little if any effect of PM2.5 exposure (McDonnell et al. 2000) compared to studies where the bulk of the population included is around 60–75 or 80 years of age at the time of death (Krewski et al. 2000, 2009; Eftim et al. 2008; Jerrett et al. 2011, 2013). The effect increases to a peak in the age range of 65–75 years, and then declines (Zeger et al. 2008). Moreover, the elderly is a subpopulation that has been shown to be sensitive to PM exposure (Simoni et al. 2015).

One air pollution cohort study of farmers and their spouses residing in Iowa and North Carolina has specifically examined mortality and chronic PM2.5 exposure in rural areas. An association was found between long-term exposure to PM2.5 and cardiovascular mortality only among males and was strongest for participants whose residential location did not change, and who therefore likely had the most precise exposure estimates. All-cause mortality was not associated with PM2.5 in the cohort as a whole. Also, an inverse relationship between all-cause mortality and PM2.5 exposure was seen among women (Weichenthal et al. 2014); however, this finding was not robust to sensitivity analyses.

In this study, we investigated impacts of long-term PM2.5 exposure in rural vs. urban areas in California on mortality from CVD, ischemic heart disease (IHD), cardiopulmonary disease (CPD), and on all-cause non-accidental (AC) mortality in an elderly population. We also assessed the effects of different methodologies used to estimate PM2.5 exposure.

Methods

Air pollution data

Ambient PM2.5 concentrations from 2000 to 2006 collected at 116 fixed monitoring sites were extracted from California’s National, State, and Local Air Monitoring Network (NAMS/SLAMS) and the Interagency Monitoring of Protected Visual Environments (IMPROVE) network (Figure 1). Seven-year (2000–2006) average concentrations of PM2.5 were calculated from monthly averages and were assigned to zip code centroids using three different exposure models: (1) the closest monitor, (2) inverse distance weighting (IDW), and (3) kriging.

Figure 1.

(Color online) California rural and urban zip codes are depicted with the locations for California’s PM2.5 fixed monitoring networks (National, State, and Local Air Monitoring Network and the IMPROVE network).

The closest monitor model assigned the seven-year average PM2.5 concentration of the closest-fixed monitor to zip code centroids as a proxy for participants’ exposures. This model was implemented using three different forms. First, the model was implemented with no restriction of the distance between the nearest fixed monitor and the zip code centroid (unrestricted). The next two forms of the model were implemented using a restriction of a radius around the fixed monitor of 50 and 10 km. Therefore, all of the zip code centroids located within these radii (either 50 or 10 km) were included in the analysis and were assigned the seven-year average PM2.5 concentration from the nearest fixed monitor.

The IDW model assigned seven-year average PM2.5 concentrations to zip code centroids by creating an interpolated pollution surface using the Geostatistical Analyst extension of ArcMap version 10.0 (ESRI, Redlands, CA). The IDW model was applied with a 50-km restriction for consistency with other studies and monitoring siting criteria (Ozkaynak et al. 2007; U.S. EPA 2009; Jerrett et al. 2011; Lipsett et al. 2011) and was set to include concentration of one monitor station regardless of the distance for areas that had no fixed monitor within a 50-km radius.

The universal kriging model, like the IDW model, assigned seven-year average PM2.5 concentrations to zip code centroids by creating an interpolated pollution surface. Universal kriging was used since in California, there is a prevailing wind trend (Chow et al. 1994; Motallebi et al. 2003; Eiguren-Fernandez et al. 2004; Rinehart et al. 2006) which can be captured by the polynomial deterministic function. Geostatistical Analyst extension of ArcMap version 10.0 (ESRI, Redlands, CA) was used to run the universal kriging model.

Study population

The study population included all those aged 65 and older (65 +) in California who died in 2006. Mortality data were obtained from the Death Statistical Master File from the California Department of Public Health. The following International Classification of Diseases (WHO 1980), 10th Revision (ICD-10) categories were used to define the mortality outcomes of interest: CVD (I00-I99), IHD (I20-I25), CPD (I00-I99 and J00-J98), and all-cause non-accidental (A00-R99, excluding V01-V99).

Zip codes were included in the analysis if they had a total of 30 or more inhabitants (Hogg et al. 1993) aged 65 +, for a total of 1535 zip codes in California. Total mortality for California at the zip code level was calculated for AC and cause-specific mortality by summing the death counts for the 65 + population. Several contextual variables of ethnicity (percentage white, black, and Hispanic) and socioeconomic status (SES) were examined as potential predictors, along with two composite SES variables (Table 1). The SES variables were selected based on the Public Health Disparities Geocoding Project Monograph guidelines for area-based (contextual) monitoring of socioeconomic inequalities in health (Krieger 1992; Krieger et al. 2002, 2003). These variables were derived from the 2000 Census for the entire population by zip code and were examined to determine their potential for confounding the PM2.5–mortality relationship. We also investigated to see if there were statistically significant differences in these variables (α = 0.05) between rural and urban California (Table 1).

Table 1.

Area-based demographic variables from the 2000 U.S. Census evaluated for the analysis of PM-associated mortality in rural vs. urban areas in California.

Demographic characteristic	California	Rural	Urban
Total N**	33,292,571	8,265,963	25,026,608
Single variables
% white population *	68.93 ± 21.07	77.83 ± 17.00	59.13 ± 20.76
% black population *	4.44 ± 8.64	1.83 ± 3.44	7.30 ± 11.31
% Hispanic population*	24.59 ± 22.44	20.52 ± 21.45	29.06 ± 22.67
% persons employed predominantly in working class-occupations i.e. nonsupervisory employees*	65.16 ± 14.85	68.04 ± 12.67	62.02 ± 16.35
% unemployment (civilian labor force aged 16 years or older)*	7.96 ± 6.24	8.84 ± 6.59	6.99 ± 5.67
% households with an income < 50 % of the U.S. median household incomes*	21.54 ± 12.18	23.47 ± 12.20	19.42 ± 11.82
% persons below the federally defined poverty level ($17,050)	14.11 ± 9.88	14.55 ± 9.67	13.62 ± 10.09
% persons aged ≥ 25 years with < high school education*	43.39 ± 18.75	45.62 ± 17.20	40.94 ± 20.04
% house households containing > 1 person per room*	49.29 ± 26.05	40.00 ± 19.66	59.48 ± 28.30
Composite variables
Factor pertaining to economic resources (Factor 1) a , c	100.00 ± 25.00	86.90 ± 16.12	114.37 ± 25.12
Townsend index b , c	0.31 ± 1.74	−0.01 ± 1.46	0.66 ± 1.93

Factor 1 contained nine variables, specifically working class, unemployment, low income, median family income, below poverty level, adults who rent their homes or apartments, adults that do not own cars, low education, and crowding.

Townsend index is a measure of economic and social deprivation consisting of a standardized Z score combining data on percent crowding, percent unemployment, percent of individuals who do not own cars, and percent of adults who rent their home or apartment.

Higher scores represent higher degrees of deprivation or economic disadvantage.

Statistically significant differences between rural and urban zip codes (p < 0.05).

Total census population from 1535 zip codes with a minimum of 30 inhabitants aged 65+.

Statistical analysis

T-tests were performed to compare estimated concentrations of PM2.5 and SES variables in rural vs. urban regions for three exposure models, i.e. closest monitor, IDW, and kriging. Poisson regression was performed to investigate relationships between mortality and chronic PM2.5 exposure using different estimation methods of exposure. The total 65 + 2006 estimated population at the zip code level was obtained from the Community Sourcebook America (ESRI 2006), and was used as the offset for the Poisson regression in the analysis. The potential confounding effects of the contextual variables were examined based on the change-in-estimate criterion using a cut-off value of 10 %. Of 11 potentially confounding SES variables, we found the percent unemployment of the civilian labor force aged 16 years or older (unemployment), and the percent of persons aged 25 and older with less than a high school education (low education), to be statistically significant predictors of mortality and therefore included them in the final regression models. To compare our results with those of other studies, rate ratios (RRs) for PM2.5 were scaled to increments of 10 μg/m3. We applied the U.S. Census Bureau’s definition of an urban census tract (i.e. having a population density equal to or greater than 500 people per square mile) to California’s zip codes, resulting in 732 urban and 803 rural zip codes (Figure 1). We conducted stratified analyses of rural vs. urban zip codes, including unemployment and low-education variables in the final models. For the statewide model, we added a rural–urban indicator variable (Crouse et al. 2012).

Results

Air pollution

This study was designed to investigate PM2.5-associated mortality in rural and urban locations. The performance of various exposure models used to estimate air pollutant concentrations was also investigated. The exposure models ranged from simple interpolation approaches, such as the closest monitor and IDW, to a more sophisticated geostatistical technique, universal kriging, which is purely spatial and accounts for patterns of space–time dependence. Depending on the exposure model used, the average concentrations of PM2.5 in urban areas ranged from 15.44 to 15.86 μg/m3 and in rural areas, the average concentrations of PM2.5 ranged from 10.16 to 10.64 μg/m3 (Table 2). PM2.5 concentrations were statistically significantly higher in urban areas than rural areas, with a mean difference of approximately 5 μg/m3 for all three exposure models. This difference was also observed for the closest monitor model restricted to a 10-km radius, which included fewer zip codes than the other models (Table 2).

Table 2.

Means and standard deviations of estimated long-term PM2.5 concentrations (2000–2006 average) by zip code for different exposure models in California.

Modeled (PM2.5 μg/m3)	California			Rural*			Urban*
	Zip code N (%)	Mean	Std Dev	Zip code N (%)	Mean	Std Dev	Zip code N (%)	Mean	Std Dev
Closest monitor	1535 (100)	12.68	5.26	803 (100)	10.16**	4.72	732 (100)	15.44**	4.35
50 km	1442 (94)	12.92	5.22	710 (88)	10.32**	4.74	732 (100)	15.44**	4.35
10 km	528 (34)	14.99	4.63	85 (11)	10.51**	3.85	443 (61)	15.86**	4.26
IDW	1535 (100)	12.94	4.70	803 (100)	10.64**	4.16	732 (100)	15.45**	3.90
Kriging	1535 (100)	13.06	5.01	803 (100)	10.58**	4.00	732 (100)	15.77**	4.58

Percentages are based on the total number of zip codes and statistically significant differences are seen between urban and rural using t-test (p < 0.05).

Differences were statistically significant for all models.

Population demographics

Selected single and composite contextual SES variables at the zip code level were analyzed to identify differences between these variables in urban vs. rural areas. Rural and urban locations had non-significant differences in the percentages of persons below the federally defined poverty level (Table 1). However, in rural locations, the percentages of the white population, those employed in working-class occupations, adults with less than a high school education, unemployed people, and lower income households were significantly higher when compared to urban locations. There were statistically significant higher percentages of blacks and Hispanics, and crowded households (more than one person per room) in urban compared to rural locations (Table 1). These observed differences between the rural and urban populations of California are comparable to differences reported for rural and urban populations in other parts of the U.S. (Probst et al. 2011).

Mortality totals by health outcome and by study area

The total study population who died in 2006 by studied area was: statewide N = 162,124, rural N = 46,753, and urban N = 115,371. The total study population who died in 2006 by health outcome and by study area was for CVD: statewide N = 69,436, rural N = 19,646, and urban N = 49,790; for IHD: statewide N = 37,781, rural N = 10,245, and urban N = 27,536; for CPD: statewide N = 89,219, rural N = 25,578, and urban N = 63,641; for AC mortality: statewide N = 161,535, rural N = 46,549, and urban N = 114,986.

Mortality RRs

The RRs for mortality in California in rural and in urban regions and statewide are shown in Figure 2(a)–(d). The magnitude of the associations between PM2.5 exposures and AC mortality and cause-specific mortality (i.e. CVD, IHD, and CPD) varied depending on the exposure model used to generate the PM2.5 estimates. From the stratified analysis, the adjusted RRs for mortality (per 10 μg/m3 increment of PM2.5) for rural areas varied from 1.09 to 1.24 for CVD, 1.16 to 1.30 for IHD, 1.08 to 1.21 for CPD, and 1.08 to 1.14 for AC. In urban areas, the adjusted RRs for mortality (per 10 μg/m3 increment of PM2.5) were somewhat lower and varied from 1.06 to 1.12 for CVD, 1.21 to 1.29 for IHD, 1.04 to 1.09 for CPD, and 0.99 to 1.02 for AC. The statewide mortality RRs after adjusting for unemployment, low education, and a rural–urban indicator ranged from 1.07 to 1.13 for CVD, 1.20 to 1.28 for IHD, 1.06 to 1.10 for CPD, and 1.01 to 1.03 for AC per 10 μg/m3 increase of PM2.5 concentration. In all three regions, i.e. rural and urban areas, and statewide, the relationships between chronic PM2.5 exposure (with all exposure models) and all mortality categories studied (CVD, IHD, CPD, and AC) were found to be statistically significant, with one exception: in urban areas, a significant association was found for AC mortality only with the closest monitor model restricted to a 10-km radius. Interestingly, the magnitude of association for all the models was higher for CVD, CPD, and AC mortality in rural California than in urban areas, but similar for IHD mortality in both rural and urban areas.

Figure 2.

(a–d) Adjusted RRs of PM2.5-associated specific-cause and all-cause mortality for all of California, and rural and urban areasa.

Discussion

This assessment found significant associations in the elderly between cause-specific mortality (CVD, IHD, and CPD) and long-term PM2.5 exposure within both rural and urban areas of California as well as statewide. The rural–urban differences that we observed for all outcomes, but less so for IHD mortality, could have been impacted by residual confounding from unmeasured factors, particularly within the urban areas. For example, it is possible that mortality in urban areas could be lower in part because of a higher quality of and better access to health care, including faster emergency response times (Stults et al. 1984; Grumbach et al. 2003; Eberhardt and Pamuk 2004; Vukmir 2004; Glover et al. 2004; Probst et al. 2004; Trivedi et al. 2013). Observed rural–urban disparities could also result from differences in individual-level health behaviors, such as physical activity, diet, or smoking (Glover et al. 2004; Patterson et al. 2004; Doescher et al. 2007; Trivedi et al. 2013). Individuals in rural and urban areas may also differ in their levels of knowledge of heart attacks and stroke (Swanoski et al. 2012) and in ways to prevent them. The accumulated health and SES disparities listed above may lead to the “non-metropolitan mortality penalty,” where higher death rates are observed in rural than in urban areas (Cosby et al. 2008; Cossman et al. 2010), hence masking the overall association between particulate pollution and mortality in urban areas. Some investigators are beginning to recognize the need to control for this penalty in epidemiological studies (Pope et al. 1995; Krewski et al. 2000; Crouse et al. 2012) by adjusting their models with a rural/urban factor. For IHD mortality, however, we found similar PM2.5-associated impacts in urban and rural regions. Hence, these unmeasured factors may not be as important for IHD mortality in either rural or urban regions. Another possibility is that the relatively smaller sample size for IHD mortality and the greater variability seen in the results may have affected our ability to detect a difference in IHD mortality between rural and urban areas.

The results for AC mortality were complex. The statewide RRs found in this assessment for AC mortality are similar to those in other studies, in that the AC mortality RRs are not as strong as those for IHD mortality (Jerrett et al. 2005; Ostro et al. 2010; Lipsett et al. 2011; Crouse et al. 2012). On the other hand, one California statewide study (Jerrett et al. 2011, 2013) found higher RRs for both AC and CVD mortality compared to those in this report, possibly because the investigators were able to control for individual as well as contextual factors. Still other studies (Lipsett et al. 2011) reported mixed findings (i.e. attenuated or no associations) in AC mortality, including the one study in rural areas in Iowa and North Carolina that did not find an AC mortality effect in male farmers, but did observe an association with PM2.5 exposure and CVD (Weichenthal et al. 2014). AC mortality includes causes that have not been linked to air pollution by plausible biological mechanisms, in contrast to CVD, CPD, or IHD, for which clear associations have been observed with chronic air pollution exposure (Pope et al. 2004; Krewski et al. 2009; Ostro et al. 2010; Lipsett et al. 2011; Crouse et al. 2012). In addition, where we observed the weakest associations for AC mortality in urban areas, subtle differences in exposure models were seen; in this case, an association with AC mortality was observed with the closest monitor model only when it was restricted to 10 km. However, significant results were found for AC mortality in rural areas with all exposure models.

Chemical and physical compositional differences have been reported between rural and urban PM, and it is possible that the long-term health impacts of ambient PM2.5 in rural and urban areas, such as the mortality effects seen in the current assessment, could be affected by these differences (Chow et al. 1994; Motallebi et al. 2003; Eiguren-Fernandez et al. 2004; Rinehart et al. 2006; Chen et al. 2007; Hu et al. 2014). Since the concentration of PM2.5 estimated for rural areas was consistently lower than that in urban areas, another possibility is that the rural–urban differences we observed may be due to a non-linear exposure–response curve for PM2.5. Various exposure–response models have been considered for relating PM2.5 to mortality (Daniels et al. 2000; Schwartz et al. 2008; Pope et al. 2009; Pope, Brook et al. 2011; Pope, Burnett et al. 2011; Cox 2012; Crouse et al. 2012), and a no-threshold model with log-linearity for PM-related mortality has been adopted by a number of investigators and policymakers (U.S. EPA 2010). However, heterogeneities in the shapes of model curves across study populations are still poorly understood (U.S. EPA 2010), particularly for different age groups, and in this regard, we cannot exclude the possibility that the rural–urban differences in PM-related mortality observed could be a result of a non-linear characteristic of exposure–response models to increasing PM concentrations.

The results of this assessment also seem to indicate that the PM2.5 concentrations estimated by simpler exposure models, such as the closest monitor and IDW, were not greatly different from the more complex universal kriging model for the statewide assessment, or for the rural and urban regions, but may have contributed to the slight differences observed in the RRs within the studied areas by health outcomes. This result may be due in part to monitor density heterogeneity shown in Figure 1 (Wong et al. 2004). Even though the density of PM2.5 monitors in California is greater than in other areas of the country, the network is not consistently distributed throughout the state, with the less-urbanized regions generally having lower monitor density. This inconsistent monitor density may influence how efficacious the IDW and kriging models are in predicting PM2.5 exposure in rural areas vs. in urban areas, as these models function best when the monitoring data adequately capture significant spatial differences. This also might be the plausible reason the IDW, which defaults to weighted distance from the available monitoring stations, provide a higher strength of association in urban than in rural area and vice versa for universal kriging. However, we feel there is no definite interpretation for the differences observed in these exposure models. Future research should focus on quantifying the variations in performance between different exposure models under different scenarios (e.g. rural vs. urban) to identify if certain models are better suited for use in specific environments, given the heterogeneity of the density of air pollution monitor networks.

As an ecological study, there are limitations to this assessment that should be considered when interpreting the results. Assigning estimated PM2.5 values to the zipcode level as a proxy measure for personal exposure likely resulted in non-differential misclassification of exposure, which can produce a bias toward null findings. However, considering the larger size of zip codes and lower monitoring station density in rural areas it is likely that non-differential misclassification bias occurred to a greater extent in the rural areas, so that the RR might have been biased toward the null more so than in the urban areas. Previous research (Lepeule et al. 2012) has shown that exposure estimates using residential address information may not reduce bias on the effects of exposure on health outcomes, but others have emphasized the need for exposure models based on residential locations (Su et al. 2009). An additional limitation is that the analyses relied on readily available data from the California death master file, and no information was available on potential individual-level confounders, such as previous occupational exposures, cigarette smoking histories, and noise exposure at the individual residence. In addition, air pollution from indoor sources, including wood smoke, which can be an indoor as well as an outdoor pollution source, was not accounted for in the study. Better methodology of estimating indoor and outdoor exposure is warranted in future research.

This study found significant associations between PM2.5 and mortality in the elderly and found that the impacts of PM2.5 on mortality were higher in rural areas for all health endpoints examined except for IHD. Though further research in this area would be helpful, significant differences were observed in the association between PM2.5 exposure and mortality in rural compared to urban environments, which suggests that residence in rural vs. urban locations should be considered in future epidemiological studies and that a rural–urban variable(s) should be used to adjust for these potential differences. Only subtle differences within the areas studied were seen in the RRs with the different models of exposure estimation used in this study. The results from this study provide evidence for the existence of rural–urban variability in the health impact of PM2.5 exposure and points to the continued need for regulations to improve air quality and health not only in urban, but also in rural areas.