Max-and-Smooth: a two-step approach for approximate Bayesian inference in latent Gaussian models

TL;DR

Max-and-Smooth is a two-step approximate Bayesian inference method for latent Gaussian models that offers efficient and accurate inference in high-dimensional settings by propagating uncertainty and leveraging Gaussian approximations.

Contribution

The paper introduces Max-and-Smooth, a novel two-step inference scheme that simplifies Bayesian inference in LGMs using Gaussian approximations, enabling faster computation.

Findings

Method is accurate on spatial and spatio-temporal data.

Significant speedup over traditional MCMC methods.

Effective in high-dimensional and large data scenarios.

Abstract

With modern high-dimensional data, complex statistical models are necessary, requiring computationally feasible inference schemes. We introduce Max-and-Smooth, an approximate Bayesian inference scheme for a flexible class of latent Gaussian models (LGMs) where one or more of the likelihood parameters are modeled by latent additive Gaussian processes. Max-and-Smooth consists of two-steps. In the first step (Max), the likelihood function is approximated by a Gaussian density with mean and covariance equal to either (a) the maximum likelihood estimate and the inverse observed information, respectively, or (b) the mean and covariance of the normalized likelihood function. In the second step (Smooth), the latent parameters and hyperparameters are inferred and smoothed with the approximated likelihood function. The proposed method ensures that the uncertainty from the first step is correctly propagated to the second step. Since the approximated likelihood function is Gaussian, the approximate posterior density of the latent parameters of the LGM (conditional on the hyperparameters) is also Gaussian, thus facilitating efficient posterior inference in high dimensions. Furthermore, the approximate marginal posterior distribution of the hyperparameters is tractable, and as a result, the hyperparameters can be sampled independently of the latent parameters. In the case of a large number of independent data replicates, sparse precision matrices, and high-dimensional latent vectors, the speedup is substantial in comparison to an MCMC scheme that infers the posterior density from the exact likelihood function. The proposed inference scheme is demonstrated on one spatially referenced real dataset and on simulated data mimicking spatial, temporal, and spatio-temporal inference problems. Our results show that Max-and-Smooth is accurate and fast.

Tables

Table 1. Table 1: Comparison of leave-one-year-out predictive performance of the four forecast schemes, CLIM, MLE, SPAT1, SPAT2 (descriptions in main text), by mean squared error of the predictive mean (MSE), continuous ranked probability score (CRPS), average width of the central 95% prediction interval (W95), and average coverage frequencies of the predictive 5, 50, and 95-percentiles (COV05, COV50, COV95, respectively).

	CLIM	MLE	SPAT1	SPAT2
MSE	0.8329	0.7843	0.7747	0.7742
CRPS	0.4952	0.4843	0.4817	0.4814
W95	3.31	3.35	3.08	3.08
COV05	0.056	0.049	0.06	0.06
COV50	0.50	0.51	0.51	0.51
COV95	0.93	0.94	0.93	0.93

Equations

$${\eta }_i={\beta }_0+\sum _{k=1}^K{\beta }_k{z}_{i,k}+\sum _{j=1}^J{u}_j{a}_{i,j}+{ϵ}_i,$$

$$\pi (\mathbf{x}|\mathbf{\theta })=\text{N}(\mathbf{x}|\mathbf{\mu }(\mathbf{\theta }),\Sigma (\mathbf{\theta })),$$

$${\mathbf{\psi }}_1=({\psi }_{1,1},\dots ,{\psi }_{1,G}{)}^{\mathsf{T}},{\mathbf{\psi }}_2=({\psi }_{2,1},\dots ,{\psi }_{2,G}{)}^{\mathsf{T}},\dots ,{\mathbf{\psi }}_M=({\psi }_{M,1},\dots ,{\psi }_{M,G}{)}^{\mathsf{T}},$$

$$\pi (\mathbf{y}|{\mathbf{\psi }}_1,{\mathbf{\psi }}_2,\dots ,{\mathbf{\psi }}_M)=\prod _{i=1}^G\pi ({\mathbf{y}}_i|{\psi }_{1,i},{\psi }_{2,i},\dots ,{\psi }_{M,i})=\prod _{i=1}^G\prod _{j\in {\mathbb{A}}_i}\pi (y_{ij}|{\psi }_{1,i},{\psi }_{2,i},\dots ,{\psi }_{M,i}),$$

$$\boldsymbol{\eta}_{1}=(\eta_{1,1},\ldots,\eta_{1,G})^{\mathsf{T}},\ \boldsymbol{\eta}_{2}=(\eta_{2,1},\ldots,\eta_{2,G})^{\mathsf{T}},\ \ldots,\ \boldsymbol{\eta}_{M}=(\eta_{M,1},\ldots,\eta_{M,G})^{\mathsf{T}},\$$

$${\mathbf{\eta }}_1$$

$${\mathbf{\eta }}_2$$

$$⋮$$

$${\mathbf{\eta }}_M$$

$$\mathbf{\eta }=({\mathbf{\eta }}_1^{\mathsf{T}},{\mathbf{\eta }}_2^{\mathsf{T}},\dots ,{\mathbf{\eta }}_M^{\mathsf{T}}{)}^{\mathsf{T}},$$

$$\mathbf{\nu }=({\mathbf{\beta }}_1^{\mathsf{T}},{\mathbf{u}}_1^{\mathsf{T}},{\mathbf{\beta }}_2^{\mathsf{T}},{\mathbf{u}}_2^{\mathsf{T}},\dots ,{\mathbf{\beta }}_M^{\mathsf{T}},{\mathbf{u}}_M^{\mathsf{T}}{)}^{\mathsf{T}},$$

$$\mathbf{ϵ}=({\mathbf{ϵ}}_1^{\mathsf{T}},{\mathbf{ϵ}}_2^{\mathsf{T}},\dots ,{\mathbf{ϵ}}_M^{\mathsf{T}}{)}^{\mathsf{T}}.$$

$${\mathbf{\mu }}_{\nu }=({\mathbf{\mu }}_{\beta ,1}^{\mathsf{T}},{\mathbf{\mu }}_{u,1}^{\mathsf{T}},{\mathbf{\mu }}_{\beta ,2}^{\mathsf{T}},{\mathbf{\mu }}_{u,2}^{\mathsf{T}},\dots ,{\mathbf{\mu }}_{\beta ,M}^{\mathsf{T}},{\mathbf{\mu }}_{u,M}^{\mathsf{T}}{)}^{\mathsf{T}},$$

$$Q_{\nu }=\text{bdiag}(Q_{\beta ,1},Q_{u,1},Q_{\beta ,2},Q_{u,2},\dots ,Q_{\beta ,M},Q_{u,M}).$$

$$Q_{ϵ}=\text{bdiag}(Q_{ϵ,1},Q_{ϵ,2},\dots ,Q_{ϵ,M}).$$

$$Z=\left(\begin{array}{cccccccccc}{X}_1& A_1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& X_2& A_2& 0& 0& 0& 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & \ddots & ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& 0& 0& 0& X_{M-1}& A_{M-1}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& X_M& A_M\end{array}\right),$$

$$\mathbf{\eta }=Z\mathbf{\nu }+\mathbf{ϵ}.$$

$$\pi (\mathbf{\eta },\mathbf{\nu },\mathbf{\theta }|\mathbf{y})\propto \pi (\mathbf{\theta })\pi (\mathbf{\eta },\mathbf{\nu }|\mathbf{\theta })\pi (\mathbf{y}|\mathbf{\eta }),$$

$$\hat{L}(\mathbf{\eta }|\mathbf{y})=\text{N}(\mathbf{\eta }|\hat{\mathbf{\eta }},{\Sigma }_{\eta y}),$$

$${\mathcal{I}}_{\eta y_i}=-{\nabla }^2\log (L({\mathbf{\eta }}_i|{\mathbf{y}}_i)){|}_{{\eta }_i={\hat{\eta }}_i},$$

$$\hat{L}({\mathbf{\eta }}_i|{\mathbf{y}}_i)=\text{N}({\mathbf{\eta }}_i|{\hat{\mathbf{\eta }}}_i,{\Sigma }_{\eta y_i})$$

$$\hat{\pi }(\mathbf{\eta },\mathbf{\nu },\mathbf{\theta }|\mathbf{y})$$

$$\propto \pi (\mathbf{\theta })\pi (\mathbf{\eta },\mathbf{\nu }|\mathbf{\theta })\text{N}(\mathbf{\eta }|\hat{\mathbf{\eta }},{\Sigma }_{\eta y})$$

$$\propto \pi (\mathbf{\theta })\pi (\mathbf{\eta },\mathbf{\nu }|\mathbf{\theta })\prod _{i=1}^G\text{N}({\mathbf{\eta }}_i|{\hat{\mathbf{\eta }}}_i,{\Sigma }_{\eta y_i}).$$

$$L({\mathbf{\eta }}_i|{\mathbf{y}}_i)\approx d_i\tilde{L}({\mathbf{\eta }}_i|{\mathbf{y}}_i)=d_i\text{N}({\mathbf{\eta }}_i|{\tilde{\mathbf{\eta }}}_i,{\Omega }_{\eta y_i}),$$

