PaperlayerView exact revision on arXiv

Source-linked AI summary

Spatial modelling with R-INLA: A review

Haakon Bakka, Håvard Rue, Geir-Arne Fuglstad, Andrea Riebler, David Bolin, Elias Krainski, Daniel Simpson, Finn Lindgren

arXiv:1802.06350v2stat.MEstat.CO

TL;DR

Spatial Bayesian models can be difficult to fit efficiently and modify for realistically sized datasets. This review explains how R-INLA and the SPDE approach address that challenge through latent Gaussian modelling, sparse precision structures, and flexible spatial extensions. It concludes that the framework supports a broad range of practical spatial and space-time models, while retaining important limitations in discretisation, point-process, and interaction settings.

  • Problem

    Spatial Gaussian random-field models are computationally difficult for large datasets, while implementing fast, modifiable inference for complex Bayesian spatial models is time-consuming.

  • Method

    The review covers R-INLA latent Gaussian models and the SPDE approach, including sparse precision representations, areal and continuously indexed models, and extensions beyond isotropy and separability.

  • Results

    R-INLA enables fast Bayesian inference for complex spatial models and supports a wide range of likelihoods, model components, and spatial or space-time structures.

  • Takeaways & Limitations

    R-INLA makes previously infeasible spatial Bayesian models practical to fit, modify, compare, and assess for sensitivity and prediction.

  • Takeaways & Limitations

    The framework retains limitations involving discretised areal models, analytically intractable spatial point-process likelihoods, and computationally costly identifiability constraints in space-time interactions.

Abstract

from arXiv · show

Coming up with Bayesian models for spatial data is easy, but performing inference with them can be challenging. Writing fast inference code for a complex spatial model with realistically-sized datasets from scratch is time-consuming, and if changes are made to the model, there is little guarantee that the code performs well. The key advantages of R-INLA are the ease with which complex models can be created and modified, without the need to write complex code, and the speed at which inference can be done even for spatial problems with hundreds of thousands of observations. R-INLA handles latent Gaussian models, where fixed effects, structured and unstructured Gaussian random effects are combined linearly in a linear predictor, and the elements of the linear predictor are observed through one or more likelihoods. The structured random effects can be both standard areal model such as the Besag and the BYM models, and geostatistical models from a subset of the Matérn Gaussian random fields. In this review, we discuss the large success of spatial modelling with R-INLA and the types of spatial models that can be fitted, we give an overview of recent developments for areal models, and we give an overview of the stochastic partial differential equation (SPDE) approach and some of the ways it can be extended beyond the assumptions of isotropy and separability. In particular, we describe how slight changes to the SPDE approach leads to straight-forward approaches for non-stationary spatial models and non-separable space-time models.

1 Introduction

The paper reviews R-INLA and the SPDE approach as computationally efficient frameworks for fitting and extending spatial models. It explains how sparse precision structures, finite-element approximations, and flexible model components support large and diverse applications.

  • Motivation: R-INLA addresses the computational difficulty of spatial modelling with Gaussian random fields, whose direct covariance-based inference becomes infeasible for very large datasets.The review positions the SPDE approach as a computationally efficient alternative implemented in R-INLA.
  • R-INLA and SPDE: R-INLA shifts computation from covariance matrices to sparse precision matrices, exploiting conditional independence among spatially close variables.Sparse precision matrices store mostly zero entries and encode the relevant dependence structure.
  • R-INLA and SPDE: The SPDE approach constructs a continuously indexed approximation of a Matérn Gaussian random field over the study area rather than a discrete covariance model at observation locations.Finite-element methods permit complex polygonal boundaries and spatially varying discretisation fidelity.
  • SPDE advantages: The finite-dimensional Gaussian approximation depends on desired resolution and is invariant to the number of observations.This property is central to scaling spatial inference to large datasets.
  • SPDE advantages: The SPDE precision matrix has a Markovian structure with non-zero entries only for close neighbours, while its sparsity pattern is invariant to the correlation range.The number of non-zero neighbours depends on the smoothness parameter ν.
  • Scope and applications: R-INLA supports complex spatial models, multiple likelihoods, and extensions beyond stationary isotropic models, including non-stationary and non-separable space-time models.The review covers areal models, continuously indexed models, and applications ranging from geostatistical models to spatial point processes.

2 Notation and background on R-INLA

