A comparison of statistical and machine learning methods for creating national daily maps of ambient PM$_{2.5}$ concentration

TL;DR

This paper compares statistical and machine learning methods for creating detailed daily maps of PM2.5 air pollution across the U.S., evaluating their accuracy and computational efficiency in various settings.

Contribution

It provides a comprehensive comparison of multiple exposure modeling approaches for PM2.5, highlighting their predictive performance and computational challenges.

Findings

Different methods vary in accuracy depending on location and season.

Machine learning approaches can offer improved predictions over traditional methods.

Computational complexity differs significantly among methods.

Abstract

A typical problem in air pollution epidemiology is exposure assessment for individuals for which health data are available. Due to the sparsity of monitoring sites and the limited temporal frequency with which measurements of air pollutants concentrations are collected (for most pollutants, once every 3 or 6 days), epidemiologists have been moving away from characterizing ambient air pollution exposure solely using measurements. In the last few years, substantial research efforts have been placed in developing statistical methods or machine learning techniques to generate estimates of air pollution at finer spatial and temporal scales (daily, usually) with complete coverage. Some of these methods include: geostatistical techniques, such as kriging; spatial statistical models that use the information contained in air quality model outputs (statistical downscaling models); linear regression modeling approaches that leverage the information in GIS covariates (land use regression); or machine learning methods that mine the information contained in relevant variables (neural network and deep learning approaches). Although some of these exposure modeling approaches have been used in several air pollution epidemiological studies, it is not clear how much the predicted exposures generated by these methods differ, and which method generates more reliable estimates. In this paper, we aim to address this gap by evaluating a variety of exposure modeling approaches, comparing their predictive performance and computational difficulty. Using PM$_{2.5}$ in year 2011 over the continental U.S. as case study, we examine the methods' performances across seasons, rural vs urban settings, and levels of PM$_{2.5}$ concentrations (low, medium, high).

Tables

Table 1. Table 1: Set of eleven predictors selected via best subset regression. In parentheses it is indicated the original database from which each predictor was obtained: North American Land Data Assimilation System (NLDAS) and North American Regional Reanalysis (NARR).

Percentage of Impervious Surface (Landsat)	Potential Evaporation (NLDAS)
Downward Longwave Radiation Flux (NLDAS)	Convective Available Potential Energy (NLDAS)
Relative Humidity (NLDAS)	V-Wind at 10m (NLDAS)
Pressure (NLDAS)	Dew-point Temperature (NARR)
Visibility (NARR)	Boundary Layer Height (NARR)
Temperature 150-180mb above ground (NARR)

Table 2. Table 2: Cross-validation results of daily PM 2.5 predictions: We compare the various statistical and machine learning methods using root mean squared error (“RMSE”), mean absolute deviation (“MAD”), correlation between predicted and observed values (“Corr”), and empirical coverage of the predictive 95% intervals. Both RMSE and MAD are in units of μ ​ g / m 3 𝜇 𝑔 superscript 𝑚 3 \mu g/m^{3} . The methods applied are ordinary least squares (“OLS”), inverse distance weighting (“IDW”), universal Kriging (“UK”), downscaler, random forests (“RF”), support vector regression (“SVR”) and Neural networks (“NN”). Methods use either CMAQ and/or other geographic covariates (“Covs”).

Method	RMSE	MAD	Corr	Coverage
OLS (CMAQ)	4.80	3.09	0.65	0.62
OLS (Covs)	4.63	2.97	0.68	0.79
OLS (CMAQ + Covs)	4.22	2.63	0.74	0.83
IDW	3.39	1.96	0.84	–
UK (CMAQ)	3.08	1.90	0.87	0.95
UK (Covs)	3.25	1.79	0.85	0.93
UK (CMAQ + Covs)	3.15	1.76	0.86	0.93
Downscaler (CMAQ)	3.10	1.70	0.87	0.94
RF (CMAQ + Covs)	3.41	2.09	0.84	0.96
SVR (CMAQ + Covs)	3.83	2.22	0.79	–
NN (CMAQ + Covs)	3.89	2.45	0.79	–

Table 3. Table 3: Cross-validation results by the number of nearby stations: Root mean squared error (RMSE) for PM 2.5 concentration predictions in μ ​ g / m 3 𝜇 𝑔 superscript 𝑚 3 \mu g/m^{3} stratified by the number of active monitoring sites within 50 miles of the prediction site for each day. Methods considered include: ordinary least squares (“OLS”), inverse distance weighting (“IDW”), universal Kriging (“UK”), downscaler, random forests (“RF”), support vector regression (“SVR”) and Neural networks (“NN”). Methods include either CMAQ and/or other geographic covariates (“Covs”).

Active nearby stations	5	5-9	10-19	$>$=20
OLS (CMAQ)	5.08	4.49	4.26	3.70
OLS (Covs)	4.88	4.22	4.34	3.84
OLS (CMAQ+Covs)	4.53	3.73	3.77	3.08
IDW	3.63	2.53	2.43	2.53
UK (CMAQ)	3.46	2.46	2.36	2.24
UK (Covs)	3.66	2.56	2.53	2.38
UK (CMAQ+Covs)	3.54	2.51	2.43	2.28
Downscaler(CMAQ)	3.51	2.42	2.31	2.22
RF (CMAQ + Covs)	4.54	3.86	3.64	3.20
SVM (CMAQ + Covs)	4.36	3.58	3.62	3.56
NN (CMAQ + covs)	4.24	3.45	3.20	3.11

Table 4. Table 4: Cross-validation results by urbanization: Root mean squared error (RMSE) for PM 2.5 concentration predictions in μ ​ g / m 3 𝜇 𝑔 superscript 𝑚 3 \mu g/m^{3} stratified by urbanization. Methods considered include: ordinary least squares (“OLS”), inverse distance weighting (“IDW”), universal Kriging (“UK”), downscaler, random forests (“RF”), support vector regression (“SVR”) and Neural networks (“NN”). Methods include either CMAQ and/or other geographic covariates (“Covs”). The last line reports the mean and standard deviation, in parenthesis, for the observed PM 2.5 concentration at urban versus non-urban monitoring sites.

