Predicting Intra-Urban Variation in Air Pollution Concentrations with Complex Spatio-Temporal Dependencies

Abstract

1 Introduction

There is a growing understanding in the literature that exposure to air pollution is associated with adverse health outcomes. The early epidemiological evidence was based on assigning exposures using area-wide monitored concentrations in different geographic regions (Dockery et al. 1993; Pope et al. 2002) or at different times within the same region (Samet et al. 2000). A weakness of area-wide monitoring approaches is that they fail to take advantage of variation between individuals living in the same geographic region. In addition, depending on the study design, there is the potential for unmeasured confounding by region or by time.

More recent cohort studies have assigned individual concentrations based on estimates of intra-urban variations in ambient concentrations. Prediction approaches have included assigning the value measured at the nearest monitor to the participant’s residential location (Miller et al. 2007; Basu et al. 2000; Ritz et al. 2006); using “land use regression” estimates based on Geographic Information System (GIS) covariates (Hoek et al. 2008; Brauer et al. 2003; Jerrett et al. 2005a); and interpolating concentrations by a geostatistical method such as kriging (Jerrett et al. 2005b; Kunzli et al. 2005) or semi-parametric smoothing (Kunzli et al. 2005). These studies and others like them have used relatively simple spatial statistical techniques for exposure assignment based on monitoring data from existing regulatory networks. Our objective is a flexible and practical methodology that accounts for the complex structure of the ambient spatio-temporal concentration field and can take full advantage of regulatory and other monitoring data to more accurately predict concentrations for individual cohort members.

The work described in this paper is motivated by the Multi-Ethnic Study of Atherosclerosis and Air Pollution (MESA Air). MESA Air is a cohort study funded by the U.S. Environmental Protection Agency (EPA) that emphasizes accurate prediction of individual exposures in order to accomplish its primary aim of assessing the relationship between chronic exposure to air pollution and sub-clinical cardiovascular disease. The MESA Air cohort is comprised of 6226 male and female subjects in six major U.S. metropolitan areas (Los Angeles, CA; New York, NY; Chicago, IL; Minneapolis-St. Paul, MN; Winston-Salem, NC; and Baltimore, MD). Although it is possible to estimate health effects based on variations in concentrations between these regions, a major thrust in MESA Air is to develop accurate exposure predictions for individuals that also incorporate intra-urban difference in ambient concentrations in order to reduce exposure misclassification, increase the study power, and obviate possible confounding by region. The primary MESA Air hypotheses relate to exposure to particulate matter of ambient origin with aerodynamic diameter less than 2.5 µm (PM_2.5). Gaseous oxides of nitrogen (NO_x) demonstrate more intra-urban heterogeneity and are also considered as a marker for traffic-related air pollution. In this paper we present examples of modeling outdoor NO_x concentrations in the Los Angeles area, but the statistical methodology is equally applicable to other regions and to PM_2.5 and will ultimately be applied in all of these settings. We also note that final exposure estimates in MESA Air will integrate predictions of outdoor concentrations with additional subject-level data, including time-activity patterns, home infiltration characteristics, address history, and employment address.

A primary source of concentration data for estimating exposures is the EPA’s regulatory Air Quality System (AQS) repository of ambient monitoring data. The AQS network includes a number of fixed site monitors in each region, each of which measures ambient air pollution levels on a regular basis, typically hourly for NO_x and less frequently for PM_2.5. Although there are some missing data and variations between sites, most AQS sites provide nearly complete NO_x concentration time series over several years at their spatial locations. MESA Air is also engaged in a supplementary measurement campaign to provide additional concentration data. The objective of the MESA Air monitoring is to more completely sample a design space that emphasizes traffic-related pollution and to capture data at actual subject home locations. For logistical reasons, the supplementary monitoring data are sampled as two-week averages based on a design with spatio-temporal misalignment that results in significant amounts of missing data at some measurement locations (Cohen et al. 2009).