$${\tilde{\mathbf{\eta }}}_i=\int {\mathbf{\eta }}_i{d}_i^{-1}L({\mathbf{\eta }}_i|{\mathbf{y}}_i)d{\mathbf{\eta }}_i,$$

$${\Omega }_{\eta y_i}=\int ({\mathbf{\eta }}_i-{\hat{\mathbf{\eta }}}_i)({\mathbf{\eta }}_i-{\hat{\mathbf{\eta }}}_i{)}^{\mathsf{T}}{d}_i^{-1}L({\mathbf{\eta }}_i|{\mathbf{y}}_i)d{\mathbf{\eta }}_i.$$

$$\tilde{\pi }(\mathbf{\eta },\mathbf{\nu },\mathbf{\theta }|\mathbf{y})$$

$$\pi (\hat{\mathbf{\eta }}|\mathbf{\eta },Q_{\eta y},\mathbf{\theta })$$

$$\pi (\mathbf{\eta }|\mathbf{\nu },\mathbf{\theta })$$

$$\pi (\mathbf{\nu }|\mathbf{\theta })$$

$$\pi (\mathbf{\eta },\mathbf{\nu },\mathbf{\theta }|\hat{\mathbf{\eta }})$$

$$\propto \pi (\mathbf{\theta })\pi (\mathbf{\eta },\mathbf{\nu }|\mathbf{\theta })\text{N}(\hat{\mathbf{\eta }}|\mathbf{\eta },Q_{\eta y}^{-1})$$

$$\propto \pi (\mathbf{\theta })\pi (\mathbf{\eta },\mathbf{\nu }|\mathbf{\theta })\hat{L}(\mathbf{\eta }|\hat{\mathbf{\eta }},{\Sigma }_{\eta y}).$$

$$\pi (\mathbf{\theta }|\hat{\mathbf{\eta }})\propto \pi (\mathbf{\theta })\frac{\pi (\hat{\mathbf{\eta }}|\mathbf{x},\mathbf{\theta })\pi (\mathbf{x}|\mathbf{\theta })}{\pi (\mathbf{x}|\hat{\mathbf{\eta }},\mathbf{\theta })}.$$

Full text

1 Introduction

Data are being generated today at an unprecedented rate. Many datasets are large and exhibit complex marginal behaviors and dependence structures. In particular, data that are indexed in space and time may indicate non-linear time trends and spatial patterns, and may be driven by complex space-time interactions. Statistical modeling and inference for high-dimensional spatio-temporal datasets becomes increasingly computationally demanding as the spatial and/or temporal dimension increases (see, e.g., Heaton et al., 2019). The same is true for, e.g., multinomial data (see, e.g., Gustafson et al., 2008) and survival data (see, e.g., Lee et al., 2011), i.e., as the number of units or individuals increases the computational cost increases.

In this paper we focus on latent Gaussian models (LGMs), which form a general and very flexible class of models that has proven to be useful in a wide range of concrete applications (see, e.g., Gelfand et al., 2007; Cooley et al., 2007; Rue et al., 2009; Margeirsson et al., 2010; Sigurdarson and Hrafnkelsson, 2016; Zinszer et al., 2017; Opitz et al., 2018; Lombardo et al., 2018, 2019). We here introduce Max-and-Smooth, a novel approximate Bayesian inference procedure for LGMs with independent data replicates that is both accurate and fast, providing significant speedups in high dimensions. Our approach has superficial similarities with the recent contribution by Risser et al. (2019), who propose a frequentist two-step inference approach and focus on the GEV distribution. In contrast, our proposed inference scheme is more general and designed for fully Bayesian inference. LGMs are Bayesian hierarchical models (Cressie and Wikle, 2011; Banerjee et al., 2014) that consist of three levels: the data level, the latent level and the hyperparameter level. Each level is specified with a probability distribution, the one at the latent level being Gaussian. For computational reasons, it is common to assume that the data are conditionally independent, given the latent process, and we also assume this here. The role of the latent process is to capture the underlying dynamics of the data (such as space-time dependence interactions). Our focus in this paper is mainly on three types of spatio-temporal LGMs that are useful in different settings, although our proposed method applies more generally, e.g., with (replicated) multinomial or survival data. The first type of LGMs that we consider assumes that the spatio-temporal dynamics of the data is described by latent parameters that vary spatially but are constant in time, and that (a potentially different number of) data time replicates are available at each spatial location. This type of LGMs focuses on capturing the data’s spatial behavior, although datasets with temporal covariates or slowly-varying temporal trends may also be modeled in the same framework. The second type of LGMs assumes that the latent parameters vary in time, and that several spatial replicates of the data are available at each time point. In this setting the latent parameters are usually constant in space, although they may also refer to the main effects of distinct regions that vary over time. The focus is therefore on capturing the data’s temporal behavior. Finally, the third type of LGMs assumes that the latent process varies in both space and time, and that several replicates for each spatio-temporal observation are available.

Several strategies have been proposed to fit LGMs. Simulation-based Markov chain Monte Carlo (MCMC) methods can be used (see Cressie and Wikle, 2011; Banerjee et al., 2014), although their application in high-dimensional settings (i.e., situations with either space-rich and/or time-rich data, or with many parameters involved at the latent level) may be limited by the computational complexity. In order to make MCMC sampling more efficient, Knorr-Held and Rue (2002) proposed a single block updating strategy for LGMs characterized by a univariate link function. Their strategy reduces the cross-correlation between the hyperparameters and the latent parameters within the posterior samples. A detailed comparison of several sampling strategies for LGMs in Filippone et al. (2013) showed that the single block updating strategy of Knorr-Held and Rue (2002) has larger effective sample size compared to sufficient augmentation, ancillary augmentation and ancillarity-sufficiency interweaving strategy (Yu and Meng, 2011), and the surrogate method (Murray and Adams, 2010). Another approach to infer LGMs was proposed by Filippone and Girolami (2014), who suggested using a pseudo-marginal sampling procedure for the marginal posterior density of the hyperparameters, which relies on the Metropolis-Hastings algorithm and importance sampling. Essentially, samples from the marginal posterior density of the hyperparameters are obtained first, and the latent parameters are sampled from the conditional posterior density of the latent parameters. Filippone and Girolami (2014) compared their pseudo-marginal approach to ancillary augmentation (Yu and Meng, 2011) and the surrogate method (Murray and Adams, 2010), and found that the effective sample size of the hyperparameters was much lower for ancillary augmentation and the surrogate method than for the pseudo-marginal approach. Filippone and Girolami (2014) concluded that this was due to ancillary augmentation and the surrogate method not being fully capable of breaking down the correlation between hyperparameters and latent parameters. The findings of Filippone et al. (2013) and Filippone and Girolami (2014) underline that sampling from the marginal posterior density of the hyperparameters in an LGM leads to more effective sampling schemes.

Alternatively, the integrated nested Laplace approximation (INLA) has proven to be very fast and accurate for approximate Bayesian inference in LGMs (Rue et al., 2009). The INLA methodology essentially bypasses MCMC sampling by performing a numerical approximation of the posterior density. Due to its computational efficiency and its convenient implementation in the package INLA for the R statistical computing environment, the INLA method has found widespread interest, and has been applied in numerous settings; see the review papers by Rue et al. (2017) and Bakka et al. (2018), and references therein. However, the current implementation in the INLA package only supports LGMs characterized by a univariate link function (i.e., with one single Gaussian linear predictor at the latent level), and with a small number (typically less than $20$) of hyperparameters. In Section 5, we discuss a linear regression model for spatio-temporal meteorological data, where the intercept, the covariate effect, and the residual variance vary spatially, thus requiring a trivariate link function in the likelihood. Generally speaking, it is common to assume that spatio-temporal data are described by LGMs of type (i), (ii) or (iii) above with multiple parameters (e.g., intercept, multiple covariate effects, scale and shape parameters, etc.) that vary spatially and/or temporally. These types of LGMs usually require multivariate link functions, and we will hereafter refer to models of this type as extended LGMs in line with Geirsson et al. (2020). Although it might be possible in principle to extend the INLA software to LGMs with multivariate link functions, this has not been implemented yet.