Method	Urban	Non-Urban
OLS (CMAQ)	4.74	5.06
OLS (Covs)	4.62	4.69
OLS (CMAQ+Covs)	4.16	4.48
IDW	3.08	3.84
UK (CMAQ)	2.96	3.60
UK (Covs)	3.16	3.69
UK (CMAQ+Covs)	3.04	3.65
Downscaler (CMAQ)	2.99	3.62
RF (CMAQ + Covs)	4.18	4.47
SVM (CMAQ + Covs)	3.78	4.06
NN (CMAQ + Covs)	3.84	4.17
PM2.5	10.16 (6.21)	9.03 (6.01)

Table 5. Table 5: Cross-validation results by level of observed PM 2.5 concentration: Root mean squared error (RMSE) for PM 2.5 concentration predictions in μ ​ g / m 3 𝜇 𝑔 superscript 𝑚 3 \mu g/m^{3} stratified by PM 2.5 concentration level. Groupings are based on: whether PM 2.5 is less than 6 μ ​ g / m 3 𝜇 𝑔 superscript 𝑚 3 \mu g/m^{3} , 6-12 μ ​ g / m 3 𝜇 𝑔 superscript 𝑚 3 \mu g/m^{3} , and greater than or equal to 12 μ ​ g / m 3 𝜇 𝑔 superscript 𝑚 3 \mu g/m^{3} (12 μ ​ g / m 3 𝜇 𝑔 superscript 𝑚 3 \mu g/m^{3} is the EPA standard). Methods considered include: ordinary least squares (“OLS”), inverse distance weighting (“IDW”), universal Kriging (“UK”), downscaler, random forests (“RF”), support vector regression (“SVR”) and Neural networks (“NN”). Methods include either CMAQ and/or other geographic covariates (“Covs”). The last line reports the mean and standard deviation, in parenthesis, for the observed PM 2.5 concentration within each category (low, medium or high).

Method	Low	Med	High
OLS (CMAQ)	2.99	2.63	7.75
OLS (Covs)	2.88	2.76	7.36
OLS (CMAQ+Covs)	2.46	2.57	6.73
IDW	2.18	2.23	4.87
UK (CMAQ)	1.72	1.90	4.93
UK (Covs)	1.92	2.24	5.03
UK (CMAQ+Covs)	1.79	1.98	5.02
Downscaler (CMAQ)	1.77	1.89	4.97
RF (CMAQ + Covs)	3.05	2.61	6.51
SVM (CMAQ + Covs)	2.04	2.00	6.36
NN (CMAQ + Covs)	2.62	2.90	5.75
PM2.5	4.14 (1.21)	8.66 (1.71)	17.48 (5.81)

Table 6. Table 6: Evaluation of downscaler methods: The table reports whether for out data analysis the methods were executed in less than a few hours on a standard PC (“Fast”), whether they ranked near the best cross-validation error (“Accurate”), whether they were fit with standard software without many tuning parameters (“Easy”) and whether predictions were paired with measures of uncertainty (“Uncertainty”).

Method	Fast	Accurate	Easy	Uncertainty
Least squares	$✓$	X	$✓$	$✓$
Geographically-weighted regression	$✓$	$✓$	$✓$	X
Universal Kriging	$✓$	$✓$	$✓$	$✓$
Downscaler	X	$✓$	X	$✓$
Random forests	$✓$	$✓$	$✓$	$✓$
Support vector regression	$✓$	X	$✓$	X
Neural networks	$✓$	X	X	X

Table 7. Table 7: RMSE Cross-validation results by distance to closest other ACTIVE monitoring station on each day. Methods considered include: ordinary least squares (“OLS”), inverse distance weighting (“IDW”), universal Kriging (“UK”), downscaler, random forests (“RF”), support vector regression (“SVR”) and Neural networks (“NN”). Methods include either CMAQ and/or other geographic covariates (“Covs”). The last line in the table reports the mean and standard deviation, in parenthesis, of PM 2.5 concentration observed at monitoring sites within each substratum.

Methods	Closest Station 50 Miles	Closest Station $>=$ 50 Miles
OLS (CMAQ)	4.67	5.50
OLS (CMAQ+Covs)	4.06	5.05
OLS (Covs)	4.48	5.46
IDW	2.95	4.55
UK (CMAQ)	2.84	4.25
UK (CMAQ+Covs)	2.92	4.33
UK (Covs)	2.95	4.72
Downscaler	2.84	4.36
RF (CMAQ + Covs)	4.09	5.03
SVM (CMAQ + Covs)	3.61	4.99
NN (CMAQ + Covs)	3.74	4.72
PM2.5	10.07 (6.12)	9.36 (6.56)

Table 8. Table 8: RMSE Cross-validation results by East/West location in the US identified by the vertical line longitude=-100. Methods considered include: ordinary least squares (“OLS”), inverse distance weighting (“IDW”), universal Kriging (“UK”), downscaler, random forests (“RF”), support vector regression (“SVR”) and Neural networks (“NN”). Methods include either CMAQ and/or other geographic covariates (“Covs”). The last line in the table provides the mean and standard deviation (in parenthesis) for PM 2.5 concentration observed at monitoring sites within each category.

Methods	West	East
OLS (CMAQ)	7.11	4.80
OLS (Covs)	6.47	4.63
OLS (CMAQ+Covs)	6.20	4.22
IDW	4.70	3.22
UK (CMAQ)	4.74	3.08
UK (Covs)	4.87	3.25
UK (CMAQ+Covs)	4.84	3.15
Downscaler (CMAQ)	4.75	3.10
RF (CMAQ + Covs)	6.07	4.23
SVM (CMAQ + Covs)	5.68	3.83
NN (CMAQ + Covs)	5.53	3.90
PM2.5	8.87 (7.69)	9.97 (6.19)

