integrated nested Laplace approximations

Let $\backslash boldsymbol\{y\}=(y\_1,\backslash dots,y\_n)$ denote the response variable (that is, the observations) which belongs to an exponential family, with the mean $\backslash mu\_i$ (of $y\_i$) being linked to a linear predictor $\backslash eta\_i$ via an appropriate link function. The linear predictor can take the form of a (Bayesian) additive model. All latent effects (the linear predictor, the intercept, coefficients of possible covariates, and so on) are collectively denoted by the vector $\backslash boldsymbol\{x\}$. The hyperparameters of the model are denoted by $\backslash boldsymbol\{\backslash theta\}$. As per Bayesian statistics, $\backslash boldsymbol\{x\}$ and $\backslash boldsymbol\{\backslash theta\}$ are random variables with prior distributions.

The observations are assumed to be conditionally independent given $\backslash boldsymbol\{x\}$ and $\backslash boldsymbol\{\backslash theta\}$:

$$\backslash pi(\backslash boldsymbol\{y\}\; |\backslash boldsymbol\{x\},\; \backslash boldsymbol\{\backslash theta\})\; =\; \backslash prod\_\{i\; \backslash in\; \backslash mathcal\{I\}\}\backslash pi(y\_i\; |\; \backslash eta\_i,\; \backslash boldsymbol\{\backslash theta\}),$$

where $\backslash mathcal\{I\}$ is the set of indices for observed elements of $\backslash boldsymbol\{y\}$ (some elements may be unobserved, and for these INLA computes a posterior predictive distribution). Note that the linear predictor $\backslash boldsymbol\{\backslash eta\}$ is part of $\backslash boldsymbol\{x\}$.

For the model to be a latent Gaussian model, it is assumed that $\backslash boldsymbol\{x\}|\backslash boldsymbol\{\backslash theta\}$ is a Gaussian Markov Random Field (GMRF) (that is, a multivariate Gaussian with additional conditional independence properties) with probability density

$$\backslash pi(\backslash boldsymbol\{x\}\; |\; \backslash boldsymbol\{\backslash theta\})\; \backslash propto\; \backslash left|\; \backslash boldsymbol\{Q\_\{\backslash theta\}\}\; \backslash right|^\{1/2\}\; \backslash exp\; \backslash left(\; -\backslash frac\{1\}\{2\}\; \backslash boldsymbol\{x\}^T\; \backslash boldsymbol\{Q\_\{\backslash theta\}\}\; \backslash boldsymbol\{x\}\; \backslash right),$$

where $\backslash boldsymbol\{Q\_\{\backslash theta\}\}$ is a $\backslash boldsymbol\{\backslash theta\}$-dependent sparse precision matrix and $\backslash left|\; \backslash boldsymbol\{Q\_\{\backslash theta\}\}\; \backslash right|$ is its determinant. The precision matrix is sparse due to the GMRF assumption. The prior distribution $\backslash pi(\backslash boldsymbol\{\backslash theta\})$ for the hyperparameters need not be Gaussian. However, the number of hyperparameters, $m=\backslash mathrm\{dim\}(\backslash boldsymbol\{\backslash theta\})$, is assumed to be small (say, less than 15).

In Bayesian inference, one wants to solve for the posterior distribution of the latent variables $\backslash boldsymbol\{x\}$ and $\backslash boldsymbol\{\backslash theta\}$. Applying Bayes' theorem

\pi(\boldsymbol{x}, \boldsymbol{\theta} | \boldsymbol{y}) = \frac{\pi(\boldsymbol{y} | \boldsymbol{x}, \boldsymbol{\theta})\pi(\boldsymbol{x} | \boldsymbol{\theta}) \pi(\boldsymbol{\theta}) }{\pi(\boldsymbol{y})},

the joint posterior distribution of $\backslash boldsymbol\{x\}$ and $\backslash boldsymbol\{\backslash theta\}$ is given by

\begin{align}

\pi(\boldsymbol{x}, \boldsymbol{\theta} | \boldsymbol{y}) & \propto \pi(\boldsymbol{\theta})\pi(\boldsymbol{x}|\boldsymbol{\theta}) \prod_i \pi(y_i | \eta_i, \boldsymbol{\theta}) \\