Posterior inference for extended LGMs in moderate or high dimensions is known to be challenging. Geirsson et al. (2020) developed an efficient block Gibbs sampling scheme, referred to as the LGM sampler, which was shown to significantly reduce the autocorrelation in posterior samples. In this paper, we propose using Max-and-Smooth, a novel two-step approximate Bayesian inference approach for extended LGMs, which borrows ideas from the INLA method and the LGM sampler, in order to be both fast and accurate in high dimensions. Essentially, our approach approximates the likelihood function by a Gaussian likelihood, similar to the Laplace approximation used in the INLA method. This allows us to perform fast inference with a correct propagation of the uncertainty. The two steps of the inference scheme are as follows: (i) In the first step (Max), we compute the maximum likelihood (ML) estimates of the latent parameters at each spatial, temporal, or spatio-temporal point (depending on the type of LGM considered), and we approximate the variance of the Gaussian approximation using the inverse observed information evaluated at the ML estimate. We also consider an alternative Gaussian approximation that uses the mean and covariance of the normalized likelihood function. (ii) In the second step (Smooth), we treat the ML estimate (or the mean of the normalized likelihood) as the observed data of the latent parameters, with a Gaussian likelihood (thus taking their estimated variance into account). We then fit the latent Gaussian model by taking advantage of the conjugacy properties of the approximate Gaussian–Gaussian model, which is hereafter referred to as the pseudo model. In other words, we essentially consider that the parameter estimates from the first step are noisy measurements of the latent field (with known noise variance) and we smooth them jointly in the second step. Notice that although our proposed approach has two consecutive steps, it properly propagates the uncertainty, and thus provides a valid approximate procedure to sample from the full posterior density. Our proposed procedure is very fast for a variety of reasons. First, the (approximate) sampling scheme is such that the unnormalized marginal posterior density of the hyperparameters can be expressed analytically, making it straightforward to sample from it. Second, the conditional density of the latent parameters given the hyperparameters is Gaussian, which is straightforward to sample from. Third, because the hyperparameters can be sampled independently from the latent parameters, similarly to Geirsson et al. (2020), their cross-correlation is reduced, which yields better MCMC mixing properties with a higher effective sample size. Fourth, as the computational cost of the second step (i.e., fitting the pseudo model) does not depend on the number of data replicates, our proposed procedure is especially well-suited for datasets with a large number of independent replicates. Finally, further speed-ups can be obtained by specifying the Gaussian prior density for the parameter vector at the latent level to be a Gaussian Markov random field (GMRF, Rue and Held, 2005) with a sparse precision matrix. When such GMRFs are used, Max-and-Smooth scales well with increasing dimension of the latent parameter vector.

Our proposed methodology involves approximations at two levels. First, the likelihood function is approximated by a Gaussian likelihood and may therefore be misspecified. Second, the variance of the Gaussian approximation to the likelihood has to be estimated from data in the first step, but is then treated as exact in the second step. Intuitively, if the shape of the “true” likelihood function is close to a Gaussian likelihood, our inference approach will be accurate. With a perfectly (or nearly) Gaussian likelihood, our inference approach will be very close to being exact provided the variance of the Gaussian approximation is properly estimated. In contrast, when the number of data replicates per latent model parameter is low, the Gaussian approximation may become a poor approximation to the likelihood, which might negatively impact the posterior inference. However, owing to the asymptotic behavior of the likelihood function in posterior inference, see, e.g., Schervish (1995), the errors of these two levels of approximation will typically become negligible as the number of data replicates per latent model parameter grows. In other words, Max-and-Smooth is expected to perform increasingly well as the number of data replicates gets larger, with a negligible effect on the overall computational time. The Gaussian approximation based on the normalized likelihood function is more suitable for highly non-Gaussian likelihood functions, because the mean and variance are propagated more precisely compared to the Gaussian approximation based on the ML estimate.

The paper is organized as follows. In Section 2, we introduce the extended latent Gaussian modeling framework, and in Section 3 we detail our proposed approximate Bayesian inference methodology and introduce Max-and-Smooth. In Section 4 we illustrate the strengths and weaknesses of our approach by simulation studies using different types of extended LGMs. We apply the proposed methodology to a real dataset in Section 5. Finally, Section 6 concludes with a discussion and directions for future research.

2 Extended latent Gaussian models

2.1 LGMs with a univariate link function

Latent Gaussian models are a subset of Bayesian hierarchical models in which parameters at the latent level have a joint Gaussian prior distribution, conditional on hyperparameters. LGMs are a subclass of structured additive regression models where the observations, $y_{i}$ ($i=1,\dots,n$), are assumed to have a density from the exponential family and the mean or a particular quantile, $\mu_{i}$, is then linked to a structured additive predictor, $\eta_{i}$, through a univariate link function $g(\cdot)$ such that $g(\mu_{i})=\eta_{i}$; see Rue et al. (2009). The structured additive predictor can then accommodate covariates and random effects in an additive way, namely,

[TABLE]

where $\beta_{0}$ is an intercept, $\{\beta_{k}\}$ are linear fixed effects of covariates $z_{i,1},\ldots,z_{i,K}$, $\{u_{j}\}$ are unknown random effects with some specified dependency structure, and with known weights $a_{i,1},\ldots,a_{i,J}$, and $\epsilon_{1},\ldots,\epsilon_{n}$ are unstructured model errors. In LGMs, the terms $\beta_{0}$, $\{\beta_{k}\}$, $\{u_{j}\}$ and $\{\epsilon_{i}\}$ all have Gaussian prior distributions. Let $\boldsymbol{x}$ contain the latent parameters, namely, $\beta_{0}$, $\{\beta_{k}\}$, $\{u_{j}\}$ and $\{\epsilon_{i}\}$. Sometimes, the vector $\boldsymbol{x}$ consists of $\beta_{0}$, $\{\beta_{k}\}$, $\{u_{j}\}$ and $\{\eta_{i}\}$, i.e., $\{\eta_{i}\}$ is included instead of $\{\epsilon_{i}\}$. Either way, the parameters in $\boldsymbol{x}$ have a joint Gaussian prior distribution, conditional on hyperparameters, $\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{p})^{\mathsf{T}}$. The hyperparameters usually do not have a Gaussian prior density. Typically, hyperparameters specify the marginal variance, correlation range, smoothness, and/or other correlation parameters of the random effects. Schematically, LGMs with a univariate link function may be represented hierarchically as follows in terms of the data level, the latent level and the hyperparameter level:

Data level:

The observations $\boldsymbol{y}=(y_{1},\ldots,y_{n})^{\mathsf{T}}$ are assumed to be dependent on the latent parameters $\boldsymbol{x}$ and have a density $\pi(\boldsymbol{y}|\boldsymbol{x},\boldsymbol{\theta})$. Often, conditional independence is assumed for simplicity, that is, $\pi(\boldsymbol{y}|\boldsymbol{x},\boldsymbol{\theta})$ factorizes as $\prod_{i}\pi(y_{i}|\boldsymbol{x},\boldsymbol{\theta})$ where $\pi(y_{i}|\boldsymbol{x},\boldsymbol{\theta})=\pi(y_{i}|\eta_{i},\boldsymbol{\theta})$ and $\mu_{i}=\textrm{E}(y_{i}|\eta_{i})=g^{-1}(\eta_{i})$ or $\mu_{i}=\textrm{Q}_{p}(y_{i}|\eta_{i})=g^{-1}(\eta_{i})$ if $\mu_{i}$ is a quantile defined in terms of an appropriate quantile function $\textrm{Q}_{p}$.

Latent level:

The latent parameters $\boldsymbol{x}$ have a Gaussian prior density and are potentially dependent on some hyperparameters $\boldsymbol{\theta}$. The density of $\boldsymbol{x}$ may be written as

[TABLE]

where the right-hand side denotes a Gaussian density with mean vector $\boldsymbol{\mu}(\boldsymbol{\theta})$ and covariance matrix $\Sigma(\boldsymbol{\theta})$. This is in line with the notation in Gelman et al. (2013).

Hyperparameter level:

The hyperparameters $\boldsymbol{\theta}$ are assigned a prior density $\pi(\boldsymbol{\theta})$.

An LGM is fully specified by the definition of these three levels. In the next section, we extend this framework to LGMs characterized by a multivariate link function.

2.2 LGMs with a multivariate link function

LGMs with a multivariate link function are referred to as extended LGMs (Geirsson et al., 2020). We assume here conditional independence at the data level for simplicity, and that the data can be lined up according to groups, e.g., sites, time points, spatio-temporal elements or categories. The models presented here have the same structure as the LGMs in Section 2.1, except that the assumption of the data density being in the exponential family is dropped and the vector $\boldsymbol{x}$ refers to several subsets of parameters found at the data level, each subset with its separate set of linear predictors at the latent level. In contrast, in the case of classical LGMs defined in Rue et al. (2009), only one single parameter of the data distribution is modeled at the latent level.

We assume that each group $i$ has $n_{i}$ observations. The total number of groups is $G$. Observations from the same group or from different groups are assumed to be conditionally independent given the latent process. The groups can represent various types of sampling setups. For example, the groups may be geological sites observed over time, i.e., each group corresponds to a site; or the groups may be time points where several observations are made at the same time point. The groups may also be spatio-temporal elements such that multiple observations are collected for each spatially and temporally-referred element. Furthermore, the groups may also represent generic categories that do not have any spatial nor temporal reference, yet, several observations are made within each category. Each group is described by $M$ parameters. The general setup is such that the $M$ subsets of parameters at the data level are mapped to $M$ subsets of parameters at the latent level through an $M$-variate link function. Each of these subsets at the latent level is modeled with a linear model of the form in (1).

The probability density function of $y_{i,j}$, the $j$-th observation from group $i$, is denoted by $\pi(y_{i,j}|\psi_{1,i},\psi_{2,i},\ldots,\psi_{M,i})$ where $\psi_{1,i},\ldots,\psi_{M,i}$ are the parameters within group $i$ such that $(\psi_{1,i},\ldots,\psi_{M,i})\in\mathcal{D}$, and $\mathcal{D}$ is a subspace of $\mathbb{R}^{M}$. Let $\boldsymbol{y}$ be the vector containing all the observations, let $\boldsymbol{y}_{i}$ be the observations from group $i$, and let