Although the primary interest in MESA Air and similar cohort studies is in predicting spatial variation of long-term average concentrations to estimate exposures, our statistical modeling approach needs to account for spatio-temporal variability and correlation structures in the data. For an overview of techniques for modeling correlated spatio-temporal data, see Banerjee et al. (2004). A recent paper by Fanshawe et al. (2008) emphasizes the role of carefully chosen covariates in obviating the need to accommodate spatio-temporal correlation in the residuals, but the model in that paper assumes a uniform time trend across locations. Paciorek et al. (2008) and Sahu et al. (2006) model particulate matter using techniques that allow for more complex spatio-temporal dependencies, however their estimation and prediction procedures are applicable only with relatively well aligned monitoring data. Smith et al. (2003) uses an expectation-maximization (EM) algorithm to allow for arbitrary missing data patterns, but their model does not accommodate complex spatio-temporal dependencies. We describe a new modeling and prediction procedure that includes sufficiently complex spatio-temporal dependencies to accurately account for variation in seasonal patterns at different locations and that naturally allows for significant amounts of missing data.

In Section 2, we describe the available AQS monitoring data as well as the sampling pattern for the MESA Air supplementary monitoring campaign. We also describe the geographic covariates that are used in this paper. In Section 3, we specify our hierarchical spatio-temporal model and discuss techniques for efficient estimation. In Section 4, we apply the model to the AQS data from the Los Angeles region and assess the quality of predictions by cross-validation. In Section 5, we conduct a simulation study to assess the added benefit of including data from the supplementary MESA Air monitoring campaign. (The MESA Air supplementary monitoring campaign and quality control process are ongoing, so these concentration measurements are not included in the present paper.) We conclude in Section 6 with a discussion.

2 Description of Data

2.1 AQS Data

The EPA manages the national AQS network of regulatory monitors. Many AQS sites report NO_x concentrations on an ongoing basis, most typically as hourly averages. For this study we are including data from 18 AQS monitors in the Los Angeles region that cover most of the area in which MESA Air cohort members reside. The monitor locations are shown on the map in Figure 1. Because the MESA Air supplementary monitoring is done at the two-week average time scale, we also aggregate the AQS monitoring data to two-week averages. There is a small amount of variability in the number of AQS measurements that contribute to each two-week average, which can result in variable amounts of measurement error. The current model assumes a common variance for the measurement error of all AQS and MESA Air two-week average concentrations. Since the data are skewed, we log-transform two-week average NO_x concentrations. Example time series for the period from July 2005 through December 2007 are shown in Figure 2. The locations of these three sites in Los Angeles, Long Beach, and Pomona are highlighted on the map in Figure 1. Notice that the time series have different mean levels as well as different patterns of seasonal variation.

Open in a separate window

Figure 1

AQS and MESA Air monitoring locations in Los Angeles, and the 200 cohort residence locations used for validation in simulation scenario. (All home locations jittered on map to protect confidentiality)

Open in a separate window

Figure 2

Example time series of log-transformed two-week average NO_x concentrations at three AQS monitors in the Los Angeles area for the period July 2005 through December 2007

3 Model and Estimation

4 Los Angeles AQS Data

We now apply the hierarchical model from Section 3 to make predictions based on two-week average log-transformed NO_x concentrations from 18 monitoring locations of AQS network in the Los Angeles area for the time period from July 2005 through December 2007. This involves estimating the model parameters and then assessing prediction accuracy for long-term average concentrations by means of cross-validation. As previously noted, the locations of the monitors we consider are shown in Figure 1, and several example time series are shown in Figure 2.

Since there are essentially complete time series of data at each location, we can conduct a multistep exploratory analysis that will validate the appropriateness of the hierarchical model described in Section 3. This multi-step analysis gives estimates for all of the parameters in the hierarchical model, and we can compare these to the estimates obtained by the unified maximum-likelihood estimation procedure of Section 3.3.

4.1 Seasonal Trends and Residual Autocorrelation