Table 9. Table 9: RMSE Cross-validation results by season. Methods considered include: ordinary least squares (“OLS”), inverse distance weighting (“IDW”), universal Kriging (“UK”), downscaler, random forests (“RF”), support vector regression (“SVR”) and Neural networks (“NN”). Methods include either CMAQ and/or other geographic covariates (“Covs”). The last line in the table provides the mean and standard deviation, in parenthesis, for monitored PM 2.5 concentration during each season.

Methods	Winter	Spring	Summer	Fall
OLS (CMAQ)	5.93	4.16	4.37	4.60
OLS (Covs)	5.58	3.89	4.44	4.50
OLS (CMAQ+Covs)	5.15	3.49	3.95	4.15
IDW	3.95	2.69	2.75	3.38
UK (CMAQ)	3.81	2.52	2.64	3.23
UK (Covs)	3.96	2.60	2.74	3.57
UK (CMAQ+Covs)	3.88	2.58	2.72	3.31
Downscaler (CMAQ)	3.82	2.55	2.64	3.27
RF (CMAQ + Covs)	5.10	3.52	4.07	4.15
SVM (CMAQ + Covs)	4.75	3.28	3.49	3.71
NN (CMAQ + Covs)	4.67	3.21	3.47	4.13
PM2.5	10.59 (7.05)	9.47 (5.81)	10.82 (5.93)	9.01 (5.74)

Equations

$$w_{0i}={\left(\frac{1}{d_{0i}}\right)}^{\varphi },$$

$${\hat{Y}}_t(\text{mbox}{\mathbf{s}}_0)=\sum _{i=1}^n{W}_{0i}{Y}_t(\text{mbox}{\mathbf{s}}_i)\text{mbox}where{W}_{0i}=\frac{w_{0i}}{{\sum }_{j=1}^n{w}_{0j}}.$$

$$\mbox{\boldmath$\lambda$}=\mbox{arg min}_{{\tilde{}\mbox{\boldmath$\lambda$}}}E\left[Y_{t}(\mbox{\bf s}_{0})-\sum_{i=1}^{n}{\tilde{\lambda}}_{i}Y_{t}(\mbox{\bf s}_{i})\right]^{2}\qquad\mbox{such that}\qquad E\left[\sum_{i=1}^{n}{\tilde{\lambda}}_{i}Y_{t}(\mbox{\bf s}_{i})\right]=E\left[Y_{t}(\mbox{\bf s}_{0})\right]$$

$${\mu }_t(s)={\gamma }_{0,t}+{\gamma }_{1,t}{Z}_t(\text{mbox}\mathbf{s})$$

$$\mu_{t}(s)=\mbox{\bf X}_{t}(\mbox{\bf s})\mbox{\boldmath$\beta$}_{t}$$

$$\mu_{t}(s)=\mbox{\bf X}_{t}(\mbox{\bf s})\mbox{\boldmath$\beta$}_{t}+Z_{t}(\mbox{\bf s})\gamma_{1,t}$$

$$Y_t(\text{mbox}\mathbf{s})={\beta }_{0,t}(\text{mbox}\mathbf{s})+{\beta }_{1,t}{Z}_t(\text{mbox}\mathbf{s})+{ϵ}_t(\text{mbox}\mathbf{s}){ϵ}_t(\text{mbox}\mathbf{s})\stackrel{iid}{\sim }N(0,{\tau }^2)$$

$$\text{mbox}Cov\left({\beta }_{0,t}(\text{mbox}{\mathbf{s}}_i),{\beta }_{0,t}(\text{mbox}{\mathbf{s}}_j)\right)={\sigma }_0^2\exp \left(-\frac{d_{ij}}{{\varphi }_0}\right)$$

$${\hat{Y}}_0=\sum _{i=1}^{n}K(\text{mbox}{\mathbf{X}}_i,\text{mbox}{\mathbf{X}}_0)({\alpha }_i-{\alpha }_i^{\star })$$

$$\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n({\alpha }_i-{\alpha }_i^{\star }{)}^T({\alpha }_j-{\alpha }_j^{\star })K(\text{mbox}{\mathbf{X}}_i,\text{mbox}{\mathbf{X}}_j)+ϵ\sum _{i=1}^n({\alpha }_i+{\alpha }_i^{\star })+\sum _{i=1}^n{Y}_i({\alpha }_i-{\alpha }_i^{\star })$$

$$\text{mbox}Outputlayer:$$

$$\text{mbox}Hiddenlayer:$$

$$\text{mbox}Inputlayer:$$

Full text

Key words: Air pollution exposure assessment; spatial modeling; kriging; Bayesian hierarchical modeling; downscaling; machine learning; random forest; support-vector regression; neural network.

1 Introduction

Accurate exposure assessment plays an essential role in the success of any environmental health study. Past air pollution epidemiological studies regularly utilize ambient air quality measurements from large monitoring networks to estimate population and individual exposures. However, measurements from these networks are spatially sparse, temporally incomplete, and preferentially located in settings with large populations and high pollution levels. There is increasing interest in developing methods to retrospectively estimate air pollution levels at fine spatial scales and with complete spatial-temporal coverage in order to minimize exposure measurement error (Alexeeff et al., 2015), support more spatially-resolved health effect analyses (Kloog et al., 2012; Hao et al., 2016), and perform impact assessments in low- and middle-income settings (Shaddick et al., 2018). Advances in Geographical Information Systems, remote sensing, and numerical model simulations have further contributed to a proliferation of modeling approaches to estimate air pollution over the past decade.