[TABLE]

denote the $M$ subsets of parameters at the data level, each vector containing only one type of parameters, e.g., all the location parameters or the regression slope coefficients for all groups. Conditional on $\boldsymbol{\psi}_{1}$, $\boldsymbol{\psi}_{2}$, …, $\boldsymbol{\psi}_{M}$, the probability density function of $\boldsymbol{y}$ is

[TABLE]

where $\mathbb{A}_{i}$ is an index set for group $i$. Let $g$ be an $M$-variate link function such that $g:\mathcal{D}\rightarrow\mathbb{R}^{M}$ with $g(\psi_{1,i},\psi_{2,i},\ldots,\psi_{M,i})=(\eta_{1,i},\eta_{2,i},\ldots,\eta_{M,i})\in\mathbb{R}^{M}$, so the domain of each $\eta_{m,i}$ is the whole real line. Linear models for the $M$ subsets of parameters at the latent level, i.e., the vectors

[TABLE]

are then specified and they may be expressed in vector notation as

[TABLE]

where $\boldsymbol{\beta}_{1}$, $\boldsymbol{\beta}_{2}$, …, $\boldsymbol{\beta}_{M}$ are fixed effects, $X_{1}$, $X_{2}$, …, $X_{M}$ are the corresponding design matrices containing covariate information, $\boldsymbol{u}_{1}$, $\boldsymbol{u}_{2}$, …, $\boldsymbol{u}_{M}$ are random effects, $A_{1}$, $A_{2}$, …, $A_{M}$ are their corresponding weight matrices, and $\boldsymbol{\epsilon}_{1}$, $\boldsymbol{\epsilon}_{2}$, …, $\boldsymbol{\epsilon}_{M}$ are independent and unstructured error terms, referred to as model errors. The terms $\boldsymbol{\beta}_{m}$, $\boldsymbol{u}_{m}$, and $\boldsymbol{\epsilon}_{m}$, $m=1,2,\ldots,M$, are assigned Gaussian prior densities and assumed to be a priori mutually independent.

In the following section, we develop an approximate inference scheme for the LGMs with a multivariate link function.

3 Approximate Bayesian inference for extended LGMs

3.1 General idea

In this section, we detail an approximation to the posterior density used to compute inference for extended LGMs (recall Section 2.2). We then introduce Max-and-Smooth, a two-step approach that is fully Bayesian and utilizes a Gaussian approximation to the likelihood. This approximation together with the conditional Gaussian prior at the latent level results in a Gaussian–Gaussian model (conditional on the hyperparameters) that is referred to here as the pseudo model. We perform inference both for the hyperparameters and for the latent parameters in the original extended LGM by exploiting this Gaussian–Gaussian model. In Section 3.2, we describe the posterior density of extended LGMs, and in Section 3.3 we detail the two consecutive steps of our inference approach. The first step relies on a Gaussian approximation to the likelihood function. We actually propose two different approximations; one is based on the ML estimates and the inverse of the observed information, while the other one is based on the mean and covariance of the normalized likelihood function. The second step consists in inferring the resulting Gaussian–Gaussian pseudo model. For more details about the approximate inference scheme, see the Supplementary Material.

3.2 The posterior density of extended LGMs

The vectors for the model in (2) are gathered as follows

[TABLE]

[TABLE]

[TABLE]

A priori, the vectors $\boldsymbol{\beta}_{1}$, $\boldsymbol{u}_{1}$, $\boldsymbol{\beta}_{2}$, $\boldsymbol{u}_{2}$,…, $\boldsymbol{\beta}_{M}$ and $\boldsymbol{u}_{M}$ are assumed to be independent. Denote the means and precision (i.e., inverse covariance) matrices of these vectors by $\boldsymbol{\mu}_{\beta,1}$, $\boldsymbol{\mu}_{u,1}$, $\boldsymbol{\mu}_{\beta,2}$, $\boldsymbol{\mu}_{u,2}$,…, $\boldsymbol{\mu}_{\beta,M}$, $\boldsymbol{\mu}_{u,M}$, $Q_{\beta,1}$, $Q_{u,1}$, $Q_{\beta,2}$, $Q_{u,2}$, …, $Q_{\beta,M}$ and $Q_{u,M}$, respectively. The prior mean of $\boldsymbol{\nu}$ is therefore

[TABLE]

while the precision matrix of $\boldsymbol{\nu}$ is a block diagonal matrix,

[TABLE]

The precision matrices of $\boldsymbol{\epsilon}_{1}$, $\boldsymbol{\epsilon}_{2}$, …, $\boldsymbol{\epsilon}_{M}$ are diagonal matrices that are denoted by $Q_{\epsilon,1}=\sigma_{\epsilon,1}^{-2}I$, $Q_{\epsilon,2}=\sigma_{\epsilon,2}^{-2}I$, …, $Q_{\epsilon,M}=\sigma_{\epsilon,M}^{-2}I$, respectively, where $I$ is the identity matrix and $\sigma_{\epsilon,1}$, $\sigma_{\epsilon,2}$, …, $\sigma_{\epsilon,M}$ are the corresponding standard deviations. The precision matrix of $\boldsymbol{\epsilon}$ is thus given by

[TABLE]

Define the matrix $Z$ based on $X_{1}$, $A_{1}$, $X_{2}$, $A_{2}$,…, $X_{M}$ and $A_{M}$ as

[TABLE]

where the zeros denote zero matrices. Simple matrix multiplication implies that (2) can be rewritten more compactly as

[TABLE]

Since the model is an extended LGM, the prior densities of $\boldsymbol{\nu}$, and $\boldsymbol{\eta}$ conditional on $\boldsymbol{\nu}$, are assumed to be Gaussian. Furthermore, the precision matrices $Q_{\beta,1}$, $Q_{\beta,2}$, …, $Q_{\beta,M}$ are assumed to be fixed, while hyperparameters govern the precision matrices for the random effects $\boldsymbol{u}_{1}$, $\boldsymbol{u}_{2}$, …, $\boldsymbol{u}_{M}$, i.e., $Q_{u,1}$, $Q_{u,2}$, …, $Q_{u,M}$. When the dimensions of $\boldsymbol{u}_{1}$, $\boldsymbol{u}_{2}$, …, $\boldsymbol{u}_{M}$ are large, and we require fast computation, then $Q_{u,1}$, $Q_{u,2}$, …, $Q_{u,M}$ need to be sparse and thus specifying them with Gaussian Markov random fields (GMRF, Rue and Held, 2005) becomes crucial.

The joint posterior density of $(\boldsymbol{\eta},\boldsymbol{\nu},\boldsymbol{\theta})$ may be expressed as

[TABLE]

where $\boldsymbol{y}$ denotes the data vector, $\pi(\boldsymbol{y}|\boldsymbol{\eta})$ is the data density defined at the data level, $\pi(\boldsymbol{\eta},\boldsymbol{\nu}|\boldsymbol{\theta})$ is the Gaussian prior density defined at the latent level, and $\pi(\boldsymbol{\theta})$ is the prior density for the hyperparameters defined at the hyperparameter level. When the data density $\pi(\boldsymbol{y}|\boldsymbol{\eta})$ is used for inference, it is referred to as the likelihood function and viewed as a function of $\boldsymbol{\eta}$. We next describe our proposed approximate inference scheme.

3.3 Max-and-Smooth: a two-step approximate inference approach

Our proposed approximate Bayesian inference scheme (Max-and-Smooth) is based on approximating the likelihood function with a Gaussian density function (Step 1, Max), and on fitting the resulting Gaussian–Gaussian pseudo model (Step 2, Smooth). We now describe each step separately.

3.3.1 Step 1 (Max): Gaussian approximation of the likelihood function

We here propose two different Gaussian likelihood approximations that we then subsequently exploit for fast fully Bayesian inference.

The first Gaussian approximation is based on the mode of the likelihood function, i.e., the ML estimate (hence the term “Max”), and the observed information evaluated at the ML estimate. Let $L(\boldsymbol{\eta}|\boldsymbol{y})$ denote the likelihood function, where $L(\boldsymbol{\eta}|\boldsymbol{y})=\pi(\boldsymbol{y}|\boldsymbol{\eta})$, and let $\hat{L}$ denote the first Gaussian approximation; then, $c\hat{L}(\boldsymbol{\eta}|\boldsymbol{y})\approx L(\boldsymbol{\eta}|\boldsymbol{y})$, where

[TABLE]

$c$ is a constant independent of $\boldsymbol{\eta}$, $\hat{\boldsymbol{\eta}}$ is the ML estimate for $\boldsymbol{\eta}$, i.e., it is the mode of $L(\boldsymbol{\eta}|\boldsymbol{y})$, and $\Sigma_{\eta y}=(-H_{\eta y})^{-1}$, where $H_{\eta y}$ denotes the Hessian matrix of $\log(L(\boldsymbol{\eta}|\boldsymbol{y}))$ evaluated at $\boldsymbol{\eta}=\hat{\boldsymbol{\eta}}$. Furthermore, the observed information evaluated at $\hat{\boldsymbol{\eta}}$ is written as $\mathcal{I}_{\eta y}=-H_{\eta y}$.

