ensembleDownscaleR: R Package for Bayesian Ensemble Averaging of PM2.5 Geostatistical Downscalers

1. Introduction

Exposure to ambient air pollution is a well established risk factors for multiple adverse health outcomes. Particulate matter of size less than 2.5 $\mu$m in diameter (${\text{PM}}_{2.5}$) has been shown to be especially harmful due to its ability to penetrate deep into the respiratory system [1,2,3,4,5,6]. ${\text{PM}}_{2.5}$ represents a chemically diverse mixture of pollutants and its major sources include electricity generation, motor vehicles, and wildland fires. Accurate estimation of ${\text{PM}}_{2.5}$ concentration is an important component of air pollution health research [7,8].

Network of monitors established for regulatory purposes provide precise, direct measurements of ${\text{PM}}_{2.5}$. However these monitors are relatively expensive to install and maintain, and therefore are both spatially sparse and preferentially located in prioritized areas. There exist two additional data sources that are highly correlated with ${\text{PM}}_{2.5}$, which enable the development of methods to predict ${\text{PM}}_{2.5}$ beyond locations with ground monitors. The first is outputs from chemical transport model (CTM) which provide simulations of air quality based on emission, physical processes and meteorological data. Gridded CTM simulations provide complete spatial coverage, but require bias-correction with observations due to uncertainties in model input, discretization, and other sources of errors [9].The second is satellite measurements of aerosol optical depth (AOD), which provide measures of aerosol for the entire atmospheric column at finer spatial resolution than CTM [10,11,12]. But AOD also requires bias-correction and is subject to high missingness due to cloud cover and retrieval error [13,14,15].

An important and active area of air pollution exposure assessment research is the development of methods that can utilize both AOD and CTM to predict ${\text{PM}}_{2.5}$ to address limitations associated with each data type. One such framework combines predictions from geostatistical models trained separately from AOD or CTM using Bayesian ensemble averaging [16]. This approach utilizes all available monitoring, AOD, and CTM data, and provides fine-scale estimates when AOD is available, and gap-filled estimates from CTM when AOD is missing. Probabilistic uncertainty quantification for each prediction is also available through prediction interval and standard deviation. This framework has been adopted by the Multi-Angle Imager for Aerosols (MAIA) project to estimate daily ${\text{PM}}_{2.5}$, ${\text{PM}}_{10}$, and speciated ${\text{PM}}_{2.5}$ major components in multiple large population centers around the world [17].

Here we introduce an R package, ensembleDownscaleR, containing a suite of functions to facilitate the adaptation of this modeling approach to other settings and air quality modeling applications. We describe the statistical method from [16], as well as model extensions implemented for MAIA. We also detail the functionality of the R package and provide an example analysis using 2018 data from the Los Angeles metropolitan area to estimate daily ${\text{PM}}_{2.5}$ at 1km spatial resolution. This work is structured as a tutorial designed to guide practitioners in the use of the ensembleDownscaleR R package. Each stage includes a methodological introduction, code examples applied to Los Angeles dataset, and corresponding results. We make the data available at Zenodo and provide all code used for model fitting, data processing, and plot/table creation at GitHub, to ensure full reproducibility of results and assist users in their own analyses. The ensembleDownscaleR R package is available for download at GitHub.

2. Case Study Data

Our analyses correspond to a spatial region overlapping the Los Angeles metropolitan area ranging from -120.50 to -115.75 longitude, and from 32.00 to 35.70 latitude, and a temporal range of January 1st, 2018, to December 31st, 2018 (Figure 1). This region contained 60 monitoring stations from the Air Quality System (AQS) that provide hourly ${\text{PM}}_{2.5}$ measurements, which we aggregate to daily averages.

AOD data were obtained from the Multi-Angle Implementation of Atmospheric Correction (MAIAC) with a spatial resolution of 1 km. Chemical transport model (CTM) data were obtained from the Community Multiscale Air Quality (CMAQ) model, which provides daily PM_2.5 simulations at a 12 km spatial resolution [18]. We spatially linked the AQS monitor, AOD and CTM data using the MAIAC 1 km grid for modeling fitting and prediction.