This paper aims to address an important gap in the current literature. When models to estimate ambient air pollution exposure are being developed, they are typically only compared to simpler models within the same modeling paradigms. Examples include comparisons between regression models with different predictors, especially in land use regression models (Tang et al., 2013; Wang et al., 2016); between geostatistical models with different spatial dependence structures (Reich et al., 2011); or between different algorithms among the machine learning toolbox (Singh et al., 2013; Reid et al., 2015). There has been very limited cross-paradigm comparisons (Adam-Poupart et al., 2014; Yu et al., 2018), likely due to the analytic effort and expertise required to carry out the different approaches. In particular, there is a lack of comparison between machine learning and advanced geostatistical approaches. As exposure modeling becomes increasingly complex and computationally demanding, there is a pressing need to better understand the advantages and limitations associated with different modeling approaches. However, synthesis of results from existing studies is challenging due to the use of different data sources and inconsistent criteria for evaluating model performance.

The main objective of this paper is to critically evaluate different exposure modeling approaches via a case study of estimating regional daily ambient fine particulate matter of aerodynamic diameter less than 2.5 $\mu m$ (PM2.5 ) across the contiguous United States. We considered major classes of recent approaches that can provide complete spatial and temporal coverage using monitoring measurements, meteorology, land use variables, and numerical model simulations. These approaches include geographically weighted regression(van Donkelaar et al., 2015) , statistical downscaling (Berrocal et al., 2010), and machine learning methods, including random forest (Hu et al., 2017), support vector machine (Liu et al., 2017), and neural networks (Di et al., 2016). Air quality estimates from these approaches have already been used in health effect and health impact analyses (Chang et al., 2011; Strickland et al., 2016; Di et al., 2017). Our evaluations focuses specifically on a model’s spatial predictive ability at locations without monitors and on computational demand.

A related issue that has hindered cross-comparison is the limited availability of data and analytic code from exposure modeling studies. This is in contrast to efforts in other fields, such as genetics, genomics and neuroimaging, where established publicly available and well-documented datasets are available for methods development and comparison (Mailman et al., 2007; Leinonen et al., 2010; Van Essen et al., 2013). Sophisticated exposure models and algorithms often require considerable efforts in input associated with data processing and parameter tuning to maximize model performance. A notable example is the air pollution modeling framework for 2-week averages from the MESA Air study (Keller et al., 2015) where an R package has been made available (Lindström et al., 2012). Hence, another contribution of this study is the documentation of all data and analytic code. We also only considered methods that can be readily implemented using packages from the statistical R software (R Core Team, 2018). This will ensure better reproducibility, allow for future evaluation of other methods, and facilitate adaptation of these modeling methods by other research groups.

2 Data sources

2.1 PM2.5 Monitoring data

For year 2011, daily PM2.5 measurements were obtained from the Air Quality System database of the US Environmental Protection Agency. A total of 829 monitors from state, local, and tribal air pollution monitors are available across our study domain of the contiguous United States. Monitors typically take measurement every 3 or 6 days, with only 15% of monitors sampling almost daily. Figure 1 shows a histogram of the number of active monitors each day: most days have less than 200 measurements across the contiguous US.

2.2 Computer model output

Computer models for air quality use information on emission sources to simulate the fate and transport of air pollutants by taking in account complex atmospheric chemistry and physics. There is increasing use of computer model outputs to estimate air pollution levels at locations without monitoring data (Berrocal et al., 2010). We obtained Community Multiscale Air Quality model (CMAQ; Byun and Schere, 2006) version 4.7 runs performed by the EPA at 12km spatial resolution across the contiguous US for year 2011. At this resolution, the contiguous US is covered by a grid of 299 by 459 cells. Each PM2.5 monitor was linked to the CMAQ grid cell that contains it. Figure 2 compares annual average PM2.5 concentrations obtained from CMAQ simulations and AQS measurements.

2.3 Meteorological and land use data

Daily average meteorological fields were obtained from the North American Land Data Assimilation System (NLDAS) and the North American Regional Reanalysis (NARR). Meteorological variables include temperature at 2m, 30m, 150-180mb above ground, downward longwave and shortwave radiation flux, relative humidity at 2m, boundary layer height, visibility, dew-point temperature, potential evaporation, potential evaporation, convective available potential energy, pressure at 2m, 10m, and 30m, and U- and V-wind speed at 10m. We also obtained the following land use variables: elevation from the US Geological Survey, major roadway lengths and percentage of forest cover from the National Land Cover database, population density from the U.S. 2010 Census Bureau at the census tract level, percentage of impervious surface from Landsat, and PM2.5 primary emission sources from the 2011 US EPA National Emissions Inventory Facility Emissions report. Although all the above mentioned variables were available at different spatial resolution (point level, in some cases, grids, etc.), they were re-gridded to the CMAQ 12km grid. Specifically, for predictors with finer spatial resolutions than CMAQ (e.g. elevation), we used the average of cell values that intersect with the CMAQ cell. For predictors with coarser spatial resolution than CMAQ (e.g. meteorology), the value of the nearest grid cell centroid was used.

2.4 Variable selection

It is well established that meteorological variables and land use characteristic are good predictors of PM2.5 concentration, as land use regression models indicate (Hoek et al., 2008). To identify which of the available predictors should be used in the exposure models considered in our subsequent model comparisons, we first performed variable selection via cross validation and best subset regression (Kutner et al., 2005). For this, we randomly split the monitors into 5 folds. For each fold, we performed best subset regression using data from the other 4 folds and predicted PM2.5 concentration at the hold-out fold. For each hold-out set, we computed the root mean square error (RMSE), comparing the predicted values with the held-out PM2.5 observations. We averaged the RMSEs across the five folds and selected the model which yielded the best predictive performance. This identified a model with 11 predictors: including an additional predictor only improved prediction by less than 0.1%. After having determined the number of predictors to include in the model, we identified the set of 11 predictors again via best subset regression. In this second step, best subset regression was performed on the full data using all meteorological variables and land use covariates. The 11 selected predictors are shown in Table 1. The set consisted mostly of meteorological variables, likely due to our interest in estimating daily PM2.5 concentrations.