Due to the assumed conditional independence, the first Gaussian approximation is straightforward to evaluate because it can be computed for each group separately. More precisely, let $L(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})=\pi(\boldsymbol{y}_{i}|\boldsymbol{\eta}_{i})$ denote the likelihood contribution of the $i$-th group such that $L(\boldsymbol{\eta}|\boldsymbol{y})=\prod_{i=1}^{G}L(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})$. The ML estimate of $\boldsymbol{\eta}_{i}$, $\hat{\boldsymbol{\eta}}_{i}$, is the mode of $L(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})$, and the inverse of the covariance matrix is equal to the observed information in $L(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})$ evaluated at $\hat{\boldsymbol{\eta}}_{i}$, i.e.,

[TABLE]

and the inverse of $\mathcal{I}_{\eta y_{i}}$ is equal to the covariance matrix $\Sigma_{\eta y_{i}}$. Now we can approximate the likelihood contribution of group $i$, $L(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})$, with $c_{i}\hat{L}(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})$, where

[TABLE]

and $c_{i}$ is a constant independent of $\boldsymbol{\eta}_{i}$.

Therefore, the approximated posterior density $\hat{\pi}(\boldsymbol{\eta},\boldsymbol{\nu},\boldsymbol{\theta}|\boldsymbol{y})$ based on the approximated full likelihood $\hat{L}(\boldsymbol{\eta}|\boldsymbol{y})=\prod_{i=1}^{G}\hat{L}(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})$, is such that $\hat{\pi}(\boldsymbol{\eta},\boldsymbol{\nu},\boldsymbol{\theta}|\boldsymbol{y})\approx\pi(\boldsymbol{\eta},\boldsymbol{\nu},\boldsymbol{\theta}|\boldsymbol{y})$, and it is given by

[TABLE]

The second Gaussian approximation relies on normalizing the likelihood function such that the function $d^{-1}L(\boldsymbol{\eta}|\boldsymbol{y})$ integrates to one (over the domain of $\boldsymbol{\eta}$), where $d\in(0,\infty)$ is an appropriate normalization constant that is independent of $\boldsymbol{\eta}$. Here we assume that $d$ is finite. If this is not the case, we may either find a more adequate model parametrization or replace the likelihood function by an alternative generalized likelihood that consists of the likelihood times an extra prior density for $\boldsymbol{\eta}$. Bayes’ Theorem ensures the finiteness of $d$ under the generalized likelihood since it is proportional to a posterior density. This second Gaussian approximation is also designed to “maximize” the match with the true likelihood function, especially in skewed scenarios.

Then, similarly to the first Gaussian approximation, the likelihood function, or the generalized likelihood function, is approximated with a Gaussian density that has mean and covariance matrix equal to those of the normalized likelihood function. If there is a need for ensuring that the mean and variance are finite then the extra prior density can be given a finite support. We exploit again the assumed conditional independence, i.e., $L(\boldsymbol{\eta}|\boldsymbol{y})=\prod_{i=1}^{G}L(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})$, now approximating the $i$-th likelihood contribution as

[TABLE]

where $d_{i}$ is a constant independent of $\boldsymbol{\eta}_{i}$, $\tilde{L}(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})$ is the Gaussian approximation, and the mean and the covariance matrix are

[TABLE]

[TABLE]

Similarly to (3), the alternative approximation to the posterior density, $\tilde{\pi}(\boldsymbol{\eta},\boldsymbol{\nu},\boldsymbol{\theta}|\boldsymbol{y})$, which is based on $\tilde{L}(\boldsymbol{\eta}|\boldsymbol{y})=\prod_{i=1}^{G}\tilde{L}(\boldsymbol{\eta}_{i}|\boldsymbol{y}_{i})$, may be expressed as

[TABLE]

Because the Gaussian density is the asymptotic form of the likelihood function under mild regularity conditions (Schervish, 1995), these two types of approximations are expected to work increasingly well when the number of data replicates grows. With a low number of replicates (i.e., less than $10$–$20$ per distinct model parameter involved at the data level and further described at the latent level), a small bias might be expected, although we have found it to be relatively negligible in the settings we have considered. More details on the quality of these Gaussian approximations are given in Section 4.

The computational benefit of the approximations in (3) and (4) lies in the fact that $\pi(\boldsymbol{\eta},\boldsymbol{\nu}|\boldsymbol{\theta})$ is Gaussian with respect to $(\boldsymbol{\eta},\boldsymbol{\nu})$ and the functional form of both $\hat{L}(\boldsymbol{\eta}|\boldsymbol{y})$ and $\tilde{L}(\boldsymbol{\eta}|\boldsymbol{y})$ with respect to $\boldsymbol{\eta}$ is proportional to a Gaussian density. As a result, the conditional posterior density of $(\boldsymbol{\eta},\boldsymbol{\nu})$ is Gaussian, and posterior samples can be obtained directly from this density and it is well known how to generate the samples from it. The information about $\boldsymbol{\eta}$ stemming from the data is quantified with reasonable accuracy in $\hat{L}(\boldsymbol{\eta}|\boldsymbol{y})$ or $\tilde{L}(\boldsymbol{\eta}|\boldsymbol{y})$ provided that at least one of these two approximations is fairly good. This information is correctly weighted against the prior information about $\boldsymbol{\eta}$ which is quantified in $\pi(\boldsymbol{\eta},\boldsymbol{\nu}|\boldsymbol{\theta})$. Since the inference scheme is Bayesian and the parameters are inferred simultaneously, then the information about $\boldsymbol{\eta}$ is correctly propagated to $\boldsymbol{\nu}$ and $\boldsymbol{\theta}$ through $\pi(\boldsymbol{\eta},\boldsymbol{\nu}|\boldsymbol{\theta})$. Notice that the likelihood approximation $\hat{L}(\boldsymbol{\eta}|\boldsymbol{y})$ may be faster to compute than $\tilde{L}(\boldsymbol{\eta}|\boldsymbol{y})$ as it is free from integrals, while $\tilde{L}(\boldsymbol{\eta}|\boldsymbol{y})$ is likely to provide a more accurate approximation in the case of a small number of replicates within each group. Hereafter, the Gaussian approximation based on the ML estimates and the inverse of the observed information will be referred to as the first Gaussian approximation, while the Gaussian approximation based on the mean and covariance of the normalized likelihood function will be referred to as the second Gaussian approximation.

3.3.2 Step 2 (Smooth): Inference for the pseudo Gaussian–Gaussian model

To infer the model presented in Section 2.2 based on the approximate posterior density in (3), we consider a pseudo model that is such that $\hat{\boldsymbol{\eta}}$ (obtained from the first Gaussian approximation) is treated as noisy measurements of the latent field. Fitting this pseudo model is equivalent to smoothing the parameters $\hat{\boldsymbol{\eta}}$ jointly (hence the term “Smooth”). A similar approach may be used based on (4) by treating $\tilde{\boldsymbol{\eta}}$ (obtained from the second Gaussian approximation) as the data. The proposed data density of the pseudo model based on (3) is $\textrm{N}(\hat{\boldsymbol{\eta}}|\boldsymbol{\eta},Q_{\eta y}^{-1})$ where $Q_{\eta y}=\Sigma_{\eta y}^{-1}$, and $Q_{\eta y}$ is known. Its numerical values are evaluated from the already observed data. This model can be written hierarchically as

[TABLE]

and $\pi(\boldsymbol{\theta})$ is the prior density for $\boldsymbol{\theta}$ as before. The posterior density for this model is given by

[TABLE]

The above posterior density stems from looking at it as a function of $\boldsymbol{\eta}$ and taking $\hat{\boldsymbol{\eta}}$ as a fixed quantity, which gives $\textrm{N}(\boldsymbol{\hat{\eta}}|\boldsymbol{\eta},Q^{-1}_{\eta y})=\textrm{N}(\boldsymbol{\eta}|\boldsymbol{\hat{\eta}},Q^{-1}_{\eta y})$. Thus, the above posterior density is exactly the same as the approximated posterior density in (3) for the extend LGM in Section 2.2. The pseudo model is a Gaussian–Gaussian model and it is convenient to approach the inference for the unknown parameters through this model.

Samples of $(\boldsymbol{x},\boldsymbol{\theta})$, where $\boldsymbol{x}=(\boldsymbol{\eta}^{\mathsf{T}},\boldsymbol{\nu}^{\mathsf{T}})^{\mathsf{T}}$, are obtained by sampling first from the marginal posterior density of $\boldsymbol{\theta}$, and then from the posterior density of $\boldsymbol{x}$ conditional on $\boldsymbol{\theta}$. The marginal posterior density of $\boldsymbol{\theta}$ given $\hat{\boldsymbol{\eta}}$ is $\pi(\boldsymbol{\theta}|\hat{\boldsymbol{\eta}})\propto\pi(\boldsymbol{\theta})\pi(\hat{\boldsymbol{\eta}}|\boldsymbol{\theta})$ and it can be represented as

