### 522 volume 117 | number 4 | April 2009 • Environmental Health Perspectives 

# Research

Chronic exposures to airborne particulate matter (PM) are related to increased mortal

ity as well as lung cancer, ischemic heart dis

ease, dysrhythmias, heart failure, and cardiac arrest (Abbey et al. 1999; Dockery et al. 1993; Finkelstein et al. 2003; Jerrett et al. 2005; Miller et al. 2007; Pope et al. 1995, 2002, 2004)

. These health effects have been shown for populations living near stationary ambient monitoring (SAM) sites and primarily and most consistently for exposure to PM with an aerodynamic diameter < 2.5 μm (PM 2.5 ). It is possible that PM 2.5 , which originates from pri

mary emissions from combustion sources and from secondary formation in the atmosphere, has different chronic health impacts than larger coarse particles (PM with aerodynamic diameters between 2.5 and 10 μm; PM 10–2.5 ) because of their different composition. In con

trast to PM 2.5 , PM 10–2.5 is typically generated from mechanical grinding or crushing, as well as from windblown dust. The chronic health effects of PM 10–2.5 , however, have been little studied (Brunekreef and Forsberg 2005). Exposure assessment in chronic PM studies has been severely limited by the lack of ambient monitoring data for PM 2.5 and especially PM 10–2.5 , which are sparse over both time and space, particularly for peri

ods before 1999. Because of this lack of data, some epidemiologic studies of PM 2.5 have restricted their health effect analyses to popu

lations near SAM sites (Dockery et al. 1993; Miller et al. 2007), whereas others have esti

mated exposures using concentrations aver

aged across metropolitan areas (Pope et al. 1995, 2002, 2004). Although this averaging may reduce instrument errors (and thereby classical exposure error), it ignores within-city spatial gradients in exposure, which can be substantial and may have important impli

cations for health. For example, in a study conducted in Los Angeles, California , Jerrett

et al. (2005) estimated chronic exposures

, using PM 2.5 data from one year (2000) and a spatial model that captured within-city spatial gradients

, and showed that within-city gra

dients in PM 2.5 were more strongly related to health effects than between-city gradients

. These results were consistent with those of Miller et al. (2007), who also found larger effects for within-city comparisons. It is possi

ble that health effects associated with chronic PM 2.5 exposure are even larger: Miller et al. (2007) and Jerrett et al. (2005) assumed that spatial gradients in PM 2.5 were constant over time, a simplifying assumption that may lead to bias (Haneuse et al. 2007). To address these limitations, we estimated monthly PM 2.5 and PM 10–2.5 concentrations using spatiotemporal models that incorporated both small

- and large- scale spatial trends, and that allowed these trends to change over time. These models were used to estimate PM 2.5 and PM 10–2.5 concentrations from 1988 through 2002 for populations living within the north

eastern and midwestern United States.


### Methods We predicted highly time- and space-resolved PM 2.5 concentrations across our study region using geographic information system (GIS)- based spatiotemporal models of monthly PM 2.5 concentrations. We did so for PM 2.5 by developing and validating two models: one for 1988

– 1998 (pre-1999) and another for 1999

– 2002 (post-1999). Both models build on the GIS-based spatiotemporal model we developed for PM with an aerodynamic diameter < 10 μm (PM 10 ) (Yanosky et al. 2008) and use publicly available air pollution, geographic, and meteorologic data

. The pre1999 model also used extinction coefficients derived from airport visibility data to predict location-specific outdoor PM 2.5 concentra

tions . For PM 10–2.5 , we estimated spatially resolved monthly values by difference using the PM 2.5 models and our previously devel

oped PM 10 model. Data. PM 2.5 and PM 10–2.5 data. Outdoor PM 2.5 concentration data from the north

eastern and midwestern United States (the study region) were included in this analy

sis (Figure 1). Data collected in states adja

cent to the study region were also included

Address correspondence to J.D. Yanosky, Harvard School of Public Health, 401 Park Dr., Boston, MA 02215 USA. Telephone: (617) 384-8836. Fax: (617) 384-8859. E-mail: jyanosky@hsph.harvard.edu

Supplemental Material is available online at http:// www.ehponline.org/members/2008/11692/suppl.pdf

This research has been funded wholly or in part by the U.S. Environmental Protection Agency (EPA) through grant 83054501-0 to Harvard University.

This research has not been subjected to the U.S. EPA’s required peer and policy review and therefore does not necessarily reflect the views of the agency, and no official endorsement should be inferred.

The authors declare they have no competing financial interests.

Received 12 May 2008; accepted 19 November 2008. 

## Predicting Chronic Fine and Coarse Particulate Exposures Using Spatiotemporal Models for the Northeastern and Midwestern United States

Jeff D. Yanosky, 1 Christopher J. Paciorek, 2 and Helen H. Suh 1

1 Exposure, Epidemiology and Risk Program, Department of Environmental Health, and 2 Department of Biostatistics, Harvard School of Public Health, Boston, Massachusetts, USA

B ackground : Chronic epidemiologic studies of particulate matter (PM) are limited by the lack of monitoring data, relying instead on citywide ambient concentrations to estimate exposures. This method ignores within-city spatial gradients and restricts studies to areas with nearby monitoring data. This lack of data is particularly restrictive for fine particles ( PM with aerodynamic diameter

< 2.5 μm; PM 2.5 ) and coarse particles ( PM with aerodynamic diameter 2.5–10 μm; PM 10–2.5 ), for which monitoring is limited before 1999. To address these limitations, we developed spatiotempo

ral models to predict monthly outdoor PM 2.5 and PM 10–2.5 concentrations for the northeastern and midwestern United States

. M ethods : For PM 2.5 , we developed models for two periods: 1988–1998 and 1999–2002. Both models included smooth spatial and regression terms of geographic information system –based and meteorologic predictors. To compensate for sparse monitoring data, the pre-1999 model also included predicted PM 10 (PM with aerodynamic diameter < 10 μm) and extinction coefficients (km –1 ). PM 10–2.5 levels were estimated as the difference in monthly predicted PM 10 and PM 2.5 , with predicted PM 10 from our previously developed PM 10 model.

R

esults : Predictive performance for PM 2.5 was strong (cross-validation R 2 = 0.77 and 0.69 for post-1999 and pre-1999 PM 2.5 models, respectively) with high precision (2.2 and 2.7 μg/m 3

, respectively). Models performed well irrespective of population density and season. Predictive per

