Fast approximate inference with INLA: the past, the present and the future

Latent Gaussian models are an extremely popular, flexible class of models. Bayesian inference for these models is, however, tricky and time consuming. Recently, Rue, Martino and Chopin introduced the Integrated Nested Laplace Approximation (INLA) method for deterministic fast approximate inference. In this paper, we outline the INLA approximation and its related R package. We will discuss the newer components of the r-INLA program as well as some possible extensions.

💡 Research Summary

**
The paper provides a comprehensive overview of the Integrated Nested Laplace Approximation (INLA) methodology for fast deterministic approximate Bayesian inference in latent Gaussian models (LGMs) and describes the associated R package, r‑INLA. LGMs are defined by three hierarchical components: an observation model (y_i|x_i\sim\pi(y_i|x_i)), a latent Gaussian field (x|\theta\sim N(\mu(\theta),Q(\theta)^{-1})), and a hyper‑parameter prior (\theta\sim\pi(\theta)). This formulation encompasses a broad spectrum of statistical models, including dynamic linear models, generalized linear and additive mixed models, spatial and spatio‑temporal models, disease mapping, stochastic volatility, and survival analysis.

Traditional Bayesian inference for LGMs relies heavily on Markov chain Monte Carlo (MCMC). While MCMC is flexible and easy to program, its computational cost scales poorly with model dimension: convergence rates are (O(N^{-1/2})) and each sample may be expensive to compute. For large spatial or hierarchical models, MCMC can require hours or days to achieve only modest accuracy.

INLA circumvents sampling by using a series of Laplace approximations. First, a Gaussian approximation (\pi_G(x|\theta,y)) is constructed that matches the true conditional posterior at its mode and curvature. This yields a Laplace approximation to the marginal posterior of the hyper‑parameters (\pi(\theta|y)). Second, for each latent component (x_i) a conditional Laplace approximation (\pi(x_i|\theta,y)) is obtained, and the marginal (\pi(x_i|y)) is computed by numerically integrating over the low‑dimensional space of (\theta). Thus, high‑dimensional integration is reduced to a series of low‑dimensional numerical quadratures, dramatically reducing computational burden.

The r‑INLA package implements this scheme in C (via the GMRF‑Lib library) and provides an R interface that mimics the familiar glm() syntax. Random effects are specified with the f() function, and the latent field is automatically treated as a Gaussian Markov random field (GMRF) when its precision matrix is sparse. The authors illustrate the workflow with a Cox proportional hazards model for childhood leukaemia in north‑west England, showing that full Bayesian inference—including spatial random effects—can be obtained in about seven seconds, and the posterior mean spatial field is visualised.

Since the original INLA paper, several extensions have been added to r‑INLA:

• Generalised observation models – an “A‑matrix” allows each observation to depend on a linear combination of latent variables, enabling area‑averaged or aggregated data.
• Enhanced f() options – replicate handles multiple independent copies of a model sharing hyper‑parameters; copy allows the same latent field to be reused with different scalings.
• Multiple likelihoods – data from different sources can be supplied as a matrix, with each column linked to a distinct likelihood function.
• Survival analysis – the inla.surv() function implements Bayesian survival models, supporting censoring and frailty terms.
• SPDE approach – stochastic partial differential equations are used to construct computationally efficient spatial GMRFs. By generating a mesh and defining an “imatern” (intrinsic Matérn) model, users can fit smooth spatial fields with non‑stationary extensions.

Future research directions outlined include:

1. Global Gaussian approximations – addressing cases where the posterior is multimodal or highly prior‑driven, by fitting a Gaussian that approximates the entire distribution rather than only around the mode.
2. Kronecker product precision matrices – enabling separable space‑time models where the precision matrix factorises as (Q = Q_{\text{time}}\otimes Q_{\text{space}}). Preliminary implementations for AR(1)⊗Matern structures are already available.
3. Non‑Gaussian random effects – extending beyond the Gaussian assumption to gamma frailty models. The authors propose an importance‑sampling decomposition “log‑Normal | {z} Prior × Correction” to approximate gamma distributions, and suggest incorporating Bayesian non‑parametric ideas for more flexible random‑effect priors.
4. Non‑separable spatio‑temporal SPDEs – developing stochastic heat‑equation type models with temporally varying anisotropy and drift, preserving physical properties such as causality and non‑reversibility.

The paper concludes that INLA offers a powerful, fast alternative to MCMC for a large class of LGMs, especially when the latent field exhibits Markov properties that can be exploited by sparse linear algebra. The r‑INLA ecosystem is actively maintained, with frequent testing releases and a community‑driven development model: users are encouraged to request new features, which are then incorporated into the code base. The authors provide URLs for the project’s source repository and documentation, emphasizing that the software is still under active development and that the “present” features will soon become “past” as newer capabilities are added.