We follow the methodology described in Section 3.4 to extract m = 2 smooth temporal basis functions that are intended to capture the range of seasonal variation in the region. After determining the first two singular vectors of the observed data matrix, we smooth them using smoothing splines as implemented in the R function smooth.spline() with four degrees of freedom per year (R Development Core Team 2008). The two seasonal trend functions are shown in Figure 3, and an additional basis function that is not shown is f₀(t) ≡ 1. For each AQS monitoring location s, we estimate values $\widehat{\beta }$_0s, $\widehat{\beta }$_1s, and $\widehat{\beta }$_2s by fitting

Y_st = β_0s + β_1sf₁(t) + β_2sf₂(t) + ν_st.

(12)

with ordinary least squares.

Open in a separate window

Figure 3

Smoothed (line) and unsmoothed (points) empirical orthogonal basis functions based on AQS NO_x two-week averages in Los Angeles area (centered and normalized to SD=0.707 for smooth version).

For each location s, the residuals from this linear model constitute a time series. Our objective in using the basis functions to model seasonal variability is to simplify the structure of the residual field ν_st, ideally allowing us to treat it as having no temporal correlation. This corresponds to the residual time series being uncorrelated. Autocorrelation plots are shown in Figure 4 for the residuals at each of the 18 AQS sites. While there is a small amount of variability between sites, these plots taken as a group validate our assumption that there is no temporal correlation at the two-week time scale.

Open in a separate window

Figure 4

Empirical autocorrelation functions for two-week average residuals after fitting to empirical orthogonal basis functions. (18 AQS monitors in Los Angeles area.)

In order to assess the Normality assumption for the residual field, we show smoothed density and Normal Q-Q plots for the combined distribution of residuals from the 18 AQS sites in Figure 5. There is evidence of skewness in the distribution and a heavy left tail, but overall the Normality assumption seems like a reasonable approximation. In principle, our spatio-temporal model could be modified to make the residual distribution closer to Normal, potentially by using a generalized logarithm transformation (Durbin and Rocke 2003) or by formulating the model on a finer time scale (e.g., daily average concentrations). However, either of these approaches would entail a significant increase in the computational burden.

Open in a separate window

Figure 5

Density plot and Normal Q-Q plot for log two-week average residuals after fitting to empirical orthogonal basis functions. (18 AQS monitors in Los Angeles area.)

4.2 $\widehat{\beta }$_is Spatial Fields

For each i = 0, 1, 2 we analyze the estimated spatial field $\widehat{\beta }$_is in terms of its relationship with geographic covariates and its residual spatial structure. The three geographic covariates considered in this analysis are: (i) distance to the coast, (ii) distance to a major road (major road is defined by census feature class code A1–A3, with distance truncated at 300 meters), and (iii) average population density in a 2000 meter buffer. These variables were extracted using GIS, as described in Section 2.3.

The results of separate linear regression model fits for each of the $\widehat{\beta }$_is fields are shown in Table 1. As expected, the estimated long-term averages $\widehat{\beta }$_0s tend to be lower for locations that are farther from roads and higher in areas of higher population density. In addition, the long-term averages tend to be higher further from the coast, which is broadly consistent with the notion that the prevailing wind concentrates pollution on the west side of the coastal mountains in the Los Angeles basin. The two sets of estimated seasonal basis function coefficients $\widehat{\beta }$_1s and $\widehat{\beta }$_2s do not have statistically significant relationships with the roadway or population variables, but both are associated with distance to the coast, indicating that the effect of meteorology varies by geography in this region.

Table 1

Linear regression parameter estimates (standard errors) for empirical orthogonal basis function coefficients of AQS monitors in Los Angeles area. Distance to nearest A1, A2, or A3 class road truncated at 300 meters. Population density based on a buffer with 2000 meter radius.

	β0		β1		β2
Intercept	3.05	(0.22)	1.10	(0.17)	−0.39	(0.14)
Distance to road (km)	−1.37	(0.57)	−0.62	(0.43)	−0.19	(0.35)
Distance to coast (km)	0.012	(0.003)	−0.009	(0.002)	0.006	(0.002)
Population density (per m2)	41.04	(10.76)	−7.89	(8.24)	9.81	(6.67)

Open in a separate window