formance for PM 10–2.5 was weaker (cross-validation R 2 = 0.39) with lower precision (5.5 μg/m 3 ). PM 10–2.5 levels exhibited greater local spatial variability than PM 10 or PM 2.5 , suggesting that PM 2.5 measurements at ambient monitoring sites are more representative for surrounding populations than for PM 10 and especially PM 10–2.5 . C onclusions : We provide semiempirical models to predict spatially and temporally resolved longterm average outdoor concentrations of PM 2.5 and PM 10–2.5 for estimating exposures of populations living in the northeastern and midwestern United States

. K

ey words : air pollution, extinction coefficient, fine particulate matter, geographic information system, generalized additive mixed models, geostatistics, spatial smoothing, spatiotemporal model

ing, visual range. Environ Health Perspect 117:522–529 (2009). doi:10.1289/ehp.11692 available
via http://dx.doi.org/ [Online 19 November 2008]Predicting chronic fine and coarse PM exposures

Environmental Health Perspectives • volume 117 | number 4 | April 2009 

### 523 in the models to minimize boundary effects (Figure 1). PM 2.5 data were obtained from the U.S. Environmental Protection Agency (EPA) Air Quality System (AQS) (U.S. EPA 2008) ; the Visibility Information Exchange Web System (2004) for Interagency Monitoring of Protected Visual Environments (IMPROVE), Stacked Filter Unit (a predecessor to IMPROVE), and Clean Air Status and Trends (CASTNet) networks ; and Harvard research studies including the Twenty-four Cities Study and Five Cities Study (Spengler et al. 1996; Suh et al. 1997). The U.S. EPA AQS network provided about 88% of the PM 2.5 concentration values and monitoring locations in the study region. Of the 546 monitoring locations, 386 were in the study region, with the remaining in adjacent states (Figure 1). Monthly average PM 2.5 concentrations were approximately log-normally distributed, with geometric mean and geometric standard deviation of 13.0 μg/m 3 and 1.5, respectively, across the study region. Measured PM 10–2.5 levels (estimated as the difference in measured monthly average PM 10 and PM 2.5 levels atco-located sites) were also approximately lognormally distributed, with geometric mean and geometric standard deviation of 6.9 μg/m 3 and 2.44, respectively. Geographic, meteorologic, and visibility data. C haracteristics of the PM monitoring sites were quantified using a GIS (ArcMap 9; ESRI, Redlands, CA ), including the following: •	Distance to nearest roadways for road classes

A1 (primary roads, typically interstates, with limited access), A2 (primary major, non

interstate roads), A3 (smaller, secondary roads, usually with more than two lanes), and A4 (two-lane, typically surface roads used for local traffic) using ESRI StreetMap data;

• Urban land use (the proportion of lowintensity residential, high-intensity residential, and industrial/commercial/transportation land uses within 1 km) using data from the U.S. Geologic al Survey (USGS) National Land Cover Dataset (2004; Homer et al. 2004); •	Block group, tract, and county population

density from the U.S. Census TIGER files from ESRI Data & Maps (ESRI); •	Point-source emissions of PM 2.5 within circu

lar buffers (1- and 10-km radii) and countylevel area-source PM 2.5 emissions from the U.S. EPA National Emissions Inventory (U.S. EPA 2005); •	Elevation from the USGS National Elevation

Dataset (USGS 2005).

Monthly average temperature, wind speed, sea level–adjusted barometric pressure , and total precipitation were obtained from the National Climatic Data Center (2005) and spatially smoothed to the monitoring locations for each month. Additional details regarding these geographic and meteorologic covariates can be found in Yanosky et al. (2008). Extinction coefficients were estimated using daily airport visual range observations (Faulke and Husar 1998; Ozkaynak et al. 1985) after correction for relative humidity and trunca

tion, spatially smoothed to monitoring and prediction locations, then averaged by month (Paciorek et al. 2008). Statistical models. To compensate for the spatially sparse PM 2.5 measurements before 1998, we constructed separate models for the pre-1999 and post-1999 time periods. Both models predicted PM 2.5 concentrations using generalized additive mixed models (GAMMs), with bivariate penalized spline terms for space and one-dimensional penalized spline terms for GIS-based and meteorologic predictors. Post-1999 PM 2.5 spatiotemporal model. The structure of the post-1999 PM 2.5 model followed that for the PM 10 model (Yanosky

et al. 2008):

y i
, t

= a + d 1 (

X

i ,1 )+...+ d Q (

X

i
, Q

)

+ g ( s i ) + f 1 ( Z i
, t

,1 )+...+ f P ( Z i
, t

, P

)

+ a t + g t ( s i ) + b i + e i
, t ; b i ~ N (0, σ 2 b ); e i
, t ~ N (0, σ 2 e,t

), [1]

where y i,t is the natural log-transformed monthly average PM 2.5 for i = 1... I

, I

= 499 sites and t = 1... T

, T = 48 monthly time peri

ods from 1999 to 2002, and s i the projected spatial coordinate pair for the i th location. The Albers Equal Area Conic USGS projec

tion was used for all geographic data (ESRI).
g

t ( s i ) accounts for residual monthly spatial variability and g ( s i ) for time-invariant spatial variability. Z i,t, 1 through Z i,t,P are time-vary

ing covariates, X i, 1 through X i,Qare timeinvariant GIS-based covariates, and a t is a monthly intercept that controls for the mean across all sites. d 1 through d Q

and f 1 through f

P are one-dimensional penalized spline smooth functions for Q time-invariant and P time-varying covariates. Note that b iis a sitespecific random effect, hence our characteriza

tion of the model as a GAMM. The post-1999 model was fit in two stages: the first to estimate site-specific terms adjust

ing for time-varying covariates and residual spatial variability, and the second to modelthe site-specific terms using site-specific, timeinvariant GIS-based predictors and residual time-invariant spatial variability. The form of the two-stage model was

y i
, t

= u i + f 1 ( Z i
, t

,1 )+...+ f P ( Z i
, t

, P

)

+ a t + g t ( s i ) + e i
, t ; e i
, t ~ N (0, σ 2

e , t

) [2]

û i

= a + d 1 (

X

i ,1 )+...+ d Q (

X

i
, Q

)

+ g ( s i ) + b i ; b i ~ N (0, σ 2 b ), [3]

