PaperlayerView exact revision on arXiv

Source-linked AI summary

spBayes for large univariate and multivariate point-referenced spatio-temporal data models

Andrew O. Finley, Sudipto Banerjee, Alan E. Gelfand

arXiv:1310.8192v1stat.COstat.ME

TL;DR

spBayes addresses computational and scalability challenges in modeling increasingly large, spatially and temporally indexed datasets. It reformulates core MCMC functions, adds predictive-process and dynamic spatio-temporal models, and reports improved computational efficiency, including a reduction from approximately 8 minutes to 0.069 minutes for 5,000 samples in one spLM comparison.

  • Problem

    Large spatial datasets make dense-matrix computations computationally demanding, while researchers increasingly need to model complex spatial and temporal data from diverse sources.

  • Method

    The paper rewrites spBayes core functions using more efficient MCMC formulations, predictive-process low-rank models, residual-variance adjustment, and dynamic models for continuous space and discrete time.

  • Results

    Approximately 8 minutes versus 0.069 minutes were required to generate 5,000 spLM samples under the previous and current implementations, respectively, while the package also improves scalability and sampler efficiency.

  • Takeaways & Limitations

    The reformulated package provides efficient estimation for univariate and multivariate point-referenced models and predictive-process options for increasingly large datasets.

  • Takeaways & Limitations

    When the spatial dependence range is short relative to knot spacing, covariance parameter estimation suffers and coarse knots can smooth spatial random-effects surfaces.

Abstract

from arXiv · show

In this paper we detail the reformulation and rewrite of core functions in the spBayes R package. These efforts have focused on improving computational efficiency, flexibility, and usability for point-referenced data models. Attention is given to algorithm and computing developments that result in improved sampler convergence rate and efficiency by reducing parameter space; decreased sampler run-time by avoiding expensive matrix computations, and; increased scalability to large datasets by implementing a class of predictive process models that attempt to overcome computational hurdles by representing spatial processes in terms of lower-dimensional realizations. Beyond these general computational improvements for existing model functions, we detail new functions for modeling data indexed in both space and time. These new functions implement a class of dynamic spatio-temporal models for settings where space is viewed as continuous and time is taken as discrete.

1. Introduction

Point-referenced spatial models face expensive matrix computations, slower convergence, and difficulties with multivariate outcomes as datasets grow. spBayes 0.3-7 reformulates core functions to improve efficiency, scalability, flexibility, and usability while adding dynamic spatio-temporal models.

  • Motivation: Hierarchical Bayesian models represent spatial and temporal dependence while accommodating multiple sources of uncertainty.Their flexibility supports richer models with less stringent assumptions than traditional approaches.
  • Motivation: Point-referenced models require expensive matrix computations that become prohibitive for large datasets.Multivariate spatial models add further matrix-computational difficulties.
  • Prior spBayes capabilities: Earlier spBayes models offered flexible Bayesian MCMC methods for Gaussian and non-Gaussian univariate and multivariate point-referenced data, but imposed substantial computational demands.Matrix decompositions scale cubically with the number of spatial locations.
  • Contributions: spBayes 0.3-7 reduces parameter space, avoids expensive matrix computations, and uses predictive process models to represent spatial processes with lower-dimensional realizations.These changes target convergence, runtime, and scalability for large datasets.
  • Contributions: The revised package adds dynamic spatio-temporal models for continuous space and discrete time.

2. Bayesian Gaussian spatial regression models

The paper formulates Bayesian Gaussian spatial regression models and develops MCMC strategies that reduce computation through marginalization, triangular solves, and low-rank representations. Full-rank estimation remains cubic in the number of observations, motivating specialized low-rank models for large datasets.

  • Model formulation: The Gaussian model combines regression effects, spatial random effects, covariance structures, and a prior on process parameters within a Bayesian hierarchy.Inference samples β, α, and θ from their posterior distribution.
  • MCMC inference: Integrating out β and α allows MCMC to sample process parameters from a marginal posterior before compositionally sampling the remaining effects.Only post-burn-in θ samples need to be stored for this procedure.
  • Computational strategy: Cholesky factorization dominates full-rank likelihood computation, giving O(n^3) work per evaluation while triangular solves and quadratic forms require lower-order operations.
  • Computational limits: Estimating spatial effects can become computationally prohibitive for datasets with thousands of observations, for which the paper recommends low-rank models.
  • Computational limits: For n = 2,000 locations and p = 2 predictors, full-rank estimation requires approximately 0.3 seconds per MCMC iteration and about 50 minutes for 10,000 iterations.
  • Low-rank models: Low-rank models set r << n and use the Sherman-Woodbury-Morrison formula to avoid repeated n × n factorizations when updating β and θ.The predictive process supplies the package’s low-rank construction.

3. Computing environment

spBayes implements its MCMC algorithms in C++ through R’s foreign-language interface and uses BLAS and LAPACK for matrix computations. Threaded numerical libraries can substantially reduce sampler runtime.

  • Implementation: The package uses C++ and R’s Foreign Language Interface to call Fortran BLAS and LAPACK routines.
  • Implementation: Table 1 lists common BLAS and LAPACK functions used in spBayes function calls and their descriptions.
  • Parallel computing: Threaded BLAS and LAPACK implementations enable multicore computation and produced near-linear MCMC speedup with the number of CPUs in the illustrative analyses.At least four CPUs were used in each function call.
  • R dependencies: The package also depends on coda, abind, magic, and Formula for posterior analysis, multivariate matrices, and symbolic model formulas.

4. Models offered by spBayes