3 Methods

3.1 Statistical methods

Let $Y_{t}(\mbox{\bf s})$ be the monitor measurement of PM2.5 at spatial location s and day $t$. Each observation is associated with land-use covariates and meteorological variables at s, $\mbox{\bf X}_{t}(\mbox{\bf s})$, and the CMAQ output at the grid cell that contains s, $Z_{t}(\mbox{\bf s})$. We assume monitoring data are available at $n$ sites $\mbox{\bf s}_{1},...,\mbox{\bf s}_{n}$ and for days $t=1,...,T$. The distance (in km) between location $\mbox{\bf s}_{i}$ and $\mbox{\bf s}_{j}$ is denoted $d_{ij}$. The objective of all methods is to make a prediction $\widehat{Y}_{t}(\mbox{\bf s}_{0})$ of PM2.5 concentration at a location $\mbox{\bf s}_{0}$ that does not have a monitor, with prediction uncertainty quantified through $\mbox{Var}[\widehat{Y}_{t}(\mbox{\bf s}_{0})]=\widehat{v}_{t}(\mbox{\bf s}_{0})$.

3.1.1 Inverse distance weighting

Inverse distance weighting (IDW) does not use the covariates and simply uses a weighted average of nearby observations as the prediction. The intuition behind IDW is that the weights for observations near the prediction location should be higher than the weights for observations far from the prediction location. The weight assigned to the observation at location $\mbox{\bf s}_{i}$ for prediction at location $\mbox{\bf s}_{0}$ decays as the distance $d_{0i}$ between $\mbox{\bf s}_{0}$ and $\mbox{\bf s}_{i}$ increases,

[TABLE]

where the rate of decay is controlled by $\phi>0$, which we select using cross validation.

The IDW prediction is

[TABLE]

The weights $W_{0i}$ are standardized to sum to one so that predictions are unbiased. For large $n$, computational speed can be increased by setting $w_{0i}=0$ for locations with large $d_{i0}$’s to reduce the number of terms in the sum in (2). Because no statistical model is assumed, there is no natural variance estimate, $v_{t}(\mbox{\bf s})$.

3.2 Geostatistical methods

In geostatistics, the problem of generating predictions of a continuous spatial process at unsampled locations is typically handled via Kriging. As with IDW, in Kriging the prediction $\hat{Y}_{t}(\mbox{\bf s}_{0})$ is a linear combination of the observations $Y_{t}(\mbox{\bf s}_{1}),\ldots,Y_{t}(\mbox{\bf s}_{n})$ with weights \mbox{\boldmath\lambda}=(\lambda_{1},\lambda_{2},\ldots,\lambda_{n}), that is, $Y_{t}(\mbox{\bf s}_{0})=\sum_{i=1}^{n}\lambda_{i}Y_{t}(\mbox{\bf s}_{i})$. The weight vector $\lambda$ is in turn determined by finding the set of $\lambda$’s that yield the Best Linear Unbiased Predictor (BLUP), e.g.

[TABLE]

The expression for the weights depends on the assumptions made on $Y_{t}(\mbox{\bf s})$, such as whether the process is second-order stationary (i.e., the covariance depends only on the separation between points) and whether its mean varies spatially and depends on covariates. For computational tractability, we assume a stationary covariance function (Cressie, 1993). In addition, given that we are generating predictions of PM2.5 concentration across the entire United States it is safe to assume that the mean concentration varies spatially as function of covariates. Thus, we will generate our predictions using Universal Kriging (Cressie, 1993). The latter, differently from Ordinary Kriging (Cressie, 1993), includes covariates in the expression of the mean function.

3.2.1 Universal Kriging

In Universal Kriging we assume that, for each day $t$, $Y_{t}(\mbox{\bf s})$ is a Gaussian process with a pre-determined covariance function and a mean $\text{E}[Y_{t}(\mbox{\bf s})]=\mu_{t}(\mbox{\bf s})$, linear function of spatially-varying covariates. In our implementation, as we are dealing with data over time, we assume that the spatial covariance does not change in time and the correlation between PM2.5 concentration at two sites decays exponentially with distance. In other words, we use a stationary exponential covariance function with covariance parameters that are constant in time (Cressie, 1993). We compare different scenarios for $\mu_{t}(\mbox{\bf s})$:

$\mu_{t}(\mbox{\bf s})$ depends only on the CMAQ output $Z_{t}(\mbox{\bf s})$, e.g.

[TABLE] 2. 2.

$\mu_{t}(\mbox{\bf s})$ depends only on meteorological and land-use covariates $\mbox{\bf X}_{t}(\mbox{\bf s})$, e.g.

[TABLE] 3. 3.

$\mu_{t}(\mbox{\bf s})$ depends on both the CMAQ output $Z_{t}(\mbox{\bf s})$ and the meteorological and land-use covariates $\mbox{\bf X}_{t}(\mbox{\bf s})$, e.g.

[TABLE]

To implement the different universal Kriging models, we use a two-stage procedure and use the gstat and geoR (Ribeiro Jr and Diggle, 2018) packages in R (R Core Team, 2018). Specifically, we first estimate the spatial covariance parameters via weighted least squares (WLS) using the gstat package, fitting an exponential semi-variogram to the empirical semi-variogram of the time-averaged residuals of the linear regressions implied by (3), (4) and (5), respectively. Using the spatial covariance parameters estimated via WLS as initial values, we used the geoR package (Ribeiro Jr and Diggle, 2018) and maximum likelihood to iteratively estimate the regression parameters parameters \mbox{\boldmath\beta}_{t}, $\gamma_{0,t}$ and $\gamma_{1,t}$ and the covariance parameters. Using such estimates, we then generate predictions of $Y_{t}(\mbox{\bf s}_{0})$ via Universal Kriging separately for each day. Because the data are assumed to be Gaussian, Kriging predictions are accompanied by predictions variance estimates ${\hat{v}}_{t}(\mbox{\bf s})$.