where û i is an estimated site-specific intercept that represents the adjusted long-term mean at each location. The first stage (Equation 2) was fit itera

tively in a back-fitting arrangement (Hastie and Tibshirani 1990) with u i + f 1 ()

... f p () esti

mated jointly and a t + g t ( s i ) estimated sepa

rately by month, such that variability in the concentrations is parsed between the covari

ates and the residual spatial terms in the first

**Figure 1.** Map of the PM 2.5 monitoring locations in the study region and adjacent states for monitoring sites reporting data from 1988 to 1998 and separately from 1999 to 2002. The locations of weather stations reporting visual range are also shown. 150 75 0 150 km Visual range stations

PM 2. 5 monitors Pre-1998 Post-1999 Study region Adjacent statesYanosky et al.



### 524 volume 117 | number 4 | April 2009 • Environmental Health Perspectives stage. The second stage (Equation 3) was fit to the estimated site-specific û i terms. Both stages were fit using calls to the gam() function (iter

atively in the first stage) in the mgcv library of R (R Development Core Team 2008). Model predictions were obtained based on generating the covariates at locations of interest for each month and, once estimated, were transformed to the original scale by exponentiation. Pre-1999 PM 2.5 spatiotemporal model. The pre-1999 model assumed a relatively sim

ple spatiotemporal structure because of the sparseness of PM 2.5 data during the earlier time period. Specifically, seasonal spatial trends were assumed to be constant across years using the natural log-transformed ratio of PM 2.5 to pre

dicted PM 10 , which was preferable to alterna

tive transformations (log

- transformed PM 2.5 , ratio of PM 2.5 to PM 10 , and logit- transformed PM 2.5 to PM 10 ratio) based on results from our exploratory models. Thus, the model described variation in the log- transformed ratio of PM 2.5 to PM 10 , accounting for spatial and temporal variability in this ratio:

y i
, t

= a + d 1 (

X

i ,1 )+...+ d Q (

X

i
, Q

)

+ g ( s i ) + c ( z ^ i
, t

) + d ( b ^

ext i
, t )

+ f

1 ( Z i
, t

,1 )+...+ f P ( Z i
, t

, P

) + a Season

+ g Season ( s i ) + h ( t ) + b i + e i
, t b i ~ N (0, σ 2 b ); e i
, t ~ N (0, σ 2 e,Season ), [4] where y i,t = 1n(PM 2.5 / P ^ M 10 ) for i = 1... I

, I = 546 sites and t = 1... T

, T = 180 monthly time periods from 1988 to 2002, z ^ i,t is the predicted log PM 10 concentration and b ^
ext i,t is the predicted log extinction coefficient after correction for relative humidity and trunca

tion. Season has four levels (winter, spring, summer, and fall); a Season is a season- specific intercept
, and g Season ( s i ) accounts for residual seasonal spatial variability. h ( t ) is a smoothly varying intercept that controls for the monthly mean across all sites. Also, c and dare onedimensional penalized spline smooth func

tions. Otherwise, the notation is the same as for the post-1999 model. The pre-1999 model was also fit in two stages:

y i
, t

= u i + c ( z ^ i
, t

) + d ( b ^

ext i
, t )

+ f

1 ( Z i
, t

,1 )+...+ f P ( Z i
, t

, P

) + [ a Season + g Season ( s i )]

+ h(t) + e i
, t ; e i
, t ~ N (0, σ 2 e,Season ) [5]

û i

= a + d 1 (

X

i ,1 )+...+ d Q (

X

i
, Q

)

+ g ( s i ) + b i ; b i ~ N (0, σ 2 b ) [6]

We obtained model predictions based on generating the covariates at locations of interest for each month, transforming the predicted PM 2.5 to PM 10 ratio back to the original scale by exponentiation, and then multiplying by predicted PM 10 . Covariate selection. Covariates that were expected a priori to have a physical influence on PM 2.5 levels were considered for inclu

sion in the post-1999 and pre-1999 models. Covariates remained in the model if their observed relationship with PM 2.5 was con

sistent with known pollutant behavior (i.e., if the smooth functions were generally in the expected direction, analogous to slopes in the expected direction for linear terms) and if their inclusion improved predictive performance. Nonlinearity in the covariate effects was accounted for using spline terms, with at most 6 degrees of freedom considered sufficient to describe the general shape of each function. Using this procedure, we found several covariates to be important predictors in the post-1999 model, including distance to near

est A1 (Figure 2A) and A2 roads, urban land use, point-source PM 2.5 emissions within 10 km, elevation (Figure 2A), wind speed, and precipitation. The covariate selection process contributed little to potential overfitting: The cross-validation R 2 was comparable among locations used for the covariate selection pro

cess and those held out from this process (0.75 vs. 0.77, respectively). Important PM 2.5 predictors for the pre-1999 model included predicted PM 10 (from Yanosky et al. 2008), extinction coefficient (Figure 2B), wind speed, temperature (Figure 2B), precipitation, monthly time trend, distance to nearest A1 road, urban land use, block group and county population density, and elevation. Model validation. We used cross-vali

dation techniques both to inform covariate selection and to compare alternative model specifications. To evaluate the performance of the post-1999 model, we randomly selected monitoring sites from the 1999

– 2002 data and assigned them exclusively to 1 of 10 sets. Data from sets 1–9 (each set contains approximately 10% of the data) were held out in turn, with predictions generated at the locations of the held-out observations. Because the covariate selection process involved fitting multiple can

didate models to the same data, we used set 10 to assess whether the covariate selection process contributed to overfitting, as reported above. Similarly, for the pre-1999 model, mon

itoring sites from the 1988

– 1998 data were assigned exclusively to one of five sets at ran

dom. Data from sets one through five were held out in turn, with predictions generated at the locations of the held-out observations. Because 1988–1998 data were not spatially rep

resentative across the region, we performed aseparate cross-validation step, where the pre1999 model was used to predict at the held-out locations during 1999. Data from 1999 for a select group of monitors were included in the model- fitting data set to ensure that the **Figure 2.** ( A ) Percent change in PM 2.5 and PM 10 as a function of distance to nearest interstate (A1 road) and elevation from the post-1999 PM 2.5 model and the PM 10 model presented by Yanosky et al. (2008), respectively. ( B ) Percent change in PM 2.5 as a function of extinction coefficient and temperature from the pre-1999 PM 2.5 model. 15

10

5

0

–5