R-INLA represents latent Gaussian models through linear predictors, sparse precision matrices, and projection matrices, then uses Laplace approximations and numerical integration for inference. Its modular structure supports varied likelihoods, spatial and temporal components, and scalable prediction workflows.

  • Model notation: R-INLA models conditionally independent observations through a linear predictor composed of fixed and Gaussian random effects.Different observations may use different likelihoods, while hyper-parameters govern likelihood and random-effect distributions.
  • Model notation: Each Gaussian random effect has a precision matrix governed by hyper-parameters and a fixed sparsity graph, with a sparse projection matrix linking latent effects to observations.The projection matrix A specifies which latent-effect entries contribute to each observation.
  • Inference: Inference becomes expensive when the hyper-parameter dimension is large because evaluating π(θ|y) requires precision-matrix factorisation and a Laplace approximation.The review recommends fewer than 20 hyper-parameters, preferably around 5.
  • Inference: Good hyper-parameter parameterisations help keep posterior distributions sufficiently regular for numerical integration, such as using log(σ) instead of σ.The posterior of σ is described as skewed and heavy-tailed, whereas log(σ) is better behaved.
  • Space-time models: Kronecker-structured space-time models allow spatial and temporal components to be mixed and matched with little additional implementation work.Temporal choices include random walks, autoregressive models, and iid replicates.
  • Prediction: Posterior samples can make large spatial prediction tasks substantially faster than adding missing-data rows at every prediction location.The review reports similar prediction results in 11 minutes using posterior samples versus 24 hours with the standard approach.
  • Inference: R-INLA uses Laplace approximations for non-Gaussian likelihoods and numerical integration over hyper-parameter posteriors to compute latent-variable marginals.The Laplace approximation targets the conditional posterior of η given y and θ rather than the likelihood itself.

3 Areal models and other discretely indexed models

R-INLA supports areal models such as Besag and BYM, while recent developments improve prior interpretability and extend models to space-time interactions. Discretely indexed models remain computationally useful but have important spatial and identifiability limitations.

  • Discretely indexed models use finite sets of indices and do not straightforwardly extend to other locations, areas, or covariate values.
  • 3.1 The Besag model: The Besag model assigns region values conditionally on neighbouring regions, with neighbourhood structure represented by a sparse precision matrix.Neighbourhoods are usually defined by shared borders, and the structure matrix has O(n) non-zero elements.
  • 3.1 The Besag model: The BYM model combines the structured Besag effect with an additional unstructured random effect for each region.
  • 3.2 Prior specification for the Besag model: The Besag precision parameter τu is difficult to interpret because its meaning depends on the underlying graph, including region count and connectivity.Scaling the precision matrix can make τu represent marginal deviation from a constant level independently of graph structure, facilitating prior specification.
  • 3.2 Prior specification for the Besag model: BYM2 reparameterizes the model with a scaled Besag component and explicitly models variance allocation between structured and unstructured effects.
  • 3.3 Discretely indexed models: Discretely indexed models are computationally efficient for smoothing but ignore sub-regional variation and can mishandle differing region sizes, shared boundaries, and changing geographic boundaries.Space-time interaction models can require many identifiability constraints, whose O(nk^2) cost may make large models computationally infeasible.

4 Continuously indexed models

The SPDE approach represents continuously indexed spatial fields through computationally efficient sparse precision structures, while FEM meshes support irregular domains and observation locations. This framework connects interpretable Matérn parameters with flexible approximations and accommodates point and areal observations.

  • Representing spatial models: Direct covariance-matrix inference becomes computationally infeasible for large spatial datasets because the precision matrix is generally not sparse.R-INLA instead targets sparse precision structures associated with conditional independence.
  • Representing spatial models: The SPDE approach constructs a continuous approximation on the study area rather than a discrete covariance model tied to locations or grid cells.This shifts the modelling focus from CAR parameters to the continuous GRF representation.
  • Constructing the approximation: Finite element methods use triangular meshes and linear hat-functions to construct continuously indexed approximations on complex domains.The mesh covers the domain and extends beyond it to account for boundary conditions.
  • Implementation and accuracy: Mesh quality controls discretisation error, with smaller and more regularly shaped triangles generally giving smaller errors.The error also depends on the random field’s correlation scale and smoothness.
  • Implementation and accuracy: The A-matrix expresses both point and areal observations as linear combinations of mesh coefficients, allowing them to share one modelling framework.This supports examples such as combining weather-station precipitation with catchment-area runoff.

5 Spatial modelling with R-INLA

