Scalable Bayesian Geostatistical Regression Model for Bias-Correcting Large-Scale Daily Satellite-Retrieved Aerosol Optical Depth and Chemical Transport Model Simulations.

Scalable Bayesian Geostatistical Regression Model for Bias-Correcting Large-Scale Daily Satellite-Retrieved Aerosol Optical Depth and Chemical Transport Model Simulations - PMC Skip to main content An official website of the United States government Here's how you know Here's how you know Search Log in Primary site navigation Search Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Atmos Environ (1994) . Author manuscript; available in PMC: 2026 Jun 5. Published in final edited form as: Atmos Environ (1994). 2026 May 19;380:122079. doi: 10.1016/j.atmosenv.2026.122079 Scalable Bayesian Geostatistical Regression Model for Bias-Correcting Large-Scale Daily Satellite-Retrieved Aerosol Optical Depth and Chemical Transport Model Simulations Wyatt G Madden 1, * , Yang Liu 2 , Howard H Chang 1 Author information Article notes Copyright and License information 1 Department of Biostatistics and Bioinformatics, Emory University, Atlanta, USA 2 Department of Environmental Health, Emory University, Atlanta, USA Author contributions Conceptualization, W.G.M. and H.H.C.; methodology, W.G.M. and H.H.C.; software, W.G.M. and H.H.C.; validation, W.G.M.; formal analysis, W.G.M. and H.H.C.; investigation, W.G.M. and H.H.C.; resources, H.H.C. and Y.L.; data curation, W.G.M., H.H.C., and Y.L.; writing—original draft preparation, W.G.M.; writing—review and editing, Y.L. and H.H.C.; visualization, W.G.M.; supervision, H.H.C.; project administration, Y.L. and H.H.C.; funding acquisition, Y.L. and H.H.C. All authors have read and agreed to the published version of the manuscript. * Author to whom any correspondence should be addressed. wyatt.madden@emory.edu Issue date 2026 Sep 1. PMC Copyright notice PMCID: PMC13225854  NIHMSID: NIHMS2178827  PMID: 42232855 The publisher's version of this article is available at Atmos Environ (1994) Abstract Ambient fine particulate matter (PM 2.5 ) is linked to numerous adverse health outcomes. Accurate large-scale PM 2.5 prediction and uncertainty quantification at fine spatial resolution are crucial for health analyses. We introduce a scalable Bayesian spatial-temporal geostatistical regression model (GRM) that performs bias-correction of gridded Chemical Transport Model (CTM) outputs or satellite-retrieved Aerosol Optical Depth (AOD) and additional spatiotemporal covariates to improve PM 2.5 predictions. This model employs nearest neighbor Gaussian process (NNGP) spatial random effects to efficiently exploit spatial correlation when analyzing large numbers of PM 2.5 monitors (7.64 times faster than a regular Gaussian process), and when predicting over large spatial grids (43.7 times faster than a regular Gaussian process). We apply the model to two case studies: an analysis using 12-km CTM outputs over the contiguous United States, and an analysis using 1-km AOD data over the state of California. For the U.S. 12-km case study, we demonstrate the spatial predictive performance and accurate uncertainty quantification of the model with several cross-validation experiments. For the California 1-km case study, in addition to the above, we also assess Random Forest (RF) and GRM-RF hybrid models which include a RF mean term in the GRM. Under the California 1-km case study the GRM-RF provides higher predictive accuracy (RMSE 7.325 μg/m 3) than both the stand-alone RF (RMSE 7.732 μg/m 3) and GRM (RMSE 8.557 μg/m 3) while maintaining accurate predictive uncertainty (95% credible interval coverage 0.957). Finally, we provide open-source R software for researchers to implement the models. Keywords: geostatistical downscaler, Bayesian, large-scale, AOD, CTM, PM 2.5 , air pollution, spatial 1. Introduction Air pollution is prevalent in both urban and natural environments, and can originate from a variety of sources such as wildfires and the burning of fossil fuels. Ambient fine particulate matter less than 2.5 μ m in aerodynamic diameter (PM 2.5 ) is a widely measured air pollutant due to its association with adverse health outcomes, including respiratory and cardiovascular disease, low birth weight, and premature mortality ( Apte et al., 2015 ; Fu et al., 2019 ; Yuan et al., 2019 ; Alexeeff et al., 2021 ; Fan et al., 2016 ; Gong et al., 2022 ). The availability of PM 2.5 concentration estimates over large geographical areas and at a fine spatial scale is important for health research to support population-level health effects studies and impact assessments. Furthermore, having complete spatial-temporal coverage of PM 2.5 concentration estimates enables investigation of how concentrations vary due to emission sources and meteorological conditions. The U.S. Environmental Protection Agency (EPA) maintains a network of monitoring stations that report daily PM 2.5 measurements, available through the Air Quality System (AQS) repository ( U.S. Environmental Protection Agency, 2025 ). However, these monitors are sparsely and preferentially located in urban locations, resulting in exposure measurement error and missing data when linked to the health outcomes of interest. To generate high-quality PM 2.5 predictions at locations without monitors from these sparse PM 2.5 observations, one approach is to bias-correct a proxy variable selected for its wider spatial availability and its high correlation with PM 2.5 . One proxy variable is PM 2.5 simulation from Chemical Transport Models (CTM) such as the EPA-developed Community Multiscale Air Quality (CMAQ) model ( US EPA Office of Research and Development, 2024 ), which uses physical governing equations and meteorological variables to simulate PM 2.5 over large grids. Another proxy variable is satellite-based Aerosol Optical Depth (AOD), a measure of the extinction of light through the atmospheric column which reflects aerosol concentration. Satellite data products have fine spatial resolution and global coverage. One approach to bias-correct gridded data is to employ a downscaler model, where PM 2.5 monitor observations are regressed on the spatially and temporally matched CTM or AOD grid cells. The fitted model is then used to predict PM 2.5 wherever the proxy data is available. A variety of regression modeling techniques have been utilized for downscaling, including both statistical and machine learning approaches ( Berrocal, Guan, et al., 2020 ). Geostatistical regression models that allow for spatially-varying additive and multiplicative biases are particularly attractive because they can capture the spatial dependence in pollution concentrations to improve spatial interpolation. Geostatistical models support uncertainty quantification under either frequentist or Bayesian inference ( Diggle and Ribeiro Jr, 2002 ). We employ Bayesian inference in our geostatistical downscaling models because it provides a coherent probabilistic framework for characterizing predictive uncertainty and propagating bias-correction error through downstream epidemiological analyses ( Berrocal, Gelfand, and Holland, 2010b ; Berrocal, Gelfand, and Holland, 2010a ; Chang, Hu, and Liu, 2014 ). Specifically, Bayesian approaches provide uncertainties for parameters that govern spatial random effects which are typically ignored in frequentist approaches. However, these models can become computationally burdensome when fitting large numbers of monitors or predicting over many fine-scale grid points. This is primarily due to reliance on Gaussian processes (GPs) for modeling spatial random effects, for which parameter estimation requires covariance inversions that scale cubically with the number of monitors ( Williams and Rasmussen, 2006 ). Machine learning approaches such as Random Forests (RFs) and Neural Networks are also effective for downscaling due to their high prediction accuracy and scalability to large datasets ( Berrocal, Guan, et al., 2020 ; Di et al., 2016 ; Hu et al., 2017 ). However, these do not provide straightforward and principled uncertainty quantification, and do not explicitly model spatial-temporal dependence in the outcome. Here we introduce a Bayesian spatial-temporal downscaler model ( Chang, Hu, and Liu, 2014 ) and the associated R package ( Madden et al., 2025 ) with two extensions that enable the model to scale to large numbers of monitors and large prediction areas, addressing key computational barriers to generating PM 2.5 estimates and uncertainty quantification for large-scale health studies. These extensions include nearest neighbor Gaussian processes (NNGPs) ( Datta et al., 2016a ) that allow covariance matrix inversions to scale linearly with respect to the number of monitors, and spatial range parameter discretization, which allows for computation of covariance matrix inverses prior to running the MCMC fitting algorithm rather than at each iteration. We demonstrate model utility and performance on two case studies for predicting daily PM 2.5 concentrations in 2018, a year with high wildfire activities: a contiguous U.S. analysis bias-correcting 12-km CMAQ simulations, demonstrating scalability for model fitting with a large number of monitors, and a California analysis bias-correcting 1-km AOD retrievals, demonstrating scalability for prediction over large fine-scale grids. As part of the California case study, we also show that GRM-RF hybrid models that include RF mean terms in the GRMs result in improved performance when compared to standalone GRM or RF models, while maintaining accurate uncertainty quantification. Finally, we detail new functionality included within the ensembleDownscaleR R package that enables researchers to employ these models on their own analyses. 2. Study Area and Data This work includes two case studies; the first is an analysis of PM 2.5 data collected over the contiguous U.S. for bias-correcting 12-km CMAQ simulations, and the second is an analysis of data collected over the state of California for bias-correcting 1-km AOD. Both case studies include daily average PM 2.5 measurements collected from AQS monitors maintained by the U.S. EPA; the California study area encompasses 158 monitors, and the U.S. study area encompasses 973 monitors, though only 861 are included in the analyses after removing monitors with less than two days of readings within each month ( figure 1 ). Monitor removal is conducted to ensure spatial and spatial-temporal random effect estimability during 10-fold cross-validation. The U.S. case study illustrates the GRM’s ability to accommodate a large number of monitors and to incorporate multiple sets of spatial random effects for different temporal periods, while the California case study demonstrates its feasibility for fine-scale prediction with multiple covariates. We include two separate case studies to demonstrate the computational advantages of using NNGP spatial random effects in GRM models, both for model fitting with large numbers of monitor locations (U.S. case study) and for prediction over large fine-scale grids (California case study). Figure 1. Open in a new tab The PM 2.5 monitors included in the (a) the California case study, and (b) the full U.S. case study. The blue points mark the PM 2.5 monitor locations in each case study. Red and black boundary lines denote the California and U.S. study areas, respectively. 2.1. Case Study 1 Data: Contiguous U.S. For the first case study, we include daily 12-km PM 2.5 simulations from the Community Multiscale Air Quality (CMAQ) model; details on the model setup and evaluation are provided in Appel et al., 2021 . There are a total of 55,508 12-km CMAQ grid cells in the contiguous U.S.. A total of 176,050 AQS PM 2.5 observations are linked to 12-km simulation grid cells for the model fitting and all CTM simulation grid cells for each day in 2018 are used for model predictions. 2.2. Case Study 2 Data: California For the second case study we include daily 1-km aerosol optical depth (AOD) retrievals obtained from the Multi-Angle Implementation of Atmospheric Correction (MAIAC) ( Lyapustin and Wang, 2018 ). There are 480,411 1-km MAIAC grid cells in the California study area. A total of 42,736 AQS PM 2.5 observations are linked to 1-km AOD retrievals for the model fitting algorithm, and all AOD retrieval grid cells for each day in July 2018 are used for model predictions. Missing AOD values were imputed as described previously ( Zhang et al., 2023 ). In addition, we match spatial land-use and spatio-temporal meteorological variables to the nearest AQS monitor and 1-km grid cell. The spatial covariates consist of population density (from Oak Ridge National Laboratory LandScan USA 2018 dataset ( Rose et al., 2019 )), elevation from the ASTER Global Digital Elevation Model ( NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team, 2019 ), while spatiotemporal covariates include cloud cover from NASA MODIS atmosphere L2 cloud product ( Platnick et al., 2015 ) and shortwave radiation, humidity at two meters above the ground, north-south and east-west wind velocities, and height of planetary boundary layer, all from the NASA NLDAS-2 Forcing Dataset ( NASA Land Data Assimilation System Team, 2025 ).Further information on covariate sources and specifications is available in the supplementary materials . 3. Methods Here we provide a description of the geostatistical regression model (GRM) employed in both case studies, with a focus on the additions made that enable scalability to large spatial regions. First we briefly describe the form and structure of the statistical model. A more complete model description and tutorial for the accompanying R package ensembleDownscaleR can be found in Madden et al . ( Madden et al., 2025 ). Next we describe the NNGP prior placed on the model’s spatial random effects, and the discretization of the spatial range parameter, which both enable scalability to the large spatial data sets used in our case studies. Finally we detail the GRM specifications and benchmark models used in each of the case studies. 3.1. Geostatistical Regression Model The Bayesian GRM regresses PM 2.5 on a primary covariate X ( s, t ) (CTM or AOD) and spatial/spatiotemporal covariates: PM 2.5 ( s , t ) = α ( s , t ) + β ( s , t ) X ( s , t ) + l ( s ) ′ γ + m ( s , t ) ′ δ + ∈ ( s , t ) where α ( s, t ) = α 0 + α ( s ) + α ( t ) and β ( s, t ) = β 0 + β ( s ) + β ( t ). Thus the intercept and slope of the primary covariate X ( s, t ) consist of fixed effects α 0 and β 0, random spatial effects α ( s ) and β ( s ), and random temporal effects α ( t ) and β ( t ). Fixed effects γ and δ are included for spatial and spatiotemporal covariates l ( s )′ and m ( s, t )′ respectively. We note that the U.S. case study does not include spatial or spatiotemporal covariates, while the California case study includes both spatial and spatiotemporal covariates, further detailed in section 2.2 . We assume ϵ ( s, t )~ N (0 , σ 2), place NNGP priors on the spatial random effects, and model the temporal random effects as first-order random walks. Models are fit using Markov chain Monte Carlo (MCMC), with the discretized NNGP range parameters sampled via a Metropolis-Hastings ( Hastings, 1970 ) step or a Gibbs update, and the remainder of the parameters sampled via Gibbs updates ( Gelfand and Smith, 1990 ). We run each MCMC algorithm for 25,000 iterations, discarding the first 5,000 as burn-in and thinning by a factor of 4 for each model fit. In our analyses all prior distributions are specified as weakly informative. 3.1.1. Nearest Neighbor Gaussian Process Gaussian processes (GPs) have a long history of being used to model spatial correlation in data ( Cressie, 2015 ; Banerjee, Carlin, and Gelfand 2003 ; Cressie and Wikle, 2011 ). For a given positive definite and symmetric covariance kernel K ( θ, d ), a continuous function of distance d and the range parameter θ , and a variance term τ 2 , a mean-zero spatial GP can be described as Y~ N ( 0 , C ( θ )), where C ( θ ) = τ 2 K ( θ, d ). Estimating the θ parameter and other model parameters typically involves inverting and calculating the determinant of C ( θ ), which costs ( n 3 ) where n is the length of Y . This is computationally infeasible for large n . Nearest neighbor Gaussian processes (NNGPs) provide a computationally tractable alternative to GPs by assuming that a given spatial location is dependent on at most m nearest neighbor locations, allowing for the full multivariate normal distribution to be decomposed into a series of n conditional distributions that depend on at most m neighbors. Evaluation of the NNGP density costs 𝒪( nm 3 ), which is linear in n and scalable to large numbers of spatial locations. Predictions are done by spatially interpolating (kriging) the spatial random effects from the nearest m monitor locations. We further describe how the NNGP is employed within the MCMC fitting algorithm, and NNGP predictions, within the supplementary materials . For the GRMs used in our analyses we employ NNGPs for all spatial random effects (the spatial intercepts and slopes for the primary covariate of interest). We use m = 10 neighbors, which has been shown to nearly match the performance of full GPs in many scenarios ( Datta et al., 2016b ); a sensitivity analysis supporting this choice is included in the supplementary materials . Distances are Euclidean, computed on coordinates projected to kilometer-scale units. 3.1.2. Spatial Range Parameter Discretization To decrease the computation time of the MCMC algorithm in both case studies, we discretize the spatial range parameter of the NNGPs which allows precomputation of kernel inverses. We first run a short GRM MCMC with continuous range parameters (5,000 burn-in iterations), then define an equally spaced grid that covers the empirical support of the posterior range samples. We then precompute and store the NNGP kernel inverses for each discrete range parameter value, which are used in density evaluations during the MCMC routine. The range parameters are inferred using a Metropolis-Hastings (MH) step with a proposal distribution that places equal probability on neighboring discretized values. For example, if θ ∈ {10, 15, 20, 25, 30}, and the current value of θ is 20, then the θ proposal will be either 15 or 25, with equal probability. We also assess the improvement in computation time of this discretization method in relation to non-discretized range parameters, as well as a similar discretization method that relies on a Gibbs update step rather than an MH step (all of which are included as fitting options in the ensembleDownscaleR package). Further details of these computations are included in the supplementary materials . 3.1.3. R Package: ensembleDownscaleR The statistical methods used in this paper are made available through the ensembleDownscaleR R package, previously introduced with a tutorial in Madden et al . ( Madden et al., 2025 ). As part of this work, we extend the ensembleDownscaleR R package with functionality for including NNGP spatial random effects, and discretizing the spatial range parameter terms. The grm() and grm_cv() functions are updated with parameters nngp, discrete.theta.gibbs, discrete.theta.alpha.values , and discrete.theta.beta.values . When nngp is set to TRUE , the GRMs are fit with NNGPs rather than GPs. When vectors of values are included for discrete.theta.alpha.values , and discrete.theta.beta.values , the spatial range parameter discretization methods are employed in the GRM MCMC algorithm using a Gibbs update. The Gibbs update can be swapped with a Metropolis-Hastings update (described in the prior section), by setting discrete.theta.gibbs to FALSE. 3.2. Model Evaluation and Comparison For both case studies, we fit a GRM with NNGP spatial intercepts and slopes for the primary covariate and generate predictions for all grid cells in each study area. For the U.S. case study, a separate set of spatial NNGP random effects (with all sets sharing the same NNGP kernel and range parameter) are used for each month in 2018, which accounts for variation in the spatial correlation present in the relationship between CTM and PM 2.5 across the year. For the California case study, we include fixed effects for the population density, elevation, cloud cover, shortwave radiation, humidity at two meters above the ground, north-south and east-west wind velocities, and height of planetary boundary layer covariates. Model performance is assessed using 10-fold cross-validation (CV) root mean squared error (RMSE) for all combinations of three NNGP kernels (Matérn with hyperparameter ν values of 1 2 , 3 2 , and 5 2 )and five CV schemes (ordinary, spatial, spatial clustered, spatial buffered with buffer size corresponding to 0.7 maximum correlation, and spatial buffered with buffer size corresponding to 0.3 maximum correlation) for which functionality is included in the ensembleDownscaleR R package. Spatial clustered CV uses k-means clustering of PM2.5 monitor locations to define folds to specifically evaluate spatial prediction at regions with no monitors. In spatial buffer CV, we drop monitors from the test-set that have spatial correlation within the threshold to ensure spatial prediction are evaluated at locations further from the training monitors. To determine the spatial buffer sizes, we run the GRM MCMC for 5,000 burn-in iterations, and compute the distance yielding the target correlation under the Matérn ( v = 1 2 ) kernel, given the maximum across the posterior mean spatial ranges. Further explanation of the CV types and spatial buffer calculation are contained within the supplementary materials . In addition, for the California case study we assess spatial 10-fold CV RMSE performance on a series of benchmark models given the larger set of land use and meteorology predictors. These include a RF model ( Breiman, 2001 ) that employs the same covariates (detailed in the supplementary materials ) as the NNGP-based GRM, a RF model that includes these covariates in addition to latitude and longitude, and two NNGP-based GRMs that include the RF models as mean terms (one for the RF fit on the original GRMs and another for the RF fit on these covariates in addition to latitude and longitude). All RF models are fit using the ranger R package ( Wright and Ziegler, 2017 ) with default hyperparameters (500 trees, minimum node size of 5, and 3 covariates randomly sampled at each split). Finally, we also assess differences in computation time between NNGP-based and GP-based GRMs, as well as between non-discretized and discretized spatial range parameter schemes, for model fitting and prediction within each case study. Model fitting is performed on all data included in each respective case study and again employs the aforementioned MCMC algorithm with 25,000 iterations, while model prediction is performed solely on data within the Los Angeles metropolitan area (defined by latitudes 32.00 to 35.70, and longitudes −120.50 to −115.75) due to infeasibility of GP predictions over large fine-scale grids, and uses all 5,000 burned and thinned MCMC samples from the model fit. All model fitting and prediction tasks are performed on a high-performance computing (HPC) cluster, with each task run as a separate job on a single compute node with 8–64 GB of RAM and one CPU core, running R version 4.2.2 ( R Core Team, 2024 ). 4. Results For both case studies, the GRMs achieve high performance and accurate prediction interval coverage (note the standard deviation of all U.S. PM 2.5 observations is 7.15 μg/m 3 and of all California PM 2.5 observations is 14.07 μg/m 3 ; the ideal coverage is 0.95). Matérn kernel with ν = 0.5 performs best (or is tied) with respect to RMSE for all CV types for both case studies, except for spatial extrapolation: spatial buffered (0.3 corr) under the U.S. case study (for which ν = 2.5 performs best) and spatial clustered and spatial buffered (0.7 corr) under the California case study (for which ν = 1.5 performs best) ( table 1 ). For the U.S. study, we generate predictions for all days in 2018 using the GRM with Matérn ν parameter set to 0.5. The mean posterior prediction for each grid cell over seasons demonstrates clear spatial-temporal patterns in PM 2.5 , with notable regions of elevated PM 2.5 across California in the October-December season which encompass multiple severe wildfires ( figure 2 ). This aligns with previous analyses that demonstrate high fire-attributed smoke PM 2.5 in California during this period ( Zhang et al., 2023 ). Table 1. PM 2.5 prediction RMSE (coverage) for each of five CV types and three Matérn ν parameter settings under both U.S. and California study areas. The spatial-buffered CV is formulated with buffer sizes corresponding to 0.3 and 0.7 spatial random effect correlation, and thus varies for each ν parameter setting and study area. The Matérn kernel with ν = 0.5 performs best (or is tied) with respect to RMSE for all CV types for both case studies, except for spatial buffered (0.3 corr) under the U.S. case study (for which ν = 2.5 performs best) and spatial clustered and spatial buffered (0.7 corr) under the California case study (for which ν = 1.5 performs best). Lowest RMSE for each study area and CV type is shown in bold. Cross-Validation Type U.S. Study Area California Study Area v = 0.5 v = 1.5 v = 2.5 v = 0.5 v = 1.5 v = 2.5 Ordinary 4.38 (0.97) 4.38 (0.97) 4.39 (0.97) 8.31 (0.96) 8.31 (0.96) 8.31 (0.96) Spatial 4.71 (0.97) 4.77 (0.97) 4.82 (0.97) 8.56 (0.97) 8.74 (0.97) 8.75 (0.97) Spatial Clustered 5.95 (0.97) 6.12 (0.97) 6.14 (0.97) 11.81 (0.96) 9.24 (0.97) 9.07 (0.97) Spatial Buffered (0.7 corr) 5.33 (0.97) 5.34 (0.97) 5.41 (0.97) 9.11 (0.97) 8.92 (0.97) 8.90 (0.97) Spatial Buffered (0.3 corr) 6.13 (0.97) 5.91 (0.97) 5.88 (0.97) 9.54 (0.97) 9.56 (0.96) 9.59 (0.96) Open in a new tab Figure 2. Open in a new tab PM 2.5 within-season (Jan-Mar, Apr-Jun, Jul-Sep, Oct-Dec) posterior mean predictions for all grid cells in the full U.S. study area, averaged over all days in 2018. For the California case study, all models assessed provide precise uncertainty quantification, with notable differences in prediction performance ( table 2 ). The two RF models (one with latitude and longitude covariates and the other without) have greater prediction accuracy (in terms of RMSE and R 2 ) than the base GRM. However, the GRMs that include each RF model as a mean term instead of separate covariate fixed effects have better prediction accuracy (RMSE/ R 2 ) than each respective RF model. These results suggest that including the RF as a mean term captures potential non-linear and interactive effects between all predictors, while the GRM excels at spatial interpolation for the spatial random effect of the primary covariate due to the explicit modeling of spatial correlation; an analysis of spatial correlation in the residuals supporting this interpretation is included in the supplementary materials . Combining these two advantages results in the most precise predictions. Table 2. California PM 2.5 prediction RMSE, R 2 , and coverage for assessed models. GRM with RF with latitude and longitude mean has the lowest RMSE and highest R 2 , while the GRM with RF mean has the best coverage. All model metrics are assessed on spatial 10-fold CV predictions. Lowest RMSE, highest R 2 , and coverage closest to 0.95, are shown in bold. Model RMSE R 2 Coverage GRM with RF + Lat/Lon Mean 7.325 0.730 0.957 RF + Lat/Lon 7.732 0.699 0.959 GRM with RF Mean 7.786 0.695 0.951 RF 8.229 0.659 0.953 GRM 8.557 0.632 0.968 Open in a new tab PM 2.5 prediction maps are also important to consider when comparing models, as similarly performant models may result in predictions with different levels of smoothness of more extreme values. We plot the California prediction map for the full California study area, and the Los Angeles metropolitan area, for the base GRM, the RF model, and the GRM with RF mean term, for July 15th, 2018 ( figure 3 ). All GRM-based models employ a Matérn kernel with ν parameter set to 0.5. The base GRM results in the smoothest PM 2.5 predictions ( figure 3 , (a) and (b) ). The GRMs with RF mean terms result in slightly smoother prediction grids and less extreme values ( figure 3 , (e) and (f) ) when compared to the RF model prediction grids ( figure 3 , (c) , (d) ). However, overall differences between these prediction grids are minor. Including latitude and longitude in RFs can sometimes introduce grid-aligned artifacts. These were not apparent here, likely because the large number of covariates included already encode the relevant spatial structure. Due to the visual similarity between the predictions from the RF model and the predictions from the RF model with latitude and longitude (and the GRMs that used these models as mean terms), we did not include these formulations here. Figure 3. Open in a new tab July 15, 2018 PM 2.5 predictions for (a) the full California study area and (b) the Los Angeles metropolitan area using the GRM, (c) the full California study area and (d) the Los Angeles metropolitan area using the RF model, and (e) the full California study area and (f) the Los Angeles metropolitan area using the GRM with RF mean term. The predictions from RF and GRM-RF models including latitude and longitude covariates are omitted, as they are visually indistinguishable from those of the respective models without these covariates. Using the NNGP prior yields substantial reductions in computation time for the U.S. case study model fitting (7.64 times faster than a regular Gaussian process when using continuous spatial range parameters) and for predictions under both case studies (26.8 times faster for the U.S. case study and 43.7 times faster for the California case study than a regular Gaussian process, when using continuous spatial range parameters), while matching the GP prior performance for the California case study model fitting ( table 3 ). This indicates that the NNGP prior can be safely adopted across a wide range of monitor counts: there is no meaningful computational penalty when fitting an NNGP-based GRM with a small number of monitors (California case study), and a substantial benefit when fitting models with a large number of monitors (U.S. case study). Furthermore the NNGP-based GRM provides a marked reduction in computation time across varying prediction grid sizes and scales. Spatial effect discretization yields computational benefits for some model fitting and prediction scenarios, though these gains are less pronounced than those achieved through use of the NNGP prior. For GP-based GRMs fit to large numbers of monitors, both the discretized MH and discretized Gibbs methods provide modest improvements in computation time. In contrast, for NNGP-based GRMs, the discretized MH method yields only marginal gains, while the discretized Gibbs method results in inferior performance. Table 3. Computation time (hours) for model fitting and prediction for all combinations of spatial random effect prior (GP or NNGP) and spatial effect fitting algorithm (continuous MH, discretized MH, or discretized Gibbs), for each case study. Fastest time within each case study and task (fit or prediction) shown in bold. Model fitting employs 25,000 MCMC iterations while predictions are computed from the 5,000 burned and thinned MCMC samples provided by each model fit. NNGP-based GRM models yield faster model fitting for scenarios with a large number of monitors (U.S. case study) and faster prediction for both case studies considered. The discretized MH method provides marginal improvements for NNGP-based GRMs, and a larger improvement for GP-based GRMs when fit on the U.S. case study data. Algorithm California Case Study U.S. Case Study GP NNGP GP NNGP Model Fit MH (Continuous) 0.98 1.06 123.27 14.27 MH (Discretized) 0.96 1.07 79.36 11.53 Gibbs (Discretized) 1.00 1.57 82.78 34.61 Prediction MH (Continuous) 74.04 1.66 3.95 0.15 MH (Discretized) 66.00 2.02 3.84 0.18 Gibbs (Discretized) 74.47 2.00 2.94 0.18 Open in a new tab 5. Discussion In this work we introduce a scalable Bayesian spatial-temporal geostatistical model that can incorporate Chemical Transport Model (CTM) outputs or satellite-retrieved Aerosol Optical Depth (AOD) readings and additional spatiotemporal covariates to improve large-scale PM 2.5 predictions. This model includes two extensions for improving scalability (when compared to the model introduced in Chang, Hu, and Liu (2014) ): NNGPs and discretization of the spatial process range parameter. NNGPs are used to model the spatial random effects, which enables key MCMC computations to scale linearly with respect to the number of PM 2.5 monitors rather than cubically, and makes prediction over large fine-scale grids feasible. The spatial range parameter discretization allows precomputation of the kernel inverses used in the NNGP MCMC algorithm and thus reduces overall computation time. Together these extensions enable Bayesian downscaler model fitting and predictions on data that were previously computationally prohibitive due to the number of PM 2.5 monitors, and prediction grid size and resolution. The first includes fitting the GRM on 861 AQS PM 2.5 monitors, and predicting over 55,508 12-km grid cells, in the contiguous U.S. The second includes fitting the GRM on 158 AQS PM 2.5 monitors, and predicting over 480,411 1-km grid cells, in the U.S. state of California. While GPs are feasible for fitting GRMs on these numbers of monitors, they are much slower for the number of monitors included in the U.S. case study, and are prohibitively computationally expensive for predicting over these sizes of prediction grids. The replacement of the GP random effects with NNGP random effects enables the large-scale predictions demonstrated in our case studies, while the range parameter discretization reduces computation for model fitting on large numbers of monitors. We note that numbers of monitors an order of magnitude greater than those here would make GPs infeasible for model fitting, in addition to prediction, which points to the utility of our NNGP-based GRM for future analyses with possibly wider scope. While both case studies demonstrate high predictive performance and accurate coverage of the GRMs, the California case study also demonstrates the utility of employing our more flexible machine learning models as a component of the GRM in a two-stage algorithm. By including the RF models as a mean term within the GRM, prediction is improved compared to each respective RF model. This suggests that the GRM is well-suited to model spatial correlation in the data that is not being fully captured by the Random Forest models, and that combining the GRM with flexible machine learning models like Random Forests is optimal with respect to prediction accuracy while precise uncertainty quantification is preserved. Here we model the mean term with Random Forests, though others may use any flexible machine learning model suited to their data (such as neural networks or ensemble methods) in conjunction with our GRM to improve predictive performance. We also note that since all RF models in this work have hyperparameters set to defaults, hyperparameter optimization has the potential to further improve the performance of both the RF and GRM-RF models. Other software is available for fitting regression models with NNGP random spatial effects. Namely, the spNNGP R package ( Finley, Datta, and Banerjee, 2022 ) provides highly efficient routines for fitting univariate Bayesian spatial regression models. However, it lacks functionality useful to Bayesian downscaling for PM 2.5 prediction. Specifically, there are no straightforward methods for including random temporal effects or multiple sets of kernel-sharing spatial random effects, and currently does not allow spatially varying coefficients. Furthermore, we elect to update the ensembleDownscaleR package with NNGP functionality rather than building a separate R package from scratch, enabling researchers to reference the tutorial contained in Madden et al . Madden et al. (2025) for combining the NNGP-based GRMs introduced in this paper with extensive CV functionality, and to use methods for fitting gap-filling ensemble models that combine GRMs based on both chemical transport model simulations and aerosol optical depth readings. While the ensembleDownscaleR package does not currently allow NNGP priors to be placed on the spatial ensemble weights employed in the namesake ensemble model, we plan to include this functionality in future versions. In this work we demonstrate the utility of Bayesian downscalers for analyses with large study areas that contain PM 2.5 monitors and fine-scale grids. We also provide accompanying software, enabling researchers to generate more precise PM 2.5 estimates and accompanying uncertainty quantification for use in downstream individual-level health studies. There are some future extensions that could further improve these developments. The California case study results showed that incorporating flexible machine learning models within the GRM resulted in higher prediction accuracy. We elected to incorporate Random Forest models as the mean term within the GRM because Random Forests are known to perform well in a wide variety of scenarios, and also because jackknife-approximated ( Tukey, 1958 ) standard errors are available in the ranger R package ( Wright and Ziegler, 2017 ). Ensemble models that incorporate both GRMs and machine learning models have the potential to further improve prediction accuracy while maintaining precise uncertainty quantification over large spatial areas. Another area for model improvement is incorporating anisotropic or nonstationary kernels that may be more suitable for data collected over periods with wildfires or other extreme weather events, such as those used in these analyses. Finally, previous work has shown that generalizing the NNGP spatial process to other Vecchia spatial process formulations, with different ordering and conditioning set schemes such as maximum-minimum distance (MMD) ordering, can result in improved spatial correlation estimation and thus improved model performance ( Katzfuss and Guinness, 2021 ; Katzfuss, Guinness, et al., 2020 ). Supplementary Material 1 NIHMS2178827-supplement-1.pdf (295.3KB, pdf) Supplementary data Supplementary materials: Document containing additional information on data and methods. Highlights. Scalable Bayesian model for daily PM2.5 prediction over large fine-scale grids Improves prediction accuracy while preserving uncertainty quantification for health analyses Employs bias correction of satellite and chemical transport simulations Open-source software supports large-scale applications Funding This work was partially supported by NASA (80NSSC19K0191, PI: Y. Liu) and the NASA MAIA investigation at the Jet Propulsion Laboratory led by D.J. Diner (Subcontract #1588347). It has not been formally reviewed by NASA. Footnotes Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. Declaration of Competing Interest I have nothing to declare. References Alexeeff Stacey E et al. (2021). “Long-term PM2. 5 exposure and risks of ischemic heart disease and stroke events: review and meta-analysis”. In: Journal of the American Heart Association 10.1, e016890. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Appel KW et al. (2021). “The Community Multiscale Air Quality (CMAQ) model versions 5.3 and 5.3.1: system updates and evaluation”. In: Geoscientific Model Development 14.5, pp. 2867–2897. doi: 10.5194/gmd-14-2867-2021. url: https://gmd.copernicus.org/articles/14/2867/2021/ . [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Apte Joshua S et al. (2015). “Addressing global mortality from ambient PM2. 5”. In: Environmental science & technology 49.13, pp. 8057–8066. [ DOI ] [ PubMed ] [ Google Scholar ] Banerjee Sudipto, Carlin Bradley P, and Gelfand Alan E (2003). Hierarchical modeling and analysis for spatial data. Chapman and Hall/CRC. [ Google Scholar ] Berrocal Veronica J, Gelfand Alan E, and Holland David M (2010a). “A spatio-temporal downscaler for output from numerical models”. In: Journal of agricultural, biological, and environmental statistics 15.2, pp. 176–197. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] — (2010b). “A bivariate space-time downscaler under space and time misalignment”. In: The annals of applied statistics 4.4, p. 1942. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Berrocal Veronica J., Guan Yawen, et al. (2020). “A comparison of statistical and machine learning methods for creating national daily maps of ambient PM2.5 concentration”. In: Atmospheric Environment 222, p. 117130. issn: 1352–2310. doi: 10.1016/j.atmosenv.2019.117130. url: https://www.sciencedirect.com/science/article/pii/S1352231019307691 . [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Breiman Leo (2001). “Random forests”. In: Machine learning 45.1, pp. 5–32. [ Google Scholar ] Chang Howard H, Hu Xuefei, and Liu Yang (2014). “Calibrating MODIS aerosol optical depth for predicting daily PM2. 5 concentrations via statistical downscaling”. In: Journal of exposure science & environmental epidemiology 24.4, pp. 398–404. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Cressie Noel (2015). Statistics for spatial data. John Wiley & Sons. [ Google Scholar ] Cressie Noel and Wikle Christopher K (2011). Statistics for spatio-temporal data. John Wiley & Sons. [ Google Scholar ] Datta Abhirup et al. (2016a). “Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets”. In: Journal of the American Statistical Association 111.514, pp. 800–812. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] — (2016b). “On nearest-neighbor Gaussian process models for massive spatial data”. In: Wiley Interdisciplinary Reviews: Computational Statistics 8.5, pp. 162–171. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Di Qian et al. (2016). “Assessing PM2. 5 exposures with high spatiotemporal resolution across the continental United States”. In: Environmental science & technology 50.9, pp. 4712–4721. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Diggle Peter J and Ribeiro Paulo J Jr (2002). “Bayesian inference in Gaussian model-based geostatistics”. In: Geographical and environmental modelling 6.2, pp. 129–146. [ Google Scholar ] Fan Jingchun et al. (2016). “The impact of PM2. 5 on asthma emergency department visits: a systematic review and meta-analysis”. In: Environmental Science and Pollution Research 23, pp. 843–850. [ DOI ] [ PubMed ] [ Google Scholar ] Finley Andrew O., Datta Abhirup, and Banerjee Sudipto (2022). “spNNGP R Package for Nearest Neighbor Gaussian Process Models”. In: Journal of Statistical Software 103.5, pp. 1–40. doi: 10.18637/jss.v103.i05. url: https://www.jstatsoft.org/index.php/jss/article/view/v103i05 . [ DOI ] [ Google Scholar ] Fu Pengfei et al. (2019). “The association between PM2. 5 exposure and neurological disorders: a systematic review and meta-analysis”. In: Science of the Total Environment 655, pp. 1240–1248. [ DOI ] [ PubMed ] [ Google Scholar ] Gelfand AE and Smith AFM (1990). “Sampling-based approaches to calculating marginal densities”. In: Journal of the American Statistical Association 85.410, pp. 398–409. doi: 10.1080/01621459.1990.10476213. [ DOI ] [ Google Scholar ] Gong Chen et al. (2022). “Maternal exposure to ambient PM2. 5 and term birth weight: a systematic review and meta-analysis of effect estimates”. In: Science of The Total Environment 807, p. 150744. [ DOI ] [ PubMed ] [ Google Scholar ] Hastings WK (1970). “Monte Carlo sampling methods using Markov chains and their applications”. In: Biometrika 57.1, pp. 97–109. doi: 10.1093/biomet/57.1.97. [ DOI ] [ Google Scholar ] Hu Xuefei et al. (2017). “Estimating PM2. 5 concentrations in the conterminous United States using the random forest approach”. In: Environmental science & technology 51.12, pp. 6936–6944. [ DOI ] [ PubMed ] [ Google Scholar ] Katzfuss Matthias and Guinness Joseph (2021). “A general framework for Vecchia approximations of Gaussian processes”. In: Statistical Science 36.1, pp. 124–141. [ Google Scholar ] Katzfuss Matthias, Guinness Joseph, et al. (2020). “Vecchia approximations of Gaussian-process predictions”. In: Journal of Agricultural, Biological and Environmental Statistics 25.3, pp. 383–414. [ Google Scholar ] Lyapustin A and Wang Y (2018). MCD19A2 MODIS/Terra+Aqua Land Aerosol Optical Depth Daily L2G Global 1 km SIN Grid V006. NASA EOSDIS Land Processes DAAC. Accessed: 2025-05-27. url: https://lpdaac.usgs.gov/products/mcd19a2v006/ . [ Google Scholar ] Madden Wyatt G. et al. (2025). “ensembleDownscaleR: R Package for Bayesian Ensemble Averaging of PM2.5 Geostatistical Downscalers”. In: Remote Sensing 17.11. issn: 2072–4292. doi: 10.3390/rs17111941. url: https://www.mdpi.com/2072-4292/17/11/1941 . [ DOI ] [ Google Scholar ] NASA Land Data Assimilation System Team (2025). NLDAS-2 Forcing Dataset Information. https://ldas.gsfc.nasa.gov/nldas/v2/forcing . NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team (Aug. 2019). ASTER Global Digital Elevation Model V003. Release date: 2019-08-05. Sioux Falls, South Dakota, USA. doi: 10.5067/ASTER/ASTGTM.003. url: https://doi.org/10.5067/ASTER/ASTGTM.003 . [ DOI ] [ Google Scholar ] Platnick S et al. (2015). MODIS atmosphere L2 cloud product (06 L2). Tech. rep. NASA MODIS Adaptive Processing System, Goddard Space Flight Center. [ Google Scholar ] R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria. url: https://www.R-project.org/ . [ Google Scholar ] Rose A et al. (2019). LandScan USA 2018. Data set. doi: 10.48690/1523375. url: https://doi.org/10.48690/1523375 . [ DOI ] [ Google Scholar ] Tukey John W. (1958). “Bias and Confidence in Not Quite Large Samples (the Jackknife, the Bootstrap and Other Methods)”. In: Annals of Mathematical Statistics 29.2, p. 614. doi: 10.1214/aoms/1177706647. url: https://doi.org/10.1214/aoms/1177706647 . [ DOI ] [ Google Scholar ] U.S. Environmental Protection Agency (2025). Air Quality System (AQS). https://www.epa.gov/aqs . Accessed September 3, 2025. US EPA Office of Research and Development (2024). CMAQ (5.5). doi: 10.5281/zenodo.13883210. url: https://doi.org/10.5281/zenodo.13883210 . [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Williams Christopher KI and Rasmussen Carl Edward (2006). Gaussian processes for machine learning. Vol. 2. 3. MIT press Cambridge, MA. [ Google Scholar ] Wright Marvin N. and Ziegler Andreas (2017). “ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R”. In: Journal of Statistical Software 77.1, pp. 1–17. doi: 10.18637/jss.v077.i01. url: https://www.jstatsoft.org/index.php/jss/article/view/v077i01 . [ DOI ] [ Google Scholar ] Yuan Sheng et al. (2019). “Long-term exposure to PM2. 5 and stroke: a systematic review and meta-analysis of cohort studies”. In: Environmental research 177, p. 108587. [ DOI ] [ PubMed ] [ Google Scholar ] Zhang Danlu et al. (2023). “Wildland fires worsened population exposure to PM2. 5 pollution in the contiguous united states”. In: Environmental Science & Technology 57.48, pp. 19990–19998. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials 1 NIHMS2178827-supplement-1.pdf (295.3KB, pdf) ACTIONS View on publisher site PDF (1.1 MB) Cite Collections Permalink RESOURCES Similar articles Cited by other articles Links to NCBI Databases Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top

Related documents

• Model-based Detection of Spatial Disease Boundaries Using Amortized Bayesian Inference#394555
• Learning De-Biased Representations for Remote-Sensing Imagery#1082918
• Satellite Remote Sensing of Global Ecosystem Water Use Efficiency: Accuracies, Patterns, and Drivers#422580
• The impact of aerosols on polarized sky radiance: model development, validation, and applications#408457
• A semi-empirical approach to model net particulate matter fluxes between vegetated surfaces and the atmosphere#11530