0 2,000 4,000 6,000**Percent change in P M** **Percent change in P M** 

0 200 400 600 800 1,000 40

30

20

10

0

–10

–30

**Distance to nearest A1 road (m) Elevation above sea-level (m**  **)**

**Extinction coefficient (km –1**  **) A ir temperature (°C)**  0.05 0.10 0.15 –10 0 1 0 2 0**Percent change in PM 2.5**  30

20

10

0

–10

–20 15

10

5

0**Percent change in PM 2.** **5 A** 

**B** PM 2. 5 PM 10Predicting chronic fine and coarse PM exposures

Environmental Health Perspectives • volume 117 | number 4 | April 2009 

### 525 number of monitoring locations and their pop

ulation density in 1999 was similar to that from 1988

– 1998. Model overfitting due to covariate selection was not examined for the pre-1999 model because of the limited amount of PM 2.5 data available between 1988 and 1998. We determined the predictive abilities of the PM models using the squared Pearson cor

relation between the held-out observations and the model predictions (cross-validation R 2 ), with both on the original rather than the log scale. We calculated prediction errors by sub

tracting held-out observations from the model predictions. Bias in model predictions was determined using the mean prediction error and the slopes from linear regression of the held-out values against the observations. We estimated the absolute precision of model pre

dictions by taking the square root of the mean of the squared prediction errors (RMSPE), and relative precision by dividing the RMSPE by the arithmetic mean of the PM measure

ments in the appropriate size category. Bias and absolute precision of predicted PM 2.5 lev

els from cross-validation were also evaluated by geographic location (state), urban land use and population density, season, monitoring network, and monitoring objective. To evaluate the ability of our models to predict PM 10–2.5 , we compared monthly pre

dicted PM 10–2.5 values with PM 10–2.5 mea

sured values at sites where PM 2.5 and PM 10 monitors were co-located. We performed sep

arate cross-validation procedures for PM 10–2.5 by again holding out co-located PM 2.5 and PM 10 sites in 10 sets. Sensitivity analyses. We compared the per

formance of our model with monthly varying spatial surfaces with alternative models with either four seasonal spatial terms or 16 sea

son by year spatial terms. A smoothly varying intercept, h ( t ), where t = 1,..., T

, T

= 180, was added to both alternative models to allow for monthly control for the mean across all sites. Additionally, we compared our post-1999 model with two simple spatial interpolation approaches, inverse distance weighting (IDW) and nearest neighbor (NN) interpolation. For the IDW approach, we used the squared inverse distance as the weighting function. For the NN approach, we excluded monitors that were not within 50 km of another monitor. To evaluate the impact of the extinction coefficient term on the predictive ability of the pre-1999 model, we fit the model withoutthis term. Additionally, we compared our pre1999 model with a simplistic ordinary least squares linear regression model that assumes a single, fixed PM 2.5 to PM 10 ratio over space and time:

y i
, t

= β 0 + β 1 exp( z ^ i
, t

)

+ e i
, t ; e i
, t ~ N (0, σ 2

e ), [7]

where y i,t is the monthly site average PM 2.5 and exp( z ^ i,t

) is the predicted PM 10 concentra

tion as in Equation 4. Data analysis. To evaluate the spatial het

erogeneity in PM levels, we considered how PM levels at unmeasured locations (as best described by model predictions) differed on average from observed values by calculating mean squared deviations (MSDs) between PM levels measured at monitoring sites and PM levels predicted at unmeasured locations

on an 8- km grid over the study region, doing this separately for PM 2.5 , PM 10 , and PM 2.5–

10 . Only monitoring sites that were within a metropolitan statistical area (MSA) and those in the AQS network with population exposure as their monitoring objective were included in this analysis. PM 10 predictions were made using our previously developed model for PM 10 (Yanosky et al. 2008). As an estimate of the amount of spatial variability at each distance relative to the total spatial vari

ability at the regional scale, we divided MSDs at distances < 400 km by the median MSD at 400 km for each pollutant. We then used the medians and 25th and 75th percentiles within 10- km bins to describe trends in spatial vari

ability as a function of distance (i.e., the vario

gram). Only monthly predictions from 1999 to 2002 for PM 2.5, PM 10 , and PM 10–2.5 were included in this analysis, because pre-1999 PM 2.5 predictions were based on predicted PM 10 levels.


### Results We first present results describing the predic

tive performance of the PM 2.5 and PM 10–2.5 models. We then present results from our sensitivity analyses, comparing these PM 2.5 models to alternative models. Finally, we pres

ent results from our evaluation of the spatial heterogeneity of PM 2.5 , PM 10 , and PM 10–2.5 levels across metropolitan areas. Model predictive performance. Model fit, cross-validation, and regression results for the post-1999 and pre-1999 PM 2.5 models are presented in Table 1. Both PM 2.5 mod

els performed well. The post-1999 PM 2.5 model explained most of the variability in monthly measured PM 2.5 levels (model fit R 2 = 0.87 and cross-validation R 2 = 0.77). The pre-1999 PM 2.5 model also performed well (cross- validation R 2 = 0.68 and 0.69 using 1988–1998 and 1999 cross-validation data, respectively). For both models, results from linear regression showed good agreement between measured values and model predic

tions, with slopes near 1 and intercepts near 0. Both the post-1999 and pre-1999 PM 2.5 models predicted PM 2.5 concentrations with lit

tle bias (–0.2 and –0.3 μg/m 3 , respectively) and high absolute and relative precision (2.2 and 2.7 μg/m 3 and 16.7% and 15.4%, respectively). Both PM 2.5 models performed well across the region and also in both rural and urban areas, al though each had slightly better precision in spring and fall seasons compared with sum

