Open access
Research Article
14 November 2014

A Unified Spatiotemporal Modeling Approach for Predicting Concentrations of Multiple Air Pollutants in the Multi-Ethnic Study of Atherosclerosis and Air Pollution

Publication: Environmental Health Perspectives
Volume 123, Issue 4
Pages 301 - 309



Cohort studies of the relationship between air pollution exposure and chronic health effects require predictions of exposure over long periods of time.


We developed a unified modeling approach for predicting fine particulate matter, nitrogen dioxide, oxides of nitrogen, and black carbon (as measured by light absorption coefficient) in six U.S. metropolitan regions from 1999 through early 2012 as part of the Multi-Ethnic Study of Atherosclerosis and Air Pollution (MESA Air).


We obtained monitoring data from regulatory networks and supplemented those data with study-specific measurements collected from MESA Air community locations and participants’ homes. In each region, we applied a spatiotemporal model that included a long-term spatial mean, time trends with spatially varying coefficients, and a spatiotemporal residual. The mean structure was derived from a large set of geographic covariates that was reduced using partial least-squares regression. We estimated time trends from observed time series and used spatial smoothing methods to borrow strength between observations.


Prediction accuracy was high for most models, with cross-validation R2 (R2CV) > 0.80 at regulatory and fixed sites for most regions and pollutants. At home sites, overall R2CV ranged from 0.45 to 0.92, and temporally adjusted R2CV ranged from 0.23 to 0.92.


This novel spatiotemporal modeling approach provides accurate fine-scale predictions in multiple regions for four pollutants. We have generated participant-specific predictions for MESA Air to investigate health effects of long-term air pollution exposures. These successes highlight modeling advances that can be adopted more widely in modern cohort studies.


Keller JP, Olives C, Kim SY, Sheppard L, Sampson PD, Szpiro AA, Oron AP, Lindström J, Vedal S, Kaufman JD. 2015. A unified spatiotemporal modeling approach for predicting concentrations of multiple air pollutants in the Multi-Ethnic Study of Atherosclerosis and Air Pollution. Environ Health Perspect 123:301–309;


Cohort studies have shown an association between long-term exposure to air pollution and cardiovascular morbidity and mortality (Brook et al. 2010; Miller et al. 2007; Pope and Dockery 2006). To estimate these associations, epidemiologic studies develop exposure prediction models to predict pollutant concentrations over long periods of time at cohort home addresses based on monitoring data from regulatory networks or study-specific monitoring campaigns. Although early models were based on region-wide averages or nearest-monitor approaches, more current methods include land-use regression (LUR) (Hoek et al. 2008; Jerrett et al. 2007), the use of satellite and remote sensing data (Kloog et al. 2011), geostatistical methods such as kriging (Beelen et al. 2009), generalized additive models (Hart et al. 2009), or a combination of these techniques (Beckerman et al. 2013; Bergen et al. 2013; Mercer et al. 2011).
The Multi-Ethnic Study of Atherosclerosis and Air Pollution (MESA Air) is investigating the association between long-term air pollution exposure and measures of cardiovascular health (Kaufman et al. 2012). MESA Air is following a cohort of > 6,000 individuals in six metropolitan regions: Baltimore, Maryland; Chicago, Illinois; Los Angeles, California; New York, New York; St. Paul, Minnesota; and Winston-Salem, North Carolina. The primary exposures of interest in MESA Air are fine particulate matter (≤ 2.5 μm; PM2.5), nitrogen dioxide (NO2), oxides of nitrogen (NOx), and black carbon, as measured by light absorption coefficient (LAC). One goal of MESA Air is the development of advanced statistical methods that incorporate extensive supplemental monitoring to improve the prediction of intra-urban pollutant variability. MESA Air health effect analyses require spatiotemporal predictions of ambient outdoor concentrations for all four pollutants in all six metropolitan regions for times ranging from 1999 through 2012.
In general, exposure prediction models developed in one city do not transfer well to another city (Allen et al. 2011), so prediction models are often study- and city-specific (e.g., Franklin et al. 2012). A challenge for multi-city studies such as MESA Air that combine data from subcohorts and include several pollutant measures is generating predictions that are of comparable quality across pollutants and cities. Here we present a unified and flexible spatiotemporal modeling framework for the four MESA Air pollutants. We apply a standardized approach to model selection for all pollutants and regions, allowing the intrinsic flexibility of our modeling framework to account for differences in the way pollution processes behave in different regions.