3.3 Downscaler

The downscaler model generates predictions of air pollution concentration at any location s by exploiting the relationship between the observed concentration measured at a monitor and the estimated air pollution concentration generated by an air quality model, CMAQ. The simplest downscaler model (Berrocal et al., 2010) relates the (log or square root of the) observed concentration at location s on day $t$, $Y_{t}(\mbox{\bf s})$, to the CMAQ output, $Z_{t}(\mbox{\bf s})$, at the grid cell that contains location s via a linear regression model with spatially and temporally-varying coefficients. Choosing the form of the downscaler model that yields the best predictive performance in various experiments, we adopt the following version of the downscaler model

[TABLE]

where $\beta_{0,t}(\mbox{\bf s})$ is a spatially and temporally-varying intercept term, while $\beta_{1,t}$ indicates a slope term, constant in space but varying in time.

As on a given day $t$, we postulate that sites located nearby have a similar intercept term, to account for spatial dependence, we model $\beta_{0,t}(\mathbf{s})$ as a stationary spatial processes with an exponential correlation function. That it, for any two pair of sites $\mbox{\bf s}_{i}$ and $\mbox{\bf s}_{j}$:

[TABLE]

In (7), $d_{ij}$ denotes the distance between sites $\mbox{\bf s}_{i}$ and $\mbox{\bf s}_{j}$, while $\sigma^{2}_{0}$ and $\phi_{0}$ represent, respectively, the spatial variability in $\beta_{0,t}(\mathbf{s})$, and the rate at which the spatial correlation vanishes.

The downscaler model in (6) is fit within a Bayesian framework (Gelman et al., 2013), and thus its specification is completed once prior distributions for all the model parameters are provided. Specifically, we assume that, for each day $t$, the spatially varying intercept $\beta_{0,t}(\mathbf{s})$ admits a constant mean $\beta_{0,t}$, which is interpreted as the overall additive calibration of the CMAQ output. In contrast, $\beta_{1,t}$ represents the overall multiplicative calibration term for CMAQ. The two global calibration parameters $\left(\beta_{0,t}\;\beta_{1,t}\right)^{\prime}$ are assumed to be independent in time and follow a bivariate normal distribution with mean $\left(0\;1\right)^{\prime}$ and with a diagonal, covariance matrix with large prior variances. The two variance parameters, $\tau^{2}$ and $\sigma^{2}_{0}$, representing, respectively, the non-spatial and spatial variability in PM2.5 concentration, are provided with vague Inverse Gamma priors, while Uniform prior on the interval $(0.0001km,0.1km)$ is placed on the decay parameters $\phi_{0}$.

Inference on model parameters is carried out through the posterior distribution, which is approximated using an MCMC algorithm (Gelman et al., 2013) we ran for 10,000 iterations, with the first 5,000 discarded for burn-in. For the case study discussed in the paper, we make a slight modification to the downscaler model, allowing the spatial dependence parameters ($\sigma^{2}$, $\phi$) and the non-spatial variance $\tau^{2}$ to vary with time. With this modification we can implement the downscaler model, using the spBayes (Finley et al., 2015) package in R.

Predictions of air pollution concentrations on day $t$ at an unsampled location $\mbox{\bf s}_{0}$ and uncertainty estimates for the predictions are obtained using the posterior predictive distribution of $Y_{t}(\mbox{\bf s}_{0})$ given the observed data. Specifically: we take as predicted concentration, the median of the posterior predictive distribution, while we characterize uncertainty in the prediction via the $95\%$ predictive interval, whose upper and lower limits are given by the 2.5-th and 97.5-th percentile of the posterior predictive distribution. Finally, we use the sample variance ${\hat{v}}_{t}(\mbox{\bf s})$ of the predicted values to quantify the variance in the predictions.

3.4 Machine learning methods

The spatial regression models in Section 3.2 and 3.3 represent the data-generating process with a small number of interpretable parameters. In contrast, machine learning algorithms have countless uninterpretable parameters and are thus essentially black-box prediction machines. However, these algorithms are flexible, general and have often outstanding predictive performance. Below we briefly describe the machine learning algorithms used in our comparison; for a detailed description of these and other machine learning algorithms see James et al. (2013).

For application of machine learning methods we do not explicitly model spatial correlation, although spatial coordinates and spatial covariates are used as predictors in these regression models. Since we do not model spatial correlation we denote the observed PM2.5 concentration and covariates relative to observation $i=1,...,n$ ($n$ is the combined sample size over space and time) as $Y_{i}$ and $\mbox{\bf X}_{i}$, respectively. The $p=15$ covariates in $\mbox{\bf X}_{i}$ are the longitude, latitude, day of year, the covariates in Table 1 and CMAQ model output.

3.4.1 Random forests

A random forest (Breiman, 2001) is an ensemble of regression trees (Morgan and Sonquist, 1963). A regression tree is grown via recursive partitioning of the covariate space called “leaves”, and then fitting separate linear models within each leaf. For each split of the covariate space, $m$ input variables are selected at random as candidates for splitting. Randomizing over the covariates decreases the correlation between trees and improves the prediction accuracy of the ensemble. The intuition is that interactions are likely present when considering the entire covariate space, but when considering small subregions of the covariate space simple linear models likely fit well within each subregion.

To construct a random forest for PM2.5 concentration, many regression trees are grown with the prediction being a weighted average of the predictions over the trees. Unlike most machine learning methods, it is possible to quantify the uncertainty in the prediction $v_{t}(\mbox{\bf s})$ using the sample variance of the trees predictions. Random forest is a desirable data mining method because it is easily understood and it is computationally efficient for very large samples. The number of trees used for prediction is selected by out-of-bag error (Hastie et al., 2001), a cross validation technique where the observations that are not selected in a tree are then used to estimate the error. For our case study, we implemented random forest using the R package randomForest (Liaw and Wiener, 2002). We used the default values $m=p/3$ at each split, where $p=15$ is the number of covariates. Based on out-of-bag error, a random forest with 500 trees is used as the final predictive model.