Additional spatial land use and spatio-temporal meteorological covariates are matched to both monitor locations and grid cells. These include elevation (elevation), population density (population), cloud cover (cloud), east-west wind velocity (u_wind), north-south wind velocity (v_wind), height of planetary boundary layer (hpbl), shortwave radiation (short_rf), and humidity at two meters above the ground (humidity_2m). Additional details on the source of each covariate and spatial resolutions are detailed in the supporting information.

3. Methods and Results

In this section we provide a case study for all stages of the Bayesian ensemble fitting and prediction process on the case study. For each stage we detail the methods used, provide code examples to illustrate the use of relevant ensembleDownscaleR package functions, and present corresponding results. We break up this complete workflow for producing Bayesian ensemble $P{M}_{2.5}$ predictions into six stages, detailed as follows.

• Fit two separate Bayesian downscaler regression models, $M_{CTM}$ and $M_{AOD}$, on the monitoring PM_2.5 data, one spatially and temporally matched with CTM data and the other matched with AOD data, following the form $p{m}_{2.5}(s,t)=\alpha (s,t)+\beta (s,t)X(s,t)+\mathbf{l}{\left(s\right)}^{\prime }\gamma +\mathbf{m}{(s,t)}^{\prime }\delta +ϵ(s,t)$, where X is either AOD or CTM, and $\mathbf{l}$ and $\mathbf{m}$ are additional spatial and spatio-temporal covariates respectively.
• Produce estimates of PM_2.5 (posterior predictive means) and variances for all locations using $M_{CTM}$ from stage 1. Produce estimates of $P{M}_{2.5}$ means and variances for all times and locations for which AOD is available, using $M_{AOD}$.
• Use cross-validation to produce two sets of out-of-sample PM_2.5 prediction means and variances using the same data and model form as in stage 1. This produces two datasets of out-of-sample prediction means and variances for each monitor observation.
• Estimate spatially varying weights from the out-of-sample prediction mean and variances from stage 2 and the monitor PM_2.5 measurements.
• Use Gaussian Process spatial interpolation (krigging) to predict weights for all grid cells in the study area.
• Use the fitted models from stage 1 and the weight estimates from stage 4 to acquire ensemble predictions of $P{M}_{2.5}$ at each grid cell in the study area.

The Los Angeles dataset included 60 monitoring stations, 15,821 observation-CTM pairs, 11,668 observation-AOD pairs and 122,735 prediction grid cells. The total computation time with 25,000 Markov chain Monte Carlo (MCMC) iterations per model fit (for two stage 1 model fits, 20 stage 3 model fits, and one stage 4 model fit), took approximately 11.29 hours on an Apple MacBook Pro laptop with an Apple M3 Max processor and 36 GB of RAM.

3.1. Stage 1: Downscaler Regression Model

This section details the model specifications available for Bayesian downscaler regression model fitting using the grm() function.

3.1.1. Model

The Bayesian downscaler regression model is formulated as a spatial-temporal regression of PM_2.5 against X, which is either AOD or CTM depending on user input.

The statistical model is as follows:

$$\begin{array}{c}\hfill p{m}_{2.5}(s,t)=\alpha (s,t)+\beta (s,t)X(s,t)+\mathbf{l}{\left(s\right)}^{\prime }\gamma +\mathbf{m}{(s,t)}^{\prime }\delta +ϵ(s,t)\end{array}$$

(1)

where $\alpha (s,t)$ and $\beta (s,t)$ are the spatial-temporal intercept and AOD/CTM slope of the regression model at location s and time t, $\gamma$ and $\delta$ are fixed effects for spatial and spatio-temporal covariates $\mathbf{l}\left(s\right)$ and $\mathbf{m}(s,t)$ respectively, and $ϵ(s,t)\sim N(0,{\sigma }^2)$. Here ${\sigma }^2$ is modeled with an inverse Gamma prior distribution, ${\sigma }^2\sim IG(a_{{\sigma }^2},b_{{\sigma }^2})$, where $a_{{\sigma }^2}$ and $b_{{\sigma }^2}$ hyperparameters are specified with the sigma.a and sigma.b arguments in the grm() function. For these and the remainder of the hyperparameter arguments, defaults are set to represent uninformative priors.

The slope and intercept parameters are composed of the following additive spatial and temporal random effects and fixed effects:

