Spatial Modeling of PM10 and NO2 in the Continental United States, 1985–2000

Background Epidemiologic studies of air pollution have demonstrated a link between long-term air pollution exposures and mortality. However, many have been limited to city-specific average pollution measures or spatial or land-use regression exposure models in small geographic areas. Objectives Our objective was to develop nationwide models of annual exposure to particulate matter < 10 μm in diameter (PM10) and nitrogen dioxide during 1985–2000. Methods We used generalized additive models (GAMs) to predict annual levels of the pollutants using smooth spatial surfaces of available monitoring data and geographic information system–derived covariates. Model performance was determined using a cross-validation (CV) procedure with 10% of the data. We also compared the results of these models with a commonly used spatial interpolation, inverse distance weighting. Results For PM10, distance to road, elevation, proportion of low-intensity residential, high-intensity residential, and industrial, commercial, or transportation land use within 1 km were all statistically significant predictors of measured PM10 (model R2 = 0.49, CV R2 = 0.55). Distance to road, population density, elevation, land use, and distance to and emissions of the nearest nitrogen oxides–emitting power plant were all statistically significant predictors of measured NO2 (model R2 = 0.88, CV R2 = 0.90). The GAMs performed better overall than the inverse distance models, with higher CV R2 and higher precision. Conclusions These models provide reasonably accurate and unbiased estimates of annual exposures for PM10 and NO2. This approach provides the spatial and temporal variability necessary to describe exposure in studies assessing the health effects of chronic air pollution.