3.4.2 Support Vector Regression (SVR)

Nonlinear support vector machines (Cortes and Vapnik, 1995) are most well known in machine learning as a classification tool, but support vectors can also be used for regression. In SVR a hyperplane is optimized to be within a certain threshold of the selected data, called the support vectors, and the hyperplane is used for regression prediction. The predicted PM2.5 concentration for an observation with covariates $\mbox{\bf X}_{0}$ is

[TABLE]

where $\alpha_{i}$ and $\alpha_{i}^{\star}$ are the support vectors and $K(\mbox{\bf X}_{i},\mbox{\bf X}_{j})=\exp(-|\mbox{\bf X}_{i}-\mbox{\bf X}_{j}|^{2}/\gamma)$ is the radial basis kernel function. The support vectors $\alpha_{i}$ and $\alpha_{i}^{*}$ are estimated as the solution to the convex optimization

[TABLE]

so that for all $i$, $0<\alpha_{i}<C$, $0<\alpha_{i}^{\star}<C$ and $\sum_{i=1}^{n}(\alpha_{i}-\alpha_{i}^{\star})=0$. The constant $C$ is the box constraint that defines the trade off between penalty for observations further than $\epsilon$ away from the hyperplane and smoothness in the prediction hyperplane.

This type of model is known to yield good predictions, especially in high dimensional space. However, we note that the calculation of the optimum solution to (8) is computationally expensive in large datasets and there is no straightforward method to estimate prediction variances. We implement support vector regression using the R package e1071 (Dimitriadou et al., 2006) and the function svm.

3.4.3 Neural networks

Neural networks are machine learning methods that have successfully been used to predict air pollution levels (Di et al., 2016). They are attractive because they can handle massive data and can model nonlinearity and interactions. We use a multilayer perceptron network (Rumelhart et al., 1986) that consists of an input layer, several hidden layers and an output layer. For an observation with covariates $\mbox{\bf X}_{0}=(X_{01},...,X_{0p})$ the predicted value with a single hidden layer is

[TABLE]

where $L_{1}$ and $L_{2}$ are the number of neurons in each layer, $f$ is the activation function, and $b_{l}^{k}$ and $w_{jl}^{k}$ are the bias and weights estimated to minimize mean squared error in the training data.

Fitting a neural network requires tuning the number of hidden layers, the numbers of neurons in each layer, and the activation function. Increasing the size of the network provides more flexibility but also may lead to overfitting. Therefore, we also tune the number of epochs, drop-out rate, learning rate, and minibatch size (Goodfellow et al., 2016). The model is fit using the deep learning package keras (Chollet et al., 2015) in R. We tried relu and sigmoid activation functions and networks with 1 to 3 hidden layers, with 500 to 2000 neurons in the the first layer, and a reduced number of neurons in each of the next layers. Further tuning with the following possible combinations were considered: number of epochs ranging from 50 to 200, drop-out rate ranging from 0.1 to 0.5, learning rate ranging from 0.0005 to 0.01, and minibatch size ranging from 128 to 1024. Based on five-fold cross-validation, the final selected model has 2 hidden layers with 2000 and 100 neurons for the first and second layer, respectively, relu activation function $f(x)=\text{max}\{0,x\}$, 200 epochs, 0.4 drop-out rate, 0.003 learning rate and 1024 minibatch size.

4 Results

We compare the statistical and machine learning methods described in Section 3 by assessing their predictive performance. We use five-fold cross validation, and randomly sample, without replacement, 20% of the sites to be in the test set in each of the five folds. Within each fold, we train each of the models on the data from the remaining 80% of the sites and we generate predictions at the test sites. We evaluate the predictive performance of the various methods in terms of root mean squared error (RMSE), mean absolute deviation (MAD), correlation between predicted and observed PM2.5 concentrations, and coverage of the 95% prediction intervals, all averaged across days and sites. Table 2 presents the results.

Predictions obtained via a linear regression model fit separately to PM2.5 concentration each day serves as a baseline for comparison, and it yields an RMSE of 4.22$\mu g/m^{3}$ when both CMAQ and the selected covariates listed in Table 1 are included in the model. Machine learning methods yield predictions that are substantially improved over those obtained via least squares, by accounting for potential non-linear and interaction effects of the covariates. Of the machine learning algorithms, Random Forest has the smallest RMSE (3.41$\mu g/m^{3}$). All of the statistical methods that explicitly model the spatial correlation in the data have smaller RMSE than the machine learning methods. Inverse distance weighting which generates prediction by simply taking a local average of the available monitoring data has an RMSE of 3.39$\mu g/m^{3}$. The best method in terms of RMSE is Universal Kriging (UK) when only CMAQ is used as a covariate. The spatial downscaler model yields a comparable RMSE (3.10$\mu g/m^{3}$ vs 3.08$\mu g/m^{3}$) but has a smaller MAD.

Table 3 examines whether the quality of the predictions is improved based on the number of nearby monitoring sites with data available on the day of the prediction. The table presents predictive performance results stratified by the number of monitoring stations within a 50 miles radius of the prediction location with same-day PM2.5 measurements. As with overall predictions, the geostatistical methods have the smallest prediction error for all numbers of active training sites. For test sites with many active neighbors, the downscaler model has slightly smaller prediction error than UK, but UK performs well in all cases. The Appendix includes analogous tables that compare the quality of the predictions generated by the various methods when we stratify by distance to nearest monitoring location, longitude, and season, while this section presents results stratified by urbanization (Table 4) and PM2.5 magnitude (Table 5). For all methods, prediction errors are smaller at urban sites, sites near other monitoring stations, sites in the Eastern US and sites with low PM2.5 concentrations. Prediction errors are also lower in the Spring and Summer than in the Fall and Winter. However, the ranking of the prediction methods remains similar across all of these strata, with Universal Kriging with CMAQ as only covariate yielding the best results. In addition, all of the spatial statistical methods generate predictions with accurate uncertainty quantification: the 95% prediction intervals constructed using spatial statistical methods have coverage at or near the nominal 95% level.