$$\begin{array}{cc}\hfill \alpha (s,t)& ={\alpha }_0+\alpha \left(s\right)+\alpha \left(t\right)\hfill \\ \hfill \beta (s,t)& ={\beta }_0+\beta \left(s\right)+\beta \left(t\right)\hfill \end{array}$$

(2)

where spatial random effects $\alpha \left(\mathbf{s}\right),\beta \left(\mathbf{s}\right)\sim GP(\mathbf{0},{\tau }^2\mathbf{K}(\theta ,D))$ follow a Gaussian Process ($GP$) depending on a user-specified kernel $\mathbf{K}$ with range parameter $\theta$ and distance matrix D. Temporal random effects $\alpha \left(\mathbf{t}\right),\beta \left(\mathbf{t}\right)$ that are set as first-order random walk. Normal priors are applied to fixed effects $\gamma$ and $\delta$ with equivalence to a ridge regression shrinkage prior.

Users can specify inclusion of any combination of additive or multiplicative spatial or temporal random effects, and can input $\mathbf{L}$ and $\mathbf{M}$ matrices for fixed effects. For example, if the user-specifies inclusion of an additive spatial effect, a multiplicative temporal effect, and no $\mathbf{L}$ and $\mathbf{M}$ matrices, the intercept/slope equations would simplify as follows:

$$\begin{array}{cc}\hfill \alpha (s,t)& ={\alpha }_0+\alpha \left(s\right)\hfill \\ \hfill \beta (s,t)& ={\beta }_0+\beta \left(t\right)\hfill \end{array}$$

(3)

The inclusion of additive or multiplicative temporal and spatial effects is specified by the user with the include.additive.temporal.effect,

include.multiplicative.temporal.effect, include.additive.spatial.effect, and include.multiplicative.spatial.effect arguments in the grm() function.

3.1.2. Spatial Random Effects

The spatial random effects are modeled using Gaussian Processes (GP), with the covariance kernel $\mathbf{K}$ specified by the user. We provide four covariance kernels for the spatial random effects: exponential, and Matérn for $\nu \in {\displaystyle \frac{1}{2}},{\displaystyle \frac{3}{2}},{\displaystyle \frac{5}{2}}$, where:

$$\begin{array}{cc}\hfill {\mathbf{K}}_{exp}=& exp\left(-{\displaystyle \frac{d}{\theta }}\right)\hfill \\ \hfill {\mathbf{K}}_{mat}=& {\displaystyle \frac{1}{\Gamma \left(\nu \right)2^{\nu -1}}}{\left({\displaystyle \frac{\sqrt{2\nu }}{\theta d}}\right)}^{\nu }{K}_{\nu }\left({\displaystyle \frac{\sqrt{2\nu }}{\theta }}d\right)\hfill \end{array}$$

(4)

We also allow for the user to specify a custom covariance kernel, which must be positive definite and symmetric and parametrized with respect to the range parameter $\theta$ and distance d such that ${\mathbf{K}}_{user}=f(\theta ,d)$.

While non-additive spatio-temporal effects are not currently implemented, we provide the option to use different sets of spatial effects for different time periods. By using different spatial effects for say seasons or months, some temporal variation in spatial effects can be accounted for. Regardless of number of spatial effect sets, Gaussian process parameters (${\tau }^2$, $\theta$) are shared across sets. For example, if spatial-time sets are specified for seasons, the additive spatial random effect is as follows:

$${[{\alpha }_n\left({\mathbf{s}}_{spr}\right),{\alpha }_n\left({\mathbf{s}}_{sum}\right),{\alpha }_n\left({\mathbf{s}}_{fal}\right),{\alpha }_n\left({\mathbf{s}}_{win}\right)]}^{\prime }\sim GP({\mathbf{1}}_4\otimes {\mathbf{0}}_n,{\tau }_{\alpha }^2{\mathbf{I}}_{4,4}\otimes \mathbf{K}{({\theta }_{\alpha },d)}_{n,n})$$

(5)

where n is the number of spatial locations and ⊗ is the Kronecker product.