Research Acute exposures to particulate and gaseous air pollutants have been associated with mor bidity and mortality in a large number of timeseries studies [Pope and Dockery 2006;U.S. Environmental Protection Agency (EPA) 1993. There are fewer cohort studies where it has been possible to examine the asso ciation of longterm exposures and mortality (Dockery et al. 1993;Finkelstein et al. 2003;Hoek et al. 2002;Jerrett et al. 2005bJerrett et al. , 2005cLaden et al. 2006;Lipfert et al. 2006;Miller et al. 2007;Nafstad et al. 2004;Nyberg et al. 2000;Pope et al. 1995Pope et al. , 2004Rosenlund et al. 2006). In most longterm studies, exposure assessment has been limited mainly to city specific average pollution measures or spatial or geographic information system (GIS)-based exposure models in small geographic areas (Adar and Kaufman 2007;Brauer et al. 2003;Briggs et al. 2000;Jerrett et al. 2005a;Liao et al. 2006;Su et al. 2008;Wheeler et al. 2008;Wong et al. 2004). One recent study has described a monthly spatiotemporal exposure model for the northeastern United States using a com bination of spatial and GISderived covari ates that outperformed models with spatial smoothing alone (Yanosky et al. 2008(Yanosky et al. , 2009). Another recent report has detailed the use of universal kriging to predict pollution levels for the European Union (Beelen et al. 2009). The purpose of this analysis is to develop nation wide models of annual exposure to particulate matter < 10 µm in diameter (PM 10 ) and nitro gen dioxide, using a combination of spatial smoothing and regression of GISderived cova riates. To date, few countrywide models have been available for these pollutants over our time scale of interest (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000). We apply the model to the addresses of the workers in the Trucking Industry Particle Study (Garshick et al. 2008;Laden et al. 2007), a retrospective cohort study of male U.S. unionized trucking company workers, to illustrate its potential use in exposure assessment for longterm epide miologic studies with members spread over the continental United States.

Methods
The Trucking Industry Particle Study. Details of the Trucking Industry Particle Study (TrIPS) are provided elsewhere (Garshick et al. 2008;Laden et al. 2007). Briefly, using personnel records from four large companies we identified 54,973 males with at least 1 day of work in 1985. Information was avail able on demographic variables, daily job and work location, and residential home address. Using an outside vendor (TeleAtlas, Lebanon, NH), we geocoded the last known residential addresses of 53,822 members living within the continental United States to at least the ZIP code level.
Pollutant data. We obtained information on annual average PM 10 (parameter codes 81102 and 85101) and NO 2 from the U.S. EPA Air Quality System (AQS). The U.S. EPA provided these annual averages on a set of DVDs compiled in 2004 for U.S. EPA Science to Achieve Results program grant 830545010. Data from 1985-2000 were used for this study if an annual mean was reported, regardless of the primary monitor ing objective of the monitor. All monitors in the continental United States were included, because excluding monitors such as those located near point or mobile sources would prevent us from incorporating all sources of spatial variability represented in the monitor ing network. Latitude and longitude of each monitor were obtained from the AQS data base and used to map the monitor locations using ArcGIS (version 9.2; ESRI, Redlands, CA). All monitors were checked for latitude/ longitude accuracy and precision to the county level before inclusion.
Modeling approach. We used generalized additive models (GAMs) to predict annual outdoor levels of PM 10 and NO 2 using smooth spatial surfaces and GISderived covariates. GAMs use semiparametric methods to model nonlinear, onedimensional, and multi dimensional functions using penalized splines (Hastie and Tibshirani 1990;Wood 2003Wood , 2004Wood , 2006. For both pollutants, models were constructed using 90% of the available monitoring locations for each calendar year. The remaining randomly selected 10% of monitors were used to perform crossvalida tion as described below.
First, the average spatial surface for each pollutant, 1985-2000, was generated in a GAM containing a bivariate thinplate spline of the projected x and ycoordinates of the monitoring locations and indicator variables for calendar year to adjust for temporal trends (Wood 2006). To obtain information on fine scale longterm spatial patterns, we included onedimensional penalized splines for a priori selected GISderived timeinvariant covariates. The covariates we considered included distance to road, population density, elevation, sur rounding land use, distance to and emission from power plants, and variables for census region of the country (northeast, west, south, and midwest) to adjust for regional patterns. These variables have previously been shown to be important predictors of ambient pollution (Adar and Kaufman 2007;Jerrett et al. 2005a;Yanosky et al. 2008Yanosky et al. , 2009. Each characteristic was assigned to the monitoring locations using ArcGIS. Information from the StreetMap data set (ESRI) was used to determine distance to the nearest road. Road segments were first classi fied by U.S. Census Feature Class Code as A1 (primary roads, typically interstates, with lim ited access), A2 (primary major, noninterstate roads), or A3 (smaller, secondary roads, usu ally with more than two lanes) (U.S. Census Bureau 1993). The distance from each location to the nearest road of each road class was then calculated in meters. Land use data were com piled from the U.S. Geological Survey (USGS) 1992 National Land Cover Dataset (USGS 2007b), which provides data on 19 categories of land use in raster image files with 1 arcsec (about 30 m) spatial resolution (Vogelmann et al. 2001). The proportion of lowintensity residential, highintensity residential, and industrial/commercial/transportation land uses within 1 km of each location was calcu lated. Population density values were assigned to each monitoring location using data from the 2000 U.S. Census at the block group level (U.S. Census Bureau 1993). Elevation data for each location were compiled from the USGS National Elevation Dataset (USGS 2007a). Information on the tons of nitrogen oxides emitted annually from all U.S. power plants in 2004 was obtained from the U.S. EPA 2006 Emissions and Generation Resource Integrated Database (U.S. EPA 2007a). The distance to and the emissions from the nearest facility were determined for each NO 2 monitoring location. Each potential covariate (or groups of covariates for distance to road, land use, and power plant distance/emissions) was first con sidered separately in models that included the bivariate spline for the 1985-2000 spa tial surface and the indicator variables for calendar year. We constructed multivariate models including all covariates that were sta tistically significant (p < 0.05) and led to a higher adjusted model R 2 . If covariates were no longer significant when included in the multivariate model, we omitted them unless they led to better model fit as determined by Akaike's information criterion (AIC) and crossvalidation testing.
To assess annual differences from the long term spatial patterns of pollution, we first cal culated the residuals from the final longterm multivariate GAM models. Then, for each calendar year, we created a bivariate smooth of the residuals using a twodimensional thin plate spline. Therefore, the annual average pollution at any location was predicted using the sum of the prediction from the longterm average surface/GISderived covariates and the prediction from the calendaryear specific residual spatial variability surface.
To perform crossvalidation, we used regression parameters from the final mod els and the annual spatial surfaces to predict annual pollutant levels at the 10% of moni toring locations that were held out from the original models. We assessed the potential bias of each final model by calculating the prediction error as the difference between the observed and predicted values at each cross validation monitoring location. We also assessed bias in the models by examining the intercept and slopes from linear regression of the predicted values on the measured val ues. The precision of the model was estimated by taking the square root of the mean of the squared prediction errors (RMSPE). In addi tion, a crossvalidation R 2 was obtained using the squared Pearson correlation between the measured values at the heldout observations and the model predictions.
For comparison, we also predicted expo sures using a simpler spatial interpolation method, inverse distance weighting (IDW), which had been frequently used in the air pollu tion literature. For the IDW models, the annual predictions for any given location (crossvalidation monitor location or cohort member address) were calculated by taking the average of the measured value at each moni tor location times the inverse of the squared distance between each location and each mon itor. IDW modeling was performed in ArcGIS (Johnston et al. 2004). The bias and precision of this simpler exposure modeling method was determined using crossvalidation.
After the final GAM models were deter mined and crossvalidated, the regression parameters were used to predict annual pollut ant levels at the 53,822 residential addresses of the TrIPS cohort members. For comparison, IDW was also used to predict annual pollut ant levels at the residential addresses. Statistical analyses were performed in PC SAS version 9.1 (SAS Institute Inc. 2006) and Unix R 2.7.0 (R Development Core Team 2006).

Results
The number of monitors used in the models and annual distributions of pollutant levels are shown in Table 1. The levels of both pol lutants decreased over time. The median value of PM 10 in 1985 was 38.2 µg/m 3 , and it fell to 23.0 µg/m 3 by 2000 (a 40% decrease). The median NO 2 level decreased 23% over the same period, from 19.0 ppb to 14.6 ppb. The distributions of the GISderived covariates at the monitor locations considered in the GAM  Table 2. The covari ate distributions were quite similar for both sets of monitors. As shown in Figure 1, the cohort participants are located throughout the continental U.S., and most live close to the monitoring locations. Specifically, the cohort members lived a median distance of 10.2 km from PM 10 monitoring sites and 16.6 km from NO 2 sites. Seventyfive percent of the cohort was no more than 21.1 km from a PM 10 moni tor included in the model and 35.6 km from an NO 2 monitor included in the model. PM 10 . The model with only the spatial spline and calendar year indicator variables had a model R 2 of 0.48. Region of the country, dis tance to all three census classes of road, block group population density, elevation, propor tion of lowintensity residential, highintensity residential, and industrial, commercial, or transportation land use within 1 km were all statistically significant independent predictors of measured PM 10 concentrations in univariate models. In a multivariate model, all predictors except population density (p = 0.15) remained statistically signifi cant predictors of meas ured PM 10 annual concentrations (Table 3). Population density was removed from the final model, because it did not increase the cross validation R 2 or model fit as determined by AIC. The final model had an R 2 of 0.49.
Increases in the proportion of surrounding land use used for highintensity residential or for industrial, commercial, or transporta tion uses were associated with increases in measured PM 10 levels. Increases in all other covariates were associated with decreases in measured PM 10 . The cross validation R 2 of the final model was 0.55. The median [and  Abbreviations: ICT, percentage of land used for industrial, commercial, or transportation; IQR, interquartile range; differ ence between the 75th and 25th percentile. Population density excluded from final PM 10 model. a All models also include indicator variables for region of the country. b R does not provide exact pvalues for those < 2 × 10 -16 . c Regression slope is linear regression of observed measurements at the holdout locations on model predictions at those locations. d Negative or positive. e Distance to and NO x from the nearest power plant were not considered for PM 10 .
interquartile range (IQR)] prediction error of the final model was 0.24 (7.0) µg/m 3 . The intercept and slope from the regression of observed and predicted meas urements were 1.49 and 0.94, respectively, and the RMSPE was 9.1 µg/m 3 . A plot of the observed versus expected values from the crossvalidation is presented in Supplemental Material, avail able online (doi:10.1289/ehp.0900840.S1 via http://dx.doi.org/). NO 2 . The model with only the spatial spline and calendar year indicators had a model R 2 of 0.73. Region of the country, distance to road, block group population den sity, elevation, surrounding land use, distance to nearest NO x emitting power plant, and the level of emissions from that power plant were all statistically significant predictors of measured NO 2 concentrations in univariate models. In a multivariate model, all predic tors remained statistically significant predic tors of measured NO 2 annual concentrations ( Table 3). The final multivariate model had an R 2 of 0.88. Increases in the block group popu lation density, NO x emissions of the nearest power plant, and the proportion of surround ing land use used for low or highintensity residential or for industrial, commercial, or transportation uses were associated with increases in measured NO 2 levels. Increases in all other covariates were associated with decreases in measured NO 2 . The cross vali dation R 2 of the final model was 0.90. The median (and IQR) prediction error of the final model was 0.10 (3.7) ppb, the inter cept and slope of the regression of observed and predicted measurements were 0.00 and 1.04, and the RMSPE was 3.5 ppb. A plot of the observed versus expected values from the crossvalidation is presented in Supplemental Material (doi:10.1289/ehp.0900840.S1).
Comparison with IDW. A summary of the crossvalidation parameters for the IDW exposure models is presented in Table 4. For both pollutants, the crossvalidation R 2 of the IDW model (R 2 = 0.44 for PM 10 and 0.67 for NO 2 ) was lower than those from the GAMs (R 2 = 0.55 for PM 10 and 0.90 for NO 2 ). For PM 10 , the slope from regression for the IDW model was 0.76 and the slope for the GAM was 0.94, indicating greater accuracy. The median prediction error for the IDW model was almost half that of the GAM, also indicating greater accuracy, but the RMSPE was higher, indicating lower pre cision. In contrast, for NO 2 the IDW predic tion error was 10fold higher than the GAM, and the RMSPE was almost twice as large.
TrIPS cohort exposures. The distribution of the GISderived variables for the residential addresses (n = 53,822) of the TrIPS cohort is presented in Table 5. The home addresses tended to be further away, on average, from each of the census road classes and from power plants than the monitors used to develop the models. The addresses were also located in areas with a lower proportion of highintensity residential or industrial, commercial, or trans portation land use, and the addresses were located further away from power plants than monitors, with lower annual emissions of NO x from the nearest plant, on average. The distri butions of the covariates tended to be tighter than those of the monitoring locations but were not significantly different. Figure 2 shows the distribution of the pollution values for each year at the cohort addresses. The mean predicted levels of the two pollutants decreased over the followup period, although there was little change in the overall spread of the distributions. The spatial   Predicted NO 2 (ppb) 1985 1987 1989 1991 1993 Year 1995 1997 1999 1985 1987 1989 1991 1993 Year 1995 1997 1999 5th 25th 50th 75th 95th volume 117 | number 11 | November 2009 • Environmental Health Perspectives distributions of the predictions for both PM 10 and NO 2 are shown in Figure 3. At all three time points shown, PM 10 values are higher in the western half of the United States than in the east. For NO 2 , however, the levels in all time periods are highest in major cities. To compare the two prediction methods, Figure 4 shows the cohort predictions for PM 10 at base line (1985), midpoint (1993), and last year of followup (2000). There is moderate cor relation between the results of the GAM and IDW PM 10 models, although the IDW mod els tend to be lower than the predictions of the GAMs (thus their lower slope of 0.76 vs. 0.94 for the GAM when both are compared with measured concentrations). The Spearman correlations between the two prediction types were 0.66 for 1985, 0.64 for 1993, and 0.77 for 2000. As shown in Figure 4, there is also moderate correlation between the GAM and IDW NO 2 models. Specifically, the Spearman correlation is 0.63 for 1985, 0.53 for 1993, and 0.51 for 2000. Overall, the IDW models tend to be lower than the GAM predictions and tend to have less variance (heterogeneity).

Discussion
Our results show that GAMs with a combina tion of spatial smoothing and GISderived covariates are a practical method for predict ing annual outdoor air pollution values for a cohort dispersed across the continental United States. The PM 10 and NO 2 GAM models were reasonably accurate and precise. The final model for NO 2 had a model R 2 of 0.88 and a crossvalidation R 2 of 0.90, whereas the final model R 2 for PM 10 was 0.49 and the cross validation R 2 was 0.55. Overall, the GAMs for both PM 10 and NO 2 outperformed the simpler IDW models, although there was a greater difference in the performance of the two modeling approaches for NO 2 . As expected, based on the growing lit erature of landuse regression models, many GISderived predictors were important in the pollution models. Distance to the nearest road of each road class, distance to and emissions from the nearest power plant, and landuse terms defining the surrounding area, variables previously shown to represent major sources of ambient NO 2 in the United States (U.S. EPA 2007b), were all statistically significant predic tors of NO 2 . In PM 10 models, distance to the nearest road of each road class was the most important class of predictors, likely representing traffic, an important local source of particulate matter (U.S. EPA 2004). These covariates did not improve the model R 2 as much for PM 10 as they did for NO 2 . It is possible that there are other important sources of PM 10 that we have not included (e.g., sea salt, crustal materials) that would improve the model R 2 more.
A growing number of studies have used spatial smoothing methods or models based on GISderived variables to predict ambient air pollution levels for use in epidemiologic studies (Adar and Kaufman 2007;Jerrett et al. 2005a;. Many of these studies have relied on proximity to spe cific pollution sources or monitoring locations to assign exposures. Others have focused on characterizing pollution from a specific source, typically onroad vehicles (Hoek et al. 2001). The most commonly used GISbased methods have used information on traffic volume and distance to roadways as surrogates of expo sure (Adar and Kaufman 2007;BayerOglesby et al. 2006;Forastiere and Galassi 2005;Garshick et al. 2003;Kan et al. 2007;Nitta et al. 1993;Oosterlee et al. 1996;Venn et al. 2005). In many of these studies, distance to road is divided into categories, or individuals are classified as exposed or not exposed, based on an a priori chosen distance. This method likely leads to exposure misclassification in many of these studies and is likely also quite sensitive to the buffer or category size selected.
Another popular GISbased exposure method is land use regression (Briggs et al. 1997;Hoek et al. 2001;Su et al. 2008). This approach is typically used in smaller areas to model local spatial variability, and roadway networks and traffic are often inputs to these mod els, although some also include information on surrounding land use, meteorology, and ambient air pollution monitoring locations. Other studies have used spatial smoothing techniques of the ambient measurements in single cities or counties (Jerrett et al. 2005b;Meng et al. 2007). Although direct compari sons are not appropriate, our NO 2 model R 2 of 0.88 is higher than those observed in many landuse regression models (0.52-0.76) (Briggs et al. 2000;Cyrys et al. 2005;Gilbert et al. 2005;Rosenlund et al. 2008) or in an EUwide model based on ordinary kriging (Beelen et al. 2009).
On a larger spatial scale, in an exposure assessment for the Women's Health Initiative, kriging in ArcGIS was used to generate daily PM 2.5 and PM 10 estimates for the entire continental United States for the year 2000 (Liao et al. 2006;Szpiro et al. 2007). For PM 10 , the authors report a median predic tion error of 0.04 µg/m 3 and an RMSPE of 19.48 µg/m 3 . In a recent exposure assessment for the Nurses' Health Study, a combination of spatial smoothing and GISderived cova riates was used to produce monthly predic tions of PM 10 1988-2002 for residences in the northeastern United States (Yanosky et al. 2008). This model has a mean prediction error of -0.4 µg/m 3 and an RMSPE of 6.4 µg/m 3 across the entire region, with no discernable differences by state or level of urbanization. Our models are similar to this modeling approach: Both include spatial smoothing and GISbased covariates to generate predictions. The Yanosky model allows the generation of monthly estimates of PM 10 through a complex spatiotemporal model and allows the inclu sion of timevarying covariates and control for seasonality. In contrast, although the model presented here also uses spatial smoothing and GISbased covariates, it is more appropriate for annual means and is less computationally intensive. Therefore, for PM 10 , the amount of bias [measured by average (mean or median) prediction error] and precision (measured by RMSPE) in our final model are comparable to that of other studies in the United States.
Our exposure model has several impor tant limitations. We rely on air pollution data from existing networks that are not uniformly distributed across the continental United States. However, the measures of precision and accuracy determined by crossvalidation for the heldout monitoring locations indicated good predictive performance of the models. Additionally, most of the members of the specific cohort we are using in this analysis live close to monitoring locations, so the mis match between monitor and subject locations is unlikely to be a large source of error in expo sure for our chosen application. For studies where the cohort is located much further from monitoring locations, this would likely be a larger source of error. In focusing our model ing on annual means, we are likely missing important seasonal and temporal variability occurring within each year. In years with fewer monitoring locations, it is possible that our model is underpowered to detect annual dif ferences from the longterm spatial trends; however, in later years, only 20-40 degrees of freedom were needed to fit these surfaces, so this may not be a large issue. Our model also does not include information on time varying covariates (such as pointsource pol lution or weather, especially wind direction and speed, mixing height, and precipitation) or interactions between our chosen covariates and calendar year. It is likely that information on these factors would improve the predic tive ability of our model; however, it would require a different modeling approach than the one we have chosen. By treating population density, distance from road, and land use as time invariant, we are assuming that these did not vary during the study period. This is not likely to be true and will lead to increased error in areas with rapidly changing infrastructure during this time period. Finally, we are using a spatial smoothing model for the entire con tinental United States. It has been suggested that regional models may be more appropri ate for the continental United States (Szpiro et al. 2007); however, it has been shown that for daily predictions, regional models do not substantially outperform a single countrywide model (Liao et al. 2006). Our models are adjusted for region of the country (using indi cator variables), and although including region did improve the fit of the models, the regional terms themselves were not significant.

Conclusions
In conclusion, our air pollution exposure model combining spatial smoothing techniques and GISbased predictors is a useful way to provide estimates of U.S.wide annual exposures for PM 10 and NO 2 . These models can be used to produce reasonably accurate and precise mea sures of pollution at the residential addresses of participants in epidemiologic studies focus ing on the adverse effects of constituents of air pollution as far back as 1985.