[TABLE]

The densities $\pi(\hat{\boldsymbol{\eta}}|\boldsymbol{x},\boldsymbol{\theta})$ and $\pi(\boldsymbol{x}|\boldsymbol{\theta})$ have precision matrices $Q_{\eta y}$ and $Q_{x}$, respectively, and if $Q_{\eta y}$ and $Q_{x}$ are sparse matrices, then the precision matrix of $\pi(\boldsymbol{x}|\hat{\boldsymbol{\eta}},\boldsymbol{\theta})$ is a sparse matrix; see details in the Supplementary Material. Samples from the marginal posterior density of $\boldsymbol{\theta}$ can be obtained by using grid sampling if the dimension of $\boldsymbol{\theta}$ is small (i.e., four or less), a Metropolis step or a Metropolis–Hastings step, or other samplers that are well suited for densities with non-tractable form. Since the conditional posterior density of $\boldsymbol{x}$ is Gaussian, samples of $\boldsymbol{x}$ are straightforward to obtain, and if the precision matrix of $\pi(\boldsymbol{x}|\hat{\boldsymbol{\eta}},\boldsymbol{\theta})$ is sparse, the computational cost is relatively low.

Further methodological and computational details on the inference scheme for the Gaussian–Gaussian pseudo model are presented in the Supplementary Material.

4 Simulation examples

4.1 Settings

In this section, we assess the accuracy of our proposed approximate Bayesian inference scheme, Max-and-Smooth, by evaluating by simulation how close the approximate posterior density of an LGM is to its exact posterior density. In particular, in Section 4.2 we consider an LGM where the data are independent mean zero Gaussian random variables on a lattice, with spatially-varying variance at each lattice point. The approximate marginal posterior densities of the latent parameters and the hyperparameters are compared to the exact posterior densities inferred by an “exact” MCMC sampler. Moreover, as the link function in this specific example is univariate, INLA can also be applied and we include it in our experiments for comparison.

In the Supplementary Material we explore three other models, namely a linear regression model on a lattice with Gaussian error terms, which is also applied to real data in Section 5, a linear regression model with $t$-distributed error terms and temporally varying coefficients, and a spatio-temporal model for Poisson counts. Our results suggest that the Gaussian likelihood approximation is accurate in finite sample sizes (even with just one single replicate in the case of the Poisson distribution with large counts), and therefore, that our Max-and-Smooth approach performs well in a rich variety of realistic settings.

4.2 Gaussian data with spatially varying log-variance

In this section, we apply Max-and-Smooth to mean-zero Gaussian data with spatially-varying variance. We first simulate a single realization of a Gaussian Markov random field (Rue and Held, 2005), $\{x_{i}\}$, on a regular lattice of size $10\times 10$ where the index $i$ corresponds to the lattice point with horizontal coordinates $i_{1}$ and vertical coordinates $i_{2}$. See further details in the Supplementary Material. With the mean fixed to zero, the conditional mean and variance of this GMRF are

[TABLE]

For the simulation, the precision parameter of the GMRF is fixed at $\tau=1$. At each of the $N=100$ lattice points, we simulate $T=10,20,50$ independent Gaussian variates with zero mean and log-variance $x_{i}$, i.e., $y_{i,t}\sim\textrm{N}(0,\exp(x_{i}))$ for $t=1,\dots,T$ and $i=1,\dots,N$. Therefore, in this simulation example, the “groups” mentioned in Section 2.2 represent the data replicates available at each lattice point, and thus, the total number of groups is equal to the number of lattice points, i.e., $G=N=100$. The goal is to exploit Max-and-Smooth to infer the latent variables $\boldsymbol{x}=\{x_{i}\}$, and the precision hyperparameter $\tau$ from the observed data $\boldsymbol{y}=\{y_{i,t}\}$. The model is simple enough to also infer the latent variables and hyperparameter using an MCMC sampler that uses the true likelihood function, to compare the approximation with an “exact” fully Bayesian procedure. The model parameters are inferred assuming that the mean at the data level is equal to zero. As the link function is here univariate (i.e., $M=1$), INLA is also used to infer the model parameters for the sake of comparison (but note that this would not be possible for more complex models with multivariate link functions). The full details of this simulation study are reported in the Supplementary Material.

Figure 1 shows the posterior densities of latent parameters $x_{i}$ at three different lattice points, inferred from datasets with $T=10,20,50$ replicates per location. We chose the three locations with the smallest, closest to zero, and largest element of $\boldsymbol{x}$. Each graph in Figure 1 shows four posterior densities, based on a fully Bayesian MCMC simulation, based on the two approximate inference schemes from Section 3.3, and based on INLA. By comparing the density curves, we see that the second Gaussian approximation (with mean and variance derived from the normalized likelihood function) is generally closer to the true posterior density than the first Gaussian approximation (with mean and variance based on the MLE and observed information). Both approximate posterior inference schemes capture the exact posterior densities very well when $T$ is sufficiently large (greater than $20$). In the case of $T=10$, the discrepancy between the approximated and exact posterior density is relatively small for the small and median values of $x_{i}$, while it is more pronounced in the case of large $x_{i}$. The reason might be that the Gaussian approximation of the likelihood function does not properly capture the right-skewness of the exact likelihood function. When the spatial prior for the vector $\boldsymbol{x}$ has a strong effect (pulling the estimate toward zero), then the right skewness of the marginal likelihood function will show up most prominently in the case of large $x_{i}$. The posterior densities estimated by INLA are slightly closer to the exact densities for these selected latent variables than our approximations.

The exact and approximate posterior densities of the precision hyperparameter $\tau$ of this model are shown in Figure 2 for sample sizes $T=10,20,50$. Figure 2 shows that when $T\geq 20$ there is a negligible difference between the densities approximated by Max-and-Smooth and the exact posterior densities. When $T=10$ the difference between the two Max-and-Smooth approximations is small and their difference with respect to the exact density appears reasonably small. For $T=10$ there is no notable difference between the quality of INLA and Max-and-Smooth, but for $T\geq 20$, the posterior density calculated by INLA is further from the exact posterior than both Max-and-Smooth approximations.

In addition to the accuracy of the approximation, we also compare the computational speed of each approach. To calculate the joint density $\pi(\boldsymbol{y},\boldsymbol{x},\tau)=\pi(\boldsymbol{y}|\boldsymbol{x},\tau)\pi(\boldsymbol{x}|\tau)\pi(\tau)$, the “exact” Bayesian inference requires $T$ evaluations of the Gaussian density per lattice point to calculate the likelihood $\pi(\boldsymbol{y}|\boldsymbol{x},\tau)$. In contrast, under the approximate inference schemes, calculating the approximation of the likelihood requires only a single evaluation of the Gaussian density per lattice point. We should thus expect the computational cost of the exact inference to scale linearly with $T$, while the computational cost of Max-and-Smooth should be constant in $T$. Figure 3 shows the time to draw 10,000 MCMC samples of $\tau$ and $\boldsymbol{x}$ for different numbers of temporal replicates $T$ and different numbers of grid points $N$. Max-and-Smooth has constant computation cost as a function of $T$, and does not suffer from the same slowdown with increasing $T$ as the exact method. As might be expected, for small sample sizes and grid sizes, the speedup is only moderate. But for large grid sizes such as $N=50\times 50$, and $T=100$ replicates, Max-and-Smooth is faster than the exact method by a factor greater than 10. Since both approximation methods use a Gaussian approximation for the likelihood function, only with a different mean and variance, there is no difference in computational speedup between them. As computation of marginal posteriors with INLA is not based on sampling, it was excluded from the comparison of computational cost in Figure 3. Inferring the marginal posterior densities of the latent field and the hyperparameters took only a few seconds with INLA, but unlike Max-and-Smooth, INLA’s run time increased with increasing number of temporal replicates.

5 Predictions of meteorological variables on a lattice

As alluded to in the introduction, a wide variety of Bayesian or frequentist approaches, involving different types of approximations or simplifications, have been proposed to deal with high-dimensional spatio-temporal data, often under the assumption of data being exactly Gaussian. These include low-rank approaches (Cressie and Johannesson, 2008), the predictive process (Banerjee et al., 2008), covariance tapering (Furrer et al., 2006; Anderes et al., 2013), multi-resolution models (Nychka et al., 2015; Katzfuss, 2017), hierarchical nearest-neighbor Gaussian processes (Datta et al., 2016), the Vecchia approximation (Vecchia, 1988; Stein et al., 2004; Katzfuss and Guinness, 2019), the integrated nested Laplace approximation (Rue et al., 2009; Bakka et al., 2018), or more recently a frequentist approach for data modeled by a generalized-extreme value (GEV) distribution (Risser et al., 2019; Russell et al., 2019). See Heaton et al. (2019) for a recent review and comparison of some of these methods. In this section, we apply our proposed Max-and-Smooth approach to analyze a moderately high-dimensional real climate dataset, assuming that the data are well described with an extended LGM.