spBayes offers hierarchical Gaussian spatial models, low-rank predictive-process variants, and computationally improved fitting and recovery tools for point-referenced data. The examples show faster sampling and useful prediction, while coarse knots can reduce fidelity and inference quality.

  • Gaussian spatial models: Gaussian spatial regression separates covariate effects, structured spatial variation, and independent measurement error, with Matérn parameters governing variance, range, smoothness, and nugget.The spatial process supplies local, dependent adjustment to the mean, while the nugget captures measurement error or microscale variation.
  • Model fitting and recovery: Posterior samples of β and spatial effects can be recovered from θ samples through posterior predictive composition sampling.The spRecover function performs this recovery from an spLM object, supporting spatial-effect summaries and mapping.
  • Low-rank predictive process models: spBayes implements full-rank and low-rank spatial models, using predictive processes to reduce computation through a lower-dimensional representation.Low-rank models specify Z(θ) with r << n; the package also provides a modified predictive process that adjusts for residual variance.
  • Predictive-process results: The modified predictive process removed the upward variance-parameter bias seen in non-modified models, with only marginal additional computational overhead.The full-rank model required 5.18 minutes for 5,000 posterior samples, while model iv produced comparable parameter estimates.
  • Limitations: Predictive-process accuracy can deteriorate when spatial dependence is short relative to knot spacing or when the knot array is coarse.The resulting smoothing can diminish predictive ability and, in some cases, parameter inference relative to a full-rank model.
  • Predictive-process results: 94.4 percent of holdout values were covered by the posterior predictive 95% confidence intervals in the synthetic-data analysis.The analysis expected approximately 950 of the true values to be covered.

5. Multivariate Gaussian spatial regression models

spBayes models multiple point-referenced outcomes jointly through multivariate spatial regression. Its linear model of coregionalization defines cross-outcome covariance, while full-rank and low-rank implementations support alternative computational scales.

  • Model specification: Multivariate spatial regression models m point-referenced outcomes at each location using outcome-specific predictors, coefficients, spatial effects, and random errors.The unstructured residuals follow a zero-centered multivariate normal distribution with dispersion matrix Ψ.
  • Cross-covariance: The linear model of coregionalization specifies cross-covariance as Cw(s,t) = AM(s,t)A′, combining a lower-triangular cross-outcome matrix with diagonal spatial correlations.This construction defines the covariance between spatial effects for different outcomes and locations.
  • Joint representation: Stacking outcomes across b locations yields a joint model with n = mb observations, block-diagonal residual covariance, and spatial covariance blocks AM(si,sj)A′.Positive-definiteness is ensured through the linear model of coregionalization.
  • Available implementations: spBayes implements full-rank multivariate Gaussian models and predictive-process and modified-predictive-process low-rank counterparts in spMvLM.The low-rank variants use fitting strategies analogous to those for univariate low-rank models.

6. Non-Gaussian models

spBayes implements spatial generalized linear models for binary and count outcomes, including predictive-process and multivariate extensions.

  • Binary responses use logit or probit regression, while count responses use Poisson regression.
  • Spatial GLMs model the transformed mean as η(s) = x(s)′β + w(s) using a suitable link function.
  • Gaussian models can marginalize spatial effects, but binary and Poisson models update them within the Gibbs sampler.
  • The spGLM function implements these models and their predictive-process counterparts, while spMvGLM extends them to multivariate settings.

7. Dynamic spatio-temporal models

The dynamic spatio-temporal models treat space as continuous and time as discrete, combining evolving spatial processes with regression and measurement error. spDynLM implements both full-rank and predictive-process versions, illustrated with New York ozone data.

  • 7. Dynamic spatio-temporal models: The framework views observations as a time series of spatial-process realizations, with continuous space and discrete time.
  • 7. Dynamic spatio-temporal models: The model combines a measurement equation with time-varying regression coefficients, spatial effects, and independent Gaussian measurement error.
  • 7.1. Model specification: The spatial process evolves through a transition equation that adds Gaussian-process innovations to the previous spatial effect.
  • 7.1. Model specification: The model assumes the same monitored locations across discrete times, while allowing missing cells in the resulting space-time data matrix.
  • 7.2. Predictive process models: Predictive-process dynamics replace the full spatial process and add an adjustment for oversmoothing and underestimated spatial variability.
  • 7.3. Example: The spDynLM function implements both dynamic formulations; moving to a low-rank representation requires specifying knot locations.
  • 7.3. Example: In the ozone example, 28 stations were observed over 62 days, with 114 missing measurements and four covariates including the intercept.
  • 7.3. Example: Three of 36 holdout measurements fell outside their 95% predicted intervals, corresponding to approximately 92% coverage.

8. Model choice

spDiag provides several measures for assessing and comparing fitted spatial models.

  • spDiag assesses model performance and comparison for spLM, spMvLM, spGLM, and spMvGLM objects.
  • Available measures include Deviance Information Criterion, posterior predictive loss, and a scoring rule.

9. Summary and future direction

The rewritten spBayes functions improve efficient MCMC estimation for univariate and multivariate point-referenced models, including predictive-process support for large datasets. Future work targets spatially misaligned, non-stationary, and more flexible spatio-temporal models with non-Gaussian and multivariate outcomes.

  • The reformulated core functions improve computational efficiency and model-specification flexibility for univariate and multivariate point-referenced MCMC models.The improvements come from reducing parameter space and avoiding expensive matrix operations.
  • Predictive-process models are provided across core functions to accommodate increasingly large datasets.
  • Future spBayes development will address multivariate spatially misaligned data by borrowing strength from covariance among outcomes within and across locations.
  • Planned extensions include non-stationary multivariate predictive processes and more flexible dynamic spatio-temporal models with non-Gaussian and multivariate outcomes.
Loading 1310.8192v1…