mer and winter [see Supplemental Material, Table 1 (online at http://www. ehponline.org/ members/2008/11692/suppl.pdf )]. Our PM 2.5 models and our previous PM 10 model (Yanosky et al. 2008) relied on similar predictors, all including elevation, urban land use, distance to nearest A1 road, wind speed, and precipitation. Further, smooth functions of these predictors common to each model were generally similar in shape (Figure 2A).

**Table 1.** Model fit, cross-validation, and regression results for the post-1999 and pre-1999 PM 2.5 models.

No. of Covariates CV results Time period Model description CV data spatial terms a included Model fit R 2 b Intercept c Slope c CV R 2
Post-1999 (1999–2002) Final model 1999–2002 48 month by year d Full set 0.85 0.8 ± 0.07 0.95 ± 0.01 0.77 Alternative season by year spatial terms 1999–2002 16 season by year d Full set 0.78 0.3 ± 0.08 1.00 ± 0.01 0.72

Alternative seasonal spatial terms 1999–2002 4 seasonal d Full set 0.76 0.5 ± 0.09 0.99 ± 0.01 0.68 IDW 1999–2002 None None

— 0.61 ± 0.11 0.92 ± 0.01 0.60 NN 1999–2002 e None None

— 3.0 ± 0.13 0.77 ± 0.01 0.61 Pre-1999 (1988–1998) Final model 1988–1998 f 4 seasonal d Full set 0.76 –0.4 ± 0.33 1.05 ± 0.02 0.68 1999 g 4 seasonal d Full set 0.76 1.09 ± 0.22 0.94 ± 0.02 0.69 Final model without extinction coefficient 1999 g 4 seasonal d Full set minus 0.76 1.19 ± 0.21 0.94 ± 0.02 0.70 extinction coefficient Alternative fixed ratio model 1999 g None None 0.61 1.5 ± 0.29 0.85 ± 0.02 0.53

CV, cross-validation.

a Corresponds to the extent of control for space–time interaction in the model. b Using data in the study region only. c Presented as parameter estimate ± SE from linear regression of heldout observations on predictions. d Number of time-varying spatial terms fit in the first stage of the model in addition to one spatial term fit in the second stage. e Only 5,210 observations available for comparison versus 10,444 observations for other models; excluded monitors did not have another monitor within 50 km. f One observation in Ohio and all observations from one site in New York excluded as outliers. g Four observations at one site in New York excluded as outliers.Yanosky et al.



### 526 volume 117 | number 4 | April 2009 • Environmental Health Perspectives For example, PM 2.5 and PM 10 levels decreased similarly from 0 to about 250 m away from the nearest A1 road, indicating that this term captured micro -scale (0–100 m) and middle

scale (100–250 m) gradients in PM 2.5 and PM 10 from traffic sources similarly. Beyond about 400 m, however, PM 2.5 and PM 10 showed a slightly different relation with dis

tance to A1 roads, as predicted PM 10 levels decreased more quickly than PM 2.5 levels. This modest difference is likely attributable to the greater relative contribution of secondary aerosols to PM 2.5 , resulting in a more homog

enous spatial distribution of PM 2.5 than PM 10 . These findings suggest that distance to A1 road described neighborhood-scale (500 m–4 km) and urban-scale (4–7 km) gradients slightly better for PM 2.5 than for PM 10 . Our models performed less well in esti

mating PM 10–2.5 than PM 2.5 . With post-1999 data, measured PM 10–2.5 values were gener

ally centered around the predicted PM 10–2.5 values (Figure 3A), with little proportional bias (slope of 1.14) and much of the vari

ability in the measurements explained by the predictions ( R 2 = 0.64). From cross-validation (Figure 3B), predicted PM 10–2.5 levels again exhibited little bias (slope = 0.89), but predic

tive performance was poorer (cross- validation R 2 of 0.39). The absolute and relative pre

cision of predicted PM 10–2.5 levels from cross-validation was relatively low (5.5 μg/m 3 and ± 61.5%, respectively). Predictive per

formance for pre-1999 data was comparable with that for post-1999 data (cross-validation R 2 = 0.33). Predictive performance of the PM models improved substantially when long-term, multi year averages (rather than monthly averages) of measured PM levels and model predictions (one mean value per site for measurements and predictions) were com

pared (cross-validation R 2 = 0.81 and 0.75 for post-1999 and pre-1999 PM 2.5 models, respectively, and 0.63 and 0.65 for post-1999 and pre-1999 PM 10–2.5 excluding one site innorthern Maine, respectively). These longterm average results are for sites with at least39 of 48 months of valid data for the post1999 models, and at least 10 of 12 months of valid data (using 1999 cross-validation data) for the pre-1999 models. Sensitivity analyses. Model fit, crossvalidation, and regression results for alterna

tive models are presented in Table 1. Our post-1999 and pre-1999 PM 2.5 models were preferable to alternative models with differ

ent spatiotemporal structures. The post-1999 model with monthly spatial terms performed better than simpler models with four seasonal or 16 season by year spatial terms (Table 1; cross-validation R 2 = 0.77 for the post-1999 model vs. 0.68 and 0.72 for the alternative models, respectively), which demonstrates the improvement in predictive performance gained by modeling space–time interaction at the monthly rather than the seasonal level. Further, the post-1999 model performed bet

ter than simple spatial interpolation methods (cross-validation R 2 = 0.60 for IDW and 0.61 for NN). With a matched data set, which addressed the fact that fewer data points were available for NN analysis (5,210 instead of 10,444), the performance of our post-1999model performance remained high (crossvalidation R 2 = 0.76) and was comparable with that across all observations. Additionally, the simple spatial interpolation methods exhibited greater proportional bias than the post-1999 model, as indicated by slopes from linear regression of 0.92 and 0.77 for IDW and NN, respectively, versus 0.95 for the post-1999 model (Table 1). Although the performance of the pre-1999 model with the extinction coefficient term was highly comparable with that without (Table 1; cross-validation R 2 = 0.69 with extinction coefficient and 0.70 without, both using 1999 cross-validation data), the extinction coeffi

cient was included in the final model

, as its smooth term indicated a nearly linear increase in PM 2.5 levels with higher extinction coeffi

cients (Figure 2B), a result consistent with light extinction theory (Ozkaynak et al. 1985). However, likely because of the truncation in airport visibility data, the distribution of the extinction coefficients was quite narrow, with 25th and 75th percentiles of 0.06 and 0.07 km –1 , respectively, leading to little predictive power in this covariate. We also noted con

siderable changes in the seasonal smooth spa

tial surfaces in pre-1999 PM 2.5 models with and without the extinction coefficient term. This may explain the lack of improvement in predictive performance when extinction coef

ficients were included in the model. The alternative pre-1999 model that assumed a fixed PM 2.5 to PM 10 ratio over time and space also exhibited lower predictive **Figure 3.** Scatter plot of monthly predicted versus measured PM 10–2.5 concentrations in the northeastern and midwestern United States from 1999 to 2002 ( A ) including all measured locations and ( B ) from cross-validation. 50

40

30

20

10

0 50

40

30

20

10

0**Measured P M** **10–2. 5** **(μg/m 3** **) Measured P** **M 10–2.** **5 (μg/m** **3 )** 

0 1 0 2 0 3 0 4 0 5 0 0 10 20 30 40 50 **Predicted PM 10–2.**  **5 (μg/m**  **3 ) P**  **redicted PM 10–2.**  **5 (μg/m**  **3 )**  1:1 line Regression line

y = –0.6 (± 0.12) + x * 1.14 (± 0.01) R 2 = 0.64 y = 1.4 (± 0.16) + x * 0.89 (± 0.02) R 2 = 0.39 **A B** 

**Figure 4.** Maps of mean predicted concentrations for ( A ) PM 2.5 across all months from 1988 to 1998 from the pre-1999 model, ( B ) PM 2.5 across all months from 1999 to 2002 from the post-1999 model, and ( C ) PM 10–2.5 across all months from 1999 to 2002, estimated as the difference in monthly predictions of PM 10 and PM 2.5 . Low and high values are 5th to 95th percentiles, respectively. **A B**  **C**

150 75 0 150 km Predicted PM 2. 5 (μg/m 3 ) High : 13.9

Low : 5.8 Predicted PM 2. 5 (μg/m 3 )

High : 16.8

Low : 6.2 Predicted PM 10–2. 5 (μg/m 3 ) High : 8.8

Low : 1.3Predicting chronic fine and coarse PM exposures

Environmental Health Perspectives • volume 117 | number 4 | April 2009 

### 527 performance than the pre-1999 model (Table 1; cross-validation R 2 = 0.53 vs. 0.69, respectively), which demonstrates the improvement in predictive performance from allowing spatial trends in the PM 2.5 to PM 10 ratio to change seasonally and from account

ing for temporal trends in this ratio. PM concentrations across the study region. Maps of the predicted PM 2.5 con

centrations averaged across all months from 1988 to 1998 are shown in Figure 4A and from 1999 to 2002 in Figure 4B. Both maps show consistent spatial patterns, with pre

dicted PM 2.5 concentrations highest in eastern Ohio, Pennsylvania, and western Maryland and in urban areas such as New York City, New York, and Detroit, Michigan. Predicted PM 2.5 levels are lowest in northern areas of Maine, New York, and Michigan. As shown in Figure 4C, mean predicted PM 10–2.5 con

centrations (1999

– 2002) have a substantially different spatial pattern, with considerably more spatial heterogeneity in PM 10–2.5 levels compared with PM 2.5 . Spatial trends in predicted PM 2.5 concen

trations were relatively similar to those found for PM 10 (Yanosky et al. 2008). However, PM 2.5 concentrations were more spatially uni

form compared with PM 10 and especially with PM 10–2.5 . Figure 5 shows trends in the propor

tion of local spatial variability relative to the total at 400 km as a function of the distance between PM monitoring locations and predic

tion grid points for PM 2.5 , PM 10 , and PM 10–
2.5 . At distances < 75 km, this proportion was generally lower and less variable for PM 2.5 than for PM 10 or PM 10–2.5 (Figure 5), consistent with the greater contribution of secondary par

ticles (such as sulfate) for PM 2.5 compared with PM 10 or PM 10–2.5 (Allen and Turner 2008). Despite this, at even moderately close dis

tances (between 0 and 25 km) relative to the scale of a typical metropolitan area, PM 2.5 , PM 10 , and especially PM 10–2.5 exhibited sub

stantial spatial heterogeneity as evidenced by the sharp increase in the 75th percentiles (Figure 5) of the MSD
x /MSD 400 , suggesting the importance of within-city gradients in PM levels. Further, these spatial gradients were larger for PM 10 and PM 10–2.5 than for PM 2.5 , as the 75th percentiles for PM 10 and PM 10–2.5 increased more quickly than those for PM 2.5 , although medians for PM 10 were only slightly higher than for PM 2.5 . Local, neighborhoodscale spatial variability (MSD 4 km ) constituted 7.8% of the total, regional scale variability (MSD 400 km ) for PM 2.5 , compared with 15.8% and 35.1% for PM 10 and PM 10–2.5 , respectively, corresponding to the values at the leftmost part of the x -axis in Figure 5. Additionally, using a standard sums of squares decomposition, the proportions of spatial, temporal, and spatiotemporal variability in predicted PM 2.5 levels were 50%, 32%, and 18%, respectively, compared with 41%, 23%, and 36% for PM 10–2.5 , respectively, indicat

ing that monthly spatiotemporal changes are more important for PM 10–2.5 than for PM 2.5 . Together, results indicate that PM 2.5 levels were more spatially homogenous than PM 10 and especially PM 10–2.5 levels. When these within-MSA population exposure monitoring sites were stratified by urbanization (above or below the median pro

portion of urban land use within 1 km and county-level population density of the moni

toring site), PM 2.5 , PM 10 , and PM 10–2.5 levels were more spatially heterogeneous in highly urbanized areas compared with less urbanized areas. Further, the amount of spatial hetero

geneity in PM levels varied by season, likely because PM 2.5 represents a smaller proportion of PM 10 in the spring and fall compared with winter and summer. In spring and fall, PM 2.5 concentrations were more spatially homoge

neous than PM 10 and especially PM 10–2.5 , as evidenced by seasonal variograms and by the number of degrees of freedom in the seasonal spatial terms of the pre-1999 PM 2.5 model (38 and 5 for spring and fall, and 355 and 303 for winter and summer, respectively). However, in winter and to a lesser extent sum

mer, PM 2.5 and PM 10 exhibited similar levels of spatial variability, with both again less vari

able than PM 10–2.5 . These results suggest that spatial variability in PM 10 in summer and especially in winter is affected primarily by variability in PM 2.5 , likely because of winter wood smoke emissions and differing meteo

rologic conditions across seasons. Results are consistent with greater mixing of the atmo

sphere in nonwinter months, with increased formation of secondary particles in summer.


### Discussion By including location-specific covariates, our spatiotemporal models provide highly spatially and temporally resolved estimates of outdoor PM concentrations that can be used to esti

mate chronic exposures to PM 2.5 and, for the first time, to PM 10–2.5 as well. Importantly, our models allow spatial trends in predicted monthly average PM 2.5 and PM 10–2.5 concen

trations to change over time, and therefore are able to provide realistic, accurate, and precise monthly location-specific predictions of PM levels that may be used as measures of chronic exposure. These values can be averaged to pro

vide exposure measures for different time peri

ods relevant to chronic health impacts (e.g., previous month, previous 3 months, previ

ous year). As demonstrated by their satisfac

tory predictive ability, the models are able to predict PM 2.5 and PM 10–2.5 levels over large regions while maintaining within-city spatial gradients using location-specific covariates and spatial smoothing. We demonstrate for the first time a methodology for estimating

**Figure 5.** The proportion of local spatial variability relative to the total (MSD x /MSD 400 ) as a function of distance for PM 2.5 , PM 10 , and PM 10–2.5 from 1999–2002 for AQS population exposure monitors within MSAs. Solid lines are medians and dotted lines are 25th and 75th percentiles. The inset shows the medians from 0 to 40 km in greater detail. 2.5

2.0

1.5

1.0

0.5

0.0**MSD x** **/MSD 400** 

50 100 150 **Distance between points (km )**  0.6

0.5

0.4

0.3

0.2

0.1

0.0**MS D** **x /MSD** **400**

10 20 30 40 **Distance between points (km )**  PM 2.5 PM 10 PM 10–2. 5Yanosky et al.



### 528 volume 117 | number 4 | April 2009 • Environmental Health Perspectives chronic exposures to ambient PM 10–2.5 over long time periods and across the northeastern and midwestern United States. Our study represents a significant improvement over previous methods of esti

mating long-term average ambient PM con

centrations, including the nearest air pollutant monitor approach and citywide or county

wide averaging. By using spatial smoothing to perform local spatial averaging, our models reduce the potential for measurement error in estimates of chronic exposure to ambient PM, as occurs with the nearest neighbor approach, while maintaining within-city spatial gradi

ents. Our post-1999 PM 2.5 model also pro

vided improved performance and exhibited less proportional bias (i.e., regression slope closer to 1) than interpolation methods. However, given their relatively strong per

formance and their simplicity, interpolation methods may be a reasonable alternative to more complicated spatiotemporal modeling of PM 2.5 when resources are limited. Our approach to model PM 2.5 concen

trations separately by time period (post1999 and pre-1999) allowed monthly PM 2.5 concentrations to be predicted for 15 years, including years before 1999 when PM 2.5 data were extremely sparse. PM 2.5 models for both time periods performed well by sea

son and population density, as demonstrated using cross-validation techniques, albeit with slightly lower precision in some states and in winter and summer seasons. Results from both models demonstrated the importance of accounting for spatial variability, temporal variability, and location-specific covariates. Model and spatiotemporal structure differed substantially between the pre- and post-1999 models, however, because of dif

ferences in the amount of available PM 2.5 monitoring data. The pre-1999 model used predicted PM 10 levels from the model described by Yanosky et al. (2008) to estimate PM 2.5 to PM 10 ratios. Although temporal trends in these ratios were adjusted for, this model assumed a fixed spatial trend in the ratio across years for each season (e.g., spa

tial trends in the ratio in winter 1988 were assumed to be the same as those in winter 1989). This innovative modeling approach leveraged the available information on PM 10 , extinction coefficients, and available PM 2.5 data while allowing for spatially and tempo

rally varying calibrations between PM 2.5 and PM 10 levels. In comparison, given the richer data set, the post-1999 PM 2.5 model could rely on measured PM 2.5 without PM 10 and was more flexible, allowing spatial trends in PM 2.5 to change monthly. Our models are similar in structure to land-use regression models in that they incorporate the effects of GIS-based predic

tors. However, our models are distinct from typical land-use regression models in that they account for residual spatial variability using spatial smoothing and allow these spa

tial surfaces to change over time. Further, because they do not involve spatial model

ing, typical land-use regression analyses are only applicable for small regions (across a city or metropolitan area), whereas our model is applicable to much larger areas (i.e., the entire northeastern and midwestern United States ). Results from our analyses showing the impor

tance of GIS-based predictors were consistent with those of Brauer et al. (2003). There the authors used GIS-based predictors in mul

tiple linear regression models (but not spatial smoothing techniques) to explain variability in annual PM 2.5 concentrations in the three areas: the Netherlands ; Munich, Germany ; and Stockholm County, Sweden. Brauer et al. reported RMSPE values between 1.1 and 1.6 μg/m 3 for annual averages. Even though our models were based on monthly averages, included spatiotemporal modeling in addition to GIS predictors, and were applied to a larger domain (the northeastern and midwestern United States as opposed to three small areas in Europe), our model performance results were similar. We found only slightly higher RMSPE values of 2.2 and 2.7 μg/m 3 for the post-1999 and pre-1999 models, respectively, for our monthly average PM 2.5 predictions. Predictive performance was lower for PM 10–2.5 than for PM 2.5 because a ) fewer monitoring sites were available (only those with co-located PM 10 and PM 2.5 provide measured values), b ) prediction errors were compounded by estimating PM 10–2.5 by dif