We expect there to be spatial correlation in the $\widehat{\beta }$_is fields. To assess the degree of spatial correlation, we examine empirical variograms for the residuals from regression on the spatial covariates. For $\widehat{\beta }$_1s and $\widehat{\beta }$_2s, only the distance to coast covariate is included in the regression since the other two covariates do not appear to be important. Empirical variogram clouds and binned values calculated using the GeoR package in R (Ribeiro and Diggle 2001) are shown in Figure 6. The variograms suggest that there is significant spatial correlation in $\widehat{\beta }$_0s, modest correlation in $\widehat{\beta }$_1s, but limited correlation in $\widehat{\beta }$_2s. We quantify this by fitting universal kriging models with exponential variograms and no nugget to each of the three fields using the likfit() function in GeoR. The resulting parameter estimates are shown in the first column of Table 2.

Open in a separate window

Figure 6

Empirical variograms for the estimated spatial fields of long-term averages ($\widehat{\beta }$_0s) and coefficients of seasonal basis functions ($\widehat{\beta }$_1s, $\widehat{\beta }$_2s). The black curve represents a classical variogram estimate, and the red curve is derived using the robust modulus method. (18 AQS monitors in Los Angeles area.)

Table 2

Hierarchical model parameter estimates for AQS data in Los Angeles area. The values in the first column are obtained by first estimating the β_is at each site s and then fitting a universal kriging model to each estimated spatial field $\widehat{\beta }$_is. An estimated field of residuals $\widehat{\nu }$_st is then derived and fit to a separate ordinary kriging model to get the parameter estimates for ν. The values in the second column are obtained by a single step maximum-likelihood estimation for the full model as described in Section 3.3.

		Multi-step estimation	Full maximum-likelihood
β0s	Intercept	3.04	3.04
	Distance to road (km)	−1.63	−1.62
	Distance to coast (km)	0.013	0.013
	Population density (m−1)	34.5	34.3
	Sill	0.049	0.049
	Range (km)	27.1	27.8
β1s	Intercept	0.96	0.97
	Distance to coast (km)	−0.009	−0.010
	Sill	0.031	0.032
	Range (km)	19.73	22.2
β2s	Intercept	−0.21	−0.21
	Distance to coast (km)	0.0052	0.0053
	Sill	0.017	0.017
	Range (km)	4.15	5.05
νst	Sill	0.031	0.033
	Range (km)	107.11	106.11
	Nugget	0.0095	0.0100

Open in a separate window

4.5 Prediction Accuracy

Using the parameter estimates derived above, it is straightforward to follow the procedure in Section 3.1 and predict long-term average concentrations at locations where measurements are not available. This amounts to predicting log-concentration values as the conditional mean of the latent Gaussian random field at locations without data, and then back-transforming to obtain concentrations as in equation (4).

We assess the accuracy of the prediction model by leave-one-out cross-validation. For each $\bar{s}$ corresponding to the spatial location of one of the 18 AQS monitoring sites, we predict the long-term average concentration Ĉs̄,lta* by applying the model from Section 3 as above, replacing Y_st by the set of observations that excludes all measurements at location $\bar{s}$

Y_st/$\bar{s}$ = {Y_st; s ≠ $\bar{s}$}.

This procedure yields a cross-validated spatial field of predicted long-term average concentrations Ĉs,lta*, where s ranges over the 18 AQS monitoring locations. A scatter-plot of cross-validated pre-dictions is shown in Figure 7. The plot suggests that the model fits well since there are no noticeable outliers. The root mean square error (RMSE) is 4.21 parts per billion (ppb), corresponding to an R² of 0.67. The formula used to compute R² is

R2=max(0,1−RMSE2/Var(Cs,lta*)),

(13)

where Cs,lta* is the spatial field of true long-term average concentrations that could be defined by replacing Ŷst* by Yst* in equation (4).

Open in a separate window

Figure 7

Cross-validated predictions of long-term average NO_x concentrations for 18 AQS monitors in Los Angeles area. The RMSE is 4.21 and the R² is 0.67. The formula used to compute R² is given in Section 4.5.

5 Simulation Study in Los Angeles