Priors are placed on the Gaussian process parameters ${\tau }_{\alpha }^2$ and ${\theta }_{\alpha }$ such that ${\tau }_{\alpha }^2\sim IG(a_{{\tau }_{\alpha }^2},b_{{\tau }_{\alpha }^2})$ and $log\left({\theta }_{\alpha }\right)\sim Gamma(a_{{\theta }_{\alpha }},b_{{\theta }_{\alpha }})$, where $a_{{\tau }_{\alpha }^2}$, $b_{{\tau }_{\alpha }^2}$, $a_{{\theta }_{\alpha }}$, and $b_{{\theta }_{\alpha }}$ are specified by the user using the tau.alpha.a, tau.alpha.b, theta.alpha.a, and theta.alpha.b arguments in the grm() function. While most parameters in the Bayesian downscaler regression model are sampled with Gibbs updates, the ${\theta }_{\alpha }$ parameter is sampled with a Metropolis-Hastings step. We employ a log-normal proposal distribution with a user-specified tuning parameter ${\theta }_{\alpha }^{*}$, such that ${\theta }_{\alpha }^{(t+1)}\sim \mathrm{Log}-\mathrm{Normal}({\theta }_{\alpha }^{\left(t\right)},{\theta }_{\alpha }^{*})$. The user can specify the tuning parameter ${\theta }_{\alpha }^{*}$ using the theta.alpha.tune argument in the grm() function, as well as the initial value for ${\theta }_{\alpha }$ using the theta.alpha.init argument. Prior specifications and hyperparameter arguments are similar for ${\tau }_{\beta }^2$ and ${\theta }_{\beta }$.

3.1.3. Temporal Random Effects

The first-order random walk temporal random effects ($\alpha \left(\mathbf{t}\right),\beta \left(\mathbf{t}\right)$) are specified such that,