ference, and c ) the PM 10 model contained little highly spatially resolved information on PM 10–2.5 sources, such as fugitive dust emis

sions and unpaved roads. In addition, relative precision estimates were worse for PM 10–2.5 than for PM 2.5 because typical PM 10–2.5 con

centrations are lower than those for PM 2.5 . However, predictive performance improved substantially when long-term, multiyear aver

ages were considered. Results from our analysis of the spatial heterogeneity of PM levels indicated that PM 2.5 had the smallest amount of local spatial variability, followed by PM 10 , and finally by PM 10–2.5 . Further, they suggest that PM 2.5 measurements at a SAM site are more repre

sentative for surrounding populations than for PM 10 , and more representative for PM 10 than for PM 10–2.5 . However, spatial variabil

ity for the three pollutants increased nearly linearly at distances from 0 to about 50 km (though slightly less steeply between 25 and 50 km), suggesting that there was no fixed distance within which PM concentrations were entirely uniform within a metropolitan area. Together, these results suggest that epi

demiologic analyses that compare the effects of PM 2.5 with PM 10–2.5 using only measure

ments at SAM sites as indicators of exposure should be interpreted with caution because of greater exposure prediction error for PM 10–2.5 than for PM 10 or PM 2.5 (Zeger et al. 2000). Our modeling approach assumes iso