In order to evaluate the added value of the MESA Air supplementary monitoring campaign, we conducted a simulation study based on sampling Y_st at all of the AQS and MESA Air sites described in Section 2 (N = 2011, n = 346). We simulated log-transformed two-week average concentrations using the hierarchical model of Section 3 with temporal basis functions and parameters as estimated from the AQS data in Section 4. We also simulated time series of concentrations at 200 additional randomly selected subject home addresses. The approximate locations are shown in Figure 1. We calculated the long-term average concentration values at each home address from July 2005 through December 2007 and regard these as a validation set Cs,lta* for evaluating prediction performance.

We simulated 48 random realizations and then estimated the parameters using maximum-likelihood as in Section 3.3. Mean estimated parameter values are shown in the second column of Table 3 along with the standard deviation of estimates and the mean standard errors (calculated based on the Hessian of the likelihood function). The estimates are generally very close to the assumed values (first column), although we note that the range and sill tend to be underestimated for the β_0s field. The sampling standard deviations are also reasonably well aligned with the standard error estimates. The standard errors are calculated by inverting the Hessian of the likelihood function at its mode. In future work, we will investigate alternative approaches such as Bayesian and empirical Bayesian calculations that may provide improved standard error estimates. The third column contains analogous values when the maximum-likelihood estimation is based only on simulated observations at the AQS and MESA Air “fixed sites” (N = 1343, n = 25). The mean values are very close to the ones obtained by using all monitoring locations, with a modest amount of additional variability. This suggests that the additional sampling at the “home outdoor” and “community snapshot” locations adds little value for parameter estimation.

However, our interest is in predicting the long-term average concentrations at subject home locations (Cs,lta*), not just in estimating the model parameters. Scatterplots of predicted values for the first twelve Monte Carlo simulations in each of two scenarios are shown: (i) using the AQS and all MESA Air monitoring sites in Figure 8, and (ii) using only the AQS and MESA Air “fixed sites” in Figure 9. The predictions obtained by using information from all of the MESA Air monitoring sites are better than those using only the MESA Air “fixed sites”, with a lower average RMSE (4.02 compared to 6.24), and a higher average R² (0.95 compare to 0.88). Since the parameter estimates are very close, this result suggests that there is significant benefit from having additional monitoring locations for the prediction step of equation (3), even though the measurements at the MESA Air “home outdoor” and “community snapshot” sites are temporally sparse.

Open in a separate window

Figure 8

Simulation study results for the first twelve Monte-Carlo realizations. Scatter plots of predicted vs. true long-term average NO_x concentrations at 200 subject homes in validation set. Results based on using all AQS and MESA Air monitoring locations for parameter estimation and prediction.

Open in a separate window

Figure 9

Simulation study results for the first twelve Monte-Carlo realizations. Scatter plots of predicted vs. true long-term average NO_x concentrations at 200 subject homes in validation set. Results based on using only AQS and MESA Air “fixed sites” for parameter estimation and prediction.

To further evaluate this potential benefit, we consider two additional prediction scenarios: (iii) all of the monitoring locations are used for parameter estimation but only the AQS and MESA Air “fixed sites” are used for prediction, and (iv) only the AQS and MESA Air “fixed sites” are used for parameter estimation but all of the monitoring locations are used for prediction. The results are shown in the red curves of Figures 10 and ​and11,11, and they validate the hypothesis that for this simulation scenario the primary benefit from including the “home outdoor” and “community snapshot” monitoring campaigns is for prediction rather than for parameter estimation. We also show analogous results for the incremental value of adding the “home outdoor” or “community snapshot” monitoring locations (green and blue curves, respectively, in Figures 10 and ​and11).11). The same pattern persists, with incremental benefit from each set of additional monitoring locations.

Open in a separate window

Figure 10

Simulation study results for 48 Monte-Carlo realizations. Average root mean squared error for predicted vs. true long-term average NO_x concentrations at 200 subject homes in validation set. Results based on using different subsets of the AQS and MESA Air monitoring locations for parameter estimation and prediction in the spatio-temporal hierarchical model.