R-INLA combines efficient spatial representations with latent Gaussian modelling to fit diverse spatial data models. Its capabilities include non-Gaussian and multi-likelihood models, continuously indexed point processes, and separable or non-separable spatio-temporal extensions.

  • Spatial GLMs and GLMMs: R-INLA’s computationally efficient spatial effects combine with its general modelling features to support many types of spatial models.The standard Gaussian geostatistical model is one basic continuously indexed example.
  • Spatial GLMs and GLMMs: Changing the observation process enables spatial GLMs for count, non-Gaussian, zero-inflated, and hurdle responses.Supported likelihoods include Binomial, Negative Binomial, Poisson, t-distribution, skew-Normal, and Gamma likelihoods.
  • Multiple likelihoods: Assigning different likelihoods to observation groups enables multivariate responses, misaligned covariates, marked point patterns, and replicated point patterns.These models can share component structures across groups.
  • Spatial point processes: Log-Gaussian Cox processes fit R-INLA’s latent Gaussian framework by modelling log-intensity with a geostatistical spatial signal.They represent doubly stochastic inhomogeneous Poisson processes suitable for aggregated point patterns.
  • Spatial point processes: Continuously indexed point-process models avoid lattice binning error but require an additional numerical integration approximation, making them more cumbersome to specify.The likelihood is analytically intractable and is approximated with a Poisson-form numerical integration.
  • Spatio-temporal models: Spatio-temporal extensions include separable Kronecker-product models and non-separable models combining the SPDE approach with temporal dynamics.The non-separable construction has precision-matrix sparsity similar to the separable model, without adding computational burden.

6 Adding complexity to the spatial effect

The SPDE approach extends spatial models beyond stationary and isotropic Matérn structure through local changes to the differential operator, while preserving valid covariance structures. The review discusses spatially varying covariance, physical barriers, and estimation challenges.

  • Stationarity and isotropy simplify covariance construction but are strong assumptions that may not hold for realistic spatial processes.
  • Local changes to the SPDE differential operator enable more complex covariance structures while ensuring global covariance validity.
  • Spatially varying covariance: Spatially varying basis-function parameters can model changing marginal variance across the study area.Figure 3 illustrates variance changing from west to east as τ varies.
  • Spatially varying covariance: Estimating covariates in both mean and covariance structures from a single realization can produce inaccurate estimates because the structures are not separately identifiable.Ingebrigtsen et al. investigate estimation using multiple realizations.
  • Physical barriers: Barrier models represent physical obstacles by reducing the operator parameter in barrier regions, forcing spatial smoothing around rather than across them.The approach is motivated by aquatic data where land obstructs correlation.

6.3 Spatially varying anisotropy

The SPDE framework supports spatially varying anisotropy, manifold domains, and non-separable space-time models. These extensions retain computationally useful sparse precision structures while representing directional, geometric, and temporal complexity.

  • Spatially varying anisotropy: Geometric anisotropy models directional dependence by applying a linear transformation through a positive-definite matrix H.
  • Spatially varying anisotropy: Allowing H to vary spatially produces continuously changing anisotropy, while varying both κ and H controls marginal deviations and correlation structure.The simulated field changes from extra horizontal dependence in the west to extra vertical dependence in the east.
  • SPDEs on manifolds: Two-dimensional meshes allow SPDEs to be solved on manifolds such as the sphere using the same approach as triangulated Euclidean domains.This supports models for global spatial data without relying on a rectangular planar subdomain.
  • Non-separable space-time model: Non-separable space-time models can represent phenomena linked to the heat equation, such as temperature evolving globally.Their precision matrices have sparsity similar to those of separable models, avoiding added computational burden.

6.6 General smoothness

The review covers extensions controlling field smoothness and allowing non-Gaussian spatial behavior. These extensions broaden the model class, but computational and inferential limitations remain for some cases.

  • General smoothness: The operator exponent α/2 determines field smoothness, but the standard SPDE approach is computable only when α is an integer.This restriction may matter when smoothness is difficult to estimate or substantively important.
  • General smoothness: A rational SPDE method is computable for any α > 0 and can be more accurate for Matérn models than the parsimonious fractional approximation.It combines finite-element spatial approximation with a rational approximation of x^−α/2.
  • Non-Gaussian spatial fields: Replacing Gaussian white noise with non-Gaussian noise produces Matérn-like fields with asymmetric or skewed marginal behavior.Normal inverse Gaussian noise adds parameters controlling skewness and marginal tails.
  • Non-Gaussian spatial fields: NIG-driven SPDE models are conditionally Gaussian and easy to simulate with R-INLA, but INLA inference is unavailable because the models are intrinsically non-Gaussian.
  • Additional work: SPDE extensions also include multivariate fields, flexible anisotropic and oscillating covariance structures, and adaptive Bayesian splines.