More precisely, the dataset used in our analysis is from a seasonal climate forecasting experiment, consisting of retrospective surface temperature forecasts, and their corresponding “verifying” observations. The forecast data were produced by a global climate model ensemble (28 forecasts started from perturbed initial conditions) downloaded from the ECMWF C3S Seasonal catalog (https://apps.ecmwf.int/data-catalogues/c3s-seasonal/) via the MARS API on 21 February 2018. Forecasts were started from perturbed May 1 initial conditions each year from 1993 to 2015, i.e., the sample size is $T=23$ in time. Each model run predicts atmospheric conditions several months into the future. The particular forecast target analyzed for this paper is surface air temperature on a 1-by-1 degree latitude-longitude grid over a rectangular region with corners 20W/40N and 40E/60N ($N=1281$ grid points in total, covering most of Europe). Here each grid point forms a group, so the total number of groups is $G=N$. The forecasts were averaged over the 28 ensemble members, and over the Boreal summer period June/July/August, yielding a single scalar prediction per grid point per year. Since these forecasts were initialized in May and predict June-August climate, the forecast lead time is 1–3 months. Verifying observations are from the ERA-Interim reanalysis dataset (Dee et al., 2011), which is available on the same 1-by-1 degree grid as the forecasts, and can also be downloaded from the ECMWF data base. Observation data were averaged over the same June–August period as the forecasts.

Due to structural errors and missing physical processes in the climate model, and due to the chaotic nature of atmospheric dynamics, numerical model forecasts have systematic biases in the forecast mean. Furthermore, for forecasts of the climate system on seasonal time scales, the correlation between forecasts and verifying observations tends to be low. The biases in the numerical model forecasts can be partly corrected through a linear regression of the observations on forecasts. The adjustment of model forecasts by linear regression, known as model output statistics (MOS), has a long tradition in the weather forecasting community (Glahn and Lowry, 1972; Glahn et al., 2009), and is part of an active area of research known as forecast recalibration or forecast post-processing (e.g., Siegert and Stephenson, 2019). In this section, we will infer the post-processing parameters on a spatial grid in an extended LGM framework, using spatial priors for the regression coefficients to reduce their estimation uncertainty, and ultimately improve the predictive skill of the recalibrated model forecasts.

The statistical model linking observed climate $y_{i,t}$ at time $t$ and location $i$ to the climate model forecast $f_{i,t}$ for the same time and location is assumed to be

[TABLE]

where $\bar{f}_{i}$ denotes the local mean forecast over the data period, and the residuals $\epsilon_{i,t}$ are independent Gaussian variates with mean zero and variance $\exp({\tau_{i}})$. It is common to estimate the regression parameters $\boldsymbol{\eta}=(\boldsymbol{\alpha}^{\mathsf{T}},\boldsymbol{\beta}^{\mathsf{T}},\boldsymbol{\tau}^{\mathsf{T}})^{\mathsf{T}}=(\boldsymbol{\eta}_{\alpha}^{\mathsf{T}},\boldsymbol{\eta}_{\beta}^{\mathsf{T}},\boldsymbol{\eta}_{\tau}^{\mathsf{T}})^{\mathsf{T}}$ individually for each grid point, using maximum likelihood estimation. However, this may lead to high estimation variability, due to the limited number of samples per grid point.

To exploit the spatial structure in the data, and borrow strength from data at neighboring grid points when estimating the regression parameters $\boldsymbol{\eta}$, we use a spatial prior distribution for the spatial fields of regression parameters as outlined in Section 2.2 of the Supplementary Material, and we exploit Max-and-Smooth for Bayesian inference. The spatial field for each regression parameter is decomposed additively into a spatially correlated component $\boldsymbol{u}$ and an unstructured component $\boldsymbol{\epsilon}$, i.e.,

[TABLE]

and respectively for $\boldsymbol{\beta}=\boldsymbol{\eta}_{\beta}$ and $\boldsymbol{\tau}=\boldsymbol{\eta}_{\tau}$. We model the structured term $\boldsymbol{u}$ as a first-order intrinsic Gaussian Markov random field on a regular lattice, (Rue and Held, 2005, Section 3.3.2, pp. 104–108) and the unstructured component $\boldsymbol{\epsilon}$ as a mean zero Gaussian process with diagonal covariance matrix. All latent processes $\boldsymbol{u}$ and $\boldsymbol{\epsilon}$ are mutually independent.

The prior model has a total of six hyperparameters: The precision parameters of the three spatially correlated fields $\boldsymbol{u}_{\alpha}$, $\boldsymbol{u}_{\beta}$ and $\boldsymbol{u}_{\tau}$, denoted $\{\tau_{u,\alpha},\tau_{u,\beta},\tau_{u,\tau}\}$, and the precision parameters of the three unstructured fields $\boldsymbol{\epsilon}_{\alpha}$, $\boldsymbol{\epsilon}_{\beta}$ and $\boldsymbol{\epsilon}_{\tau}$, denoted $\{\tau_{\epsilon,\alpha},\tau_{\epsilon,\beta},\tau_{\epsilon,\tau}\}$. The hyperparameters $\boldsymbol{\theta}=(\tau_{u,\alpha},\tau_{\epsilon,\alpha},\tau_{u,\beta},\tau_{\epsilon,\beta},\tau_{u,\tau},\tau_{\epsilon,\tau})^{\mathsf{T}}$ are a priori independent, and spatially homogeneous. We used independent penalized complexity (PC) priors (Simpson et al., 2017) for all hyperparameters. Specifically, the PC priors were specified as exponential distributions with rate parameter 1 for the standard deviations, which leads to a prior density for the precision parameter $\tau$ proportional to $\tau^{-1.5}\exp(-1/\sqrt{\tau})$.

We begin by exploring the marginal posterior distributions of the precision hyperparameters $\boldsymbol{\theta}$ obtained using Max-and-Smooth (recall Section 3). Preliminary studies suggested that the hyperparameters are nearly uncorrelated under their joint posterior, which allows us to explore their posterior distributions individually. We evaluated each unnormalized marginal posterior at 41 points, that are equidistant on the log-scale and centered around the posterior mode, spanning $\pm 4$ posterior standard deviations (estimated via a Laplace approximation from the numerical second derivative).

The (normalized) marginal posteriors of the hyperparameters are shown in Figure 4. For better interpretability, the hyperparameters were transformed from precision $\tau_{v,l}$ to standard deviation $\sigma_{v,l}=\tfrac{1}{\sqrt{\tau_{v,l}}}$ for $v\in\{u,\epsilon\}$ and $l\in\{\alpha,\beta,\tau\}$. The standard deviations of the unstructured components $\sigma_{\epsilon,l}$ are small compared to the spatial variability seen in the ML estimates of the regression coefficients, and compared to the standard deviation parameters of the spatial components $\sigma_{u,l}$. This suggests that a simpler model might be fitted to the data that does not include an unstructured component. The standard deviation corresponding to the spatial effect of the intercept, $\sigma_{u,\alpha}$, is large compared to $\sigma_{u,\beta}$ and $\sigma_{u,\tau}$. Thus, there is a greater spatial variability in the intercept than in the slope and log-variance. By comparing the posterior densities with their PC prior distributions (not shown), we found no evidence that the posterior is unduly influenced by the prior, indicating that the amount of data available is sufficient to properly constrain the model’s hyperparameters.

Figure 5 compares the posterior means of the latent fields $\boldsymbol{\eta}_{l}=\boldsymbol{u}_{l}+\boldsymbol{\epsilon}_{l}$ for $l\in\{\alpha,\beta,\tau\}$ (conditional on the posterior mode of the hyperparameters) with the corresponding maximum likelihood estimates. We used the first Gaussian approximation (ML estimates and observed information) for Figure 5, noting that the corresponding plots for the second Gaussian approximation (mean and variance of normalized likelihood) are almost indistinguishable. The intercept parameters of the regression model in the vector $\boldsymbol{\alpha}$ are not smoothed very much compared to the corresponding ML estimates. This is because the maximum likelihood estimates of $\boldsymbol{\alpha}$ are generally well-constrained by the data, as indicated by the low sampling uncertainty derived from the observed information matrix. In contrast, $\boldsymbol{\beta}$ and $\boldsymbol{\tau}$ exhibit considerable smoothing. The pointwise posterior standard deviations of the latent fields (not shown) are smaller than the sampling standard deviations of the maximum likelihood estimates: they are on average about 5% smaller for $\alpha$ and around 60% smaller for $\beta$ and $\tau$.

We then test the performance of the spatial regression model using cross-validated predictions of surface temperatures calculated using equation (6), integrated over the posterior distributions of $\boldsymbol{\alpha}$, $\boldsymbol{\beta}$ and $\boldsymbol{\tau}$. In practice, forecast and observation data are typically available at all locations in the historical dataset, but no data are available for the present year, for which bias-corrected forecasts are required. We mimic this setting by applying "leave-one-year-out" cross-validation: We leave out all forecasts and observations from one year, fitting the regression parameters using data from the remaining years, and then using the fitted parameters to forecast the temperature at all grid points for the left-out year. For each year and grid point, we generated $N_{s}=1000$ samples from the posterior predictive distribution by repeating the following algorithm $N_{s}$ times:

• •

Sample hyperparameters $\boldsymbol{\theta}^{(s)}$ independently from their marginal posteriors, i.e., $\pi(\tau_{u,l}|\boldsymbol{y}_{-t})$, and $\pi(\tau_{\epsilon,l}|\boldsymbol{y}_{-t})$, where $l\in\{\alpha,\beta,\tau\}$ and $\boldsymbol{y}_{-t}$ is forecast and observation data with year $t$ being left out of the training dataset. Sampling of hyperparameters is done by grid sampling, i.e., after evaluating the marginal density of each hyperparameter at 21 equally spaced points, one of the 21 values is sampled with probability proportional to the marginal density.

• •

Draw an independent sample $\boldsymbol{\eta}^{(s)}$ of each of the latent fields $\boldsymbol{\eta}_{\alpha}$, $\boldsymbol{\eta}_{\beta}$ and $\boldsymbol{\eta}_{\tau}$ from their conditional posterior distribution $\pi(\boldsymbol{\eta}|\hat{\boldsymbol{\eta}},\boldsymbol{\theta}^{(s)})$.

• •

Using the sampled regression parameter in $\boldsymbol{\eta}^{(s)}$, simulate a spatial field of responses $\boldsymbol{y}_{t}$ from the regression model conditional on $\boldsymbol{\eta}^{(s)}_{\alpha}$, $\boldsymbol{\eta}^{(s)}_{\beta}$, and $\boldsymbol{\eta}^{(s)}_{\tau}$, and using the covariates $\boldsymbol{f}_{t}$ from the left-out year $t$.

The procedure is repeated leaving each year out in turn, ultimately resulting in $N_{s}$ out-of-sample posterior predictive samples per grid point for each year. The predictive samples are compared to the verifying observations by the mean squared error (MSE) of the posterior predictive mean, by the continuous ranked probability score (CRPS, see Winkler (1969); Gneiting and Raftery (2007)) of the posterior predictive distribution, by the average widths of the central 95% posterior prediction interval (W95), and by the average coverage frequencies of the posterior predictive $5$, $50$ and $95$-percentiles (COV05, COV50, COV95, respectively). The CRPS is preferred here because it can be evaluated solely from samples drawn from the posterior predictive distribution, without having to know the predictive distribution in closed form, using the following relation:

[TABLE]

where $x_{1},\dots,x_{N_{s}}$ are samples from the posterior predictive distribution, and $y$ is the corresponding observed temperature value (see Gneiting and Raftery, 2007, Equation (20)). Lower CRPS values indicate better forecasts. The different measures of sharpness, reliability and accuracy of the posterior predictive distribution are compared to two benchmark predictions. We draw the same number of samples from the classical predictive distributions using the pointwise ML estimates for the linear regression models. This forecast scheme is denoted MLE and is used to characterize the performance of the forecast post-processing method achieved when spatial correlation in the regression parameters is ignored. An even simpler benchmark model is constructed by sampling $1000$ times from a Gaussian distribution which has the mean and standard deviation (taken over time) of the observations $\boldsymbol{y}_{i}$, calculated separately at each grid point $i$. This benchmark prediction is almost constant between the years, with slight differences only due to the leave-one-out procedure. This forecast scheme is denoted CLIM (for climatology) and quantifies the average predictability of atmospheric surface temperature if no further forecast information is available. The prediction scheme CLIM is thus used to quantify the merit of the forecast information available from the numerical model. Finally, we use the acronyms SPAT1 and SPAT2 to denote the posterior predictive distribution derived from the extended LGM with a spatial prior, inferred using Max-and-Smooth based on the first and second Gaussian approximations, respectively. Due to computational limitations we have not implemented an exact Bayesian inference using an MCMC sampler.

Table 1 shows that spatially smoothing the regression estimates using a hierarchical prior leads to improvements in mean squared error of the predictive mean, and CRPS of the predictive distribution. Simultaneously, forecast uncertainty, as indicated by the narrower average prediction intervals, is reduced by our method. Coverage frequencies are quite good overall.

The seemingly small improvement in MSE and CRPS of SPAT1/2 compared to MLE may suggest a rather negligible improvement of our approach. However, the magnitude of this improvement has to be compared to the improvement of MLE over CLIM. Informally speaking, the improvement of MLE over CLIM quantifies the advantage of a forecaster who has access to the numerical model forecast compared to a forecaster who only knows the climatological distribution. Due to the inherent unpredictability of the atmosphere on long time scales, this advantage tends to be small in seasonal climate prediction. The improvement of SPAT1/2 over MLE thus has to be judged relative to the improvement of MLE over CLIM. The improvement of SPAT1 and SPAT2 over MLE relative to the improvement of MLE over CLIM is about 20% for the MSE, and about 25% for the CRPS. Both improvements are substantial.

To further assess statistical significance of the improvement of SPAT1 over MLE, we approximated the sampling distribution of the average CRPS difference between the MLE and SPAT1 predictions using block bootstrapping. We treat all data as independent in time because the time increment between successive samples is one year, long enough to ignore the temporal correlation in small scale meteorological settings. To account for spatial correlation, we divided the spatial domain into a number of $S$ rectangular non-overlapping blocks. This yields a total of $23\times S$ blocks. We resample these blocks $23\times S$ times with replacement, and then create a bootstrap sample of the average CRPS difference. The bootstrap distribution was estimated from 500 replicates. When using $S=3$ spatial blocks (i.e. assuming only 3 effective spatial degrees of freedom) we obtain a bootstrap mean of $0.0025$ and boostrap standard deviation of $0.0013$ for the CRPS difference, and a bootstrap p-value of $0.032$. For $S=10$ spatial degrees of freedom, we obtain a bootstrap mean of $0.0025$, bootstrap standard deviation of $0.0012$, and a bootstrap p-value of $0.006$. For $S>10$ the bootstrap p-value is essentially zero. Based on visual inspection of maps of CRPS differences, a conservative estimate of the spatial degrees of freedom is at least 20, so the improvement of SPAT over MLE, albeit small, can thus be considered statistically significant with high confidence.

Finally, it should be noted that evaluating the marginal posterior distributions of the 6 hyperparameters at 21 values each, drawing 1000 samples of the spatial fields of the regression parameters from their posterior distributions, and calculating 1000 posterior predictive samples of atmospheric temperature over the entire grid takes less than 1 minute on a standard laptop.

6 Discussion

In this paper, we have introduced Max-and-Smooth, a new two-step approximate posterior sampling scheme for extended LGMs with independent data replicates. Extended LGMs include a diverse range of important applications, such as regression models with spatial or temporal effects in regression coefficients and error variances. Our proposed Max-and-Smooth approach is fast and well-suited for this class of LGMs with a complex and high-dimensional latent structure and multivariate link functions. The first step of the inference scheme involves approximating the likelihood function around its mode by its asymptotic Gaussian density, and we also explored another Gaussian approximation which uses the mean and variance of the normalized likelihood. We suggest that the adequacy of the Gaussian approximation be explored case by case. The second step involves Bayesian inference for the latent parameters and the hyperparameters of the model such that the uncertainty from the first step is correctly propagated into the uncertainty of the posterior distributions. Max-and-Smooth contrasts with the INLA software, which is designed for LGMs with a univariate link function but not for extended LGMs.

Our approach scales well with the dimensions of the multivariate link function and the latent parameter vector given that the precision matrices of the Gaussian priors at the latent level are sparse. This is important for high-dimensional applications, e.g., models with high-dimensional temporal and spatial effects. Additionally, exploiting conjugacy of the Gaussian–Gaussian pseudo model (similar to the INLA method), posterior samples from the marginal density of the hyperparameters can be generated efficiently with low sample autocorrelation. For our approach to be successful, we need the Gaussian likelihood approximation in the first step of our two-step inference approach to be accurate. It turns out that models with non-Gaussian likelihood functions can often be aggregated over groups of independently replicated observations, such that the joint likelihood function of the group is close to a multivariate Gaussian likelihood. In our paper, we have provided evidence that Max-and-Smooth indeed provides accurate and reliable inference in various settings. The computational cost of our proposed inference scheme is close to being insensitive to the number of independent data replicates. In the case of a high-dimensional latent vector modeled with sparse precision matrices and a large number of independent data replicates, the speedup is substantial in comparison to a Markov chain Monte Carlo inference scheme which samples from the exact posterior density. A case study comparing Max-and-Smooth with INLA (in a simplified scenario with a univariate link function) showed that INLA seems to do better at approximating marginal posteriors of latent variables, whereas Max-and-Smooth produces better approximations of marginal posteriors of hyperparameters. The observed differences between posterior approximations are minor, and the extent to which the results generalize to different applications will be subject of future research.

Max-and-Smooth is designed for a rich class of flexible models. It is straightforward to implement, and when sparse precision matrices are used at the latent level, high-dimensional latent vectors can be handled. This paves the way for having an additional set of complex models that are feasible for the analysis of complex and high-dimensional datasets. Future research involves exploring ways to drop the conditional independence assumption at the observation level, and computing the Gaussian approximation based on the normalized likelihood when analytical results for the mean and variance are not available, and when the dimension of the parameter vector within groups is high.