Open in a separate window

Figure 11

Simulation study results for 48 Monte-Carlo realizations. Average R² for predicted vs. true long-term average NO_x concentrations at 200 subject homes in validation set. Results based on using different subsets of the AQS and MESA Air monitoring locations for parameter estimation and prediction in the spatio-temporal hierarchical model. The formula used to compute R² is given in Section 4.5.

We note that this simulation study demonstrates the benefit of having additional sampling from the “home outdoor” and “community snapshot” campaigns in the prediction step. In addition to the simulation study findings, there may turn out to be additional benefit for parameter estimation in the actual MESA Air study. This will be determined by the final choice of geographic covariates and how well the various sampling campaigns span the range of values for those covariates. The “community snapshot” campaign was specifically designed to capture near-road effects, and as we develop more refined traffic-related covariates we expect the data from this sub-campaign to be particularly important for estimating the relevant regression coefficients (Wilton et al. 2008).

In addition to the simulation scenario described above, we have also considered a second set of simulations in which we assume that there is greater variability in the temporal trend at a fine spatial scale. This additional variability is induced by allowing the β_1s and β_2s fields to be dependent on the distance to a major road covariate. The results are very similar, so we do not show them here.

6 Discussion

The methodology described in this paper has two features that make it attractive for exposure prediction in environmental epidemiology, and more generally for applications that benefit from accurate prediction using spatio-temporal data. First, our model has a very flexible correlation structure that allows for non-separability of space and time by modeling seasonal and long-term trends using empirical orthogonal basis functions with spatially correlated random fields of coefficients. Second, the unified estimation approach can be implemented with standard software and accommodates arbitrary missing data patterns, as long as there are sufficiently rich time series at a subset of locations to derive temporal basis functions.

An important consideration in implementing this model is determining what covariates are helpful for modeling the spatial fields of temporal basis function coefficients (recall that one of the temporal basis functions is the constant function, representing the long-term average). In this paper we have restricted our attention to a relatively small number of covariates that are easy to calculate using GIS. However, as part of our work on the MESA Air project, we are developing a more comprehensive set of spatial covariates, and we expect that these will be very valuable for making predictions at subject homes. In particular, we are investigating more complex covariates to account for the influence of local traffic patterns. This includes estimating actual traffic densities on the road network and incorporating the results of meteorology through physics-based plume modeling. Preliminary exploratory analysis in the Los Angeles region indicates that accounting for meteorology and traffic patterns will contribute significantly to improved predictions (Wilton et al. 2008).

One of our key findings in the simulation study is that the temporally sparse components of the MESA Air monitoring campaign (“home outdoor” and “community snapshot”) contribute to improved predictions, but that this improvement is primarily through better interpolation in the prediction step and not through improved parameter estimation. The details of this finding are likely connected to the choice of geographic covariates. For example, when we use more refined covariates to represent traffic density, we expect that the “community snapshot” monitoring will prove important for accurate estimation of the regression coefficients for these covariates. This is because the “community snapshot” monitoring includes extensive sampling at gradients within a few hundred meters of major roads in multiple directions corresponding to up- and down-wind locations. Thus, our expectation is that the benefit of MESA Air monitoring data in the prediction step will persist when we use additional covariates, and that additional gains in prediction accuracy will be realized through more accurate parameter estimation based on these data.

We have focused on making point predictions and evaluating the accuracy of these predictions compared to true values. It would be straightforward to also calculate a prediction variance at each location, taking a similar approach to Fanshawe et al. (2008). We do not pursue this here, however, because separate uncertainty estimates for each location are not helpful if the objective is to use the predicted concentrations to estimate the health effect in an environmental epidemiology study. To obtain valid health effect standard errors, we need to properly account for the uncertainty in the exposure prediction surface taken as a whole, including correlations between locations. In a future paper we will adapt recently developed measurement error correction procedures to accomplish this objective (Szpiro et al. 2008). Standard errors for the spatio-temporal model parameter estimates like those reported in Table 3 play an important role in the applicable measurement error correction procedures.

Footnotes

References