& \propto \pi(\boldsymbol{\theta}) \left| \boldsymbol{Q_{\theta}} \right|^{1/2} \exp \left( -\frac{1}{2} \boldsymbol{x}^T \boldsymbol{Q_{\theta}} \boldsymbol{x} + \sum_i \log \left[\pi(y_i | \eta_i, \boldsymbol{\theta}) \right] \right).

\end{align}

Obtaining the exact posterior is generally a very difficult problem. In INLA, the main aim is to approximate the posterior marginals

\begin{array}{rcl}

\pi(x_i | \boldsymbol{y}) &=& \int \pi(x_i | \boldsymbol{\theta}, \boldsymbol{y}) \pi(\boldsymbol{\theta} | \boldsymbol{y}) d\boldsymbol{\theta} \\

\pi(\theta_j | \boldsymbol{y}) &=& \int \pi(\boldsymbol{\theta} | \boldsymbol{y}) d \boldsymbol{\theta}_{-j} ,

\end{array}

where $\backslash boldsymbol\{\backslash theta\}\_\{-j\}\; =\; \backslash left(\backslash theta\_1,\; \backslash dots,\; \backslash theta\_\{j-1\},\; \backslash theta\_\{j+1\},\; \backslash dots,\; \backslash theta\_m\; \backslash right)$.

A key idea of INLA is to construct nested approximations given by

\begin{array}{rcl}

\widetilde{\pi}(x_i | \boldsymbol{y}) &=& \int \widetilde{\pi}(x_i | \boldsymbol{\theta}, \boldsymbol{y}) \widetilde{\pi}(\boldsymbol{\theta} | \boldsymbol{y}) d\boldsymbol{\theta} \\

\widetilde{\pi}(\theta_j | \boldsymbol{y}) &=& \int \widetilde{\pi}(\boldsymbol{\theta} | \boldsymbol{y}) d \boldsymbol{\theta}_{-j} ,

\end{array}

where $\backslash widetilde\{\backslash pi\}(\backslash cdot\; |\; \backslash cdot)$ is an approximated posterior density. The approximation to the marginal density $\backslash pi(x\_i\; |\; \backslash boldsymbol\{y\})$ is obtained in a nested fashion by first approximating $\backslash pi(\backslash boldsymbol\{\backslash theta\}\; |\; \backslash boldsymbol\{y\})$ and $\backslash pi(x\_i\; |\; \backslash boldsymbol\{\backslash theta\},\; \backslash boldsymbol\{y\})$, and then numerically integrating out $\backslash boldsymbol\{\backslash theta\}$ as

\begin{align}

\widetilde{\pi}(x_i | \boldsymbol{y}) = \sum_k \widetilde{\pi}\left( x_i | \boldsymbol{\theta}_k, \boldsymbol{y} \right) \times \widetilde{\pi}( \boldsymbol{\theta}_k | \boldsymbol{y}) \times \Delta_k,

\end{align}

where the summation is over the values of $\backslash boldsymbol\{\backslash theta\}$, with integration weights given by $\backslash Delta\_k$. The approximation of $\backslash pi(\backslash theta\_j\; |\; \backslash boldsymbol\{y\})$ is computed by numerically integrating $\backslash boldsymbol\{\backslash theta\}\_\{-j\}$ out from $\backslash widetilde\{\backslash pi\}(\backslash boldsymbol\{\backslash theta\}\; |\; \backslash boldsymbol\{y\})$.

To get the approximate distribution $\backslash widetilde\{\backslash pi\}(\backslash boldsymbol\{\backslash theta\}\; |\; \backslash boldsymbol\{y\})$, one can use the relation

\begin{align}

{\pi}( \boldsymbol{\theta} | \boldsymbol{y}) = \frac{\pi\left(\boldsymbol{x}, \boldsymbol{\theta}, \boldsymbol{y} \right)}{\pi\left(\boldsymbol{x} | \boldsymbol{\theta}, \boldsymbol{y} \right) \pi(\boldsymbol{y})},

\end{align}