tropy and stationarity of the spatial surfaces, normality and homoscedasticity of residu

als, independence of the time-varying spa

tial terms, uniform effects of GIS-based and meteorologic covariates across space, and no interactions of these covariates; see Paciorek

et al. (2008) for evidence supporting these simplifying assumptions. The satisfactory performance of our two models demonstrates their suitability for esti

mating long-term average population exposures to ambient PM concentrations. These models have the potential to reduce exposure measure

ment error in air pollution health effect studies. We are currently using these models to esti

mate chronic PM 2.5 and PM 10–2.5 exposures in the Nurses' Health Study, a large prospective cohort study of U.S. women, to examine the chronic health impacts of these pollutants.

**R efe**  **r ences** 

Abbey D, Nishino N, McDonnell W, Burchette R, Knutsen S, Beeson W, et al. 1999. Long-term inhalable particles and other air pollutants related to mortality in nonsmokers. Am J Respir Crit Care Med 159:373–382. Allen D, Turner J. 2008. Transport of atmospheric fine particu

late matter: Part 1. Findings from recent field programs on the extent of regional transport within North America. J Air Waste Manag Assoc 58:254–264. Brauer M, Hoek G, van Vliet P, Meliefste K, Fischer P, Gehring U, et al. 2003. Estimating long-term average particulate air pollution concentrations: application of traffic indica