Spatial predictions of log PM2.5 over the entire contiguous United States can be seen in Figures 3 and 4 for the dates of January 1, 2011 and August 1, 2011, respectively. As the figures indicate, PM2.5 concentration tends to be higher in the winter (January 1) than in the summer (August 1), with areas in the North East and the North West, especially along the coasts, having higher PM2.5 levels. In contrast, on August 1, it is the interior region in the Eastern United States that experienced higher levels of PM2.5 compared to the rest of the country.

The predictive surfaces of PM2.5 concentration yielded by Universal Kriging and the downscaler model present similar spatial gradients on both days. In both cases, the surfaces are generally smooth, but presents fine spatial variation particularly in areas where more monitoring data is available. Although geographical space is utilized as another independent variable in Random Forest (through the use of latitude and longitude as covariates), the spatial maps of PM2.5 concentration yielded by Random Forest are typically smooth and tend to display similar patterns to those observed in the predictive surfaces generated by the spatial statistical methods (UK and the downscaler). However, Random Forest’s predictive surfaces generally have less variability than the predictive surfaces obtained using spatial statistics methods, and tend to generate predictive distributions that do not have long tails as the predicted PM2.5 values. As a result, predictions tend to shrink closer to the mean. Finally, while the IDW method seems to yield a spatial surface that is closer to that obtained using the other methods on August 1, the predictive surface is overly smooth and does not show fine-scale spatial variability. In particular, since IDW generates predictions using only the information available from nearby monitors on each day, on days where few observations are available the predictions tend to be the same across large subregions. This leads to maps with a blocky appearance, and sudden, unrealistic, sharp changes in regions where more monitoring data is available. As such, predictive PM2.5 maps generated via IDW are not recommended as maps capturing the spatial distribution of PM2.5 concentration in the contiguous United States.

5 Discussion

In this paper, we present a review of commonly used statistical and machine learning methods for air pollution exposure assessment with the goal of improving our understanding of the fundamental differences among these methods. The contribution of this paper is to provide a rigorous comparison of these exposure modeling methods and make analytic code available to the broader scientific community to implement these techniques through a GitHub repository https://github.com/yawenguan/DataFusion. In our case study on PM2.5 concentration in the contiguous United States, we have found that the downscaler model and Universal Kriging yield better predictive performance than machine learning algorithms, potentially due to the fact that these methods explicitly account for the spatial dependence in PM2.5 concentration. This may hold true for other environmental exposure (e.g. ozone, or other pollutants) where spatial dependence is not negligible.

Machine learning algorithms are attractive methods for analyzing large data sets due to their computational speed and easy implementation for massive data, partly driven by the recent availability of highly optimized computing software. In this review paper, we have chosen Random Forest, Support Vector Regression and Neural Network for comparison, because these methods have already been used for exposure modeling (Hu et al., 2017; Liu et al., 2017; Reid et al., 2015) and software within R is readily available.

In our case study, we have found that the machine learning methods yield worse predictive performance than the statistical methods that account for spatial dependence. An explanation for this finding can be due to the number of variables used in our case study, which is smaller than the number of variables typically used in applications of these methods. Another reason for the underperformance of machine learning techniques might be due to the fact that these methods do not account explicitly for spatial dependence, which we have shown to be important for exposure modeling. For instance, Table 7 compares RMSE cross-validation results of daily PM2.5 predictions by distance to the closest active monitoring sites. When there are monitoring stations nearby the prediction site ($<$50 miles), the RMSE from the best machine learning method is worse than the best spatial statistical model: 3.61$\mu g/m^{3}$ for the best machine learning method versus 2.84$\mu g/m^{3}$ for the best spatial statistical model (ratio=1.27); on the other hand, when the monitoring stations are far ($\geq$50 miles) from the prediction site, the RMSE from the best machine learning method is 4.72$\mu g/m^{3}$, which is comparable to the RMSE of 4.25$\mu g/m^{3}$ for the best geostatistical model (ratio=1.08). The better performance by spatial statistical models in the first case is likely due to the borrowing of strength from neighboring observations. This also indicates that to improve the predictive performance of machine learning methods, algorithms that leverage values from neighboring sites for prediction, such as it is done in convolution neural net, might be useful for this application.

Acknowledgments

The authors thank Drs. Xuefie Hu and Yang Liu for providing the land use and meteorology data, and James Mulholland and Niru Senthilkumar for providing the CMAQ and air quality data. This material was based upon work partially supported by the National Science Foundation under Grant DMS-1638521 to the Statistical and Applied Mathematical Sciences Institute. BJR is also supported by NIH R01ES027892, DOI-JFSP 14-1-04-9 and KAUST 3800.2; VJB is supported by NIH P30ES017885. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Intitutes of Health and the National Science Foundation, U.S.

Appendix

This Appendix presents maps with spatial predictions of PM2.5 concentrations in $\mu g/m^{3}$ over the entire contiguous United States on January 1, 2011 and August 1, 2011, respectively. In addition, it presents a comparison of the various methods for PM2.5 concentration estimation, similar to those discussed in Section 4.

Specifically, Table 6 provides a simplified view on characteristics of the various methods compared in this case study with respect to computational speed, prediction accuracy, implementation difficulty and whether prediction uncertainty can be obtained for the various statistical and machine learning methods. We note that the computational speed varies across machines and it depends heavily on the software and code optimization. Hence, the computational speed here only serves as a rough guideline of the computational complexity and the comparisons are made based on the authors’ experience.

Additional tables compare the quality of the predictions by the various methods with stratification by distance to nearest monitoring location (Table 7), longitude (Table 8), and season (Table 9).