$$\begin{array}{cc}\hfill E\left[\alpha \right(t\left)\right]& =\left\{\begin{array}{cc}{\rho }_{\alpha }\alpha (t+1)\hfill & \mathrm{if}t=1\hfill \\ {\rho }_{\alpha }{\displaystyle \frac{\alpha (t-1)+\alpha (t+1)}{2}}\hfill & \mathrm{if}1<t<T\hfill \\ {\rho }_{\alpha }\alpha (t-1)\hfill & \mathrm{if}t=T\hfill \end{array}\right.\hfill \\ \hfill Var\left[\alpha \right(t\left)\right]& =\left\{\begin{array}{cc}{\omega }_{\alpha }^2\hfill & \mathrm{if}t=1\hfill \\ {\displaystyle \frac{{\omega }_{\alpha }^2}{2}}\hfill & \mathrm{if}1<t<T\hfill \\ {\omega }_{\alpha }^2\hfill & \mathrm{if}t=T\hfill \end{array}\right.\hfill \end{array}$$

(6)

with $\beta \left(\mathbf{t}\right)$ similarly specified. To reduce computation burden, each $\rho$ is discretized as $2,000$ evenly spaced values between 0 and 1, and each $\omega$ determines the temporal smoothness level. Initial values for ${\rho }_{\alpha }$ and ${\rho }_{\beta }$ can be specified by the user using the rho.alpha.init and rho.beta.init arguments in the grm() function. An inverse gamma prior is place on ${\omega }_{\alpha }^2$ and ${\omega }_{\beta }^2$ such that ${\omega }_{\alpha }^2\sim IG(a_{{\omega }_{\alpha }^2},b_{{\omega }_{\alpha }^2})$ and ${\omega }_{\beta }^2\sim IG(a_{{\omega }_{\beta }^2},b_{{\omega }_{\beta }^2})$. The user can specify the hyperparameters $a_{{\omega }_{\alpha }^2}$, $b_{{\omega }_{\alpha }^2}$, $a_{{\omega }_{\beta }^2}$, and $b_{{\omega }_{\beta }^2}$ using the omega.alpha.a, omega.alpha.b, omega.beta.a, and omega.beta.b arguments in the grm() function.

3.1.4. Fixed Effects

The user is able to specify inclusion of fixed effects $\gamma$ and $\delta$ for spatial and spatio-temporal covariates, $\mathbf{L}$ and $\mathbf{M}$, respectively. The fixed effects are modeled with normal priors, $\gamma \sim N(0,\left({\lambda }_{\gamma }\mathbf{I}\right))$ and $\delta \sim N(0,\left({\lambda }_{\delta }\mathbf{I}\right))$ where inverse gamma priors are placed on the ${\lambda }_{\gamma }$ and ${\lambda }_{\delta }$ parameters such that ${\lambda }_{\gamma }\sim IG(a_{{\sigma }^2},b_{{\sigma }^2})$ and ${\lambda }_{\delta }\sim IG(a_{{\sigma }^2},b_{{\sigma }^2})$. Thus the $\lambda$ parameters share the same hyperparameter settings as those for ${\sigma }^2$, specified by the user using the $\mathtt{sigma}.\mathtt{a}$ and $\mathtt{sigma}.\mathtt{b}$ arguments in the grm() function.

3.1.5. Stage 1 Code Example

We load the ensembleDownscaleR package and fit the Bayesian downscaler regression models for CTM using the previously described Los Angeles $P{M}_{2.5}$ dataset. The code for fitting the AOD model is omitted for brevity, but is similar to the CTM model fit and is included in the tutorial code at GitHub.

3.2. Stage 2: Produce PM_2.5 Estimates and Predictions with Available CTM and AOD Data

CTM data is available for all times and locations in the study area while AOD data availability depends on time period. In this stage we use the grm_pred() to produce posterior predictive means and variances for all CTM and AOD data variables. For example, we first input the fitted $M_{CTM}$ model and CTM data for all times and locations in the study area, to produce ${\mu }_{st}^{CTM}$ and ${\sigma }_{st}^{2,CTM}$ for all locations s and times t. We then input the fitted $M_{AOD}$ model and the sparser AOD data, to produce ${\mu }_{st}^{AOD}$ and ${\sigma }_{st}^{2,AOD}$ for all times and locations for which AOD data is available. We note that grm_pred() outputs NA values for prediction data with locations identical to monitor locations because in this case the observed concentrations can be used.

3.2.1. Stage 2 code example

Using the fitted downscaler regression models from stage 1, we now produce full PM_2.5 predictions for all locations and times in the study area using the CTM-based fitted model, and for all locations and times for which AOD is available using the AOD-based fitted model (Figure 2). The code for producing AOD predictions is omitted here, but is included in the full tutorial code at GitHub.

3.3. Stage 3: Use Cross-Validation to Produce out-of-Sample Prediction Means and Variances

3.3.1. Cross-Validation Details

K-fold cross-validation prevents overfitting by separating the data set into k number of folds, iteratively fitting the model to $k-1$ folds and predicting the remaining fold. We provide two functions to perform k-fold cross-validation with the geostatistical regression model. The first is the function create_cv() which creates cross-validation indices according to a user specified sampling scheme. The second function, grm_cv(), returns the out-of-sample PM_2.5 predictions, calculated according to user inputted cross-validation indices (either obtained from the create_cv() function or created by the user), and arguments similar to those used for the grm() function to specify the downscaler regression model. The out-of-sample predictions are stacked into a dataset of the same length and order as the original dataset on which the cross-validation is applied.

The create_cv() allows specification of the following types of cross-validation (Figure 3):

• Ordinary: Folds are randomly assigned across all observations
• Spatial: Folds are randomly assigned across all spatial locations.
• Spatial Clustered: K spatial clusters are estimated using k-means clustering on spatial locations. These clusters determine the folds.
• Spatial Buffered: Folds are randomly assigned across all spatial locations. For each fold, observations are dropped from the training set if they are within a user-specified distance from the nearest test set point.

In Figure 3 we visually detail how the folds are assigned for each type of cross-validation, using the monitor locations and times in our study area as an example and assuming 5 cross-validation folds. We plot all fold assignments for four randomly chosen days, for ordinary, spatial, spatial clustered cross-validation, with color representing fold assignment. For the spatial buffered cross-validation we plot only one fold assignment, with color representing if a location is in the first test fold, the first training fold, or dropped due to being within a 30km buffer of a location in the first test fold.

When assigning folds we enforce that each fold contains at least one observation from each spatial location and spacetime indicator. Locations for which there are fewer observations than folds should be filtered out of the dataset prior to analysis. Data from the first and last time point are left unassigned, and out-of-sample predictions for these data are output as missing data. If the out-of-sample dataset has a larger temporal range than the in-sample dataset for a given fold, the out-of-sample predictions for the extra time points are also output as missing data.

3.3.2. Producing out-of-Sample Prediction Means and Variance

The grm_cv() function uses the previously detailed cross validation indices and model specifications to produce estimates of $f_{CTM}\left(y_{st}\right|M_{CTM})$ and $f_{CTM}\left(y_{st}\right|M_{AOD})$, where $y_{st}$ is the is $P{M}_{2.5}$ value at location s, time t, and $f_{CTM}$ and $f_{AOD}$ are the posterior predictive distributions based on models $M_{CTM}$ and $M_{AOC}$ respectively, the downscaler regression models. Specifically, grm_cv() outputs posterior predictive means ${\mu }_{st}$ and variances ${\sigma }_{st}^2$ that are used in stage 4 to fit the full ensemble model.

3.3.3. Stage 3 Code Example

We now create the cross-validations indices with the create_cv() function for both AOD and CTM linked monitors, and then use the grm_cv() function to produce out-of-sample PM_2.5 predictions for all monitor observations, for both the CTM and AOD data. The code for producing AOD out-of-sample predictions is omitted here, but is included in the tutorial code at GitHub.

3.4. Stage 4: Estimate Spatially Varying Weights

At this stage we use the ensemble_spatial() function to fit the ensemble model $p\left(y_{st}\right|M_{CTM},M_{AOD})=w_s{f}_{CTM}\left(y_{st}\right|M_{CTM})+(1-w_s)f_{AOD}\left(y_{st}\right|M_{AOD})$, where $w_s$ are spatially varying weights. We estimate the weights $w_s$ by fitting the ensemble model $p\left(y_{st}\right|M_{CTM},M_{AOD})$ on the out-of-sample predictions produced during stage 3,

$\widehat{f_{CTM}\left(y_{st}\right|M_{CTM})}$ and $\widehat{f_{CTM}\left(y_{st}\right|M_{CTM})}$, and the original ${\text{PM}}_{2.5}$ data at all times and a monitor locations. We place a Gaussian Process prior on the weights, $q_s=logi{t}^{-1}\left(w_s\right)\sim GP(0,{\tau }_w^2\mathbf{K}({\theta }_w,d))$, where $\mathbf{K}$ is an exponential kernel, and ${\tau }_w^2$ and ${\theta }_w$ are the distance and range parameters, respectively. Similar to the spatial processes in stage 1, we place an inverse gamma prior on ${\tau }_w^2$ and a gamma prior on ${\theta }_w$, such that ${\tau }_w^2\sim IG(a_{{\tau }_w^2},b_{{\tau }_w^2})$ and $log\left({\theta }_w\right)\sim Gamma(a_{{\theta }_w},b_{{\theta }_w})$. The user can specify the hyperparameters $a_{{\tau }_w^2}$, $b_{{\tau }_w^2}$, $a_{{\theta }_w}$, and $b_{{\theta }_w}$ using the tau.a, tau.b, theta.a, and theta.b arguments in the ensemble_spatial() function.

The ensemble_spatial() function accepts the output from grm_cv() in stage 3 as input, and outputs the full posterior distribution samples of $q_s$ where $q_s=logi{t}^{-1}\left(w_s\right)$.

3.4.1. Stage 4 Code Example

At this stage we estimate the spatially varying weights $w_s$ for the ensemble model, using the out-of-sample predictions produced in stage 3, and the original monitor PM_2.5 measurements.

We display the out-of-sample predictions produced in stage 3 to these ensemble weight estimates and the original monitor $P{M}_{2.5}$ measurements, while including the weight estimates for reference (Figure 4). If desired, one can additionally use these weight estimates with the gap_fill() function to produce ensemble-based predictions for each location at which $P{M}_{2.5}$ is observed, though these outputs are not used in later stages. We do this here to compare the ensemble model performance to the CTM-based and AOD-based models on the times and locations for which both CTM and AOD data are observed, employing ordinary, spatial, spatial-clustered, and spatial-buffered cross-validation for each model (Table 1, see supplementary materials for more details). For all cross-validation formulations, the ensemble model outperforms both the CTM-based and AOD-based models in terms of RMSE and $R^2$, while maintaining accurate 95% prediction interval coverage. Note that gap_fill() function uses the weight estimates to produce ensemble-based predictions for times and locations at which both CTM and AOD data are observed, and fills in the remaining times and locations with the CTM-based predictions.

3.5. Stage 5: Predict Weights for All Locations

At this stage we use the weight_pred() function to spatially interpolate the posterior samples of $w_s$ garnered in stage 4, across all 1km x 1km grid cells in the study area. Specifically $w_{predictions}=logi{t}^{-1}\left(q_{predictions}\right)$ where $q_{predictions}=\mathrm{krig}\left(q_s\right)$, with $q_s$ output from stage 4. These weights are used in the final stage to produce ensemble-based $P{M}_{2.5}$ predictions for all locations in the study area.

3.5.1. Stage 5 Code Example

We spatially interpolate the posterior samples of $w_s$ from stage 4 across all locations in the study area using the weight_pred() function. (Figure 5). This provides the weights used in the final stage to produce full ensemble estimates for all times and locations in the study area.

3.6. Stage 6: Compute Ensemble Predictions for All Locations

The last stage comprises using the posterior means and variances for all CTM and AOD data produced in stage 2, and the spatially interpolated weights from stage 5, to compute ${\text{PM}}_{2.5}$ posterior predictive means ${\widehat{y}}_{st}$ and variances ${\widehat{\sigma }}_{y_{st}}^2$ for all times t and locations s in the study area, done with the gap_fill() function.

For times and locations for which both CTM and AOD are observed, gap_fill() outputs ensemble-based estimates, where ${\widehat{y}}_{st}=w_s{\mu }_{st}^{CTM}+(1-w_s){\mu }_{st}^{AOD}$ and ${\widehat{\sigma }}_{y_{st}}^2=w_s^2{\sigma }_{st}^{2,CTM}+{(1-w_s)}^2{\sigma }_{st}^{2,AOD}$. For times and locations for which solely CTM is available, gap_fill() outputs posterior predictive means and variances identical to those produced in stage 2 from $M_{CTM}$.

3.6.1. Stage 6 Code Example

Here we input the posterior means and variances for all CTM and AOD data produced in stage 2, and the spatially interpolated weights from stage 5, into the gap_fill() function, which outputs $P{M}_{2.5}$ posterior predictive means ${\widehat{y}}_{st}$ and variances ${\widehat{\sigma }}_{y_{st}}^2$ for all times t and locations s in the study area (Figure 5).

4. Discussion

In this work we introduce the ensembleDownscaleR package for fitting Bayesian geostatistical regression and ensemble models, designed for predicting PM_2.5 using CTM simulations and AOD measurements. We also provide a code tutorial based on a case study of Los Angeles metropolitan area data from 2018. The purpose of this work is guide practitioners in generating robust PM_2.5 predictions and uncertainty quantification with a flexible and well-documented software workflow. The framework can also be applied to other air pollutants and data integration problems.

There are areas for future improvement that are worth noting. The Gaussian Process spatial random effects employed by our model are appropriate for the size of the Los Angeles case study data used here. For much larger datasets (with many more monitors and/or much larger prediction grids) the inference of the Gaussian Process parameters will be prohibitively slow. Incorporation of scalable random processes, such as Nearest Neighbor Gaussian Processes [19], could make these methods feasible for much larger datasets than those assessed here. Furthermore, spatial covariance specifications are currently limited to isotropic, stationary kernels. There are cases where PM_2.5 data may exhibit correlation that suggests anisotropic or nonstationary kernels would be more appropriate, such as data that includes periods with high wind or localized wildfires. Inclusion of covariates such as wind information can often resolve this. Inspecting residual correlation and covariance parameter can ensure that the covariance is reasonably specified for a given dataset. Finally, the current software does not support integrated parallelization. For example, the grm_cv function fits a model on each cross-validation fold sequentially rather than exploiting multiple cores or compute nodes to fit each model concurrently. Model fitting time could be substantially lowered without sacrificing software ease-of-use by incorporating parallelization specifications directly in the ensembleDownscaleR package functions.