tors and geographic information systems. Epidemiology 14:228–239. Brunekreef B, Forsberg B. 2005. Epidemiological evidence of effects of coarse airborne particles on health. Eur Respir J 26:309–318. Dockery D, Pope CI, Xu X, Spengler J, Ware J, Fay M, et al. 1993. An association between air pollution and mortality in six US cities. N Engl J Med 329:1753–1759. Faulke S, Husar R. 1998. Maps of PM 2.5 over the U.S. derived from regional PM 2.5 and surrogate visibility and PM 10 monitoring data. Report number 98-A918. Available: http:// capita.wustl.edu/CAPITA/DataSets/FM_VISIB/FMVISIB. HTML [accessed 3 March 2006]. Finkelstein M, Jerrett M, DeLuca P, Finkelstein N, Verma D, Chapman K, et al. 2003. Relation between income, air pollution and mortality: a cohort study. Can Med Assoc J 169:397–402. Haneuse S, Wakefield J, Sheppard L. 2007. The interpretation of exposure effect estimates in chronic air pollution stud

ies. Stat Med 26:3172–3187. Hastie T, Tibshirani R. 1990. Generalized Additive Models. New York:Chapman and Hall. Homer C, Huang C, Yang L, Wylie B, Coan M. 2004. Development of a 2001 national landcover database for the United States. Photogrammetric Engineering and Remote Sensing 70(7):829-840. Jerrett M, Burnett R, Ma R, Pope CI, Krewski D, Newbold K,
et al. 2005. Spatial analysis of air pollutions and mortality in Los Angeles. Epidemiology 16:727–736. Miller K, Siscovick D, Sheppard L, Shepherd K, Sullivan J, Anderson G, et al. 2007. Long-term exposure to air pollu

tion and incidence of cardiovascular events in women. N Engl J Med 356:447–458. National Climatic Data Center. 2005. Home page. Available: http:// www.ncdc.noaa.gov/oa/ncdc.html [accessed 1 June 2005]. Ozkaynak H, Schatz A, Thurston G, Isaacs R, Husar R. 1985.Predicting chronic fine and coarse PM exposures

Environmental Health Perspectives • volume 117 | number 4 | April 2009 

### 529 Relationships between aerosol extinction coefficients derived from airport visual range observations and alter

native measures of airborne particle mass. J Air Pollut Control Assoc 35:1176–1185. Paciorek C, Yanosky J, Suh H. 2008. Practical large-scale spatio-temporal modeling of particulate matter con

centrations. Ann Appl Stat. Available: http://www.imstat. org/aoas/next_issue.html [accessed 24 February 2009 ]. Pope CA III, Burnett R, Thun M, Calle E, Krewski D, Ito K, et al.2002. Lung cancer, cardiopulmonary mortality, and longterm exposure to fine particulate air pollution. JAMA 287:1132–1141. Pope CA III, Burnett R, Thurston G, Thun M, Calle E, Krewski D,
et al. 2004. Cardiovascular mortality and long term expo

sure to particulate air pollution. Circulation 109:71–77. Pope CA III, Thun M, Namboodiri M, Dockery D, Evans J, Speizer F, et al. 1995. Particulate air pollution as a predic

tor of mortality in a prospective study of U.S. adults. Am J Respir Crit Care Med 151:669–674. R Development Core Team. 2008. R: A language and environ

ment for statistical computing. Vienna:R Foundation for Statistical Computing. Available: http://www.R-project.org [accessed 26 February 2009]. Spengler J, Koutrakis P, Dockery D, Raizenne M, Speizer F. 1996. Health effects of acid aerosols on North American children: air pollution exposures. Environ Health Perspect 104:492–499. Suh H, Nishioka Y, Allen G, Koutrakis P, Burton R. 1997. The Metropolitan Acid Aerosol Characterization Study: results from the summer 1994 Washington, D.C. field study. Environ Health Perspect 105:826–834. U.S. EPA. 2005. U.S. Environmental Protection Agency National Emissions Inventory. Available: http://www.epa.gov/ttn/ chief/trends/ [accessed 6 October 2005]. U.S. EPA. 2009. U.S. Environmental Protection Agency Air Quality System. Available: http://www.epa.gov/ttn/airs/ airsaqs/ [accessed 24 February 2009]. USGS. 2004. U.S. Geological Survey National Land Cover Dataset. Available: http://nlcd.usgs.gov/ [accessed 27 May 2004]. USGS. 2005. U.S. Geological Survey National Elevation Dataset. Available: http://ned.usgs.gov/ [accessed 2 May 2005]. Visibility Information Exchange Web System. 2004. Home page. Available: http://vista.cira.colostate.edu/views/ [accessed 27 May 2004]. Wood S. 2006. Generalized Additive Models: An Introduction with R. Boca Raton, FL:Chapman and Hall/CRC Press. Yanosky J, Paciorek C, Schwartz J, Laden F, Puett R, Suh H. 2008. Spatio-temporal modeling of chronic PM 10 exposure for the Nurses' Health Study. Atmos Environ 42:4047–4062. Zeger S, Thomas D, Dominici F, Samet J, Schwartz J, Dockery D, et al. 2000. Exposure measurement error in time-series studies of air pollution: concepts and consequences. Environ Health Perspect 108:419–426.