Environmental Science & Technology Temporal CV R2 values were similarly calculated by ﬁrst averaging weekly predictions across locations (one mean per week) and then comparing week-speciﬁc average measurements and predictions. We assessed bias in model predictions using the slope from linear regression of the left-out measurements against model predictions, and examined the precision of model predictions by calculating the square root of the mean of the squared prediction errors (RMSPE), where prediction errors were deﬁned as predicted minus measured values. We also calculated the mean fractional bias (MFB) and mean fractional error (MFE) of model predictions according to the formulas provided in the Supporting Information. We used a generalized additive mixed model (GAMM) to describe spatial and temporal variation in weekly average GSH OP levels measured across greater London from 2002 to 2006. GAMMs were ﬁt using the gamm () function in the mgcv package32−34 for R35 with smoothing parameter selection performed using the maximum likelihood method. We also used linear mixed eﬀects models with a compound symmetric covariance structure for preliminary identiﬁcation of GIS covariates using PROC MIXED in SAS v9.2 (Cary, NC), selecting those that had the lowest (best) AIC or most signiﬁcant (smallest) p-value in separate models. Further model selection was performed in several phases, described below. First, we evaluated an initial space-time GAMM, GSH OPit = α + g(si) + h(t) + eit (1) eit ∼ N(0, σ 2e), Cov(eit , eit′) = σ 2 eλitt′; i = 1, ..., n where GSH OPit is the weekly average at location i (i = 1...I; I = 34) and week t, and g(si) is a two-dimensional spatial smooth function of the projected coordinates (easting and northing) at location i speciﬁed using thin-plate penalized splines (basis dimension k = 0.9I). The week index t ranged from 1 to T = 247 (corresponding to the weeks beginning December 31, 2001 and September 18, 2006, respectively). Though weekly averages were available at ≥4 sites for most weeks (73%), for two weeks no data were available. Because of this missing data in time, we did not specify an intercept for each time period αt. Instead, we speciﬁed h(t) as a thin-plate spline (basis dimension k = 0.9T) which represents the adjusted smoothed mean value across all sites for each week. Our choice of basis dimension was as large as practical, to allow the data to determine the complexity of the ﬁtted functions32−34. eit is assumed a mean-zero, Gaussian- distributed error term with constant variance σ e. σ2eλitt′ speciﬁes 2 covariances between errors at location i, with λitt′ parametrized as to induce either a compound symmetric (ﬁxed correlation between measurements at the same site and no correlation between measurements at diﬀerent sites: λitt′ = σ1/σ2e) or an autoregressive (AR(1); exponentially decreasing correlation in time between measurements at the same site and no correlation between measurements at diﬀerent sites: λitt′ = ρ|t−t′|) covariance structure. Next, we removed the spatial and temporal terms in separate 2 models. To the “best” (highest CV R ) of these three models (later referred to as the “referent” model), we then added GIS covariates one at a time in separate models (“univariate” models). In these univariate models, we evaluated both linear terms (βzi) and nonlinear smooth functions ( f(zi); using a thin-plate spline basis with basis dimension k = 4) of each GIS covariate zi. From these models, we obtained a set of GIS covariates that improved predictive accuracy compared to the referent model, and only Article these GIS covariates were considered in later stages of model selection. Next, we included these selected GIS covariates in ‘multi- variable’ GAMMs, removing any nonsigniﬁcant (i.e., p > 0.05) covariates one at a time, beginning with the least signiﬁcant. We also evaluated these models for parsimony, removing each remaining term one at a time, and keeping only those that improved predictive accuracy. We then evaluated whether the eﬀects of the remaining GIS covariates varied seasonally (i.e., winter, spring, summer, autumn) by including seasonal interaction terms, again removing nonsigniﬁcant terms. Finally, we evaluated whether any spatiotemporal interaction remained in the data, which could potentially explain additional variability in the outcome, by adding seasonal spatial terms, gseason(si), and weekly spatial terms, gt(si), separately, to the model. Details regarding space-time interaction and diﬀerent spatial models (thin-plate spline smoothing vs ordinary and simple kriging) are given in the Supporting Information. To obtain predicted weekly average GSH OP at unmeasured locations, we ﬁrst generated the selected GIS covariates at prediction locations and then used the ﬁtted GAMMs. GIS and meteorological predictors were truncated to their range among the monitoring locations to avoid extrapolation. Prediction involved only the ﬁxed-eﬀects of the model; the nonzero covariance assumed among the within-site errors provided control for overﬁtting of the smooth terms and other ﬁxed- eﬀects. PM10. We also developed a geostatistical model to predict weekly average PM10 concentrations at unmeasured locations anywhere within the study domain of central London, U.K. throughout 2002−2006. Because high-resolution spatial in- formation on annual-average PM10 levels was available from the hybrid emissions-dispersion/regression model developed for greater London,30,31 the goal of the weekly model was to leverage this rich spatial data and explicitly model weekly temporal variation. To do this, we used a similar approach as for GSH OP, but considered for inclusion as a covariate only the three-year mean of annual-average PM10 predictions from the hybrid emissions-dispersion/regression model (M0306HED). Cross- validation, model speciﬁcation with respect to spatiotemporal structure, and model prediction were performed for the PM10 model in the same way as for GSH OP. Exposure Metrics. Thus, two exposure metrics can be calculated from the predictive models: (1) PM10 and (2) GSH OP-weighted PM10 . This formulation allows direct comparison of the eﬀects of ÊPM10 and ÊGSH OP in health eﬀect analyses because the PM10 mass concentration component is the same in both. Modeling (GSH OP)(PM10) in units of μg m−3 would not aﬀord this direct comparison. ■ RESULTS GSH OP data were available at 34 monitoring sites (2118 weekly averages), and PM10 data at 66 sites (12 043 weekly averages), from 2002 to 2006. Weekly average GSH OP values in units of OP μg−1 were approximately normally distributed, with 0%, 25%, 50%, mean, 75%, and 100% values of 0, 0.45, 0.72, 0.76, 1.05, and 7614 dx.doi.org/10.1021/es3010305 | Environ. Sci. Technol. 2012, 46, 7612−7620