To predict outdoor concentration of pollutants at MESA Air participant residences, we fit a separate spatiotemporal exposure prediction model for each pollutant (PM2.5, NO2, NOx, and LAC) in each metropolitan region. Briefly, our model decomposes the space–time field of concentrations into spatially varying long-term (i.e., duration of study period) averages, spatially varying seasonal and long-term trends, and spatially correlated but temporally independent residuals, and accommodates data from the complex monitoring design described below. We modeled on a 2-week time scale because of the 2-week sampling of MESA Air supplementary monitoring instruments; this allows for flexible aggregation of predictions to time scales of interest for health effects analyses (e.g., 12 months before clinic visit).
Monitoring data. We used three sources of outdoor air monitoring data. We obtained hourly NO2 and NOx and daily PM2.5 concentration measurements in each metropolitan region from 1 January 1999 through 31 March 2012 from the U.S. Environmental Protection Agency (EPA) Air Quality System (AQS) (, including data from the Interagency Monitoring of Protected Visual Environments (IMPROVE) network ( No AQS data were used for black carbon because of differences in collection methods from the MESA Air LAC measurement methods described below. We aggregated hourly data into daily averages and subsequently averaged daily values to the 2-week scale. AQS monitors that had < 2 years of data or had irregular temporal coverage (e.g., operated only in the summer) were not used.
In each metropolitan region, we defined the modeling area to be locations within approximately 75 km of each metropolitan center (Figure 1). AQS monitors within the modeling regions were considered for inclusion in the model, and predictions at participant residences were restricted to locations within these modeling regions. In New York, MESA Air participants were recruited from both New York City and Rockland County, so the modeling region included locations near both areas. In Winston-Salem, only one AQS monitoring location for NO2 and NOx met inclusion criteria. To have a complete time series for the 14-year modeling period, an AQS monitor in Charlotte, North Carolina, was included for estimating time trends. In Chicago, the modeling region was further restricted to locations west of –87.5°W longitude because some covariates were unavailable east of that meridian. In Los Angeles, only locations south and/or west of the San Gabriel Mountains were included.
Figure 1 Maps of the modeling areas (denoted by dashed black line) in the six metropolitan regions, including monitor and subject locations. Abbreviations: Fixed, MESA Air fixed monitoring sites; Home, MESA Air home monitoring sites; Snapshot, MESA Air snapshot monitoring sites; Participant, MESA Air participant residence location (moved slightly to protect confidentiality).
To better capture the within-city variability of pollutant concentrations, MESA Air conducted a supplementary monitoring campaign targeting the study cohort from July 2005 through August 2009. The MESA Air measurements were 2-week cumulative measurements that began and ended on Wednesdays. Measurements of NO2 and NOx were made using Ogawa passive samplers, and PM2.5 mass was measured on Harvard Personal Environmental Monitor impactors using Teflon filters. LAC was computed from the Teflon filters via reflectance. A detailed description of the data collection and site selection procedures has been previously published (Cohen et al. 2009).
The MESA Air monitoring campaign included three types of monitoring sites: fixed, home, and snapshot. Fixed sites were operated for the duration of the 4-year MESA Air sampling period to provide long time series of measurements, with one fixed site collocated with an AQS monitor in each region. Samples of participant residences in each metropolitan region were selected for monitoring as home outdoor sites on a rotating basis, with most locations monitored one to three times in different seasons. Snapshot sites, which measured only NO2 and NOx, were located in clusters to capture gradients near sources (e.g., primary roadways) and monitored for three 2-week periods, one each in winter, summer, and either spring or fall.
In New York City, data from the New York City Community Air Survey (NYCCAS) were used to supplement the AQS and MESA Air data (Matte et al. 2013; NYC Department of Health 2014). The NYCCAS data consist of 2-week measurements of PM2.5, NO2, NOx, and LAC collected during December 2008–December 2010 in a manner consistent with MESA Air sampling protocols. Five NYCCAS reference sites (one in each borough) collected measurements throughout the sampling period, and 150 NYCCAS distributed sites were monitored once per season during this time. Because of the similarity in monitoring scheme, we treated NYCCAS reference sites in the same manner as MESA Air fixed sites, and NYCCAS distributed sites in the same manner as MESA Air home sites, in our models. The NYCCAS data and a small subset of the MESA Air 2-week data were centered on different weeks than most of the MESA Air measurements. To align these measurements with the rest of the MESA Air data, we treated these measurements as if they were made 1 week earlier or later, as appropriate.
Between 0.4% (LAC) and 1.6% (NO2) of the pollutant measurements were below the limit of detection (LOD) and were replaced with the value LOD/2. The number of each type of monitoring site by region and pollutant is provided in Table 1. The range of the number of PM2.5 observations at each monitoring site during the study period is provided in Table 2, along with summary statistics for the site means. Corresponding statistics for NO2, NOx, and LAC observations are provided in Supplemental Material, Tables S1–S3.
Table 1 Number of monitors by site type, region, and pollutant.
Site typePM2.5NO2NOxLAC
Baltimore, MD
MESA fixed5555
MESA home86878786
MESA snapshot 104104
Chicago, IL
MESA fixed6666
MESA home136113113136
MESA snapshot 129129
Los Angeles, CA
MESA fixed7777
MESA home113120120113
MESA snapshot 252250
New York, NY
MESA fixed3333
MESA home107119118107
MESA snapshot 157157
NYCCAS reference5555
NYCCAS distributed150150150150
St. Paul, MN
MESA fixed3443
MESA home126132132129
MESA snapshot 107107
Winston-Salem, NC
MESA fixed4444
MESA home114117117114
MESA snapshot 121121
Table 2 Summary of PM2.5 monitoring data.
Site typeNo. of observations per siteSite means (μg/m3)
MinimumMaximumMinimumMaximumMean ± SD
Baltimore, MD
AQS6434510.916.913.4 ± 1.4
MESA fixed189212.115.413.7 ± 1.25
MESA home137.322.714.3 ± 3.1
Chicago, IL
AQS7132011.716.414.0 ± 1.3
MESA fixed68712.214.013.1 ± 0.75
MESA home145.219.511.5 ± 3.2
Los Angeles, CA
AQS8234510.722.816.2 ± 3.5
MESA fixed768513.719.316.2 ± 2.0
MESA home120.742.616.9 ± 6.1
New York, NY
AQS513429.317.112.5 ± 1.8
MESA fixed498311.515.713.7 ± 2.1
MESA home133.541.615.1 ± 4.9
NYCCAS reference51528.89.99.4 ± 0.42
NYCCAS distributed686.819.811.0 ± 11.0
St. Paul, MN
AQS553057.911.610.0 ± 0.91
MESA fixed81899.610.510.0 ± 0.46
MESA home155.027.410.3 ± 3.8
Winston-Salem, NC
AQS8634610.315.913.4 ± 1.5
MESA fixed809313.013.813.4 ± 0.35
MESA home149.022.814.3 ± 2.6
Geographic covariates. We compiled > 300 geographic covariates for use in the model (see Supplemental Material, Table S4). These covariates included proximity measures (distance to nearest major road, intersection, truck route, railway, railyard, coastline, airport, and port) and buffer measures (major road length, truck route length, land-use category, long-term vegetation index, population density, and emission sources). We included a long-term average of the dispersion model output from a modified implementation of the Caline3QHCR line-source model (Eckhoff and Braverman 1995). The Caline3QHCR model incorporates distance, traffic volume, meteorology, and diurnal traffic patterns in each region.
Geographic covariates with minimal variation or potentially highly influential values were excluded from the modeling process. Specifically, variables were removed if a) > 80% of monitoring sites had the same value, b) > 2% of observations were more than 5 SDs away from mean, c) the standard deviation of the distribution of values at participant residences was more than five times the standard deviation of the distribution of values at monitoring locations, or d) the maximum value was 10% among all monitoring sites (for land-use variables only). These filters were applied separately for each pollutant and region.
Spatiotemporal model. The monitoring data were highly unbalanced, with a small number of locations providing long time series of several years’ duration and a larger number of locations providing broader spatial coverage, but at a relatively small number of time points. A hierarchical spatiotemporal model had been previously developed to accommodate the unbalanced nature of the MESA Air data (Lindström et al. 2014; Sampson et al. 2011; Szpiro et al. 2010). This model can be written as
C(s,t) = μ(s,t) + v(s,t), [1]
where C(s,t) represents the log-transformed 2-week average pollutant concentration at location s and time t. The μ(s,t) term is the spatiotemporal mean surface, and the v(s,t) term represents spatiotemporal residual variation.
We break down the spatiotemporal mean into components
where β0(s) is the long-term mean at location s, fi(t) are smooth time trends, and βi(s) are spatially varying coefficients for the time trends.
The time trends are estimated from AQS and MESA Air fixed sites (and NYCCAS reference sites in New York) using a procedure developed by Fuentes et al. (2007) and previously described in detail by Sampson et al. (2011). In brief, we applied an expectation-maximization procedure to fill in missing values in the time series and derived the trends from a singular value decomposition. We smoothed the trends using splines, controlling the smoothness with the degrees of freedom (df) parameter. The model assumes the time trends account for enough of the temporal structure that the residuals are uncorrelated in time.
The long-term averages β0(s) and time trend coefficients βi(s) are modeled as spatial random fields with a spatial mean, distributed as
βi ~ N[Xi(si, Σii, σi, τi)], i = 0, 1, …, m. [3]
Here, Xi(s) are reduced-dimension summaries of the geographic covariates (described in detail below) at location s, and αi are vectors of coefficients to be estimated. The covariance structure for βi, denoted by Σi, is either an independence model with variance τi or a spatial smoothing model with exponential covariance function parameterized by range φi, partial sill σi, and nugget τi (Cressie 1993).
The zero-mean spatiotemporal residual term v(s,t) in Equation 1 has a spatial correlation structure and is assumed independent at each time point. It includes a random effect for each time point to model short-term variations that affect an entire region, such as large-scale meteorological events.
Partial least squares (PLS) scores. Rather than include each of the hundreds of geographic covariates directly in the model or use variable selection methods, we reduced the dimensionality of the covariates using PLS. In a manner similar to principal components analysis (PCA), PLS computes linear combinations, called scores, of the columns of a data matrix. Unlike PCA, the PLS procedure constructs scores that maximize the covariance between the scores and an outcome rather than the variance between the scores. A technical explanation of the PLS algorithm is provided by Abdi (2010). Sampson et al. (2013) described the application of PLS for spatial models, and here we describe how we applied the method to spatiotemporal data.
PLS regression requires a single outcome value for each location. Because the MESA Air data are unbalanced time series, we first derived values that could be used as outcomes in PLS regression. For each AQS, fixed, and NYCCAS reference site s, we regressed the time series of observations C(s,t) on the smoothed time trends using ordinary least squares regression with mean function E[C(s,t)] = γs0 f0(t) + … + γsm fm(t) to get estimates (^γs0,…, ^γsm) for each location. For each time trend, PLS regression was performed separately with the ^γsi as the outcomes and the matrix of geographic covariates as the predictors. This gave a set of PLS scores for each location that was different for each time trend. PLS scores at home and snapshot monitoring sites were predicted using the geographic variables at those locations and the score definitions defined from the regression at fixed sites. PLS regression was performed using the pls package (Mevik et al. 2011), in R (R Core Team; Scores were included in the model as the Xi(s) in Equation 3.
Parameter estimation and model selection. Once the PLS scores Xi(s) and time trends fi(t) were computed, the remaining parameters were calculated via maximum likelihood using the SpatioTemporal package, version 1.1.7 (Lindström et al. 2012), in R. We varied several model parameters and used cross-validation to find the best-fitting model in each metropolitan region, as described below. We considered different values for the number of time trends (either 1 or 2), the df for smoothing time trends (either 4 or 8 per year), the number of PLS scores per time trend (2 or 3), and the covariance structure of the βi fields (spatial smoothing or no spatial smoothing).
Cross-validation procedure. The primary interest of MESA Air is in long-term average exposures, so we assessed model performance using cross-validation of long-term averages (LTAs). Because the highly unbalanced structure of the monitoring data means that LTAs at home sites are computed from a handful of observations of a few weeks’ duration, whereas LTAs at fixed sites are computed from long time series, we performed cross-validation separately for each site type. For home sites and NYCCAS distributed sites, we used 10-fold cross-validation, which leaves out one-tenth of the data in turn. For AQS, fixed, and NYCCAS references sites, we used leave-one-out cross-validation because the total number of monitors was relatively small. For snapshot sites, we used 10-fold cross-validation, with monitors in the same cluster left out together. For all three schemes, the covariance parameters (but not the time trends or PLS scores) were re-estimated using all but the left-out sites. Pollutant concentrations at the left-out sites were predicted using the parameters estimated from the remaining data.
We assessed cross-validation performance using two measures: root mean-squared error (RMSE) and cross-validation R2 (denoted by R2CV). Letting yj denote the mean observed value and ^yj the mean of the predicted values for the observed time points at monitoring site j, RMSE and R2CV were computed on the original scale of the data according to the formulas
R2CV = max(0,1 – RMSE2/MSEobs), [5]
where MSEobs = 1/n Σnj = 1 (yjy)2 is the mean-squared error of the observed values. R2CV provides a measure of fit to the 1-1 line, in contrast to the typical regression-based R2 (R2CVreg), which measures fit to the regression line and is computed as the square of the correlation coefficient between the cross-validation predictions and the observed values. R2CV reflects the contrast of interest because our goal is accurate prediction at unmeasured locations, and it is typically lower than R2CVreg. Although R2CV was the primary metric for our model evaluation, we present R2CVreg for comparison with published results from other authors. Because we are most interested in spatial contrasts between individual exposures within each region, we prioritized home-site over fixed-site R2CV and RMSE in the model selection process.
In the context of the hierarchical model (Equation 1), it is challenging to separate the spatial and temporal contributions to R2CV for cross-validation of temporally sparse data sets such as the home sites and NYCCAS distributed sites. Lindström et al. (2014) proposed three temporally adjusted adaptations of R2CV that use data from the AQS and fixed sites as the reference MSE instead of MSEobs in Equation 5 in order to focus on spatial prediction accuracy. R2Avg uses the average values at AQS and fixed sites within that region. R2Close uses the closest (in absolute distance) AQS or fixed site. R2Smooth uses the smoothed time trend at the closest AQS or fixed site.
Prediction at participant locations. Using the best models from each metropolitan region, predictions of pollutant log-concentrations at participant residences were made on a 2-week scale from January 1999 through March 2012. We back-transformed these predictions using exponentiation to return them to the original scale of concentration measurements and computed averages of 2-week predictions over the study period.


Model structure. Table 3 provides an overview of the model structure selected for each metropolitan region and pollutant. Most models have only one time trend, although all the New York models have two. The two smoothed trends for the NO2 model in Los Angeles are shown in Figure 2. Figure 2 also includes plots of the fitted trends for a selected AQS site and fixed site.
Table 3 Model structure for the best model (selected by cross-validation) for each pollutant and metropolitan region.
ModelNo. of time trendsaNo. of PLS scoresbdf/year in time trendcSpatial smoothingd
Long-term average (β0)Time trend coefficients (βi)
Baltimore, MD
Chicago, IL
Los Angeles, CA
New York, NY
St. Paul, MN
Winston‑Salem, NC
aSelected from either 1 or 2 time trends. bSelected from either 2 or 3 PLS scores; scores were covariates in the mean component of the long-term average (β0) and time trend (βi) fields [denoted by Xi(s) in Equation 3]. cSelected from either 4 or 8 degrees of freedom (df) per year; controls smoothness of estimated time trends. dYes, exponential covariance structure. No, independent covariance structure.
Figure 2 Time trends for the NO2 model in Los Angeles. The top panel shows the smoothed time trends calculated from AQS and fixed sites. The middle and bottom panels show the observed data and fitted trends at an AQS site and fixed site, respectively.
For PM2.5, there was noticeable heterogeneity of the best models across metropolitan regions (Table 3). New York and Winston-Salem had no spatial smoothing in the long-term PM2.5 average [β0(s)], and no model had spatial smoothing in the PM2.5 time trend coefficients [βi(s)]. Half of the regions had two time trends, whereas the other half had only a single time trend for PM2.5. For NO2, all of the models had spatial smoothing for the long-term average, and the same was true for NOx except in New York.
The relative contribution of geographic covariates to the PLS scores varied by region and pollutant. In the Supplemental Material, Figure S1 shows the correlations between covariates and PLS scores for the NO2 model in Chicago, which are representative of the general patterns in the other regions (data not shown). Overall, the distance-to-feature covariates and vegetation measures [Normalized Difference Vegetation Index (NDVI) and low development, open development, forest, and wetland land use] tended to have the opposite correlation from emissions and traffic measures (A1, A2/A3, and truck route lengths and intersection counts) within buffers.
Model results. Table 4 shows the cross-validation metrics for all models, broken down by pollutant and region. These metrics assess how well the site means are modeled, incorporating both the spatial and temporal components of the predictions. Scatter plots of predictions and observed values are provided in Figure 3 for AQS and fixed sites and in Supplemental Material, Figure S2, for home sites. Metrics for the snapshot sites (for NO2 and NOx) and for AQS and fixed sites (all four pollutants) on the 2-week scale are reported in the Supplemental Material, Tables S5 and S6.
Table 4 Cross-validation measures of predictive accuracy for site means at monitoring locations.
RegionAQS and MESA fixed sitesMESA home sites
Los Angeles1.280.830.842.920.770.78
New Yorkb0.590.910.912.800.540.56
St. Paul0.600.450.841.780.780.79
Los Angeles2.230.880.893.130.770.78
New Yorkb1.860.920.933.820.780.78
St. Paul1.270.930.941.240.870.87
Los Angeles6.740.870.875.690.910.92
New Yorkb8.850.610.8916.660.500.50
St. Paul1.690.980.983.580.830.84
Los Angeles0.1140.700.930.2660.690.71
New Yorkb0.1470.750.790.3290.510.52
St. Paul0.0430.910.920.0740.690.69
aUnits for RMSE are μg/m3 (PM2.5), ppb (NO2 and NOx), and 10–5/m (LAC). bNew York models include NYCCAS reference sites with AQS and fixed sites, and NYCCAS distributed sites with home sites.
Figure 3 Long-term averages of cross-validated predictions and observations for AQS and fixed monitoring locations for each pollutant.
Predictive accuracy was generally good (R2CV > 0.6) to excellent (R2CV > 0.8) in almost all regions for each pollutant. NOx models in Baltimore and Los Angeles had the best performance at MESA home sites (R2CV of 0.92 and 0.91, respectively) (Table 4). The lowest R2CV for home sites was in the Chicago NO2 model (0.45), but its RMSE (3.31 ppb) was comparable with those in New York and Los Angeles (3.82 and 3.13 ppb, respectively). At AQS and fixed sites, R2CV was very good (0.70 for LAC in Los Angeles) to excellent (0.98 for NOx in St. Paul), with two notable exceptions: Winston-Salem NOx (0.00) and St. Paul PM2.5 (0.45). However, in both cases the RMSE was comparable with the corresponding RMSE for models in other regions. The small range of observed data (9.8–22.4 ppb) (see Supplemental Material, Table S2) and the small number of monitors in Winston-Salem (Table 1) explain the low R2CV for NOx in that city.
Table 5 provides three versions of temporally adjusted R2CV at home sites (R2Avg, R2Close, and R2Smooth). For NO2 and NOx, these temporally adjusted R2CV are fairly similar to the unadjusted R2CV reported in Table 4, suggesting that we are predicting spatial differences well. For PM2.5, however, the temporally adjusted R2CV are consistently lower than the unadjusted R2CV.
Table 5 Temporally adjusted cross-validation measures of predictive accuracy for home site means.
Los Angeles0.400.230.43
New Yorka0.480.360.38
St. Paul0.230.290.62
Los Angeles0.630.660.66
New Yorka0.890.780.64
St. Paul0.770.890.90
Los Angeles0.810.850.88
New Yorka0.720.640.52
St. Paul0.790.880.85
Los Angeles0.280.340.48
New Yorka0.590.650.53
St. Paul0.670.800.84
General formula for R2 measures is R2 = max (0,1 – RMSE2/MSEobs). R2Avg uses the mean-squared error of the average observed values at AQS and fixed sites within the region as MSEobs. R2Close uses the mean-squared error of the observed values at the closest AQS or Fixed site as MSEobs. R2Smooth uses the mean-squared error of the values of the smoothed time trend at the nearest AQS or Fixed site as MSEobs. aIncludes NYCCAS distributed sites.
Box plots of the long-term averages of predictions at participant residences are provided in Figure 4. On average, predicted concentrations tended to be higher in New York and Los Angeles, consistent with the higher observed monitoring values in those regions. Variability in predictions is also greatest in these two cities, especially in the tails of the distributions.
Figure 4 Pollutant- and region-specific box plots of long-term averages of predictions from 1999 through early 2012 at participant residence locations. Metropolitan region abbreviations: Bal, Baltimore; Chi, Chicago; LA, Los Angeles; NY, New York; SP, St. Paul; W-S, Winston-Salem. Boxes extend from the 25th to the 75th percentile, horizontal bars represent the median, whiskers extend 1.5 times the length of the interquartile range above and below the 75th and 25th percentiles, respectively, and outliers are presented as points.
Supplemental Material, Table S7, provides performance metrics for the New York models when the NYCCAS data were excluded from the modeling process. Without the NYCCAS data, R2CV was noticeably lower for PM2.5 (0.79 vs. 0.91 at AQS and fixed sites, and 0.36 vs. 0.54 at home sites) and LAC (0.55 vs. 0.75 at AQS and fixed sites, and 0.43 vs. 0.51 at home sites).


We present here a complex and successful approach to predicting long-term air pollution concentrations for application in a cohort study. Although this approach was tailored to this particular well-characterized cohort study—taking advantage of cohort-specific monitoring, for example—the success of the approach demonstrates modeling improvements that can be adopted for application in future population-based research on spatially varying pollutants. We believe that this approach to capturing variation in within-region pollution highlights advances that should be adopted in the next generation of air pollution cohort studies, both for understanding contrasts at relatively low concentrations in the United States and at the higher concentrations experienced globally.
We describe a unified framework for implementing exposure prediction models of four air pollutants in six metropolitan regions that easily incorporates spatially and temporally unbalanced monitoring data. The application of a consistent modeling framework to all regions and pollutants is important for studies such as MESA Air that use exposure estimates from multiple subcohorts together in health analyses. Although we applied the same approach in all regions, we varied the model structure to best fit the data for each region and pollutant. This unified modeling approach was shown to have very good model performance (R2CV > 0.70) for almost all of the pollutants and regions. The architecture for this modeling approach is publicly available through the SpatioTemporal R package (Lindström et al. 2012).
As a result of the success of our spatiotemporal modeling approaches, we are confident in using these approaches to model outdoor pollutant concentrations in epidemiological analyses in this cohort and in other populations residing in these same communities. We have also found that implementation of portions of this approach, such as PLS regression of a large set of geographic covariates combined with spatial smoothing via universal kriging, can be used with good success in other regions to predict pollutant concentrations without the same level of small-area monitoring (Sampson et al. 2013). The NYCCAS data increased the spatial density of the monitoring data in New York, which was likely one reason for the improved model performance when the data were included. For LAC, the NYCCAS data provided particular benefit because they allowed the model to be extended through 2010, which would not have been possible with only the MESA Air data.
A majority of models included spatial smoothing in the long-term average. This suggests that although PLS scores derived from geographic covariates can predict much of the spatial variation in the data, benefit is gained from borrowing strength across observations nearby in space.
Differences in the underlying pollutant variability likely caused some of the differences seen in temporally adjusted R2CV. PM2.5 tended to exhibit less small-scale spatial variation and greater temporal variability, leading to temporally adjusted R2CV that are noticeably lower than the unadjusted measures. The NO2 and NOx data tended to exhibit greater spatial variability, and the similarity of the unadjusted and temporally adjusted R2CV values suggests that the unadjusted R2CV are not overly inflated by well-predicted temporal variation.
The modeling approach presented here does have several limitations. First, we used geographic covariates that were constant in time [although the modeling framework readily extends to spatiotemporal covariates (Lindström et al. 2014)]. Changes in these variables likely occurred during the study decade, but we nonetheless believe that the time-constant geographic variables still provided a useful means to predict long-term pollutant concentrations. Second, the calculation of PLS scores was limited to AQS and fixed sites because they had long time series. For LAC in particular, this means that the scores were based on a very small number of locations because the LAC model relied only on MESA Air data (plus NYCCAS data for New York). Third, the cross-validation model selection procedure conditioned on the time trends and PLS scores. Overfitting may have occurred in the cross-validation of the AQS and fixed sites, because the left-out observations were used to estimate the time trends and PLS scores. However, because the home sites were not used in estimating time trend or in defining the PLS scores, any overfitting was restricted to the AQS and fixed-site cross-validation. This provides further motivation for prioritization of cross-validation metrics from home sites when selecting the best models.


Our unified spatiotemporal modeling method successfully characterized outdoor concentrations of multiple air pollutants at the homes of cohort members in multiple metropolitan regions. This flexible and powerful modeling approach can incorporate an unbalanced monitoring data structure, leveraging data from supplemental monitoring campaigns that increase the spatial coverage of monitoring data. The method was easily transferred between regions and pollutants, allowing for straightforward comparison between model fits across regions. Although aspects of our techniques are particularly tailored to the unique data and resources of MESA Air, lessons learned here can be applied to understand the spatial and temporal variation of pollutants in future cohort studies. Advances in fine-scale modeling resolved in both space and time are important for the next generation of cohort studies assessing health effects of environmental agents.

Supplemental Material

(1.7 MB) PDF
Click here for additional data file.


Abdi H. 2010. Partial least squares regression and projection on latent structure regression (PLS Regression). Wiley Interdiscip Rev Comput Stat 2:97-106;
Allen RW, Amram O, Wheeler AJ, Brauer M. 2011. The transferability of NO and NO2 land use regression models between cities and pollutants. Atmos Environ 45:369-378;
Beckerman BS, Jerrett M, Serre M, Martin RV, Lee SJ, van Donkelaar Aet al. 2013. A hybrid approach to estimating national scale spatiotemporal variability of PM2.5 in the contiguous United States. Environ Sci Technol 47:7233-7241;
Beelen R, Hoek G, Pebesma E, Vienneau D, de Hoogh K, Briggs DJ. 2009. Mapping of background air pollution at a fine spatial scale across the European Union. Sci Total Environ 407:1852-1867;
Bergen S, Sheppard L, Sampson PD, Kim SY, Richards M, Vedal Set al. 2013. A national prediction model for PM2.5 component exposures and measurement error-corrected health effect inference. Environ Health Perspect 121:1017-1025;
Brook RD, Rajagopalan S, Pope CA, Brook JR, Bhatnagar A, Diez-Roux AVet al. 2010. Particulate matter air pollution and cardiovascular disease: an update to the scientific statement from the American Heart Association. Circulation 121:2331-2378;
Cohen MA, Adar SD, Allen RW, Avol E, Curl CL, Gould Tet al. 2009. Approach to estimating participant pollutant exposures in the Multi-Ethnic Study of Atherosclerosis and Air Pollution (MESA Air). Environ Sci Technol 43:4687-4693
Cressie NAC. 1993. Statistics for Spatial Data, Revised Edition. Hoboken NJ John Wiley & Sons.
Eckhoff P, Braverman T. 1995. Addendum to the User’s Guide to CAL3QHC Version 2.0 (CAL3QHCR user’s guide). Available: [accessed 4 March 2015].
Franklin M, Vora H, Avol E, McConnell R, Lurmann F, Liu Fet al. 2012. Predictors of intra-community variation in air quality. J Expo Sci Environ Epidemiol 22:135-147;
Fuentes M, Guttorp P, Sampson PD. 2007. Using transforms to analyze space-time processes. In Statistical Methods for Spatio-Temporal Systems (Finkenstädt B, Held L, Isham V, eds.). Boca Raton, FL:Chapman Hall/CRC, 77–150.
Hart JE, Yanosky JD, Puett RC, Ryan L, Dockery DW, Smith TJet al. 2009. Spatial modeling of PM10 and NO2 in the continental United States, 1985–2000. Environ Health Perspect 117:1690-1696;
Hoek G, Beelen R, de Hoogh K, Vienneau D, Gulliver J, Fischer Pet al. 2008. A review of land-use regression models to assess spatial variation of outdoor air pollution. Atmos Environ 42:7561-7578;
Jerrett M, Arain M, Kanaroglou P, Beckerman B, Crouse D, Gilbert NLet al. 2007. Modeling the intraurban variability of ambient traffic pollution in Toronto, Canada. J Toxicol Environ Health A 70:200-212;
Kaufman JD, Adar SD, Allen RW, Barr RG, Budoff MJ, Burke GLet al. 2012. Prospective study of particulate air pollution exposures, subclinical atherosclerosis, and clinical cardiovascular disease: the Multi-Ethnic Study of Atherosclerosis and Air Pollution (MESA Air). Am J Epidemiol 176:825-837;
Kloog I, Koutrakis P, Coull BA, Lee HJ, Schwartz J. 2011. Assessing temporally and spatially resolved PM2.5 exposures for epidemiological studies using satellite aerosol optical depth measurements. Atmos Environ 45:6267-6275;
Lindström J, Szpiro AA, Sampson PD, Bergen S, Oron AP. 2012. SpatioTemporal: Spatio-Temporal Model Estimation. R Package Version 1.1.7. Available: [accessed 1 September 2013].
Lindström J, Szpiro AA, Sampson PD, Oron AP, Richards M, Larson TVet al. 2014. A flexible spatio-temporal model for air pollution with spatio-temporal covariates. Environ Ecol Stat 21:411-433;
Matte TD, Ross Z, Kheirbek I, Eisl H, Johnson S, Gorczynski JEet al. 2013. Monitoring intraurban spatial patterns of multiple combustion air pollutants in New York City: design and implementation. J Expo Sci Environ Epidemiol 23:223-231;
Mercer LD, Szpiro AA, Sheppard L, Lindström J, Adar SD, Allen RWet al. 2011. Comparing universal kriging and land-use regression for predicting concentrations of gaseous oxides of nitrogen (NOx) for the Multi-Ethnic Study of Atherosclerosis and Air Pollution (MESA Air). Atmos Environ 45:4412-4420;
Mevik BH, Wehrens R, Liland KH. 2011. PLS: Partial Least Squares and Principal Component regression. R Package Version 2.3-0. Available: [accessed 1 October 2012].
Miller KA, Siscovick DS, Sheppard L, Shepherd K, Sullivan JH, Anderson GLet al. 2007. Long-term exposure to air pollution and incidence of cardiovascular events in women. N Engl J Med 356:447-458
NYC Department of Health. 2014. NYC Community Air Survey. Available: [accessed 13 January 2014].
Pope CA, Dockery DW. 2006. Health effects of fine particulate air pollution: lines that connect. J Air Waste Manag Assoc 56:709-742;
Sampson PD, Richards M, Szpiro AA, Bergen S, Sheppard L, Larson TVet al. 2013. A regionalized national universal kriging model using partial least squares regression for estimating annual PM2.5 concentrations in epidemiology. Atmos Environ 75:383-392;
Sampson PD, Szpiro AA, Sheppard L, Lindström J, Kaufman JD. 2011. Pragmatic estimation of a spatio-temporal air quality model with irregular monitoring data. Atmos Environ 45:6593-6606;
Szpiro AA, Sampson PD, Sheppard L, Lumley T, Adar SD, Kaufman JD. 2010. Predicting intra-urban variation in air pollution concentrations with complex spatio-temporal dependencies. Environmetrics 21:606-631;

Information & Authors


Published In

Environmental Health Perspectives
Volume 123Issue 4April 2015
Pages: 301 - 309
PubMed: 25398188


Received: 17 January 2014
Accepted: 11 November 2014
Published online: 14 November 2014



Joshua P. Keller
Department of Biostatistics, and
Casey Olives
Department of Environmental and Occupational Health Sciences, University of Washington, Seattle, Washington, USA
Sun-Young Kim
Department of Environmental and Occupational Health Sciences, University of Washington, Seattle, Washington, USA
Institute of Health and Environment, Seoul National University, Seoul, Korea
Lianne Sheppard
Department of Biostatistics, and
Department of Environmental and Occupational Health Sciences, University of Washington, Seattle, Washington, USA
Paul D. Sampson
Department of Statistics, University of Washington, Seattle, Washington, USA
Adam A. Szpiro
Department of Biostatistics, and
Assaf P. Oron
Core for Biomedical Statistics, Seattle Children’s Hospital, Seattle, Washington, USA
Johan Lindström
Centre for Mathematical Science, Lund University, Lund, Sweden
Sverre Vedal
Department of Environmental and Occupational Health Sciences, University of Washington, Seattle, Washington, USA
Joel D. Kaufman
Department of Environmental and Occupational Health Sciences, University of Washington, Seattle, Washington, USA


Address correspondence to J.P. Keller, Department of Biostatistics, University of Washington, Health Sciences Building, Box 357232, 1705 NE Pacific St., Seattle, WA 98195 USA. Telephone: (206) 543-1044. E-mail: [email protected]

Competing Interests

The views expressed in this document are solely those of the authors, and the U.S. EPA does not endorse any products or commercial services mentioned in this publication.

Competing Interests

The authors declare they have no actual or potential competing financial interests.

Funding Information

This research was developed under a STAR research assistance agreement, RD831697 [The Multi-Ethnic Study of Atherosclerosis and Air Pollution (MESA Air)], awarded by the U.S. Environmental Protection Agency (EPA). It has not been formally reviewed by the U.S. EPA. Additional support was provided by U.S. EPA grants RD-83479601-0 (J.P.K., C.O., S.Y.K., L.S., P.D.S., A.A.S., A.P.O., S.V., J.D.K.) and CR-834077101-0 (S.Y.K., L.S., A.A.S., J.L.), National Institute of Environmental Health Sciences (NIEHS) grants K24ES013195 (C.O., L.S., J.D.K.) and P30ES007033 (J.D.K.), and the Biostatistics, Epidemiologic, and Bioinformatic Training in Environmental Health Training Grant from the NIEHS (T32ES015459) (J.P.K.).

Metrics & Citations


About Article Metrics


Download citation

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click DOWNLOAD.

Cited by

  • Air Pollutants and Risk of Parkinson’s Disease among Women in the Sister Study, Environmental Health Perspectives, 10.1289/EHP13009, 132, 1, (2024).
  • Comparison of Air Pollution Exposures and Health Effects Associations Using 11 Different Modeling Approaches in the Women’s Health Initiative Memory Study (WHIMS), Environmental Health Perspectives, 10.1289/EHP12995, 132, 1, (2024).
  • Coarse Particulate Matter and Markers of Inflammation and Coagulation in the Multi-Ethnic Study of Atherosclerosis (MESA) Population: A Repeat Measures Analysis, Environmental Health Perspectives, 10.1289/EHP12972, 132, 2, (2024).
  • Outdoor Ultrafine Particulate Matter and Risk of Lung Cancer in Southern California, American Journal of Respiratory and Critical Care Medicine, 10.1164/rccm.202305-0902OC, 209, 3, (307-315), (2024).
  • Representativeness of the US EPA PM monitoring site locations to the US population: implications for air pollution prediction modeling, Journal of Exposure Science & Environmental Epidemiology, 10.1038/s41370-024-00644-3, (2024).
  • Ambient air pollution and rate of spontaneous abortion, Environmental Research, 10.1016/j.envres.2023.118067, 246, (118067), (2024).
  • Long-term ozone exposure and lung function in middle childhood, Environmental Research, 10.1016/j.envres.2023.117632, 241, (117632), (2024).
  • Evaluating low-cost monitoring designs for PM2.5 exposure assessment with a spatiotemporal modeling approach, Environmental Pollution, 10.1016/j.envpol.2023.123227, 343, (123227), (2024).
  • A comprehensive review of the development of land use regression approaches for modeling spatiotemporal variations of ambient air pollution: A perspective from 2011 to 2023, Environment International, 10.1016/j.envint.2024.108430, 183, (108430), (2024).
  • Traffic-related air pollution and dementia incidence in the Adult Changes in Thought Study, Environment International, 10.1016/j.envint.2024.108418, 183, (108418), (2024).
  • See more

View Options

View options


View PDF

Get Access

Restore your content access

Enter your email address to restore your content access:

Note: This functionality works only for purchases done as a guest. If you already have an account, log in to access the content to which you are entitled.







Copy the content Link

Share on social media