- Email: [email protected]
Spatiotemporal modeling with temporal-invariant variogram subgroups to estimate fine particulate matter PM2.5 concentrations
Atmospheric Environment 54 (2012) 1e8
Contents lists available at SciVerse ScienceDirect
Atmospheric Environment journal homepage: www.elsevier.com/locate/atmosenv
Spatiotemporal modeling with temporal-invariant variogram subgroups to estimate fine particulate matter PM2.5 concentrations Chu-Chih Chen a, *, Chang-Fu Wu b,1, Hwa-Lung Yu c, Chang-Chuan Chan b, Tsun-Jen Cheng b a
Division of Biostatistics and Bioinformatics, Institute of Population Health Sciences, National Health Research Institutes, Taiwan Institute of Occupational Medicine and Industrial Hygiene, College of Public Health, National Taiwan University, Taiwan c Department of Bioenvironmental Systems Engineering, National Taiwan University, Taiwan b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 7 July 2011 Received in revised form 12 December 2011 Accepted 6 February 2012
Short-term exposure estimation of daily air pollution levels incorporating geographic information system (GIS) into spatiotemporal modeling remains a great challenge for assessing corresponding acute adverse health effects. Due to daily meteorological effects on the dispersion of pollutants, explanatory spatial covariables and their coefficients may not be the same as in classical land-use regression (LUR) modeling for long-term exposure. In this paper, we propose a two-stage spatiotemporal model for daily fine particulate matter (PM2.5) concentration prediction: first, daily nonlinear temporal trends are estimated through a generalized additive model, and second, GIS covariates are used to predict spatial variation in the temporal trend-removed residuals. To account for spatial dependence on meteorological conditions, the dates of the study period are divided by the sill of the daily empirical variogram into approximately temporal-invariant subgroups. Within each subgroup, daily PM2.5 estimations are obtained by combining the temporal and spatial parts of the estimations from the two stages. The proposed method is applied to the modeling of spatiotemporal PM2.5 concentrations observed at 18 ambient air monitoring stations in Taipei metropolitan area during 2006e2008. The results showed that the PM2.5 concentrations decreased whereas the relative humidity and wind speed increased with the sill subgroups, which may be due to the effects of daily meteorological conditions on the dispersions of the particles. Also, the covariates and their coefficients of the LUR models varied with subgroups and had in general higher adjusted R-squares and smaller root mean square errors in prediction than those of a single overall LUR model. 2012 Elsevier Ltd. All rights reserved.
Keywords: Empirical variogram Geographic information system Kriging Land-use regression Temporal trend
1. Introduction Numerous epidemiological and toxicological studies have shown that increased particulate matter (PM) is associated with increased cardiovascular mortality and morbidity and that fine PM particles with aerodynamic diameter 2.5 mm (PM2.5) are more hazardous than coarse particles, such as PM10 (Cheng et al., 2003; Chuang et al., 2005; Schwartz, 2004; Zanobetti and Schwartz, 2009). Many of the studies relate mortality or hospital emergency visit data to mean air pollution level measured at a central site or from several monitoring stations. However, PM2.5 concentrations are subject to within-city or intra-urban variations (Jerrett et al., 2005; Wilson et al., 2005), and individuals may be exposed to Abbreviations: GIS, geographic information system; LUR, land-use regression; OK, ordinary kriging; PM, particulate matter; RMSE, root mean square error; UK, universal kriging. * Corresponding author. E-mail addresses: [email protected] (C.-C. Chen), [email protected] (C.-F. Wu). 1 Equal contribution with the correspondence. 1352-2310/$ e see front matter 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.atmosenv.2012.02.015
different levels depending on their residential area. To estimate small-area air pollution concentrations, in addition to using geostatistical interpolation based on sampled site measurements, land-use regression (LUR) models using data on nearby traffic and land-use covariates have been successfully applied to explain spatial variation (Clougherty et al., 2008; Henderson et al., 2007; Hoek et al., 2008; Ross et al., 2007; Ryan and LeMasters, 2007). Typically such models are built by regressing annual averages of pollutant concentrations over geographical information system (GIS) variables that are invariant to temporal variations. As a result, standard LUR models are incapable of studying acute health effects due to short-term exposures. Furthermore, temporal variations in PM concentrations due to meteorological conditions on dispersion effects or weekdays also need to be considered to relate to individual acute adverse effects. Therefore, a spatiotemporal model to predict PM concentrations at unsampled locations on a specific date is required for individual exposure estimation. Using a semiparametric latent variable regression model, Gryparis et al. (2007) extend geoadditive models (Kammann and
2
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
Wand, 2003) by incorporating temporal covariables for spatiotemporal modeling of mobile particle concentrations. Maynard et al. (2007) use corresponding short-term exposure modeling to study associated mortality risks considering GIS-based measures of cumulative traffic density. Similarly, Alexeeff et al. (2011) use the same modeling to estimate subjects’ daily black carbon exposure to relate to their repeated inflammation markers. However, the approach specifies a lack of spaceetime interaction over time across locations, which may not be valid due to the fact that spatial dispersion of air pollutants is temporal-dependent (Bogaert et al., 2009; Yu and Wang, 2010; Yu et al., 2011). To separate temporal and spatial variations, modeling based on two-stage approach by separating temporal trend from withinregion spatial variation in the mean structure have been proposed (Dadvand et al., 2011; Fanshawe et al., 2008; Yanosky et al., 2009). However, additional information about spatiotemporal variability and correlation structures in the data need to be considered (Banerjee et al., 2004; Bogaert et al., 2009; Yu and Wang, 2010). To deal with the complex spatial and temporal covariance dependence, Bogaert et al. (2009) divide and estimate the spatial variogram of ozone distribution by month. Alternatively, Szpiro et al. (2010) decompose a hierarchical model incorporating dependence on GIS covariates with temporal trends and a residual with spatial correlation structure. However, it may be computationally difficult to extend both of these spatiotemporal models to predict daily air pollution. In this study, we adopt a two-stage approach: (1) in the first stage, the daily mean temporal trend of the (log-transformed) PM2.5 concentrations of the study area is estimated by using a generalized additive model (GAM); (2) in the second stage, LUR covariables are incorporated for spatial prediction of the (temporal trend removed) residuals from the first stage. Because temporal dependence may still be present, in stage 2, the dates of the study period were classified into subgroups based on the sill of the daily empirical semivariogram of the residuals from stage 1. Dates falling within the same group were then aggregated together to obtain the time-averaged residuals and a LUR model was applied separately, instead of an overall LUR model. Kriging method was further applied to the residuals of the LUR models in the presence of spatial autocorrelations. We then predicted the final concentration level at a specific date by restoring the temporal trend to the spatial prediction. Daily averaged PM2.5 concentrations from 18 ambient air monitoring stations in the Taipei metropolitan area between 2006 and 2008 were analyzed as an example.
2. Materials 2.1. PM2.5 and meteorological measurements Air pollution levels were obtained from the Air Quality Monitoring (AQM) Network of the Environmental Protection Administration (EPA) of Taiwan, from the 18 monitoring stations that are located within the Taipei metropolitan area, including Taipei City and New Taipei City within the Taipei Basin (Fig. 1). Daily averages of PM2.5, wind speed, temperature, and relative humidity were calculated based on the available monitoring data from January 1, 2006, to December 31, 2008. 2.2. GIS data GIS-based data for the 18 AQM stations fall into three categories: (1) land use (residential, commercial, urban green, forest, industry, and open space) within buffers of several radii (100, 300, 500, and 1000 m); (2) types and lengths of roads (all roads, highways, and major roads) within buffer zone of 25, 50, 100, 300, 500, and 1000 m; (3) population size within radii of 100, 300, 500, and 1000 m. Also, distance to the nearest major roads is included as a covariate. The land-use database was obtained from the National Land Survey and Mapping Center of Taiwan. Road network information was acquired from the Institute of Transportation, Ministry of Transportation and Communications, and census population data were accessed from the Civil Affairs Bureau of Taipei City and New Taipei City, Taiwan. ArcGIS 9.3 was employed to calculate the lengths of roads and population size within each buffer zone of the 18 AQM stations. 3. Statistical method 3.1. Two-stage spatiotemporal modeling While PM2.5 concentration levels may change nonlinearly with temperature and relative humidity, the meteorological effect on the dispersion of particles is expected to be uniform across a mediumsized study region such as the Taipei metropolitan area (approximately 2325 km2). This implicit assumption suggests that for assessing pollutant levels of a metropolitan area, spatial and temporal variations in concentrations may be decomposed separately. Suppose the log-transformed measurement of monitoring
Fig. 1. Map of the Taipei metropolitan area and the locations of the 18 AQM stations.
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
station i at location si and time t is y(si,t), i ¼ 1, ., K; t ¼ 1, ., T, a two-stage spatiotemporal modeling method is proposed as follows. 3.1.1. Stage 1. Modeling for daily mean temporal trend Let 0
YðtÞ ¼ x ðtÞb þ
H X
3.2. Grouping of dates based on the daily empirical semivariogram For each day, an empirical semivariogram was constructed based on the residuals of the log-transformed daily average PM2.5 concentrations of the 18 monitoring stations (Eq. (2)) as follows:
g^ k ðtÞ ¼ fh ðqh ðtÞÞ þ ZðtÞ
(1)
h¼1
be the mean temporal GAM model on day t, where YðtÞ is the average of the K monitoring station measurements Y(i,t) on day t, x(t) is a vector of the linear or categorical temporal predictors such as weekdays and months, b is the corresponding vector of coefficients, and fh(qh(t)), h ¼ 1, ., H is a nonlinear function of the predictors q1, ., qH on day t. The correlations of model residuals Z(t) may be assumed to follow an autoregressive process AR(p) of order p. Let the trend-removed residuals at station si be
^ dðsi ; tÞ ¼ Yðsi ; tÞ YðtÞ;
(2)
^ where YðtÞ is the mean estimate on day t. 3.1.2. Stage 2. Land-use regression with spatial interpolation Given a fixed date t, a spatial model is applied to the residuals d(si,t) obtained from stage 1 with GIS information as explanatory covariates. Because the spatial covariance of the residuals d(si,t) may still change with time, spatial interpolation using the kriging method may vary with different dates. Also, while GIS-based measures are invariant with time, the significant LUR model covariates and the corresponding coefficients may change with daily dispersion conditions. It is thus reasonable to aggregate the residuals across different days with similar spatial correlations as a group. The procedure is accomplished by dividing the sills of the daily empirical variograms into approximately homogeneous groups as described in the next subsection. Within each group, LUR covariables are then incorporated (with kriging in the presence of spatial autocorrelation) for spatial prediction of the subgroupaveraged residuals. The LUR covariables are selected following the procedure given in Vienneau et al. (2010): 1) the increase in adjusted R2 of the predictor variable greater than 1% was selected; 2) the coefficient conformed to the pre-specified direction of effect; 3) the direction of effect for predictors already included in the model did not change; and 4) the p-value for all variables already in the model remained below 0.1. When the spatial correlations among the residuals are negligible, the LUR model may be written as
dl ðsi Þ ¼ w0l ðsi Þ4l þ eðsi Þ;
3
1 2NBk
X
dðsi ; tÞ d sj ; t
2
;
(5)
fði;jÞ:ksi sj k˛Bk g
where NBk is the number of pairs of sites with distance fall within the bin Bk. Ten bins were divided with 15 pairs of sites each for each day. Because the daily meteorological conditions may vary substantially from day to day, the effects on dispersions of the PM2.5 particles should also be different. The dispersions may further be affected by the differential PM2.5 concentrations across the monitoring sites (e.g., traffic sites versus rural sites). For example, spatial correlation of PM2.5 concentrations is expected to be smaller and the corresponding variogram is thus larger in rainy days because of cleaner air. Similarly, higher wind speed could keep the atmosphere mixed better. As a result, spatial correlation also tends to be smaller. ^ k ðtÞ may Therefore, the daily empirical semivariogram values g change accordingly, which may be assessed by the sill estimated by the daily maximum empirical semivariogram, i.e.,
g^ max ðtÞ ¼ max g^ k ðtÞ: k
(6)
The corresponding dates were then divided into subgroups by ^ max ðtÞ from the smallest to the largest. For the dates sorting g belonging to the same group, a LUR model was fitted separately. As a comparison, UK with a variogram model fitting within each of the subgroups was also applied for prediction. 4. Results The averaged daily PM2.5 concentrations across the 18 AQM stations of the Taipei metropolitan area had a mean 28.35 mg m3, ranged from 0.25 to 128.20 mg m3 with an inter-quartile range of 20.04 mg m3 during the study period 2006e2008. Fig. 2 shows the monthly PM2.5 concentrations of the AQM stations at the general sites (n ¼ 13), traffic sites (n ¼ 3) and the background sites (n ¼ 2), respectively. The results showed that the monthly mean PM2.5 concentrations were in general higher during the cold months (JanuaryeApril) compared to those of the warm months (JuneeOctober). Site effects of the AQM stations were also apparent, with the mean concentrations of the traffic sites the highest, followed by those of the general sites. Although the background stations had the lowest mean PM2.5 concentrations, the values
(3)
where dl ðsi Þ is the averaged residual of d(si,t) with t ˛ Tl, Tl is the l-th set of dates with similar variograms, l ¼ 1, ., L, wl(si) is a vector of the LUR covariables at site si, 4l is the corresponding vector of coefficients, and e(si) are assumed to be independent, zero-mean, and a common variance of s2e . In the presence of spatial correlation, the resultant model for universal kriging (UK, LUR þ kriging) becomes
dl ðsi Þ ¼ w0l ðsi Þ4l þ Sðsi Þ þ eðsi Þ;
(4)
where {S(s):s ˛ R2} is a stationary zero-mean stochastic process with spatial covariance matrix C, and e(si) is the nugget effect with spatially-independent variance s2e . Finally, given a day t ˛ Tl, a prediction of the log-PM2.5 concentration at location s0 can be obtained from the results of Eqs. (1)e(4).
Fig. 2. Monthly mean PM2.5 concentrations of the general, traffic, and background stations of the Taipei metropolitan area from January 2006 to December 2008.
4
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
could also reach 25 mg m3 for certain months, which could be due to the long-range Asian sand storms transported from outside of Taiwan (Chan et al. 2008). 4.1. Site-specific annual-average PM2.5 predictions using simple LUR and kriging method By aggregating the daily average PM2.5 concentrations of the 18 AQM stations from 2006 to 2008, a simple annual LUR model was obtained using the GIS-based measures. Of the 48 explanatory variables, only the three covariables: all roads within 50 m, commercial area within 1000 m, and distance to the nearest major roads remained in the final model. The resultant LUR model and ordinary kriging (OK) were applied separately to predict the mean PM2.5 concentrations of the 18 AQM stations for cross-validation. Because the number of monitoring stations was barely enough to yield stable variogram model parameter estimates, the empirical variogram was generated by overlaying daily estimates from the 1096 days together, assuming no spatialetemporal interaction. The cross-validation adjusted R-square of OK was substantially lower than that of LUR (0.43 versus 0.75). Also, the adjusted R-square 0.68 of UK was lower than that of LUR alone, which could be due to additional variation from the variogram model fitting (Eq. (4)) with non-significant residual spatial autocorrelation (Moran’s I statistic p-value 0.28). The results are summarized together with those of the two-stage models in Table 2 below. 4.2. Temporal trend estimation and grouping of the dates A GAM model was fitted to the daily average log-PM2.5 concentrations across the 18 AQM stations using the SAS procedure PROC GAM. Daily mean temperature, relative humidity, and wind speed all had significant nonlinear effects. Month also had a significant effect on the PM concentrations. An autocorrelation model of lag 2 (AR(2)) was fitted to the residuals based on the results of a DubineWatson test. Fig. 3 displays the model residuals of the average log-PM2.5 concentrations along the dates of the study period. To illustrate the differences across various sites, mean residuals of the general, traffic, and background monitoring stations were plotted separately for every 7 days. As shown in the figure, the mean residuals of the general stations clustered around the line 0, whereas the mean residuals of the traffic sites and the background sites tended to be positive (underestimated) and negative (overestimated), respectively. The pattern is expected because of traffic-related PM2.5 concentration variations across different sites. However, the mean residuals of different sites appeared to be temporally invariant across the dates. Therefore, the assumption that the temporal and spatial variations may be decomposed separately at two stages held approximately legitimate. After the residuals were obtained following Eq. (2), daily ^ max ðtÞ for the sill estimate maximal empirical semivariograms g were calculated. A separate GAM model was applied to the daily g^ max ðtÞ to examine the temporal effects on the spatial variances. The result showed that the daily mean temperature, relative humidity, and wind speed all had significant effects on the daily g^ max ðtÞ (Fig. 4). Because of the dependence of the residuals on meteorological ^ max ðtÞ were divided into 5 subgroups so that conditions, the sorted g the variograms within each subgroup were approximately homo^ max ðtÞ was much geneous. Due to the fact that the distribution of g skewed to the left, we chose to classify the subgroups with the cutoffs of the percentiles at 60%, 70%, 80%, and 90%, respectively. Corresponding dates within the same subgroup were aggregated for separate LUR model fitting at the second stage. Fig. 5 displays
Fig. 3. Mean temporal trend-removed residuals of the log-PM2.5 concentrations across the dates of the study period. The mean residuals of the general, traffic, and background stations are plotted for every 7 days and are denoted by “”, “D”, and “”, respectively.
^ max ðtÞ. The medians of the nonthe box plots of the 5 subgroups of g ^ max ðtÞ indicated inter-group heterogeneity overlapping subgroup g of spatial variances. Notice that there were some estimates of g^ max ðtÞ in the fifth subgroup approximately 10 times greater than those of the first subgroup, which indicated that the daily empirical variograms were quite heterogeneous. 4.3. Subgroup LUR model fitting and cross-validations of different prediction methods Table 1 summarizes the means and standard deviations of PM2.5 concentration, atmospheric pressure, temperature, relative humidity, and wind speed of the dates corresponding to the 5 subgroups. Group differences of these variables showed that, except for the atmospheric pressure, all the other variables were highly significant with a p-value < 0.0001 using an analysis of variance (ANOVA) test. The PM2.5 concentration and temperature decreased with the subgroups, whereas the relative humidity, wind speed increased. Tests for trend of these variables were also all highly significant (p-value < 0.0001). Therefore, the grouping of the ^ max ðtÞ effectively classified the interrelationship between dates by g the daily mean PM2.5 concentrations and the mean temperature, relative humidity, and wind speed. A LUR model was fitted to each of the sill subgroups separately. The resultant models are summarized in Table 2. Except for group 2, all the other groups had covariates included different from those of the overall LUR model. Similar to the overall model, the covariables of the length of all roads within 50 m and distance to the nearest major roads were included in all subgroup LUR models, which indicated that traffic activity was the major predictor. Commercial area within 1000 m was included in all except for subgroup 1. Length of major roads within 300 m was included in subgroups 2 and 3, and open space within 1000 m was included in subgroup 5 only. All of the other covariates, such as the portion of green space and population size at various buffers, were non-significant. Adjacent subgroups tended to have similar covariables because of similarities in the empirical sill. Notice that for the same covariable included in the subgroup LUR models, the coefficients tended to increase with the groups (for example, the length of all types of
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
5
Fig. 4. Semiparametric spline curves of the sill of the daily empirical variograms of the temporal trend-removed residuals with the predictors (a) temperature, (b) relative humidity, and (c) wind speed.
roads within 50 m, commercial area within 1000 m, and distance to the nearest major roads), which was consistent with the scales of g^ max ðtÞ of the subgroups. In addition to the LUR model fittings for each of the subgroups, Table 2 lists the adjusted R-squares and RMSEs of the LUR model building and cross-validations. Because of much higher variation in estimating daily concentrations than the average 3-year level, the subgroup adjusted R-squares were lower
Fig. 5. Box plots of the sills of the daily empirical semivariograms of the temporal trend-removed residuals for each of the 5 subgroups.
than those of the overall model for annual predictions. In comparison with the results using the same overall LUR model throughout the 5 subgroups, the subgroup LUR models in general had higher adjusted R-squares and lower RMSEs. Together with the outcomes shown in Table 1, these results support the hypothesis that meteorological conditions had an effect on daily spatial associations of PM2.5 concentrations as stated in Section 3.2. Also, the adopted procedure by fitting a subgroup LUR model separately effectively improved over a single overall LUR model, especially for the dates with larger scale in the sill of the daily empirical variogram. Similar as before, UK, LUR, and OK were compared for their performance in daily PM2.5 concentration predictions. The adjusted R-squares ranged from 0.39 to 0.59 for UK, 0.39 to 0.59 for LUR only, and 0.09 to 0.34 for OK, respectively. Notice that the adjusted Rsquares increased from group 1 to group 5 with increasing RMSEs for all three methods because the subgroups were divided by increasing spatial variance, resulting in increasing biases. Similar to the overall model, the OK method had the lowest R-squares and highest RMSE for all 5 subgroups. Possibly due to extra variation in variogram model fitting with non-significant spatial associations (Moran’s I statistic p-values were 0.82, 0.82, 0.81, 0.69, and 0.34, respectively), the results of UK were essentially the same as those of LUR within each subgroup. However, all the three methods at stage 2 had higher R-squares compared with those methods using estimates based only on mean temporal trend (the last column of Table 2), which indicated the existence of spatial variations in the temporal trend-removed residuals. To assess the prediction performance at different sites, Table 3 further compared the AQM stations at the general, traffic, and background sites. The results show that the proposed procedure had the best performance at the traffic sites, followed by the
6
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
Table 1 Characteristics of the 5 subgroups: number of days, PM2.5, atmospheric pressure, temperature, relative humidity, and wind speed. Group
N (days)
PM2.5 (mg m3) Mean
std
Mean
std
Mean
std
Mean
std
Mean
std
1 2 3 4 5
656 110 109 109 110
32.9 25.2 23.0 19.9 18.0
14.2 13.1 13.7 12.1 11.3
1012.5 1011.2 1012.0 1013.1 1013.6
6.9 7.3 6.3 6.7 6.5
23.3 22.1 22.3 21.8 22.1
4.9 5.0 5.2 5.0 5.3
71.8 77.3 79.7 80.7 80.6
6.9 6.6 6.4 6.7 8.2
2.24 2.25 2.20 2.45 2.49
0.59 0.71 0.76 0.81 1.06
Temperature ( C)
Atmospheric pressure (hPa)
general sites. Sites with background concentrations, or national parks in general, had relatively lower adjusted R-squares than those of the other sites. 5. Discussions We have created a two-stage spatiotemporal model to predict daily PM2.5 concentrations by first removing the mean temporal trend and then accounting for spatial variations in the residuals. Because of the dependence of pollutant dispersion on daily meteorological conditions, the dates of the study period were divided by ^ max ðtÞ of daily empirical semivariogram into the sill estimate g subgroups that were approximately temporally invariant. This approach allows us to incorporate GIS covariables for spatial associations of the temporal trend-removed residuals at the second stage without an additional assumption of temporal independence. In contrast with a constant LUR model for long-term exposure estimation, subgroup LUR models with different covariates and coefficient estimates were obtained as a result, which had in general better performance in adjusted R-squares and RMSE for prediction compared to those of a single overall model. ^ max ðtÞ of the temporal trend-removed Our results showed that g residuals depended on mean temperature, wind speed, and relative humidity (Fig. 5). The results given in Table 1 also demonstrate that ^ max ðtÞ effectively classified the interrelationship the grouping of g between the PM2.5 concentrations and the daily meteorological conditions. Thus, the assumption of lacking a spaceetime
Relative humidity (%)
Wind speed (m s1)
interaction with log-scale daily air pollution concentrations (Gryparis et al., 2007; Maynard et al., 2007) does not hold. Although similar two-stage modeling has been proposed previously (Dadvand et al., 2011; Fanshawe et al., 2008; Yanosky et al., 2009), the ambient pollutant concentration predictions of these works are mainly based on weekly or monthly time scales. Therefore, there is essentially no need to adjust for significant terms and corresponding coefficients of GIS-based covariates on daily meteorological effects. The model proposed by Szpiro et al. (2010), though appealing with temporal basis functions approach, lacks an intuitive interpretation, and the computational burden becomes cumbersome when the time scale is for daily mean concentrations. We have shown in Table 2 that the overall and subgroup crossvalidation R-squares using LUR only were slightly higher than those of UK and substantially higher than those of OK. This indicates that the spatial associations among the temporal trend-removed residuals at different sites could be essentially explained by the GIS covariates without further spatial interpolations, as was supported by the non-significant Moran’s I test statistics for each of the subgroups. The results given in Table 3 also indicate that the PM2.5 concentration predictions may perform best at sites near heavy traffic. However, for rural or background sites with little to no traffic, the predictions may not be accurate. Possible explanations of this phenomenon are that the covariates of the subgroup LUR models were all essentially traffic-based and that there were only two background sites in our study. Therefore, estimates that were largely based on the general or the traffic sites tended to
Table 2 LUR models for the full study period and each of the 5 subgroups, performance statistics for model building, and cross-validation with different methods. The rows below each of the subgroup LUR models are the corresponding results by fitting the overall LUR model. Group
LUR covariables
Model building 2
Adj R
Overall
1
2
3
4
5
a
Annual All roads 50 m, commercial 1000 m, distance to the nearest major roads (1.69, 0.00083, 0.048)a Empirical variogram subgroups All road 50 m, distance to the nearest major roads (1.16, 0.037) Overall All roads 50 m, commercial 1000 m, distance to the nearest major roads (1.85, 0.00098, 0.054) Overall All roads 50 m, commercial 1000 m, distance to the nearest major roads, major roads 300 m (2.44, 0.0012, 0.049, 0.11) Overall All roads 50 m, commercial 1000 m, distance to the nearest major roads, major roads 300 m (2.73, 0.0013, 0.055, 0.13) Overall All roads 50 m, commercial 1000 m, distance to the nearest major roads, open space 1000 m (3.00, 0.0014, 0.075, 0.00028) Overall
RMSE
Cross-validation LUR þ Kriging
LUR
Adj R2
RMSE
Adj R2
RMSE
Adj R2
RMSE
OK
1st stage model Adj R2
0.72
0.13
0.68
0.14
0.75
0.12
0.43
0.19
e
0.66
0.10
0.39
0.37
0.39
0.37
0.31
0.39
0.28
e 0.75
e 0.14
0.27 0.48
0.34 0.42
0.40 0.49
0.37 0.41
e 0.28
e 0.49
0.30
e 0.80
e 0.15
0.40 0.59
0.33 0.43
0.49 0.59
0.42 0.43
e 0.34
e 0.54
0.37
e 0.79
e 0.17
0.42 0.59
0.39 0.46
0.57 0.59
0.44 0.45
e 0.09
e 0.68
0.33
e 0.74
e 0.24
0.52 0.53
0.44 0.60
0.56 0.53
0.47 0.59
e 0.16
e 0.79
0.26
e
e
0.48
0.57
0.49
0.62
e
e
The numbers in the parenthesis were the coefficients corresponding to the LUR covariables, which were multiplied by 1000 for ease of interpretation.
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
7
Table 3 Subgroup cross-validation of daily PM2.5 concentration classified by AQM station characteristics based on the established LUR models. Group
Background/National Park (n ¼ 2) LUR þ kriging 2
1 2 3 4 5
LUR
General (n ¼ 13) LUR þ kriging
OK 2
2
2
Traffic (n ¼ 3) LUR
LUR þ kriging
OK 2
2
2
LUR
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
0.32 0.36 0.45 0.43 0.20
0.41 0.45 0.46 0.53 0.83
0.32 0.36 0.45 0.43 0.20
0.41 0.45 0.46 0.53 0.83
0.27 0.25 0.32 0.25 0.19
0.43 0.49 0.52 0.61 0.83
0.33 0.40 0.53 0.54 0.44
0.38 0.43 0.44 0.46 0.60
0.34 0.42 0.53 0.54 0.45
0.38 0.42 0.44 0.46 0.59
0.34 0.38 0.47 0.28 0.32
0.38 0.43 0.47 0.57 0.66
0.38 0.41 0.55 0.48 0.57
0.30 0.34 0.33 0.35 0.35
0.38 0.41 0.55 0.49 0.59
overestimate the PM2.5 concentrations. Estimates based on kriging methods or simply average measurements for the rural sites themselves may be sufficient in this case. Tables 2 and 3 show that the adjusted R-squares increased with ^ max ðtÞ. An intuitive explanation of this groups with increasing g outcome is that the LUR models reflected better spatial variation ^ max ðtÞ, and because of decreasing spatial covariance divided by g thus made better predictions in latter subgroups. Note that the mean and standard deviation of PM2.5 concentration decreased with the subgroups (Table 1), which also explained the order of the adjusted R-squares. The meteorological covariates temperature, relative humidity, wind speed, and atmospheric pressure further justified the grouping of the dates to apply separate LUR models accordingly. Based on the results of Tables 1 and 2, extra attention to the PM2.5 concentration prediction may need to be paid under certain meteorological conditions, especially for subgroup 5 with ^ max ðtÞ. the most heterogeneity in g We estimated and removed mean temporal trends because temporal covariates are expected to be spatially homogeneous across the whole Taipei metropolitan area. For a larger study region with temporal covariates that are spatially dependent, a trend estimate and residuals (Eqs. (1) and (2)) may also be modified to be site-specific. The rest of the procedure remains the same. The GIS covariates were coded based on classifications of road type and land-use information obtained from the National Land Survey and Mapping Center of Taiwan. The measures could vary with different coding criteria, and different LUR models may be generated. However, spatial predictions should be approximately the same, with minor changes in adjusted R-squares. Also, the primary purpose of this study is to illustrate temporal-dependent LUR models rather than geographical variable coding. There were some limitations with the study. First, the spatiotemporal model was developed based on the air monitoring data and GIS information of only 18 AQM stations. Therefore, the estimated variogram model may not be stable for spatial extrapolations in kriging. Second, there was a lack of traffic density in our data. We used the lengths of all roads, major roads, highway, and distance to the nearest major roads as surrogate measures. However, even with only 18 sites, we were able to identify subgroups that were approximately temporally invariant with different LUR models classified by the sill of daily empirical variogram. The resultant overall and subgroup LUR models also demonstrated that the road length covariates worked very well as surrogates of traffic-related PM2.5 concentrations. Third, the grouping of the dates was somewhat arbitrary. However, similar results were obtained with other grouping schemes (data not shown). Therefore, it appears that the grouping scheme does not greatly affect the predictions, and the number of subgroups depends on the scope and the length of the study period. Further exploration is needed of how to achieve optimal prediction with approximately temporally invariant subgroups.
OK 2
RMSE
Adj R2
RMSE
0.30 0.34 0.33 0.35 0.35
0.32 0.36 0.51 0.41 0.51
0.31 0.36 0.34 0.37 0.38
6. Conclusion In summary, we have demonstrated that the proposed method successfully divides the dates of the study period into approximately temporal-invariant homogeneous subgroups with separate LUR models for daily PM2.5 concentration predictions, which had in general better performance in comparison with a single LUR model. Depending on mean PM2.5 concentration, traffic density, and the meteorological covariates of the date used for prediction, the crossvalidated R-squared value varied from 0.39 to 0.59 across different subgroups using the LUR model alone. Predictions using UK could outperform the LUR model if the prediction region and the number of AQM stations were larger than those used in the present study. Further investigation of this possibility is required. Acknowledgments The authors thank Ms. You-Sin Liang for calculating the GISbased measures and Mr. Jih-Shin Liu for the statistical computations and programming. This research was supported by projects PH-99-PP-09 of the National Health Research Institutes and NSC 98-EPA-M-002-001 of the National Science Council of Taiwan. References Alexeeff, S.E., Coull, B.A., Gryparis, A., Suh, H., Sparrow, D., Vokonas, P.S., Schwartz, J., 2011. Medium-term exposure to traffic-related air pollution and markers of inflammation and endothelial function. Environmental Health Perspectives 119 (4), 481e486. Banerjee, S., Carlin, B.P., Gelfand, A.E., 2004. Hierarchical Modeling and Analysis for Spatial Data. Chapman and Hall, CRC, Boca Raton. Bogaert, P., Christakos, G., Jerrett, M., Yu, H.-L., 2009. Spatiotemporal modeling of ozone distribution in the state of California. Atmospheric Environment 43, 2471e2480. Chan, C.C., Chuang, K.J., Chen, W.J., Chang, W.T., Lee, C.T., Peng, C.M., 2008. Increasing cardiopulmonary emergency visits by long-range transported Asian dust storms in Taiwan. Environmental Research 106, 393e400. Cheng, T.J., Hwang, J.S., Wang, P.Y., Tsai, C.F., Chen, C.Y., Lin, S.H., et al., 2003. Effects of concentrated ambient particles on heart rate and blood pressure in pulmonary hypertensive rats. Environmental Health Perspectives 111 (2), 147e150. Chuang, K.J., Chan, C.C., Chen, N.T., Su, T.C., Lin, L.Y., 2005. Effects of particle size fractions on reducing heart rate variability in cardiac and hypertensive patients. Environmental Health Perspectives 113, 1693e1697. Clougherty, J.E., Wright, R.J., Baxter, L.K., Levy, J.I., 2008. Land use regression modeling of intra-urban residential variability in multiple traffic-related air pollutants. Environmental Health 7, 17. Dadvand, P., Rushton, S., Diggle, P.J., Goffe, L., Rankin, J., Pless-Mulloli, T., 2011. Using spatio-temporal modeling to predict long-term exposure to black smoke at fine spatial and temporal scale. Atmospheric Environment 45, 659e664. Fanshawe, T.R., Diggle, P.J., Rushton, S., Sanderson, R., Lurz, P.W.W., Glinianaia, S.V., Pearce, S., Parker, L., Charlton, M., Pless-Mulloli, T., 2008. Modeling spatiotemporal variation in exposure to particulate matter: a two-stage approach. Environmetrics 19, 549e566. Gryparis, A., Coull, B.A., Schwartz, J., Suh, H.H., 2007. Semiparametric latent variable regression models for spatiotemporal modeling of mobile source particles in the greater Boston area. Journal of the Royal Statistical Society e Series C: Applied Statistics 56, 183e209. Henderson, S.B., Beckerman, B., Jerrett, M., Brauer, M., 2007. Application of land use regression to estimate long-term concentrations of traffic-related nitrogen oxides and fine particulate matter. Environmental Science and Technology 41, 2422e2428.
8
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
Hoek, G., Beelen, R., de Hoogh, K., Vienneau, D., Gulliver, J., Fischer, P., Briggs, D., 2008. A review of land-use regression models to assess spatial variation of outdoor air pollution. Atmospheric Environment 42, 7561e7578. Jerrett, M., Arain, A., Kanaroglou, P., Beckerman, B., Potoglou, D., Sahsuvaroglu, T., Morrison, J., Giovis, C., 2005. A review and evaluation of intraurban air pollution exposure models. Journal of Exposure Analysis and Environmental Epidemiology 15, 185e204. Kammann, E.E., Wand, M.P., 2003. Geoadditive models. Applied Statistics 52, 1e18. Maynard, D., Coull, B.A., Gryparis, A., Schwartz, J., 2007. Mortality risk associated with short-term exposure to traffic particles and sulfates. Environmental Health Perspectives 115, 751e755. Ross, Z., Jerrett, M., Ito, K., Tempalski, B., Thurston, G.D., 2007. A land use regression for predicting fine particulate matter concentrations in the New York City region. Atmospheric Environment 41, 2255e2269. Ryan, P.H., LeMasters, G.K., 2007. A review of land-use regression models for characterizing intraurban air pollution exposure. Inhalation Toxicology 19 (Suppl. 1), 127e133. Schwartz, J., 2004. The effects of particulate air pollution on daily deaths: a multicity case crossover analysis. Occupational Environmental Medicine 61, 956e961.
Szpiro, A.A., Sampson, P.D., Sheppard, L., Lumle, T., Adar, S.D., Kaufman, J.D., 2010. Predicting intra-urban variation in air pollution concentrations with complex spatio-temporal dependencies. Environmetrics 21, 606e631. Vienneau, D., de Hoogh, K., Beelen, R., Fischer, P., Hoek, G., Briggs, D., 2010. Comparison of land-use regression models between Great Britain and the Netherlands. Atmospheric Environment 44, 688e696. Wilson, J.G., Kingham, S., Pearce, J., Andrew, P., Sturmana, B., 2005. A review of intraurban variation in particulate air pollution: Implication for epidemiological research. Atmospheric Environment 39, 6444e6462. Yanosky, J.D., Paciorek, C.J., Suh, H.H., 2009. Predicting chronic fine and coarse particulate exposures using spatiotemporal models for the northeastern and midwestern United States. Environmental Health Perspectives 117 (4), 522e529. Yu, H.-L., Wang, C.-H., 2010. Retrospective prediction of intraurban spatiotemporal distribution of PM10 in Taipei. Atmospheric Environment 44, 3053e3065. Yu, H.-L., Wang, C.-H., Liu, M.-C., Kuo, Y.-M., 2011. Estimation of fine particulate matter in Taipei using land use regression and bayesian maximum entropy methods. International Journal of Environmental Research and Public Health 8, 2153e2169. Zanobetti, A., Schwartz, J., 2009. The effect of fine and coarse particulate air pollution on mortality: a national analysis. Environmental Health Perspectives 117 (6), 898e903.
Contents lists available at SciVerse ScienceDirect
Atmospheric Environment journal homepage: www.elsevier.com/locate/atmosenv
Spatiotemporal modeling with temporal-invariant variogram subgroups to estimate fine particulate matter PM2.5 concentrations Chu-Chih Chen a, *, Chang-Fu Wu b,1, Hwa-Lung Yu c, Chang-Chuan Chan b, Tsun-Jen Cheng b a
Division of Biostatistics and Bioinformatics, Institute of Population Health Sciences, National Health Research Institutes, Taiwan Institute of Occupational Medicine and Industrial Hygiene, College of Public Health, National Taiwan University, Taiwan c Department of Bioenvironmental Systems Engineering, National Taiwan University, Taiwan b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 7 July 2011 Received in revised form 12 December 2011 Accepted 6 February 2012
Short-term exposure estimation of daily air pollution levels incorporating geographic information system (GIS) into spatiotemporal modeling remains a great challenge for assessing corresponding acute adverse health effects. Due to daily meteorological effects on the dispersion of pollutants, explanatory spatial covariables and their coefficients may not be the same as in classical land-use regression (LUR) modeling for long-term exposure. In this paper, we propose a two-stage spatiotemporal model for daily fine particulate matter (PM2.5) concentration prediction: first, daily nonlinear temporal trends are estimated through a generalized additive model, and second, GIS covariates are used to predict spatial variation in the temporal trend-removed residuals. To account for spatial dependence on meteorological conditions, the dates of the study period are divided by the sill of the daily empirical variogram into approximately temporal-invariant subgroups. Within each subgroup, daily PM2.5 estimations are obtained by combining the temporal and spatial parts of the estimations from the two stages. The proposed method is applied to the modeling of spatiotemporal PM2.5 concentrations observed at 18 ambient air monitoring stations in Taipei metropolitan area during 2006e2008. The results showed that the PM2.5 concentrations decreased whereas the relative humidity and wind speed increased with the sill subgroups, which may be due to the effects of daily meteorological conditions on the dispersions of the particles. Also, the covariates and their coefficients of the LUR models varied with subgroups and had in general higher adjusted R-squares and smaller root mean square errors in prediction than those of a single overall LUR model. 2012 Elsevier Ltd. All rights reserved.
Keywords: Empirical variogram Geographic information system Kriging Land-use regression Temporal trend
1. Introduction Numerous epidemiological and toxicological studies have shown that increased particulate matter (PM) is associated with increased cardiovascular mortality and morbidity and that fine PM particles with aerodynamic diameter 2.5 mm (PM2.5) are more hazardous than coarse particles, such as PM10 (Cheng et al., 2003; Chuang et al., 2005; Schwartz, 2004; Zanobetti and Schwartz, 2009). Many of the studies relate mortality or hospital emergency visit data to mean air pollution level measured at a central site or from several monitoring stations. However, PM2.5 concentrations are subject to within-city or intra-urban variations (Jerrett et al., 2005; Wilson et al., 2005), and individuals may be exposed to Abbreviations: GIS, geographic information system; LUR, land-use regression; OK, ordinary kriging; PM, particulate matter; RMSE, root mean square error; UK, universal kriging. * Corresponding author. E-mail addresses: [email protected] (C.-C. Chen), [email protected] (C.-F. Wu). 1 Equal contribution with the correspondence. 1352-2310/$ e see front matter 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.atmosenv.2012.02.015
different levels depending on their residential area. To estimate small-area air pollution concentrations, in addition to using geostatistical interpolation based on sampled site measurements, land-use regression (LUR) models using data on nearby traffic and land-use covariates have been successfully applied to explain spatial variation (Clougherty et al., 2008; Henderson et al., 2007; Hoek et al., 2008; Ross et al., 2007; Ryan and LeMasters, 2007). Typically such models are built by regressing annual averages of pollutant concentrations over geographical information system (GIS) variables that are invariant to temporal variations. As a result, standard LUR models are incapable of studying acute health effects due to short-term exposures. Furthermore, temporal variations in PM concentrations due to meteorological conditions on dispersion effects or weekdays also need to be considered to relate to individual acute adverse effects. Therefore, a spatiotemporal model to predict PM concentrations at unsampled locations on a specific date is required for individual exposure estimation. Using a semiparametric latent variable regression model, Gryparis et al. (2007) extend geoadditive models (Kammann and
2
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
Wand, 2003) by incorporating temporal covariables for spatiotemporal modeling of mobile particle concentrations. Maynard et al. (2007) use corresponding short-term exposure modeling to study associated mortality risks considering GIS-based measures of cumulative traffic density. Similarly, Alexeeff et al. (2011) use the same modeling to estimate subjects’ daily black carbon exposure to relate to their repeated inflammation markers. However, the approach specifies a lack of spaceetime interaction over time across locations, which may not be valid due to the fact that spatial dispersion of air pollutants is temporal-dependent (Bogaert et al., 2009; Yu and Wang, 2010; Yu et al., 2011). To separate temporal and spatial variations, modeling based on two-stage approach by separating temporal trend from withinregion spatial variation in the mean structure have been proposed (Dadvand et al., 2011; Fanshawe et al., 2008; Yanosky et al., 2009). However, additional information about spatiotemporal variability and correlation structures in the data need to be considered (Banerjee et al., 2004; Bogaert et al., 2009; Yu and Wang, 2010). To deal with the complex spatial and temporal covariance dependence, Bogaert et al. (2009) divide and estimate the spatial variogram of ozone distribution by month. Alternatively, Szpiro et al. (2010) decompose a hierarchical model incorporating dependence on GIS covariates with temporal trends and a residual with spatial correlation structure. However, it may be computationally difficult to extend both of these spatiotemporal models to predict daily air pollution. In this study, we adopt a two-stage approach: (1) in the first stage, the daily mean temporal trend of the (log-transformed) PM2.5 concentrations of the study area is estimated by using a generalized additive model (GAM); (2) in the second stage, LUR covariables are incorporated for spatial prediction of the (temporal trend removed) residuals from the first stage. Because temporal dependence may still be present, in stage 2, the dates of the study period were classified into subgroups based on the sill of the daily empirical semivariogram of the residuals from stage 1. Dates falling within the same group were then aggregated together to obtain the time-averaged residuals and a LUR model was applied separately, instead of an overall LUR model. Kriging method was further applied to the residuals of the LUR models in the presence of spatial autocorrelations. We then predicted the final concentration level at a specific date by restoring the temporal trend to the spatial prediction. Daily averaged PM2.5 concentrations from 18 ambient air monitoring stations in the Taipei metropolitan area between 2006 and 2008 were analyzed as an example.
2. Materials 2.1. PM2.5 and meteorological measurements Air pollution levels were obtained from the Air Quality Monitoring (AQM) Network of the Environmental Protection Administration (EPA) of Taiwan, from the 18 monitoring stations that are located within the Taipei metropolitan area, including Taipei City and New Taipei City within the Taipei Basin (Fig. 1). Daily averages of PM2.5, wind speed, temperature, and relative humidity were calculated based on the available monitoring data from January 1, 2006, to December 31, 2008. 2.2. GIS data GIS-based data for the 18 AQM stations fall into three categories: (1) land use (residential, commercial, urban green, forest, industry, and open space) within buffers of several radii (100, 300, 500, and 1000 m); (2) types and lengths of roads (all roads, highways, and major roads) within buffer zone of 25, 50, 100, 300, 500, and 1000 m; (3) population size within radii of 100, 300, 500, and 1000 m. Also, distance to the nearest major roads is included as a covariate. The land-use database was obtained from the National Land Survey and Mapping Center of Taiwan. Road network information was acquired from the Institute of Transportation, Ministry of Transportation and Communications, and census population data were accessed from the Civil Affairs Bureau of Taipei City and New Taipei City, Taiwan. ArcGIS 9.3 was employed to calculate the lengths of roads and population size within each buffer zone of the 18 AQM stations. 3. Statistical method 3.1. Two-stage spatiotemporal modeling While PM2.5 concentration levels may change nonlinearly with temperature and relative humidity, the meteorological effect on the dispersion of particles is expected to be uniform across a mediumsized study region such as the Taipei metropolitan area (approximately 2325 km2). This implicit assumption suggests that for assessing pollutant levels of a metropolitan area, spatial and temporal variations in concentrations may be decomposed separately. Suppose the log-transformed measurement of monitoring
Fig. 1. Map of the Taipei metropolitan area and the locations of the 18 AQM stations.
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
station i at location si and time t is y(si,t), i ¼ 1, ., K; t ¼ 1, ., T, a two-stage spatiotemporal modeling method is proposed as follows. 3.1.1. Stage 1. Modeling for daily mean temporal trend Let 0
YðtÞ ¼ x ðtÞb þ
H X
3.2. Grouping of dates based on the daily empirical semivariogram For each day, an empirical semivariogram was constructed based on the residuals of the log-transformed daily average PM2.5 concentrations of the 18 monitoring stations (Eq. (2)) as follows:
g^ k ðtÞ ¼ fh ðqh ðtÞÞ þ ZðtÞ
(1)
h¼1
be the mean temporal GAM model on day t, where YðtÞ is the average of the K monitoring station measurements Y(i,t) on day t, x(t) is a vector of the linear or categorical temporal predictors such as weekdays and months, b is the corresponding vector of coefficients, and fh(qh(t)), h ¼ 1, ., H is a nonlinear function of the predictors q1, ., qH on day t. The correlations of model residuals Z(t) may be assumed to follow an autoregressive process AR(p) of order p. Let the trend-removed residuals at station si be
^ dðsi ; tÞ ¼ Yðsi ; tÞ YðtÞ;
(2)
^ where YðtÞ is the mean estimate on day t. 3.1.2. Stage 2. Land-use regression with spatial interpolation Given a fixed date t, a spatial model is applied to the residuals d(si,t) obtained from stage 1 with GIS information as explanatory covariates. Because the spatial covariance of the residuals d(si,t) may still change with time, spatial interpolation using the kriging method may vary with different dates. Also, while GIS-based measures are invariant with time, the significant LUR model covariates and the corresponding coefficients may change with daily dispersion conditions. It is thus reasonable to aggregate the residuals across different days with similar spatial correlations as a group. The procedure is accomplished by dividing the sills of the daily empirical variograms into approximately homogeneous groups as described in the next subsection. Within each group, LUR covariables are then incorporated (with kriging in the presence of spatial autocorrelation) for spatial prediction of the subgroupaveraged residuals. The LUR covariables are selected following the procedure given in Vienneau et al. (2010): 1) the increase in adjusted R2 of the predictor variable greater than 1% was selected; 2) the coefficient conformed to the pre-specified direction of effect; 3) the direction of effect for predictors already included in the model did not change; and 4) the p-value for all variables already in the model remained below 0.1. When the spatial correlations among the residuals are negligible, the LUR model may be written as
dl ðsi Þ ¼ w0l ðsi Þ4l þ eðsi Þ;
3
1 2NBk
X
dðsi ; tÞ d sj ; t
2
;
(5)
fði;jÞ:ksi sj k˛Bk g
where NBk is the number of pairs of sites with distance fall within the bin Bk. Ten bins were divided with 15 pairs of sites each for each day. Because the daily meteorological conditions may vary substantially from day to day, the effects on dispersions of the PM2.5 particles should also be different. The dispersions may further be affected by the differential PM2.5 concentrations across the monitoring sites (e.g., traffic sites versus rural sites). For example, spatial correlation of PM2.5 concentrations is expected to be smaller and the corresponding variogram is thus larger in rainy days because of cleaner air. Similarly, higher wind speed could keep the atmosphere mixed better. As a result, spatial correlation also tends to be smaller. ^ k ðtÞ may Therefore, the daily empirical semivariogram values g change accordingly, which may be assessed by the sill estimated by the daily maximum empirical semivariogram, i.e.,
g^ max ðtÞ ¼ max g^ k ðtÞ: k
(6)
The corresponding dates were then divided into subgroups by ^ max ðtÞ from the smallest to the largest. For the dates sorting g belonging to the same group, a LUR model was fitted separately. As a comparison, UK with a variogram model fitting within each of the subgroups was also applied for prediction. 4. Results The averaged daily PM2.5 concentrations across the 18 AQM stations of the Taipei metropolitan area had a mean 28.35 mg m3, ranged from 0.25 to 128.20 mg m3 with an inter-quartile range of 20.04 mg m3 during the study period 2006e2008. Fig. 2 shows the monthly PM2.5 concentrations of the AQM stations at the general sites (n ¼ 13), traffic sites (n ¼ 3) and the background sites (n ¼ 2), respectively. The results showed that the monthly mean PM2.5 concentrations were in general higher during the cold months (JanuaryeApril) compared to those of the warm months (JuneeOctober). Site effects of the AQM stations were also apparent, with the mean concentrations of the traffic sites the highest, followed by those of the general sites. Although the background stations had the lowest mean PM2.5 concentrations, the values
(3)
where dl ðsi Þ is the averaged residual of d(si,t) with t ˛ Tl, Tl is the l-th set of dates with similar variograms, l ¼ 1, ., L, wl(si) is a vector of the LUR covariables at site si, 4l is the corresponding vector of coefficients, and e(si) are assumed to be independent, zero-mean, and a common variance of s2e . In the presence of spatial correlation, the resultant model for universal kriging (UK, LUR þ kriging) becomes
dl ðsi Þ ¼ w0l ðsi Þ4l þ Sðsi Þ þ eðsi Þ;
(4)
where {S(s):s ˛ R2} is a stationary zero-mean stochastic process with spatial covariance matrix C, and e(si) is the nugget effect with spatially-independent variance s2e . Finally, given a day t ˛ Tl, a prediction of the log-PM2.5 concentration at location s0 can be obtained from the results of Eqs. (1)e(4).
Fig. 2. Monthly mean PM2.5 concentrations of the general, traffic, and background stations of the Taipei metropolitan area from January 2006 to December 2008.
4
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
could also reach 25 mg m3 for certain months, which could be due to the long-range Asian sand storms transported from outside of Taiwan (Chan et al. 2008). 4.1. Site-specific annual-average PM2.5 predictions using simple LUR and kriging method By aggregating the daily average PM2.5 concentrations of the 18 AQM stations from 2006 to 2008, a simple annual LUR model was obtained using the GIS-based measures. Of the 48 explanatory variables, only the three covariables: all roads within 50 m, commercial area within 1000 m, and distance to the nearest major roads remained in the final model. The resultant LUR model and ordinary kriging (OK) were applied separately to predict the mean PM2.5 concentrations of the 18 AQM stations for cross-validation. Because the number of monitoring stations was barely enough to yield stable variogram model parameter estimates, the empirical variogram was generated by overlaying daily estimates from the 1096 days together, assuming no spatialetemporal interaction. The cross-validation adjusted R-square of OK was substantially lower than that of LUR (0.43 versus 0.75). Also, the adjusted R-square 0.68 of UK was lower than that of LUR alone, which could be due to additional variation from the variogram model fitting (Eq. (4)) with non-significant residual spatial autocorrelation (Moran’s I statistic p-value 0.28). The results are summarized together with those of the two-stage models in Table 2 below. 4.2. Temporal trend estimation and grouping of the dates A GAM model was fitted to the daily average log-PM2.5 concentrations across the 18 AQM stations using the SAS procedure PROC GAM. Daily mean temperature, relative humidity, and wind speed all had significant nonlinear effects. Month also had a significant effect on the PM concentrations. An autocorrelation model of lag 2 (AR(2)) was fitted to the residuals based on the results of a DubineWatson test. Fig. 3 displays the model residuals of the average log-PM2.5 concentrations along the dates of the study period. To illustrate the differences across various sites, mean residuals of the general, traffic, and background monitoring stations were plotted separately for every 7 days. As shown in the figure, the mean residuals of the general stations clustered around the line 0, whereas the mean residuals of the traffic sites and the background sites tended to be positive (underestimated) and negative (overestimated), respectively. The pattern is expected because of traffic-related PM2.5 concentration variations across different sites. However, the mean residuals of different sites appeared to be temporally invariant across the dates. Therefore, the assumption that the temporal and spatial variations may be decomposed separately at two stages held approximately legitimate. After the residuals were obtained following Eq. (2), daily ^ max ðtÞ for the sill estimate maximal empirical semivariograms g were calculated. A separate GAM model was applied to the daily g^ max ðtÞ to examine the temporal effects on the spatial variances. The result showed that the daily mean temperature, relative humidity, and wind speed all had significant effects on the daily g^ max ðtÞ (Fig. 4). Because of the dependence of the residuals on meteorological ^ max ðtÞ were divided into 5 subgroups so that conditions, the sorted g the variograms within each subgroup were approximately homo^ max ðtÞ was much geneous. Due to the fact that the distribution of g skewed to the left, we chose to classify the subgroups with the cutoffs of the percentiles at 60%, 70%, 80%, and 90%, respectively. Corresponding dates within the same subgroup were aggregated for separate LUR model fitting at the second stage. Fig. 5 displays
Fig. 3. Mean temporal trend-removed residuals of the log-PM2.5 concentrations across the dates of the study period. The mean residuals of the general, traffic, and background stations are plotted for every 7 days and are denoted by “”, “D”, and “”, respectively.
^ max ðtÞ. The medians of the nonthe box plots of the 5 subgroups of g ^ max ðtÞ indicated inter-group heterogeneity overlapping subgroup g of spatial variances. Notice that there were some estimates of g^ max ðtÞ in the fifth subgroup approximately 10 times greater than those of the first subgroup, which indicated that the daily empirical variograms were quite heterogeneous. 4.3. Subgroup LUR model fitting and cross-validations of different prediction methods Table 1 summarizes the means and standard deviations of PM2.5 concentration, atmospheric pressure, temperature, relative humidity, and wind speed of the dates corresponding to the 5 subgroups. Group differences of these variables showed that, except for the atmospheric pressure, all the other variables were highly significant with a p-value < 0.0001 using an analysis of variance (ANOVA) test. The PM2.5 concentration and temperature decreased with the subgroups, whereas the relative humidity, wind speed increased. Tests for trend of these variables were also all highly significant (p-value < 0.0001). Therefore, the grouping of the ^ max ðtÞ effectively classified the interrelationship between dates by g the daily mean PM2.5 concentrations and the mean temperature, relative humidity, and wind speed. A LUR model was fitted to each of the sill subgroups separately. The resultant models are summarized in Table 2. Except for group 2, all the other groups had covariates included different from those of the overall LUR model. Similar to the overall model, the covariables of the length of all roads within 50 m and distance to the nearest major roads were included in all subgroup LUR models, which indicated that traffic activity was the major predictor. Commercial area within 1000 m was included in all except for subgroup 1. Length of major roads within 300 m was included in subgroups 2 and 3, and open space within 1000 m was included in subgroup 5 only. All of the other covariates, such as the portion of green space and population size at various buffers, were non-significant. Adjacent subgroups tended to have similar covariables because of similarities in the empirical sill. Notice that for the same covariable included in the subgroup LUR models, the coefficients tended to increase with the groups (for example, the length of all types of
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
5
Fig. 4. Semiparametric spline curves of the sill of the daily empirical variograms of the temporal trend-removed residuals with the predictors (a) temperature, (b) relative humidity, and (c) wind speed.
roads within 50 m, commercial area within 1000 m, and distance to the nearest major roads), which was consistent with the scales of g^ max ðtÞ of the subgroups. In addition to the LUR model fittings for each of the subgroups, Table 2 lists the adjusted R-squares and RMSEs of the LUR model building and cross-validations. Because of much higher variation in estimating daily concentrations than the average 3-year level, the subgroup adjusted R-squares were lower
Fig. 5. Box plots of the sills of the daily empirical semivariograms of the temporal trend-removed residuals for each of the 5 subgroups.
than those of the overall model for annual predictions. In comparison with the results using the same overall LUR model throughout the 5 subgroups, the subgroup LUR models in general had higher adjusted R-squares and lower RMSEs. Together with the outcomes shown in Table 1, these results support the hypothesis that meteorological conditions had an effect on daily spatial associations of PM2.5 concentrations as stated in Section 3.2. Also, the adopted procedure by fitting a subgroup LUR model separately effectively improved over a single overall LUR model, especially for the dates with larger scale in the sill of the daily empirical variogram. Similar as before, UK, LUR, and OK were compared for their performance in daily PM2.5 concentration predictions. The adjusted R-squares ranged from 0.39 to 0.59 for UK, 0.39 to 0.59 for LUR only, and 0.09 to 0.34 for OK, respectively. Notice that the adjusted Rsquares increased from group 1 to group 5 with increasing RMSEs for all three methods because the subgroups were divided by increasing spatial variance, resulting in increasing biases. Similar to the overall model, the OK method had the lowest R-squares and highest RMSE for all 5 subgroups. Possibly due to extra variation in variogram model fitting with non-significant spatial associations (Moran’s I statistic p-values were 0.82, 0.82, 0.81, 0.69, and 0.34, respectively), the results of UK were essentially the same as those of LUR within each subgroup. However, all the three methods at stage 2 had higher R-squares compared with those methods using estimates based only on mean temporal trend (the last column of Table 2), which indicated the existence of spatial variations in the temporal trend-removed residuals. To assess the prediction performance at different sites, Table 3 further compared the AQM stations at the general, traffic, and background sites. The results show that the proposed procedure had the best performance at the traffic sites, followed by the
6
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
Table 1 Characteristics of the 5 subgroups: number of days, PM2.5, atmospheric pressure, temperature, relative humidity, and wind speed. Group
N (days)
PM2.5 (mg m3) Mean
std
Mean
std
Mean
std
Mean
std
Mean
std
1 2 3 4 5
656 110 109 109 110
32.9 25.2 23.0 19.9 18.0
14.2 13.1 13.7 12.1 11.3
1012.5 1011.2 1012.0 1013.1 1013.6
6.9 7.3 6.3 6.7 6.5
23.3 22.1 22.3 21.8 22.1
4.9 5.0 5.2 5.0 5.3
71.8 77.3 79.7 80.7 80.6
6.9 6.6 6.4 6.7 8.2
2.24 2.25 2.20 2.45 2.49
0.59 0.71 0.76 0.81 1.06
Temperature ( C)
Atmospheric pressure (hPa)
general sites. Sites with background concentrations, or national parks in general, had relatively lower adjusted R-squares than those of the other sites. 5. Discussions We have created a two-stage spatiotemporal model to predict daily PM2.5 concentrations by first removing the mean temporal trend and then accounting for spatial variations in the residuals. Because of the dependence of pollutant dispersion on daily meteorological conditions, the dates of the study period were divided by ^ max ðtÞ of daily empirical semivariogram into the sill estimate g subgroups that were approximately temporally invariant. This approach allows us to incorporate GIS covariables for spatial associations of the temporal trend-removed residuals at the second stage without an additional assumption of temporal independence. In contrast with a constant LUR model for long-term exposure estimation, subgroup LUR models with different covariates and coefficient estimates were obtained as a result, which had in general better performance in adjusted R-squares and RMSE for prediction compared to those of a single overall model. ^ max ðtÞ of the temporal trend-removed Our results showed that g residuals depended on mean temperature, wind speed, and relative humidity (Fig. 5). The results given in Table 1 also demonstrate that ^ max ðtÞ effectively classified the interrelationship the grouping of g between the PM2.5 concentrations and the daily meteorological conditions. Thus, the assumption of lacking a spaceetime
Relative humidity (%)
Wind speed (m s1)
interaction with log-scale daily air pollution concentrations (Gryparis et al., 2007; Maynard et al., 2007) does not hold. Although similar two-stage modeling has been proposed previously (Dadvand et al., 2011; Fanshawe et al., 2008; Yanosky et al., 2009), the ambient pollutant concentration predictions of these works are mainly based on weekly or monthly time scales. Therefore, there is essentially no need to adjust for significant terms and corresponding coefficients of GIS-based covariates on daily meteorological effects. The model proposed by Szpiro et al. (2010), though appealing with temporal basis functions approach, lacks an intuitive interpretation, and the computational burden becomes cumbersome when the time scale is for daily mean concentrations. We have shown in Table 2 that the overall and subgroup crossvalidation R-squares using LUR only were slightly higher than those of UK and substantially higher than those of OK. This indicates that the spatial associations among the temporal trend-removed residuals at different sites could be essentially explained by the GIS covariates without further spatial interpolations, as was supported by the non-significant Moran’s I test statistics for each of the subgroups. The results given in Table 3 also indicate that the PM2.5 concentration predictions may perform best at sites near heavy traffic. However, for rural or background sites with little to no traffic, the predictions may not be accurate. Possible explanations of this phenomenon are that the covariates of the subgroup LUR models were all essentially traffic-based and that there were only two background sites in our study. Therefore, estimates that were largely based on the general or the traffic sites tended to
Table 2 LUR models for the full study period and each of the 5 subgroups, performance statistics for model building, and cross-validation with different methods. The rows below each of the subgroup LUR models are the corresponding results by fitting the overall LUR model. Group
LUR covariables
Model building 2
Adj R
Overall
1
2
3
4
5
a
Annual All roads 50 m, commercial 1000 m, distance to the nearest major roads (1.69, 0.00083, 0.048)a Empirical variogram subgroups All road 50 m, distance to the nearest major roads (1.16, 0.037) Overall All roads 50 m, commercial 1000 m, distance to the nearest major roads (1.85, 0.00098, 0.054) Overall All roads 50 m, commercial 1000 m, distance to the nearest major roads, major roads 300 m (2.44, 0.0012, 0.049, 0.11) Overall All roads 50 m, commercial 1000 m, distance to the nearest major roads, major roads 300 m (2.73, 0.0013, 0.055, 0.13) Overall All roads 50 m, commercial 1000 m, distance to the nearest major roads, open space 1000 m (3.00, 0.0014, 0.075, 0.00028) Overall
RMSE
Cross-validation LUR þ Kriging
LUR
Adj R2
RMSE
Adj R2
RMSE
Adj R2
RMSE
OK
1st stage model Adj R2
0.72
0.13
0.68
0.14
0.75
0.12
0.43
0.19
e
0.66
0.10
0.39
0.37
0.39
0.37
0.31
0.39
0.28
e 0.75
e 0.14
0.27 0.48
0.34 0.42
0.40 0.49
0.37 0.41
e 0.28
e 0.49
0.30
e 0.80
e 0.15
0.40 0.59
0.33 0.43
0.49 0.59
0.42 0.43
e 0.34
e 0.54
0.37
e 0.79
e 0.17
0.42 0.59
0.39 0.46
0.57 0.59
0.44 0.45
e 0.09
e 0.68
0.33
e 0.74
e 0.24
0.52 0.53
0.44 0.60
0.56 0.53
0.47 0.59
e 0.16
e 0.79
0.26
e
e
0.48
0.57
0.49
0.62
e
e
The numbers in the parenthesis were the coefficients corresponding to the LUR covariables, which were multiplied by 1000 for ease of interpretation.
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
7
Table 3 Subgroup cross-validation of daily PM2.5 concentration classified by AQM station characteristics based on the established LUR models. Group
Background/National Park (n ¼ 2) LUR þ kriging 2
1 2 3 4 5
LUR
General (n ¼ 13) LUR þ kriging
OK 2
2
2
Traffic (n ¼ 3) LUR
LUR þ kriging
OK 2
2
2
LUR
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
RMSE
Adj R
0.32 0.36 0.45 0.43 0.20
0.41 0.45 0.46 0.53 0.83
0.32 0.36 0.45 0.43 0.20
0.41 0.45 0.46 0.53 0.83
0.27 0.25 0.32 0.25 0.19
0.43 0.49 0.52 0.61 0.83
0.33 0.40 0.53 0.54 0.44
0.38 0.43 0.44 0.46 0.60
0.34 0.42 0.53 0.54 0.45
0.38 0.42 0.44 0.46 0.59
0.34 0.38 0.47 0.28 0.32
0.38 0.43 0.47 0.57 0.66
0.38 0.41 0.55 0.48 0.57
0.30 0.34 0.33 0.35 0.35
0.38 0.41 0.55 0.49 0.59
overestimate the PM2.5 concentrations. Estimates based on kriging methods or simply average measurements for the rural sites themselves may be sufficient in this case. Tables 2 and 3 show that the adjusted R-squares increased with ^ max ðtÞ. An intuitive explanation of this groups with increasing g outcome is that the LUR models reflected better spatial variation ^ max ðtÞ, and because of decreasing spatial covariance divided by g thus made better predictions in latter subgroups. Note that the mean and standard deviation of PM2.5 concentration decreased with the subgroups (Table 1), which also explained the order of the adjusted R-squares. The meteorological covariates temperature, relative humidity, wind speed, and atmospheric pressure further justified the grouping of the dates to apply separate LUR models accordingly. Based on the results of Tables 1 and 2, extra attention to the PM2.5 concentration prediction may need to be paid under certain meteorological conditions, especially for subgroup 5 with ^ max ðtÞ. the most heterogeneity in g We estimated and removed mean temporal trends because temporal covariates are expected to be spatially homogeneous across the whole Taipei metropolitan area. For a larger study region with temporal covariates that are spatially dependent, a trend estimate and residuals (Eqs. (1) and (2)) may also be modified to be site-specific. The rest of the procedure remains the same. The GIS covariates were coded based on classifications of road type and land-use information obtained from the National Land Survey and Mapping Center of Taiwan. The measures could vary with different coding criteria, and different LUR models may be generated. However, spatial predictions should be approximately the same, with minor changes in adjusted R-squares. Also, the primary purpose of this study is to illustrate temporal-dependent LUR models rather than geographical variable coding. There were some limitations with the study. First, the spatiotemporal model was developed based on the air monitoring data and GIS information of only 18 AQM stations. Therefore, the estimated variogram model may not be stable for spatial extrapolations in kriging. Second, there was a lack of traffic density in our data. We used the lengths of all roads, major roads, highway, and distance to the nearest major roads as surrogate measures. However, even with only 18 sites, we were able to identify subgroups that were approximately temporally invariant with different LUR models classified by the sill of daily empirical variogram. The resultant overall and subgroup LUR models also demonstrated that the road length covariates worked very well as surrogates of traffic-related PM2.5 concentrations. Third, the grouping of the dates was somewhat arbitrary. However, similar results were obtained with other grouping schemes (data not shown). Therefore, it appears that the grouping scheme does not greatly affect the predictions, and the number of subgroups depends on the scope and the length of the study period. Further exploration is needed of how to achieve optimal prediction with approximately temporally invariant subgroups.
OK 2
RMSE
Adj R2
RMSE
0.30 0.34 0.33 0.35 0.35
0.32 0.36 0.51 0.41 0.51
0.31 0.36 0.34 0.37 0.38
6. Conclusion In summary, we have demonstrated that the proposed method successfully divides the dates of the study period into approximately temporal-invariant homogeneous subgroups with separate LUR models for daily PM2.5 concentration predictions, which had in general better performance in comparison with a single LUR model. Depending on mean PM2.5 concentration, traffic density, and the meteorological covariates of the date used for prediction, the crossvalidated R-squared value varied from 0.39 to 0.59 across different subgroups using the LUR model alone. Predictions using UK could outperform the LUR model if the prediction region and the number of AQM stations were larger than those used in the present study. Further investigation of this possibility is required. Acknowledgments The authors thank Ms. You-Sin Liang for calculating the GISbased measures and Mr. Jih-Shin Liu for the statistical computations and programming. This research was supported by projects PH-99-PP-09 of the National Health Research Institutes and NSC 98-EPA-M-002-001 of the National Science Council of Taiwan. References Alexeeff, S.E., Coull, B.A., Gryparis, A., Suh, H., Sparrow, D., Vokonas, P.S., Schwartz, J., 2011. Medium-term exposure to traffic-related air pollution and markers of inflammation and endothelial function. Environmental Health Perspectives 119 (4), 481e486. Banerjee, S., Carlin, B.P., Gelfand, A.E., 2004. Hierarchical Modeling and Analysis for Spatial Data. Chapman and Hall, CRC, Boca Raton. Bogaert, P., Christakos, G., Jerrett, M., Yu, H.-L., 2009. Spatiotemporal modeling of ozone distribution in the state of California. Atmospheric Environment 43, 2471e2480. Chan, C.C., Chuang, K.J., Chen, W.J., Chang, W.T., Lee, C.T., Peng, C.M., 2008. Increasing cardiopulmonary emergency visits by long-range transported Asian dust storms in Taiwan. Environmental Research 106, 393e400. Cheng, T.J., Hwang, J.S., Wang, P.Y., Tsai, C.F., Chen, C.Y., Lin, S.H., et al., 2003. Effects of concentrated ambient particles on heart rate and blood pressure in pulmonary hypertensive rats. Environmental Health Perspectives 111 (2), 147e150. Chuang, K.J., Chan, C.C., Chen, N.T., Su, T.C., Lin, L.Y., 2005. Effects of particle size fractions on reducing heart rate variability in cardiac and hypertensive patients. Environmental Health Perspectives 113, 1693e1697. Clougherty, J.E., Wright, R.J., Baxter, L.K., Levy, J.I., 2008. Land use regression modeling of intra-urban residential variability in multiple traffic-related air pollutants. Environmental Health 7, 17. Dadvand, P., Rushton, S., Diggle, P.J., Goffe, L., Rankin, J., Pless-Mulloli, T., 2011. Using spatio-temporal modeling to predict long-term exposure to black smoke at fine spatial and temporal scale. Atmospheric Environment 45, 659e664. Fanshawe, T.R., Diggle, P.J., Rushton, S., Sanderson, R., Lurz, P.W.W., Glinianaia, S.V., Pearce, S., Parker, L., Charlton, M., Pless-Mulloli, T., 2008. Modeling spatiotemporal variation in exposure to particulate matter: a two-stage approach. Environmetrics 19, 549e566. Gryparis, A., Coull, B.A., Schwartz, J., Suh, H.H., 2007. Semiparametric latent variable regression models for spatiotemporal modeling of mobile source particles in the greater Boston area. Journal of the Royal Statistical Society e Series C: Applied Statistics 56, 183e209. Henderson, S.B., Beckerman, B., Jerrett, M., Brauer, M., 2007. Application of land use regression to estimate long-term concentrations of traffic-related nitrogen oxides and fine particulate matter. Environmental Science and Technology 41, 2422e2428.
8
C.-C. Chen et al. / Atmospheric Environment 54 (2012) 1e8
Hoek, G., Beelen, R., de Hoogh, K., Vienneau, D., Gulliver, J., Fischer, P., Briggs, D., 2008. A review of land-use regression models to assess spatial variation of outdoor air pollution. Atmospheric Environment 42, 7561e7578. Jerrett, M., Arain, A., Kanaroglou, P., Beckerman, B., Potoglou, D., Sahsuvaroglu, T., Morrison, J., Giovis, C., 2005. A review and evaluation of intraurban air pollution exposure models. Journal of Exposure Analysis and Environmental Epidemiology 15, 185e204. Kammann, E.E., Wand, M.P., 2003. Geoadditive models. Applied Statistics 52, 1e18. Maynard, D., Coull, B.A., Gryparis, A., Schwartz, J., 2007. Mortality risk associated with short-term exposure to traffic particles and sulfates. Environmental Health Perspectives 115, 751e755. Ross, Z., Jerrett, M., Ito, K., Tempalski, B., Thurston, G.D., 2007. A land use regression for predicting fine particulate matter concentrations in the New York City region. Atmospheric Environment 41, 2255e2269. Ryan, P.H., LeMasters, G.K., 2007. A review of land-use regression models for characterizing intraurban air pollution exposure. Inhalation Toxicology 19 (Suppl. 1), 127e133. Schwartz, J., 2004. The effects of particulate air pollution on daily deaths: a multicity case crossover analysis. Occupational Environmental Medicine 61, 956e961.
Szpiro, A.A., Sampson, P.D., Sheppard, L., Lumle, T., Adar, S.D., Kaufman, J.D., 2010. Predicting intra-urban variation in air pollution concentrations with complex spatio-temporal dependencies. Environmetrics 21, 606e631. Vienneau, D., de Hoogh, K., Beelen, R., Fischer, P., Hoek, G., Briggs, D., 2010. Comparison of land-use regression models between Great Britain and the Netherlands. Atmospheric Environment 44, 688e696. Wilson, J.G., Kingham, S., Pearce, J., Andrew, P., Sturmana, B., 2005. A review of intraurban variation in particulate air pollution: Implication for epidemiological research. Atmospheric Environment 39, 6444e6462. Yanosky, J.D., Paciorek, C.J., Suh, H.H., 2009. Predicting chronic fine and coarse particulate exposures using spatiotemporal models for the northeastern and midwestern United States. Environmental Health Perspectives 117 (4), 522e529. Yu, H.-L., Wang, C.-H., 2010. Retrospective prediction of intraurban spatiotemporal distribution of PM10 in Taipei. Atmospheric Environment 44, 3053e3065. Yu, H.-L., Wang, C.-H., Liu, M.-C., Kuo, Y.-M., 2011. Estimation of fine particulate matter in Taipei using land use regression and bayesian maximum entropy methods. International Journal of Environmental Research and Public Health 8, 2153e2169. Zanobetti, A., Schwartz, J., 2009. The effect of fine and coarse particulate air pollution on mortality: a national analysis. Environmental Health Perspectives 117 (6), 898e903.