as the starting point. Then $\backslash widetilde\{\backslash pi\}(\; \backslash boldsymbol\{\backslash theta\}\; |\; \backslash boldsymbol\{y\})$ is obtained at a specific value of the hyperparameters $\backslash boldsymbol\{\backslash theta\}\; =\; \backslash boldsymbol\{\backslash theta\}\_k$ with Laplace's approximation

\begin{align}

\widetilde{\pi}( \boldsymbol{\theta}_k | \boldsymbol{y}) &\propto \left . \frac{\pi\left(\boldsymbol{x}, \boldsymbol{\theta}_k, \boldsymbol{y} \right)}{\widetilde{\pi}_G\left(\boldsymbol{x} | \boldsymbol{\theta}_k, \boldsymbol{y} \right)} \right \vert_{\boldsymbol{x} = \boldsymbol{x}^{*}(\boldsymbol{\theta}_k)}, \\

& \propto \left . \frac{\pi(\boldsymbol{y} | \boldsymbol{x}, \boldsymbol{\theta}_k)\pi(\boldsymbol{x} | \boldsymbol{\theta}_k) \pi(\boldsymbol{\theta}_k)}{\widetilde{\pi}_G\left(\boldsymbol{x} | \boldsymbol{\theta}_k, \boldsymbol{y} \right)} \right \vert_{\boldsymbol{x} = \boldsymbol{x}^{*}(\boldsymbol{\theta}_k)},

\end{align}

where $\backslash widetilde\{\backslash pi\}\_G\backslash left(\backslash boldsymbol\{x\}\; |\; \backslash boldsymbol\{\backslash theta\}\_k,\; \backslash boldsymbol\{y\}\; \backslash right)$ is the Gaussian approximation to $\{\backslash pi\}\backslash left(\backslash boldsymbol\{x\}\; |\; \backslash boldsymbol\{\backslash theta\}\_k,\; \backslash boldsymbol\{y\}\; \backslash right)$ whose mode at a given $\backslash boldsymbol\{\backslash theta\}\_k$ is $\backslash boldsymbol\{x\}^\{*\}(\backslash boldsymbol\{\backslash theta\}\_k)$. The mode can be found numerically for example with the Newton-Raphson method.

The trick in the Laplace approximation above is the fact that the Gaussian approximation is applied on the full conditional of $\backslash boldsymbol\{x\}$ in the denominator since it is usually close to a Gaussian due to the GMRF property of $\backslash boldsymbol\{x\}$. Applying the approximation here improves the accuracy of the method, since the posterior $\{\backslash pi\}(\; \backslash boldsymbol\{\backslash theta\}\; |\; \backslash boldsymbol\{y\})$ itself need not be close to a Gaussian, and so the Gaussian approximation is not directly applied on $\{\backslash pi\}(\; \backslash boldsymbol\{\backslash theta\}\; |\; \backslash boldsymbol\{y\})$. The second important property of a GMRF, the sparsity of the precision matrix $\backslash boldsymbol\{Q\}\_\{\backslash boldsymbol\{\backslash theta\}\_k\}$, is required for efficient computation of $\backslash widetilde\{\backslash pi\}(\; \backslash boldsymbol\{\backslash theta\}\_k\; |\; \backslash boldsymbol\{y\})$ for each value $\{\backslash boldsymbol\{\backslash theta\}\_k\}$.

Obtaining the approximate distribution $\backslash widetilde\{\backslash pi\}\backslash left(\; x\_i\; |\; \backslash boldsymbol\{\backslash theta\}\_k,\; \backslash boldsymbol\{y\}\; \backslash right)$ is more involved, and the INLA method provides three options for this: Gaussian approximation, Laplace approximation, or the simplified Laplace approximation. For the numerical integration to obtain $\backslash widetilde\{\backslash pi\}(x\_i\; |\; \backslash boldsymbol\{y\})$, also three options are available: grid search, central composite design, or empirical Bayes.

{{reflist}}

• {{cite book |first=Virgilio |last=Gomez-Rubio |title=Bayesian inference with INLA |location= |publisher=Chapman and Hall/CRC |year=2021 |isbn=978-1-03-217453-2 }}

Category:Computational statistics

Category:Bayesian inference