7 R-packages building on R-INLA

Several R-packages build on R-INLA to simplify spatial modelling, add model components, support post-processing, or provide user-friendly interfaces for applied researchers.

  • Diseasemapping provides Poisson regression with fixed effects, BYM or BYM2 models, and PC priors.
  • SSTCDapp offers a Shiny interface for estimating discrete space and space-time models without writing R code.It also provides descriptive statistics and supports several spatial and temporal components.
  • Geostatsp simplifies SPDE-based geostatistical modelling, while excursions calculates credible regions and excursion sets from R-INLA output.
  • Inlabru simplifies spatial and spatiotemporal model specification by hiding inla.stack code and extends models to mildly nonlinear predictors.It also supports nonlinear posterior prediction from independent posterior samples.

8 Discussion

R-INLA makes complex spatial and spatio-temporal Bayesian models computationally feasible and comparatively easy to extend, though model reformulation and preprocessing require substantial effort. Important frontiers include prior specification, non-stationary and non-separable models, and scalability for high-dimensional space-time problems.

  • R-INLA combines spatial effects with fixed effects, other random effects, and complex or multiple likelihoods with little additional implementation effort.The INLA method also enables fast Bayesian inference and repeated runs such as cross-validation.
  • The SPDE approach converts continuously indexed spatial models into sparse structures suitable for R-INLA.This extends precision-based methods beyond directly specified CAR models while retaining computational advantages.
  • SPDE-based models incur implementation and preprocessing costs from constructing SPDEs, triangulations, and GMRF representations.The authors regard these costs as unavoidable for efficient computation.
  • Using the FEM or SPDE framework requires learning to construct a reasonable mesh and understand its connection to spatial resolution, especially for statistical researchers.The resulting models are described as generally straightforward and stable once the approach is understood.
  • Current development focuses on implementing and documenting non-stationary and non-separable models before establishing when each model should be used.Non-stationary models often require replicates, placing them in competition with non-separable models for space-time data.
  • R-INLA handles space-time representations with a hundred thousand dimensions, but a 100 × 100 spatial-temporal grid producing one million dimensions exceeds its capability at the time of review.Parallel computing and improved approximations are being investigated as possible solutions.
  • The review frames model choice as finding the best implementable option within available time, especially for applications requiring advanced models.The authors report success across many applications and continued development to extend applied statistical modelling.

A.1 Details on the BYM2 model

The BYM2 model uses penalised complexity priors to control its marginal precision and structured-component weight, shrinking flexibly toward simpler spatial models.

  • The BYM2 model assigns a penalised complexity prior to τ through an exponential prior on its standard deviation 1/√τ.The parameter λ is elicited from the probability statement π((1/√τ) > U) = α, yielding λ = −log(α)/U.
  • The marginal standard deviation assumption 0.92/2.58 ≈0.36 gives U = 1.16 and α = 1% for the PC prior on τ.
  • The PC prior for w represents the proportion of variability contributed by the structured component and is specified through Prob(w < U) = α.Its parameter has no closed-form expression.
  • PC priors shrink the flexible BYM2 model first toward no spatial variation and then toward a model with only unstructured variability.This sequential shrinkage is reported as unavailable with other commonly used priors.

A.2 Details on the approximation of the Mat´ern field on a grid

The grid approximation replaces a second-order differential operator with a sparse central-difference matrix whose discrete output closely approximates the operator applied to the continuous function.

  • Computational stencils encode matrix entries according to the central grid cell and its neighbours, with h−2 multiplying the entire matrix.Stencils may be rotationally symmetric or use larger neighbourhoods.
  • Applying the operator to f(x) = 5x2 + sin(15x) produces the second-order derivative of f.
  • A second-order differential operator is discretised on an equally spaced grid using a central-difference computational stencil.For grid spacing h, the stencil represents the standard central-difference approximation to the second derivative.
  • Figure 8 compares the continuous function and its vector discretisation with the discrete operator output and the discretised transformed function.The blue-circle operator values are close to the discretisation of g = L1Df.
  • The resulting matrix has −2h−2 on the diagonal, 1h−2 on the upper and lower diagonals, and zero elsewhere.This structure corresponds to first-order neighbouring grid cells, with periodic boundary conditions used to avoid boundary corrections.
Loading 1802.06350v2…