# Fast, Scalable Approximations to Posterior Distributions in Extended Latent
Gaussian Models

Alex Stringer
Department of Statistics and Actuarial Science, University of Waterloo
and
Patrick Brown
Department of Statistical Sciences, University of Toronto
Centre for Global Health Research
and
Jamie Stafford
Department of Statistical Sciences, University of Toronto  The authors
gratefully acknowledge funding from the Natural Sciences and Engineering
Research Council of Canada.

###### Abstract

We define a novel class of additive models, called Extended Latent Gaussian
Models, that allow for a wide range of response distributions and flexible
relationships between the additive predictor and mean response. The new class
covers a broad range of interesting models including multi-resolution spatial
processes, partial likelihood-based survival models, and multivariate
measurement error models. Because computation of the exact posterior
distribution is infeasible, we develop a fast, scalable approximate Bayesian
inference methodology for this class based on nested Gaussian, Laplace, and
adaptive quadrature approximations. We prove that the error in these
approximate posteriors is $o_{p}(1)$ under standard conditions, and provide
numerical evidence suggesting that our method runs faster and scales to larger
datasets than methods based on Integrated Nested Laplace Approximations and
Markov Chain Monte Carlo, with comparable accuracy. We apply the new method to
the mapping of malaria incidence rates in continuous space using aggregated
data, mapping leukaemia survival hazards using a Cox Proportional-Hazards
model with a continuously-varying spatial process, and estimating the mass of
the Milky Way Galaxy using noisy multivariate measurements of the positions
and velocities of star clusters in its orbit.

Keywords: Additive model; Approximate inference; Bayesian methods; Laplace
approximation; Spatial statistics.

##  1 Introduction

###  1.1 Latent Gaussian and Extended Latent Gaussian Models

Latent Gaussian Models (LGMs; Rue et al. 2009) are a class of additive
regression models that are of broad interest in modern practice. The class of
LGMs includes generalized linear models (Fing et al., 2010), generalized
additive models (Wood, 2011), semi-parametric (Martino et al., 2011) and joint
(van Niekerk et al., 2021) models for survival data, and spatial models
(Lindgren et al., 2011; Lindgren and Rue, 2015). However, a central assumption
of the LGM (1) is that the response variables $y_{i}$ are conditionally
independent, given a single natural parameter $\eta_{i}$ each, or
equivalently, that for each $i\in[n]$, $\mu_{i}=\mathbb{E}(y_{i}|\eta_{i})$ is
a function of only $\eta_{i}$. This limits the breadth of models that belong
to the class of LGMs.

The contributions of this paper are: (a) a class of models is introduced,
called Extended Latent Gaussian Models (ELGMs), that permit much more general
dependence between the $\mu_{i}$ and the natural parameters (§3.1), (b) a fast
and scalable algorithm is developed for making approximate Bayesian inferences
in complex ELGMs and with large data sets (§3.2), and (c) convergence of the
approximate posteriors is established (Theorem 1 in §4). The new class of
ELGMs includes multi-resolution spatial point processes (Nandi et al., 2020;
Lee et al., 2017; Diggle et al., 2013), partial likelihood-based survival
(Diggle et al., 2010) and multinomial logistic (Stringer et al., 2021) models,
zero-inflated models (Hall, 2000) including those for spatially-correlated
data (Diggle and Giorgi, 2016), and models incorporating multivariate
measurement errors (Eadie et al., 2017).

###  1.2 Existing Methods for Latent Gaussian Models

Bayesian inferences for LGMs are based upon a posterior distribution which is
intractable, and approximate inference algorithms mitigate this challenge. A
number of Markov Chain Monte Carlo (MCMC) algorithms are available, including
those based on block-updating (Knorr-Held and Rue, 2002) and elliptical slice
sampling (Murray et al., 2010), as well as the general algorithms available in
the STAN language (Carpenter et al., 2017). Beyond MCMC, deterministic
approximation methods exist that offer speed and scalability to large data.
The Integrated Nested Laplace Approximation (INLA; Rue et al. (2009); Lindgren
and Rue (2015); Martins et al. (2013); Rue et al. (2017)) is one such method
based on Gaussian, Laplace, and quadrature approximations. Variational
inference (Blei et al., 2017) is another such method, and has been applied in
models for which INLA cannot be used (Ding et al., 2018; Law et al., 2018).

The use of Laplace and quadrature approximations in fitting generalized linear
mixed models is well precedented (Breslow and Clayton, 1993; Pinheiro and
Bates, 1995). Recent approximation algorithms for generalized additive models
make use of Gaussian approximations (Wood, 2016). Taken separately, these
posterior approximations each have established convergence theory (Walker,
1969; Shun and McCullagh, 1995; Jin and Andersson, 2020). However, as
described in §2 of this paper, INLA is based on a nested application of these
three types of approximations. Consequently, no formal convergence results
have been derived for INLA. The method we propose in §3 is based on a nested
approximation strategy similar to INLA, but applies to a richer class of
models (§3), scales well with increasing sample size (§5), and has
approximation error shown to be $o_{p}(1)$ under standard conditions (Theorem
1, §4; Online Supplement A).

###  1.3 Computation For Extended Latent Gaussian Models

An LGM (Equation 1 in §2) depends on an additive predictor or natural
parameter $\boldsymbol{\eta}=\boldsymbol{Z}\boldsymbol{w}\in\mathbb{R}^{n}$,
where $\boldsymbol{w}\in\mathbb{R}^{m}$ are parameters, and
$\boldsymbol{Z}\in\mathbb{R}^{n\times m}$ is a known design matrix. The
likelihood may depend on additional parameters $\boldsymbol{\theta}_{1}$ (§2).
Inferences for $\boldsymbol{w}$ and $\boldsymbol{\theta}_{1}$ based on
Gaussian and Laplace approximations require evaluation of the (negative)
Hessian matrix,
$\boldsymbol{C}_{\boldsymbol{w}}(\boldsymbol{w},\boldsymbol{\theta}_{1})\in\mathbb{R}^{m\times
m}$ of the log-likelihood
$\pi(\boldsymbol{Y}|\boldsymbol{\eta},\boldsymbol{\theta}_{1})=\prod_{i=1}^{n}\pi(Y_{i}|\eta_{i},\boldsymbol{\theta}_{1})$.
However, INLA employs a further modification of the LGM (Rue et al., 2009),
basing inference off of the noisy parameter
$\boldsymbol{w}^{*}=(\boldsymbol{\eta}^{*},\boldsymbol{w})$, where
$\boldsymbol{\eta}^{*}\sim\text{N}(\boldsymbol{\eta},\tau^{-1}\boldsymbol{I}_{n})$
for some large, fixed $\tau>0$. Consequently, the Hessian matrix that INLA
must evaluate to make approximate inferences is
$\boldsymbol{C}_{\boldsymbol{w}^{*}}(\boldsymbol{w}^{*},\boldsymbol{\theta}_{1})\in\mathbb{R}^{(m+n)\times(m+n)}$,
a larger matrix whose size increases linearly with the sample size, $n$,
irrespective of the properties of the model being fit.

The conditional independence of $y_{i}$ given $\eta_{i}$ that is required in
an LGM forces the upper left $n\times n$ block of
$\boldsymbol{C}_{\boldsymbol{w}^{*}}(\boldsymbol{w}^{*},\boldsymbol{\theta}_{1})$
to be diagonal, and the use of algorithms for efficient decomposition of
sparse matrices (Rue, 2001) make computations involving this matrix feasible
for LGMs fit to moderately sized data. In contrast, in the ELGMs we introduce
in §3, this conditional independence is relaxed, and the upper left $n\times
n$ block of
$\boldsymbol{C}_{\boldsymbol{w}^{*}}(\boldsymbol{w}^{*},\boldsymbol{\theta}_{1})$
is no longer diagonal and may, in fact, be dense. In §3.2 we present a
computational strategy for approximate inference in ELGMs that does not make
use of INLA’s noisy additive predictor. This leads to matrices which may be
denser in the case of LGMs, but, cruicially, are always smaller by $n$
dimensions than the corresponding matrices that would be required if
implementing INLA’s approximations for ELGMs. The substantial performance
gains of our strategy, especially at large sample sizes, are demonstrated in
§5.

###  1.4 Plan of Paper

This paper is organized as follows. In §2 we describe existing approaches to
approximate inference in LGMs. In §3 we introduce the class of ELGMs and an
algorithm (Algorithm 1) for making approximate Bayesian inferences based on
them. In §4 we present a formal convergence result for our posterior
approximation in the form of Theorem 1, and discuss model-specific theoretical
considerations. In §5 we provide numerical evidence suggesting that our
approximations are faster and scale to larger datasets than other methods,
while remaining accurate.

In §6, we apply our new methodology to fit three challenging models which are
beyond the class of LGMs but belong to the new class of ELGMs. We fit a fully
Bayesian analysis of aggregated spatial point process data, a challenge which
Li et al. (2012) and others address using data augmentation techniques and
MCMC; the existence of software for doing so in this example enables direct
comparison of our method with MCMC techniques where we make similar inferences
at substantial reduction in computational cost. We fit a Cox proportional
hazards model with partial likelihood and spatially-varying hazard; previous
analyses of these data employ parametric (Henderson et al., 2002; Lindgren et
al., 2011) and semi-parametric (Martino et al., 2011) models for the hazard
due to the difficulty of working with both the partial likelihood and the
latent spatial process. Finally, we fit the Galactic Mass Estimation model of
Eadie et al. (2017) and Eadie and Juric (2019) for estimating the mass of the
Milky Way galaxy in the presence of multivariate measurement errors.
Implementation makes use of the publicly available aghq package in the R
language, and code for all examples is made available at
https://github.com/awstringer1/elgm-paper-code.

##  2 Preliminaries

###  2.1 Notation and Definitions

In what follows, we use the notation $[n]=\left\\{1,\ldots,n\right\\}$, denote
(column) vectors by lower-case
$\boldsymbol{v}=(v_{i})_{i\in[n]}\in\mathbb{R}^{n}$, and sub-vectors by
$\boldsymbol{v}_{S}=(v_{i})_{i\in S}$ for any $S\subseteq[n]$. We denote the
concatenation of vectors $\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{p}$,
having dimensions $d_{1},\ldots,d_{p}$, by
$\boldsymbol{v}=(\boldsymbol{v}_{j})_{j\in[p]}=(v_{11},\ldots,v_{pd_{p}})\in\mathbb{R}^{d_{1}+\cdots+d_{p}}$.
We denote matrices by upper-case
$\boldsymbol{Z}=(Z_{ij})_{i\in[n],j\in[d]}\in\mathbb{R}^{n\times d}$. An LGM
takes the following form:

$\displaystyle y_{i}|\eta_{i},\boldsymbol{\theta}_{1}$
$\displaystyle\overset{ind}{\sim}\pi(y_{i}|\eta_{i},\boldsymbol{\theta}_{1}),\
\mu_{i}=\mathbb{E}(Y_{i}|\eta_{i})=g(\eta_{i}),\
\eta_{i}=\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}+\sum_{q=1}^{r}u_{q}(z_{iq}),$
(1)

The response variable is $\boldsymbol{y}=(y_{i})_{i\in[n]}$ having likelihood
$\pi(\boldsymbol{y}|\boldsymbol{\eta},\boldsymbol{\theta}_{1})=\prod_{i=1}^{n}\pi(y_{i}|\eta_{i},\boldsymbol{\theta}_{1})$.
Denote the conditional mean of each response by
$\mu_{i}=\mathbb{E}(y_{i}|\eta_{i}),i\in[n]$, and define an inverse link
function $g:\mathbb{R}\to\mathbb{R}$ such that $\mu_{i}=g(\eta_{i})$. The
parameters $\boldsymbol{\theta}_{1}\in\mathbb{R}^{s_{1}}$, where $s_{1}$ is
assumed small, represent additional parameters of the likelihood, such as
dispersion parameters. We denote the observed covariates, assumed fixed and
known, by $\boldsymbol{x}_{i}\in\mathbb{R}^{p}$ and
$\boldsymbol{z}_{i}\in\mathbb{R}^{r}$.

Our treatment of the unknown functions $u_{q}(\cdot),q\in[r]$ follows that of
Wood (2016). We define a basis expansion
$u_{q}(z)=\sum_{j=1}^{d_{q}}\phi_{qj}(z)u_{qj}$, where the $\phi_{qj}(\cdot)$
are known basis functions (subject to identifiability constraints) and the
$\boldsymbol{u}_{q}=(u_{qj})_{j\in[d_{q}]}$ are unknown weights to be
inferred. Different choices of $\phi_{qj}(\cdot)$ lead to different models,
including: random effects models (§5), semi-parametric regression/generalized
additive models, spatial models (§6.1,§6.2), and measurement error models
(§6.3). The full vector of parameters (basis function weights) relating to
$u_{q}(\cdot),q\in[r]$ is
$\boldsymbol{u}=(\boldsymbol{u}_{q})_{q\in[r]}\in\mathbb{R}^{d}$, with
$d=\sum_{q=1}^{r}d_{q}$.

Define
$\boldsymbol{w}=(\boldsymbol{u},\boldsymbol{\beta})\in\mathbb{R}^{m},m=d+p$,
and write $\boldsymbol{\eta}=\boldsymbol{Z}\boldsymbol{w}$ as in §1, for
design matrix $\boldsymbol{Z}\in\mathbb{R}^{n\times m}$ depending on
$\phi_{qj},\boldsymbol{z}_{i},\boldsymbol{x}_{i},i\in[n],j\in[d],q\in[r]$. For
the remainder, we will use both
$\pi(\boldsymbol{y}|\boldsymbol{\eta},\boldsymbol{\theta}_{1})$ and
$\pi(\boldsymbol{y}|\boldsymbol{w},\boldsymbol{\theta}_{1})$ to denote the
likelihood, as context dictates. A Gaussian prior is placed on
$\boldsymbol{u}|\boldsymbol{\theta}_{2}\sim\text{N}[\boldsymbol{0},\boldsymbol{Q}_{\boldsymbol{u}}^{-1}(\boldsymbol{\theta}_{2})]$,
where the known precision (inverse covariance) matrix
$\boldsymbol{Q}_{\boldsymbol{u}}(\boldsymbol{\theta}_{2})$ depends upon a
further parameter $\boldsymbol{\theta}_{2}\in\mathbb{R}^{s_{2}}$, where
$s_{2}$ is assumed small. A Gaussian prior
$\boldsymbol{\beta}\sim\text{N}\left(\boldsymbol{0},\boldsymbol{\Sigma}_{\boldsymbol{\beta}}\right)$
with fixed covariance $\boldsymbol{\Sigma}_{\boldsymbol{\beta}}$, usually
taken to be diagonal, is placed on $\boldsymbol{\beta}$. Let
$\boldsymbol{Q}(\boldsymbol{\theta}_{2})=\text{diag}\left\\{\boldsymbol{Q}_{\boldsymbol{u}}(\boldsymbol{\theta}_{2}),\boldsymbol{\Sigma}_{\boldsymbol{\beta}}^{-1}\right\\}$,
so
$\boldsymbol{w}|\boldsymbol{\theta}_{2}\sim\text{N}[\boldsymbol{0},\boldsymbol{Q}^{-1}(\boldsymbol{\theta}_{2})]$.
Further define
$\boldsymbol{\theta}=(\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2})\in\mathbb{R}^{s},s=s_{1}+s_{2}$,
with prior $\pi(\boldsymbol{\theta})$.

###  2.2 Posteriors and Approximate Inference

Inferences for the parameters of interest, $\boldsymbol{w}$ and
$\boldsymbol{\theta}$, are based, in principle, on the posterior distributions

$\displaystyle\pi(\boldsymbol{w}|\boldsymbol{y})$
$\displaystyle=\int\pi(\boldsymbol{w}|\boldsymbol{y},\boldsymbol{\theta})\pi(\boldsymbol{\theta}|\boldsymbol{y})d\boldsymbol{\theta},$
(2) $\displaystyle\pi(\boldsymbol{\theta}|\boldsymbol{y})$
$\displaystyle=\int\pi(\boldsymbol{w},\boldsymbol{\theta}|\boldsymbol{y})d\boldsymbol{w},$

where

$\pi(\boldsymbol{w}|\boldsymbol{y},\boldsymbol{\theta})\propto\exp\left\\{-\frac{1}{2}\boldsymbol{w}^{\texttt{T}}\boldsymbol{Q}(\boldsymbol{\theta}_{2})\boldsymbol{w}+\log\pi(\boldsymbol{y}|\boldsymbol{w},\boldsymbol{\theta}_{1})\right\\}.$
(3)

All of these posteriors are intractable, and in practice, inferences are based
on approximations. For fixed $\boldsymbol{\theta}$, define
$\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}=\text{argmax}_{\boldsymbol{w}}\log\pi(\boldsymbol{w},\boldsymbol{y},\boldsymbol{\theta})$
and
$\boldsymbol{H}_{\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}}(\boldsymbol{\theta})=-\partial^{2}_{\boldsymbol{w}}\log\pi(\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}},\boldsymbol{y},\boldsymbol{\theta})=\boldsymbol{Q}(\boldsymbol{\theta}_{2})+\boldsymbol{C}_{\boldsymbol{w}}(\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}},\boldsymbol{\theta}_{1})$.
A Gaussian approximation to
$\pi(\boldsymbol{w}|\boldsymbol{y},\boldsymbol{\theta})$ is given by

$\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\boldsymbol{w}|\boldsymbol{y},\boldsymbol{\theta})=(2\pi)^{-m/2}\left|\boldsymbol{H}_{\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}}(\boldsymbol{\theta})\right|^{1/2}\exp\left\\{-\frac{1}{2}(\boldsymbol{w}-\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}})^{\texttt{T}}\boldsymbol{H}_{\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}}(\boldsymbol{\theta})(\boldsymbol{w}-\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}})\right\\}.$
(4)

Define $\mathcal{Q}(s,k)$ to be the set of points, and
$\omega_{k}:\mathcal{Q}(s,k)\to\mathbb{R}$ to be the weights, from a
quadrature rule (to be specified) in $s$ dimensions with $k$ points per
dimension. Let $\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})=\pi(\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}},\boldsymbol{\theta},\boldsymbol{y})/\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}|\boldsymbol{y},\boldsymbol{\theta})$,
$\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt
LA}}}=\text{argmax}_{\boldsymbol{\theta}}\log\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})$,
$\widehat{\boldsymbol{H}}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta})=-\partial^{2}\log\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})$, and
$\widehat{\boldsymbol{H}}_{\text{\tiny{\tt
LA}}}^{-1}(\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt
LA}}})=\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\widehat{\boldsymbol{L}}_{\text{\tiny{\tt LA}}}^{\texttt{T}}$ where
$\widehat{\boldsymbol{L}}_{\text{\tiny{\tt LA}}}$ is the lower Cholesky
triangle. A Laplace approximation to $\pi(\boldsymbol{\theta}|\boldsymbol{y})$
is given by:

$\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta}|\boldsymbol{y})=\frac{\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})}{\left|\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\right|\sum_{\boldsymbol{z}\in\mathcal{Q}(s,k)}\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\boldsymbol{z}+\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt
LA}}},\boldsymbol{y})\omega_{k}(\boldsymbol{z})}.$ (5)

Finally, define the nested approximation

$\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})=\left|\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\right|\sum_{\boldsymbol{z}\in\mathcal{Q}(s,k)}\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\boldsymbol{w}|\boldsymbol{y},\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\boldsymbol{z}+\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt
LA}}})\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\boldsymbol{z}+\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt
LA}}}|\boldsymbol{y})\omega_{k}(\boldsymbol{z}),$ (6)

where the same adapted points and weights have been used to approximate this
second, different integral.

Rue et al. (2009) make approximate inferences starting from (6), with three
important differences. First, their INLA method focusses on approximating
marginal posteriors $\widetilde{\pi}(w_{j}|\boldsymbol{y}),j\in[m]$, for which
an additional Laplace approximation is used, using a computational strategy
recently expanded upon by Wood (2020). Our focus in the present work is on the
joint posterior $\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})$, as evidenced
in the examples of §5 and §6, so we do not pursue this here for ELGMs. Second
is the choice of quadrature rule, for which they consider the empirical
performance of several different options. As we demonstrate in the proof of
Theorem 1 in Online Supplement A, the specific properties of adaptive (Gauss-
Hermite) quadrature are used when proving convergence of the overall
approximation (6), and it is not immediate that similar properties hold for
other rules.

The third difference is that INLA bases inferences off of a modified model,
with noisy parameter
$\boldsymbol{w}^{*}=(\boldsymbol{\eta}^{*},\boldsymbol{w})$. Evaluating the
Gaussian approximation $\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\boldsymbol{w}^{*}|\boldsymbol{y},\boldsymbol{\theta})$ required to
compute (4), (5), and (6) under the modification employed by INLA therefore
requires computing, storing, and decomposing the matrix
$\boldsymbol{H}^{*}_{\widehat{\boldsymbol{w}}^{*}_{\boldsymbol{\theta}}}(\boldsymbol{\theta})=\boldsymbol{Q}^{*}(\boldsymbol{\theta}_{2})+\boldsymbol{C}_{\boldsymbol{w}^{*}}(\widehat{\boldsymbol{w}}^{*}_{\boldsymbol{\theta}},\boldsymbol{\theta}_{1})\in\mathbb{R}^{(n+m)\times(n+m)}$,
where

$\displaystyle\boldsymbol{Q}^{*}(\boldsymbol{\theta}_{2})=\tau\begin{pmatrix}\boldsymbol{I}_{n}&-\boldsymbol{Z}\\\
-\boldsymbol{Z}^{\texttt{T}}&\tau^{-1}\boldsymbol{Q}(\boldsymbol{\theta}_{2})+\boldsymbol{Z}^{\texttt{T}}\boldsymbol{Z}\end{pmatrix},\qquad\boldsymbol{C}_{\boldsymbol{w}^{*}}(\widehat{\boldsymbol{w}}^{*}_{\boldsymbol{\theta}},\boldsymbol{\theta}_{1})=\begin{pmatrix}\boldsymbol{C}_{\boldsymbol{\eta}}(\widehat{\boldsymbol{\eta}}^{*}_{\boldsymbol{\theta}},\boldsymbol{\theta}_{1})&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}\end{pmatrix}.$

Here
$\boldsymbol{w}^{*}\sim\text{N}\left\\{(\boldsymbol{\eta},\boldsymbol{0}),\boldsymbol{Q}^{*}(\boldsymbol{\theta}_{2})^{-1}\right\\}$,
$\boldsymbol{C}_{\boldsymbol{\eta}}(\boldsymbol{\eta},\boldsymbol{\theta}_{1})=-\partial^{2}_{\boldsymbol{\eta}}\log\pi(\boldsymbol{y}|\boldsymbol{\eta},\boldsymbol{\theta}_{1})\in\mathbb{R}^{n\times
n}$, and hence
$\boldsymbol{C}_{\boldsymbol{w}^{*}}(\widehat{\boldsymbol{w}}^{*}_{\boldsymbol{\theta}},\boldsymbol{\theta}_{1})$
depends on $\widehat{\boldsymbol{w}}^{*}_{\boldsymbol{\theta}}$ through
$\widehat{\boldsymbol{\eta}}^{*}_{\boldsymbol{\theta}}$. Crucially, in an LGM,
the matrix
$\boldsymbol{C}_{\boldsymbol{\eta}}(\widehat{\boldsymbol{\eta}}^{*}_{\boldsymbol{\theta}},\boldsymbol{\theta}_{1})$
is diagonal, owing to the conditional independence of $y_{i}$ given
$\eta_{i}$. The feasibility of computations under the modified model therefore
depends on the conditional independence constraint enforced in an LGM.

In ELGMs, however, this constraint is relaxed, and as a result, the $n\times
n$ matrix
$\boldsymbol{C}_{\boldsymbol{\eta}}(\widehat{\boldsymbol{\eta}}^{*}_{\boldsymbol{\theta}},\boldsymbol{\theta}_{1})$
is no longer diagonal, and in fact may be dense. To summarize, basing
inferences off of the modified model has several implications:

  1. 1.

It becomes challenging to fit models with more complex dependence structures
than LGMs allow, because of the need to store and manipulate large,
potentially dense matrices,

  2. 2.

It becomes challenging to fit LGMs to large datasets, because the size of the
sparse matrices involved scales linearly with $n$ for every such model fit,

  3. 3.

It becomes challenging to study the theoretical properties of the
approximations, because convergence results pertaining to the marginal Laplace
approximation (Tierney and Kadane, 1986) do not apply in this context (Shun
and McCullagh, 1995).

We now define the class of ELGMs and describe our approach to inference for
them which mitigates these challenges.

##  3 Extended Latent Gaussian Models

###  3.1 Definition of ELGMs

In ELGMs, the vector of additive predictors is redefined as
$\boldsymbol{\eta}=(\eta_{j})_{j\in[N_{n}]}$, with dimension
$N_{n}\in\mathbb{N}$ depending on, but not necessarily equal to, $n$. This
dependence is unrestricted, and we may have in any particular example
$N_{n}<n,N_{n}>n$, or $N_{n}=n$. Each mean response $\mu_{i}$ will now depend
on some _subset_ , $\mathcal{J}_{i}\subseteq[N_{n}]$, of indices of
$\boldsymbol{\eta}$. We require that
$\bigcup_{i=1}^{n}\mathcal{J}_{i}=[N_{n}]$ and
$1\leq\left|\mathcal{J}_{i}\right|\leq N_{n}$, but otherwise enforce no
restrictions on the $\mathcal{J}_{i}$. This means, informally, that
observation means may depend on multiple, shared additive predictors, in a
very general manner. In contrast, LGMs always have
$\mathcal{J}_{i}=\left\\{i\right\\},i\in[n]$, and $N_{n}=n$. The inverse link
function $g(\cdot)$ is redefined for each observation to be a many-to-one
mapping $g_{i}:\mathbb{R}^{\left|\mathcal{J}_{i}\right|}\to\mathbb{R}$, such
that $\mu_{i}=g_{i}(\boldsymbol{\eta}_{\mathcal{J}_{i}})$.

Extended Latent Gaussian Models take the following form:

$\displaystyle
y_{i}|\boldsymbol{\eta}_{\mathcal{J}_{i}},\boldsymbol{\theta}_{1}$
$\displaystyle\overset{ind}{\sim}\pi(y_{i}|\boldsymbol{\eta}_{\mathcal{J}_{i}},\boldsymbol{\theta}_{1}),i\in[n],\
\mu_{i}=\mathbb{E}(y_{i}|\boldsymbol{\eta}_{\mathcal{J}_{i}})=g_{i}(\boldsymbol{\eta}_{\mathcal{J}_{i}}),\
\eta_{j}=\boldsymbol{x}_{j}^{\texttt{T}}\boldsymbol{\beta}+\sum_{q=1}^{r}u_{q}(z_{jq}),j\in[N_{n}],$
(7)

which may be contrasted with the LGM (1).

The potentially complex dependence between $\mu_{i}$ and $\boldsymbol{\eta}$
introduces computational burden depending directly on the sets
$\mathcal{J}_{i}$. Specifically, the sparsity pattern of
$\boldsymbol{C}_{\boldsymbol{\eta}}(\boldsymbol{\eta},\boldsymbol{\theta}_{1})$
is determined as follows: for any $k,l\in[N_{n}],$

$\exists
i\in[n]:k,l\in\mathcal{J}_{i}\implies\boldsymbol{C}_{\boldsymbol{\eta}}(\boldsymbol{\eta},\boldsymbol{\theta}_{1})_{kl}\neq
0.$

Because of this relationship, the $N_{n}$-dimensional
$\boldsymbol{C}_{\boldsymbol{\eta}}(\boldsymbol{\eta},\boldsymbol{\theta}_{1})$
matrix becomes denser as the model becomes more complex, as specifically
measured by the complexity of the index sets $\mathcal{J}_{i}$. This in turn
affects the computations required for approximate Bayesian inference in ELGMs.

###  3.2 Computational Considerations

Algorithm 1 shows the full approximation procedure for ELGMs. We highlight two
important points here. First, we base inferences off of $\boldsymbol{w}$
instead of $\boldsymbol{w}^{*}$, and hence the sparsity structure of the
Hessian
$\boldsymbol{H}_{\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}}(\boldsymbol{\theta})=-\partial^{2}_{\boldsymbol{w}}\log\pi(\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}|\boldsymbol{y},\boldsymbol{\theta})=\boldsymbol{Q}(\boldsymbol{\theta}_{2})+\boldsymbol{C}_{\boldsymbol{w}}(\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}},\boldsymbol{\theta}_{1})$
is determined both by that of $\boldsymbol{Q}(\boldsymbol{\theta}_{2})$ and

$\displaystyle\boldsymbol{C}_{\boldsymbol{w}}(\boldsymbol{w},\boldsymbol{\theta}_{1})=\boldsymbol{Z}^{\texttt{T}}\boldsymbol{C}_{\boldsymbol{\eta}}(\boldsymbol{\eta},\boldsymbol{\theta}_{1})\boldsymbol{Z},$

which is (now) of dimension $m\times m$ where $m=\text{dim}(\boldsymbol{w})$,
for any $n$. While
$\boldsymbol{C}_{\boldsymbol{w}}(\boldsymbol{w},\boldsymbol{\theta}_{1})$ may
be less sparse than
$\boldsymbol{C}_{\boldsymbol{w}^{*}}(\boldsymbol{w}^{*},\boldsymbol{\theta}_{1})$
in general, and hence less amenable to the use of efficient sparse matrix
algorithms (Rue, 2001), it is also much smaller, with dimension $n$ less than
that of
$\boldsymbol{C}_{\boldsymbol{w}^{*}}(\boldsymbol{w}^{*},\boldsymbol{\theta}_{1})$.
This has the following implications:

  1. 1.

It becomes feasible to fit ELGMs, because the dimension of the large,
potentially dense matrices involved in fitting is reduced by $n$ compared to
if inferences were based on $\boldsymbol{w}^{*}$,

  2. 2.

It becomes more efficient to fit LGMs to larger datasets (Table 1 in §5), for
this same reason,

  3. 3.

The theoretical properties of the approximation error are more feasible to
study (§4) because the convergence of the marginal Laplace approximation is
established (Tierney and Kadane, 1986).

We still exploit sparsity when available on a problem-specific basis, using
the efficient routines in the Matrix package for R (Bates and Maechler, 2019).
Further, the reduced size of matrices involved is expected to improve
computational efficiency of our approximations at larger sample sizes, where
this difference is more pronounced. We investigate these notions empirically
in §5.

Second, we choose Gauss-Hermite quadrature for constructing the rule
$\mathcal{Q}(s,k),\omega_{k}(\cdot)$ that is adapted to the renormalization of
$\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\boldsymbol{\theta},\boldsymbol{y})$.
A number of suitable extensions $\mathcal{Q}(s,k)$ of the univariate rule
$\mathcal{Q}(1,k)$ to $s$ dimensions are permitted, including product (Jin and
Andersson, 2020), sparse (Heiss and Winschel, 2008), and nested (Petras, 2003)
rules. Beyond being readily available in standard software (Weiser, 2016), the
asymptotic properties of the adaptive GHQ rule were recently established by
Bilodeau et al. (2021), and we use these results in the proof of Theorem 1
(§4; Online Supplement A).

###  3.3 Implementation

Algorithm 1 is implemented in the aghq package for the R language. The core of
the implementation requires computing
$\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt LA}}}$ and
$\widehat{\boldsymbol{H}}_{\text{\tiny{\tt
LA}}}(\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt LA}}})$, which requires
optimization and two numerical derivatives of
$\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\boldsymbol{\theta},\boldsymbol{y})$.
This is complicated by the fact that each evaluation of
$\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\boldsymbol{\theta},\boldsymbol{y})$
requires evaluating $\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\boldsymbol{w}|\boldsymbol{y},\boldsymbol{\theta})$, which requires
further optimization and derivatives of
$\pi(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{y})$ with respect to
$\boldsymbol{w}$.

Our recommendation is to implement
$-\log\pi(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{y})$ using the TMB
package for automatic differentiation in R (Kristensen et al., 2016), which
provides two automatic derivatives of
$\log\pi(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{y})$, automatic
evaluation of $\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})$, and—critically—an automatic
derivative of $\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})$ that avoids repeatedly computing
$\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}$. When a TMB template is not
available, we use trust region optimization to compute
$\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}$, with sparse (Braun, 2014) or
dense (Geyer, 2020) matrix algebra depending on the application, and numerical
derivatives of $\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})$. Assumption 2 in Online Supplement
A states that the objective function
$\pi(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{y})$ is log-concave
(locally, with high probability as $n\to\infty$), making trust region
optimization stable and efficient for this application.

Input:

  * Likelihood $\pi(\boldsymbol{y}|\boldsymbol{w},\boldsymbol{\theta})$, precision matrix $\boldsymbol{Q}(\boldsymbol{\theta}_{2})\in\mathbb{R}^{m\times m}$, prior $\pi(\boldsymbol{\theta})$, Hessian $\boldsymbol{H}_{\boldsymbol{w}}(\boldsymbol{\theta})=-\partial^{2}_{\boldsymbol{w}}\log\pi(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{y})$, number of samples $B\in\mathbb{N}$, number of quadrature points $k\in\mathbb{N}$ and base quadrature rule $\mathcal{Q}(s,k)$.

Compute:

  1. 1.

$\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt
LA}}}=\text{argmax}_{\boldsymbol{\theta}}\ \widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})$ (Equation 5) using numerical
optimization. For each evaluation of $\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})$ at point $\boldsymbol{\theta}$,
compute:

    1. (a)

$\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}=\text{argmax}_{\boldsymbol{w}}\
\pi(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{y})$ using trust region
optimization (§3.3),

    2. (b)

Hessian
$\boldsymbol{H}_{\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}}(\boldsymbol{\theta})$
analytically or using automatic differentiation,

    3. (c)

$\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}|\boldsymbol{y},\boldsymbol{\theta})=(2\pi)^{-m/2}\left|\boldsymbol{H}_{\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}}(\boldsymbol{\theta})\right|^{1/2}$,
and hence $\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})=\pi(\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}},\boldsymbol{\theta},\boldsymbol{y})/\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}}|\boldsymbol{y},\boldsymbol{\theta})$.

  2. 2.

$\widehat{\boldsymbol{H}}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta})=-\partial^{2}_{\boldsymbol{\theta}}\log\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta},\boldsymbol{y})$ using numeric or automatic
differentiation,

  3. 3.

The lower Cholesky triangle $\widehat{\boldsymbol{L}}_{\text{\tiny{\tt LA}}}$
of $\widehat{\boldsymbol{H}}_{\text{\tiny{\tt
LA}}}(\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt
LA}}})^{-1}=\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\widehat{\boldsymbol{L}}_{\text{\tiny{\tt LA}}}^{\texttt{T}}$, directly.

Output:

  1. 1.

Posterior of summaries of $\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta}|\boldsymbol{y})$, using further numerical
integration,

  2. 2.

Independent samples
$(\boldsymbol{w}_{b})_{b\in[B]}\overset{iid}{\sim}\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})$
(Equation 6):

  3. Input: $B\in\mathbb{N}$, $\lambda(\boldsymbol{z})=\left|\widehat{\boldsymbol{L}}_{\text{\tiny{\tt LA}}}\right|\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\widehat{\boldsymbol{L}}_{\text{\tiny{\tt LA}}}\boldsymbol{z}+\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt LA}}}|\boldsymbol{y})\omega_{k}(\boldsymbol{z}),\boldsymbol{z}\in\mathcal{Q}(s,k)$.

  4. For: $b=1,\ldots,B$:

    1. a.

Draw $\boldsymbol{z}_{b}$ from $\mathcal{Q}(s,k)$ with probabilities
$\mathbb{P}(\boldsymbol{z}_{b}=\boldsymbol{z})=\lambda(\boldsymbol{z}),\boldsymbol{z}\in\mathcal{Q}(s,k)$,

    2. b.

Draw $\boldsymbol{w}_{b}$ from $\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\cdot|\boldsymbol{y},\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\boldsymbol{z}_{b}+\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt LA}}})$
(Rue, 2001),

Algorithm 1 Fitting Extended Latent Gaussian Models

##  4 Theoretical Considerations

###  4.1 Approximation error

Inferences made from Algorithm 1 are based on the approximations
$\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\boldsymbol{\theta}|\boldsymbol{y})$
and $\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})$. Convergence results
pertaining to $\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta}|\boldsymbol{y})$ are standard, see Tierney and
Kadane (1986); Kass et al. (1990); Shun and McCullagh (1995); Ogden (2021,
2017) as well as Equations (3) and (5) in Online Supplement A, but note that
these results do not immediately apply when inferences are based on
$\boldsymbol{w}^{*}$ instead of $\boldsymbol{w}$, since
$\text{dim}(\boldsymbol{w}^{*})=O(n)$ (Shun and McCullagh, 1995; Ogden, 2021).
Even when basing inferences off of $\boldsymbol{w}$ as we do, convergence of
$\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})$ is non-trivial.

Define the total variation error
$\left\lVert\pi(\boldsymbol{w}|\boldsymbol{y})-\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})\right\rVert_{\text{{TV}}}=\sup_{\mathcal{K}\in\mathcal{B}(\mathbb{R}^{m})}\left|\int_{\mathcal{K}}\pi(\boldsymbol{w}|\boldsymbol{y})-\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})d\boldsymbol{w}\right|$,
where $\mathcal{B}(\mathbb{R}^{m})$ are the Borel sets on $\mathbb{R}^{m}$. We
have the following result:

######  Theorem 1.

Under Assumptions 1 – 7 (Online Supplement A),
$\left\lVert\pi(\boldsymbol{w}|\boldsymbol{y})-\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})\right\rVert_{\text{{TV}}}=o_{p}(1).$

The proof of Theorem 1 is given in Online Supplement A. Practically, Theorem 1
guarantees convergence in probability of the approximate coverage
$\widetilde{\alpha}=\int_{\mathcal{K}}\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})d\boldsymbol{w}$
of any credible set $\mathcal{K}\in\mathcal{B}(\mathbb{R}^{m})$ to its nominal
level
$\alpha=\int_{\mathcal{K}}\pi(\boldsymbol{w}|\boldsymbol{y})d\boldsymbol{w}$.

Our nested approximation strategy (Equation 6) assumes that each individual
Gaussian (4) and Laplace (5) approximation converges in probability to their
exact counterpart, which is the case when the likelihood satisfies,
respectively, Assumption 6 and Assumptions 1–5, 7 in Online Supplement A. In
that case, Theorem 1 guarantees that our nested approximation (6) preserves
this property, showing the convergence to zero in probability of the total
variation measure between the approximate and exact posterior. This is a non-
trivial result, for two reasons, both of which are also present in the method
of, but not discussed by, Rue et al. (2009, §4.1).

First, the Gaussian approximation $\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\boldsymbol{w}|\boldsymbol{y},\boldsymbol{\theta})$ is evaluated at
$\boldsymbol{\theta}=\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\boldsymbol{z}+\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt
LA}}},\boldsymbol{z}\in\mathcal{Q}(s,k)$, which depend on $\boldsymbol{y}$,
and hence is an approximation to the posterior density
$\pi(\boldsymbol{w}|\boldsymbol{y},\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\boldsymbol{z}+\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt LA}}})$
based on the likelihood
$\pi(\boldsymbol{y}|\boldsymbol{w},\widehat{\boldsymbol{L}}_{\text{\tiny{\tt
LA}}}\boldsymbol{z}+\widehat{\boldsymbol{\theta}}_{\text{\tiny{\tt LA}}})$.
Even under the strong assumption (which we do not make, see Assumptions 1 – 7
in Online Supplement A) that the model is well-specified, in the sense that
$\boldsymbol{y}$ has joint density
$\pi(\boldsymbol{y}|\boldsymbol{w}_{0},\boldsymbol{\theta}_{0})$ for some
$\boldsymbol{w}_{0},\boldsymbol{\theta}_{0}$, this likelihood is misspecified,
because it is not evaluated at $\boldsymbol{\theta}_{0}$, and standard
Bernstein von-Mises theory (van der Vaart, 1998, Ch. 10.2) does not
immediately apply. This notion also applies to the Gaussian and Laplace
approximations employed in this manner by Rue et al. (2009) and Wood (2020),
although this is not discussed.

Second, the quadrature points and weights are adapted to integrate
$\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\boldsymbol{\theta},\boldsymbol{y})$,
but are then used a second time to integrate $\widetilde{\pi}_{\text{\tiny{\tt
G}}}(\boldsymbol{w}|\boldsymbol{y},\boldsymbol{\theta})\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta}|\boldsymbol{y})$ (for any $\boldsymbol{w})$. This
strategy was initially proposed by Rue et al. (2009) and has obvious
computational benefits, but requires theoretical justification. To our
knowledge, the theoretical implications of this have not been investigated in
any other work. In proving Theorem 1, we explicitly handle this theoretical
complexity.

###  4.2 Limitations and practical considerations

Beyond the challenges present in establishing Theorem 1, as noted by Rue et
al. (2009, §4.1), there are a number of practical issues with many
applications of LGMs that make study of the approximation error challenging.
Similar issues also apply to ELGMs, and here we briefly expand on this
discussion. The two main challenges are:

  1. 1.

Despite basing inferences on $\boldsymbol{w}$ instead of $\boldsymbol{w}^{*}$,
many specific ELGMs have $\text{dim}(\boldsymbol{w})$ increasing with $n$ by
construction, and the convergence of the Laplace approximation is not
immediately implied by standard results (Shun and McCullagh, 1995),

  2. 2.

Even under the assumption that the model (7) is well-specified, the basis
function expansion $u_{q}(z)=\sum_{j=1}^{d_{q}}\phi_{qj}(z)u_{qj}$ induces
misspecification into the model, and model-specific work is required to
validate Assumptions 1 – 7.

Theorem 1 can be expected to hold for well-specified models with fixed
parameter dimension. However, it is not the case that Theorem 1 is false for
models exhibiting one or both of the above challenges. Rather, efforts must be
made on a model-specific basis to establish whether the assumptions underlying
Theorem 1 hold. The need to verify complex analytical criteria on a model-
specific basis in order to apply general convergence thoery is not
unprecedented; for example, Ogden (2017); Kleijn and van der Vaart (2012).
Asymptotics in spatial problems are especially challenging (Zhang, 2004;
Baghishania and Mohammadzadeh, 2012; Simpson et al., 2016; Fuglstad et al.,
2019).

Model-specific challenges are not addressed by Theorem 1. We emphasize that
(a) Theorem 1 reduces the task of establishing posterior convergence to that
of verifying whether the three _individual_ approximations converge for a
given model, by addressing all the complexities associated with applying these
approximations in a non-standard manner, and (b) the problems with nested
approximations that Theorem 1 addresses are all present when fitting LGMs
using INLA, and remain unaddressed in that context.

Some results exist that are relevant to the examples given in §5 and §6. Ogden
(2017) shows that Assumption 7 holds for a Bernoulli Generalized Linear Mixed
Model of the type we consider in §5, and reinforces the claim that the model-
specific calculations required are often non-trivial. Ogden (2021) expands on
the general application of Laplace approximations in models where
$\text{dim}(\boldsymbol{w})$ increases with $n$, which is also relevant to the
examples of §6.2 and §6.3. Lindgren et al. (2011) analyze the error in the
basis function approximation used in §6.1, and Simpson et al. (2016) discuss
convergence of the approximate likelihood used in that example.

##  5 Computational Comparisons

In this section we compare our procedure (Algorithm 1) to the R-INLA
implementation of the INLA method of Rue et al. (2009) for a model which both
can fit, and compare both to MCMC using rstan (Carpenter et al., 2017). We
find that our procedure, which does not make inferences based on the modified
model (§1.3, §2.2), runs faster than R-INLA (Table 1), and that the difference
becomes more pronounced at larger sample sizes.

We use the testing version of R-INLA accessed on 2021/10/22 with the PARDISO
sparse matrix library. Computations were done on a private virtual server with
180 Gb of memory and 16 cores. We parallelize our likelihood computations
across $16$ cores. R-INLA and PARDISO were run in parallel using the
num.threads and num.blas.threads options using $16$ cores as well, except when
$n>1,000,000$ when a single core was used to avoid crashing. We set
strategy=’gaussian’ in R-INLA, yielding an approximation most similar to
$\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})$ out of the available
approximations, and hence the fairest comparison of computation time between
the two approaches. We implement MCMC using the self-tuning No-U-Turn sampler
under the default control parameters (Hoffman and Gelman, 2014) in the rstan
language (Carpenter et al., 2017), through the tmbstan package (Monnahan and
Kristensen, 2018), using the same TMB template as we use to implement
Algorithm 1.

We perform our comparisons by fitting a Bernoulli Generalized Linear Mixed
Model to a large dataset. A complete record of every person who was discharged
from a publicly funded drug or alcohol treatment facility in the United States
between 2006 – 2011 and was listed as having used opioids was obtained from
the Inter-University Consortium for Political and Social Research through the
University of Michigan (ICPSR, 2015). The dataset contains person-level
records of the year of admission, reason for discharge (including successful
completion of the program), and gender, race, and living arrangement
(homeless, independent or dependent), for each such person. There are
$n=7,283,575$ records after filtering out missing data. We quantify the
association between the probability of a subject successfully completing their
program and their gender, race, and living arrangements, using discrete,
nested random effects to account for correlation between outcomes from
subjects living in the same state and the same town.

Denote by $y_{i}$ the observed indicator of whether the $i^{th}$ subject
completed their treatment, $i=1\ldots n$, and the vector of covariates for
each subject by $\boldsymbol{x}_{i}$. There are $d_{1}=47$ states and
$d_{2}=262$ towns with reported data, and we denote by $z_{i1}\in[d_{1}]$ the
state and by $z_{i2}\in[d_{2}]$ the town in which the $i^{th}$ subject lives.
We consider the following model:

$\displaystyle Y_{i}|\eta_{i}$
$\displaystyle\overset{ind}{\sim}\text{Bernoulli}(p_{i}),i\in[n],$ (8)
$\displaystyle\log\left(\frac{p_{i}}{1-p_{i}}\right)$
$\displaystyle=\eta_{i}=\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}+u_{1}(z_{i1})+u_{2}(z_{i2}),$
$\displaystyle u_{1j}|\sigma_{1}$
$\displaystyle\overset{iid}{\sim}\text{N}(0,\sigma_{1}^{2}),j\in[d_{1}],$
$\displaystyle u_{2t}|\sigma_{2}$
$\displaystyle\overset{iid}{\sim}\text{N}(0,\sigma_{2}^{2}),t\in[d_{2}],$

where $u_{1}(z_{i1})=u_{1j}$ whenever $z_{i1}=j,j\in[d_{1}]$, and
$u_{2}(z_{i2})=u_{2t}$ whenever $z_{i2}=t,t\in[d_{2}]$, so that
$\boldsymbol{u}_{1}=(u_{1j})_{j\in[d_{1}]},\boldsymbol{u}_{2}=(u_{1t})_{t\in[d_{2}]}$,
and $\boldsymbol{u}=(\boldsymbol{u}_{1},\boldsymbol{u}_{2})\in\mathbb{R}^{d}$
where $d=d_{1}+d_{2}$. The regression coefficients all correspond to
categorical predictors, with a total of $p=7+1$ non-reference levels and an
intercept. The index sets are $\mathcal{J}_{i}=\left\\{i\right\\},i\in[n]$ and
hence this model is a LGM, compatible with the methodology of Rue et al.
(2009). The full parameter vector is
$\boldsymbol{w}=(\boldsymbol{u},\boldsymbol{\beta})$, and
$\text{dim}(\boldsymbol{w})=d_{1}+d_{2}+p=47+262+8=317$. In contrast, when
running R-INLA, inferences are made based on the modified model
$\boldsymbol{w}^{*}=(\boldsymbol{\eta}^{*},\boldsymbol{w})$, having
$\text{dim}(\boldsymbol{w}^{*})=317+n=7,283,892$ for the full dataset,
resulting in much larger matrices to be stored and decomposed during fitting.

We fit the model to subsamples of the dataset ranging from a size of $n=1,000$
to $n=1,000,000$. Ten different subsamples were taken of each size, and the
computing times and posterior distributions were compared. The full dataset
could not be fit with rstan, and caused R-INLA to crash except when restricted
to a single core. Table 1 shows the mean and standard deviation of computing
times for each scenario. The ELGM procedure (Algorithm 1) is faster than
R-INLA, with the difference increasing as the size of the dataset increases.
For $n=1,000,000$, the ELGM is roughly six times faster, and the full dataset
with $n=7,283,575$ is $15$ times faster with ELGM, with the caveat that R-INLA
ran in single core mode.

Table 2 compares the similarity of the posterior distributions for the three
inference methods, as measured by the KS statistic (the maximum distance
between empirical CDFs). The ELGM is, for the most part, closer to rstan than
R-INLA for the $\sigma$ parameters, and the posterior for the $\beta$
parameters are broadly comparable. For $n=1,000,000$ observations, ELGM is
perceptible closer to rstan than R-INLA for the $\sigma$ parameters, and the
discrepancy between R-INLA and ELGM is substantial when run on the full
dataset. The posterior for $\sigma_{1}$ appears to be over-estimated by R-INLA
when using the full dataset, as the posterior median of $3.19$ and $95\%$
credible interval of $(2.86,3.46)$ is well above the R-INLA estimate from
$1,000,000$ observations of $2.14$ $(1.35,3.29)$, the ELGM estimates of $1.88$
$(1.17,2.88)$ and $2.14$ $(1.33,3.29)$ for the full and reduced datasets
respectively, and rstan’s estimates from the reduced dataset of $2.14$
$(1.35,3.30)$.

| Time, seconds | Iterations
---|---|---
$n$ | ELGM | R-INLA | ELGM | R-INLA
1,000 | 1.25 (0.235) | 2.57 (0.230) | 308 | 630
10,000 | 2.17 (0.193) | 4.75 (0.157) | 176 | 385
100,000 | 9.44 (0.645) | 27.7 (2.47) | 36.3 | 107
1,000,000 | 100 (4.75) | 641 (18.5) | 10.2 | 65.3
7,283,575 | 926 (5.49) | 14,169 | – | –
Table 1: Mean (std. dev) run times in seconds for $10$ runs each of ELGM and R-INLA, run in parallel using $16$ cores, on subsets of the opioid treatment centre data of §5. For $n=7,283,575$, R-INLA was run once, single-threaded (to avoid crashing), and rstan was not run, to avoid excessive computation times. Also shown are the number of iterations of rstan that could have been run in the same time it took to fit ELGM and R-INLA, based on a single run of $1,000$ iterations of rstan at each sample size. $n$ | Methods | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ | $\beta_{5}$ | $\beta_{6}$ | $\beta_{7}$ | $\sigma_{1}$ | $\sigma_{2}$
---|---|---|---|---|---|---|---|---|---|---|---
$1,000$ |  ELGM/rstan | 0.036 | 0.011 | 0.021 | 0.026 | 0.022 | 0.048 | 0.033 | 0.102 | 0.046 | 0.053
R-INLA/rstan | 0.067 | 0.012 | 0.043 | 0.015 | 0.036 | 0.014 | 0.028 | 0.063 | 0.051 | 0.093
ELGM/R-INLA | 0.099 | 0.017 | 0.050 | 0.021 | 0.021 | 0.058 | 0.056 | 0.048 | 0.019 | 0.047
$10,000$ |  ELGM/rstan | 0.019 | 0.016 | 0.031 | 0.035 | 0.009 | 0.040 | 0.023 | 0.031 | 0.017 | 0.069
R-INLA/rstan | 0.020 | 0.022 | 0.010 | 0.012 | 0.020 | 0.009 | 0.017 | 0.031 | 0.011 | 0.101
ELGM/R-INLA | 0.027 | 0.018 | 0.031 | 0.030 | 0.023 | 0.046 | 0.012 | 0.056 | 0.012 | 0.043
$100,000$ |  ELGM/rstan | 0.049 | 0.013 | 0.023 | 0.007 | 0.025 | 0.024 | 0.012 | 0.028 | 0.013 | 0.016
R-INLA/rstan | 0.026 | 0.024 | 0.017 | 0.011 | 0.018 | 0.014 | 0.023 | 0.017 | 0.013 | 0.031
ELGM/R-INLA | 0.037 | 0.021 | 0.014 | 0.013 | 0.025 | 0.019 | 0.016 | 0.031 | 0.015 | 0.031
$1,000,000$ |  ELGM/rstan | 0.017 | 0.009 | 0.016 | 0.009 | 0.011 | 0.010 | 0.010 | 0.017 | 0.018 | 0.009
R-INLA/rstan | 0.040 | 0.011 | 0.012 | 0.012 | 0.009 | 0.008 | 0.013 | 0.010 | 0.016 | 0.035
ELGM/R-INLA | 0.041 | 0.009 | 0.010 | 0.008 | 0.015 | 0.012 | 0.014 | 0.025 | 0.030 | 0.038
$7,283,575$ |  ELGM/R-INLA | 0.329 | 0.019 | 0.016 | 0.012 | 0.014 | 0.022 | 0.011 | 0.011 | 0.957 | 0.329
Table 2: KS statistic (maximum distance between empirical CDFs) for ELGM,
R-INLA, and rstan, for the ten marginal posterior distributions at each sample
size, for the Opioid example of §5. The ELGM results are generally as close to
rstan as R-INLA, but computed in a fraction of the time (Table 1).

##  6 Examples of Extended Latent Gaussian Models

In this section we demonstrate the breadth of the ELGM class and practical
utility of our approximate inference method for it through three challenging
examples. We fit a point process model to spatially-aggregated data (Li et
al., 2012; Nandi et al., 2020) and compare the accuracy of our approach to
MCMC. We fit a Cox Proportional Hazards model with partial likelihood for
mapping spatial variation in Leukaemia survival times, which is an example of
a model in which the Hessian of the log-likelihood is fully dense and hence
the most computationally-intensive type of ELGM. We then fit the Galactic Mass
Estimator model of Eadie et al. (2017) and Eadie and Juric (2019) for
estimating the mass of the Milky Way galaxy in the presence of multivariate
measurement uncertainties, a challenging model. All three examples are beyond
the class of LGMs but belong to the class of ELGMs. The computations make use
of the aghq package in the R language, and code for all examples is available
from https://github.com/awstringer1/elgm-paper-code.

###  6.1 Spatially Aggregated Point Process Data

A spatial point process is a stochastic process which generates random points
within a fixed study area. Often, for reasons of privacy or lack of available
resources, the exact point locations are not recorded and are instead
aggregated to counts within predefined regions within this area. In such
cases, covariate information may be available at a different spatial
resolution than the response, complicating inference. Spatial _downscaling_ or
_disaggregation_ is a challenging task where an aggregated response is
combined with fine-scale covariate information to infer spatial variation in a
phenomenon of interest at a higher resolution than the observed counts. This
task has been approached using frequentist methods (Lee et al., 2017; Brown
and Stafford, 2021) and Bayesian inference using specialized MCMC algorithms
(Li et al., 2012; Diggle et al., 2013).

More recently, Nandi et al. (2020) develop software for approximate Bayesian
inference in these models where a single Laplace approximation (corresponding
to adaptive Gauss-Hermite quadrature with $k=1$ point) is used for the joint
posterior of all model parameters, making use of the TMB package (Kristensen
et al., 2016). They describe an example of inferring Malaria incidence in
Madagascar using aggregated case counts and fine-scale environmental
covariates. We fit this example here within the ELGM framework, using more
than one point for the adaptive quadrature, and obtaining more accurate
approximations to the posteriors when compared to a long MCMC run.

Denote the study region of Madagascar by $\mathcal{M}\subset\mathbb{R}^{2}$
and define an inhomogenous Poisson point process
$X=\left\\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathcal{M}\right\\}$ with rate
function $\lambda(\cdot)$. We do not observe point locations
$X(\boldsymbol{s})$, and instead observe case counts $\boldsymbol{y}$ with
$y_{i}=\left\lVert\boldsymbol{s}\in\mathcal{M}:X(\boldsymbol{s})\in
S_{i}\right\rVert$ for predefined regions $S_{i}\subset\mathcal{M}$ with
$\cup_{i=1}^{n}S_{i}=\mathcal{M}$ and $S_{i}\cap S_{j}=\emptyset$ for $i\neq
j\in[n]$. For any $\boldsymbol{s}\in\mathcal{M}$ define spatially-varying
covariates $\boldsymbol{x}(\boldsymbol{s}):\mathbb{R}^{2}\to\mathbb{R}^{p}$
corresponding to elevation, vegetation index and land surface temperature,
with $p=4$ including the intercept, and known population offset function
$P(\boldsymbol{s}):\mathbb{R}^{2}\to\mathbb{R}$ such that the population in
any fixed area $S\subseteq\mathcal{M}$ is
$\int_{S}P(\boldsymbol{s})d\boldsymbol{s}$. We wish to infer incidence risk
$\lambda(\cdot)$ using the spatially-aggregated log-Gaussian Cox process:

$\displaystyle Y_{i}|\lambda(\cdot),v_{i}$
$\displaystyle\overset{ind}{\sim}\text{Poisson}\left[\exp(v_{i})\int_{S_{i}}P(\boldsymbol{s})\lambda(\boldsymbol{s})d\boldsymbol{s}\right],\
\log\lambda(\boldsymbol{s})=\beta_{0}+\boldsymbol{x}(\boldsymbol{s})^{\texttt{T}}\boldsymbol{\beta}+u(\boldsymbol{s}),$
$\displaystyle u(\boldsymbol{s})$
$\displaystyle\sim\mathcal{GP}\left\\{0,\text{M}_{\nu}(\cdot;\sigma,\rho)\right\\},\
\text{Cov}\\{u(\boldsymbol{s}+\boldsymbol{h}),u(\boldsymbol{s})\\}=\text{M}_{\nu}(\left\lVert\boldsymbol{h}\right\rVert;\sigma,\rho),\boldsymbol{s}\in\mathcal{M},\boldsymbol{h}\in\mathcal{M},$
$\displaystyle v_{i}$ $\displaystyle\overset{iid}{\sim}\text{N}(0,\tau^{-1}),\
\tau^{-1/2}\sim\text{Exp}(\xi_{\tau}),\sigma\sim\text{Exp}(\xi_{\sigma}),\rho^{-1}\sim\text{Exp}(\xi_{\rho}).$

Here $\text{M}_{\nu}(\cdot;\sigma,\rho)$ is a Matern covariance function with
fixed shape parameter $\nu=1$ in the parametrization of Brown (2015). The
polygon-level effects $v_{i}$ are included to account for overdispersion in
the polygon counts.

The risk surface $\lambda(\boldsymbol{s})$ must be discretely inferred, and
this requires approximations to the continuously-defined covariates
$\boldsymbol{x}(\boldsymbol{s})$ and spatial process $u(\boldsymbol{s})$. For
the spatial process we use a basis-function representation (Lindgren et al.,
2011):

$u(\boldsymbol{s})\approx\sum_{j=1}^{d}\phi_{j}(\boldsymbol{s})u_{j},$ (9)

where $\phi_{j}(\cdot)$ are fixed, known piecewise-linear functions. The
unknown parameter is $\boldsymbol{u}=(u_{j})_{j\in[d]}$. This approximation is
constructed on a _mesh_ with triangular cells $\mathcal{M}_{j},j\in[d]$ such
that $\mathcal{M}\subseteq\cup_{j=1}^{d}\mathcal{M}_{j}$ and
$\mathcal{M}_{j}\cap\mathcal{M}_{k}=\emptyset,j\neq k\in[d]$ (Lindgren and
Rue, 2015). See Lindgren et al. (2011) and Lindgren and Rue (2015) for further
details.

The covariates are observed on a fine grid with cells
$\mathcal{Q}_{t},t\in[T]$ such that
$\cup_{t=1}^{T}\mathcal{Q}_{t}=\mathcal{M}$ and
$\mathcal{Q}_{t}\cap\mathcal{Q}_{l}=\emptyset$ for $t\neq l\in[T]$. The
population counts in each grid cell
$P_{it}\equiv\int_{S_{i}\cap\mathcal{Q}_{t}}P(\boldsymbol{s})d\boldsymbol{s}$
are available (Figure 1(c)). Further approximate
$x(\boldsymbol{s})\approx\boldsymbol{x}(\mathcal{Q}_{t})$ for
$\boldsymbol{s}\in\mathcal{Q}_{t}$, where $\boldsymbol{x}(\mathcal{Q}_{t})$ is
shorthand for taking $\boldsymbol{s}$ to be the centre of the cell
$\mathcal{Q}_{t}$. Inferences for $\lambda(\cdot)$ are made on the grid
$\mathcal{Q}_{t}$ by defining
$\lambda(\boldsymbol{s})\approx\lambda(\mathcal{Q}_{t})$ for each
$\boldsymbol{s}\in\mathcal{Q}_{t}$. Finally, define
$\boldsymbol{\lambda}=\left\\{\lambda(\mathcal{Q}_{t}):t\in[T]\right\\}$.
Figure 1 shows the observed data at the polygon level and on the grid, as
appropriate. The discrete model is then:

$\displaystyle y_{i}|\boldsymbol{\lambda}_{\mathcal{J}_{i}},v_{i}$
$\displaystyle\overset{ind}{\sim}\text{Poisson}\left(\exp(v_{i})\sum_{t\in\mathcal{J}_{i}}P_{it}\lambda(\mathcal{Q}_{t})\right),\
\mathcal{J}_{i}=\left\\{t\in[T]:S_{i}\cap\mathcal{Q}_{t}\neq\emptyset\right\\},$
$\displaystyle\eta_{t}$
$\displaystyle=\log\lambda(\mathcal{Q}_{t})=\beta_{0}+\boldsymbol{x}(\mathcal{Q}_{t})^{\texttt{T}}\boldsymbol{\beta}+u(\mathcal{Q}_{t}),\
u(\mathcal{Q}_{t})=\sum_{j=1}^{d}\phi_{j}(\mathcal{Q}_{t})u_{j},$
$\displaystyle\boldsymbol{U}|\boldsymbol{\theta}_{2}$
$\displaystyle\sim\text{Normal}\left[\boldsymbol{0},\boldsymbol{\Sigma}(\boldsymbol{\theta}_{2})\right],\
v_{i}\overset{iid}{\sim}\text{N}(0,\tau^{-1}),\
\tau^{-1/2}\sim\text{Exp}(\xi_{\tau}),\sigma\sim\text{Exp}(\xi_{\sigma}),\rho^{-1}\sim\text{Exp}(\xi_{\rho})$

where $\boldsymbol{\Sigma}^{-1}(\boldsymbol{\theta}_{2})$ is the sparse
approximation to the precision matrix of the weights $\boldsymbol{u}$ on the
mesh (Lindgren et al., 2011). Following Nandi et al. (2020) the nonlinear
parameters $\boldsymbol{\theta}=(\log\tau,\log\sigma,\log\rho)$ are assigned
priors such that $\tau^{-1/2},\sigma$ and $\rho^{-1}$ follow independent
Exponential distributions with parameters chosen to satisfy
$\mathbb{P}(\tau>0.1)=\mathbb{P}(\sigma>0.1)=\mathbb{P}(\rho<3)=0.01$ (which
are Penalized Complexity priors of Simpson et al. 2017 and Fuglstad et al.
2019), as well as $\beta_{0}\sim\text{N}(0,2^{2})$ and
$\beta_{j}\sim\text{N}(0,0.4^{2}),j\in[p]$ independently.

We fit this model using our ELGM approach (Algorithm 1) with $k=7$ quadrature
points, adapting the excellent software of Nandi et al. (2020) for use with
the aghq package with a computation time of about $27$ minutes. We also fit
the model using an MCMC run of $16$ chains with $5,000$ iterations (including
a warmup of $1,000$) each for a total of $64,000$ usable iterations, using the
tmbstan package (Monnahan and Kristensen, 2018) at a wall computation time of
approximately $29$ hours when running the chains in parallel. Figure 2 shows
the posterior mean incidence rates $\lambda(\cdot)$ obtained from $100$
independent posterior samples from the approximate posterior and from MCMC, as
well as mean excess spatial variation $u(\cdot)$ and probabilty of exceeding
$20\%$ risk, $\mathbb{P}[\lambda(\cdot)>0.2|\boldsymbol{Y}]$. Table 3 shows
the posterior mean and standard deviation for the regression coefficients and
nonlinear parameters computed using the ELGM procedure and MCMC. Finally,
Figure 3 shows the posteriors for $(\tau^{-1/2},\sigma,\rho)$ computed using
our approach, the full Laplace approximation, and MCMC. The ELGM approach
generally matches the MCMC output closely, and in the case of the posteriors
for the nonlinear parameters (Figure 3), appears more accurate than the full
Laplace approximation of Nandi et al. (2020).

(a) Malaria case counts $Y_{i}$

(b) Population counts $P_{i}$

(c) Population raster $P(\mathcal{Q}_{t})$

(d) Elevation $x_{1}(\mathcal{Q}_{t})$

(e) Vegetation index $x_{2}(\mathcal{Q}_{t})$

(f) Surface temp. $x_{3}(\mathcal{Q}_{t})$

Figure 1: Case counts $Y_{i}$, population counts within each polygon
$P_{i}=\int_{S_{i}}P(\boldsymbol{s})d\boldsymbol{s}$ and at the pixel level
$P(\mathcal{Q}_{t})$, and covariate rasters
$\boldsymbol{x}(\mathcal{Q}_{t})=\left\\{1,x_{1}(\mathcal{Q}_{t}),x_{2}(\mathcal{Q}_{t}),x_{3}(\mathcal{Q}_{t})\right\\},t\in[T]$
for Malaria incidence in Madagascar, §6.1.

(a) $\mathbb{E}\left\\{\lambda(\mathcal{Q}_{t})|\boldsymbol{Y}\right\\}$, ELGM

(b) $\mathbb{E}\left\\{\lambda(\mathcal{Q}_{t})|\boldsymbol{Y}\right\\}$, MCMC

(c) $\mathbb{E}\left\\{u(\mathcal{Q}_{t})|\boldsymbol{Y}\right\\}$, ELGM

(d) $\mathbb{E}\left\\{u(\mathcal{Q}_{t})|\boldsymbol{Y}\right\\}$, MCMC

(e) $\mathbb{P}\left\\{\lambda(\mathcal{Q}_{t})>0.2|\boldsymbol{Y}\right\\}$,
ELGM

(f) $\mathbb{P}\left\\{\lambda(\mathcal{Q}_{t})>0.2|\boldsymbol{Y}\right\\}$,
MCMC

Figure 2: Posterior summaries from the ELGM procedure and MCMC on the grid
$\mathcal{Q}_{t},t\in[T]$, for the Malaria example of §6.1.

(a) $\pi(\tau^{-1/2}|\boldsymbol{Y})$

(b) $\pi(\sigma|\boldsymbol{Y})$

(c) $\pi(\rho|\boldsymbol{Y})$

Figure 3: Approximate posteriors for $\tau^{-1/2},\sigma,\rho$ using the ELGM procedure (—) and the Laplace approximation of Nandi et al. (2020) (- - -), along with posterior samples from MCMC (${\color[rgb]{.75,.75,.75}\blacksquare}$) and prior ($-\cdot-$), for the Malaria example of §6.1. | Mean | SD | $2.5\%$ | $97.5\%$
---|---|---|---|---
Param. | ELGM | MCMC | ELGM | MCMC | ELGM | MCMC | ELGM | MCMC
$\beta_{0}$ | -3.054 | -3.15 | 0.363 | 0.348 | -3.75 | -3.84 | -2.29 | -2.48
$\beta_{\texttt{elev}}$ | -0.465 | -0.464 | 0.187 | 0.188 | -0.825 | -0.831 | -0.096 | -0.091
$\beta_{\texttt{vege}}$ | 0.385 | 0.384 | 0.204 | 0.206 | -0.017 | -0.016 | 0.784 | 0.791
$\beta_{\texttt{lst}}$ | 0.203 | 0.185 | 0.263 | 0.265 | -0.318 | -0.334 | 0.713 | 0.703
$\tau^{-1/2}$ | 0.660 | 0.652 | 0.088 | 0.0848 | 0.502 | 0.493 | 0.847 | 0.827
$\sigma$ | 1.02 | 1.03 | 0.198 | 0.194 | 0.668 | 0.693 | 1.45 | 1.46
$\rho$ | 2.57 | 2.47 | 0.984 | 0.833 | 1.32 | 1.31 | 5.00 | 4.54
Table 3: Coefficient and nonlinear parameter posterior summaries from the ELGM
procedure and MCMC, for the Malaria example of §6.1. Computation time for the
ELGM was $27$ minutes, and for MCMC $29$ hours.

###  6.2 Spatial survival regression with partial likelihood

The Cox Proportional Hazards model with partial likelihood is a standard model
for survival analysis problems. The partial (log) likelihood has a Hessian
matrix which is completely dense, and hence is not compatible with the LGM
framework, and is in fact the most computationally-intensive type of ELGM. We
are unaware of any other approach that has been used to make Bayesian
inferences based a Cox Proportional Hazards model with both a partial
likelihood and a continuously-indexed latent spatial process, which we do in
this example, and therefore this example serves to illustrate the breadth of
our proposed method.

We consider a classic dataset of Leukaemia survival times in northern England
where the goal is to infer the spatial variation in survival time. These data
have been analyzed using survival models with parametric hazard functions
(Henderson et al., 2002; Lindgren et al., 2011), resulting in simpler
computations than when using partial likelihood. More recently Martino et al.
(2011) adapt the methodology of Rue et al. (2009) to analyze these data using
a semi-parametric hazard model which is still less computationally intensive
than using partial likelihood, and where spatial variation is restricted to
occur between predefined geographical regions. Unlike these previous analyses,
we fit a Cox Proportional Hazards model using partial likelihood and a
continuously-varying spatial model using the exact observed point locations,
which is feasible in the ELGM framework introduced in this paper.

The data consist of the survival times of $n=1043$ Leukaemia patients in
northern England. Of these $n_{0}=879$ were observed to die during the study
period and $n-n_{0}=164$ were right-censored. Denote the study region by
$\mathcal{M}\subset\mathbb{R}^{2}$. Each subject has a point location of
residence $\boldsymbol{s}_{i}\in\mathcal{M}$, and covariate vector
$\boldsymbol{x}_{i}\in\mathbb{R}^{p}$ with $p=4$ containing age, sex, white
blood cell count (WBC), and the Townsend Deprivation Index (TPI), a measure of
social deprivation. Let $y_{i},i\in[n]$ denote the survival time of the
$i^{th}$ subject. For convenience suppose $0<y_{1}<\cdots<y_{n_{0}}$ are
observed and $(y_{i})_{i=n_{0}+1,\ldots,n}$ are censored. We use the following
hierarchical model:

$\displaystyle\pi\left(\boldsymbol{y}|\boldsymbol{\lambda}\right)$
$\displaystyle=\prod_{i=1}^{n_{0}}\frac{\lambda_{i}}{\sum_{j=i}^{n}\lambda_{j}},\
\eta_{i}=\log\lambda_{i}=\beta_{0}+\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}+u(\boldsymbol{s}_{i}),i\in[n],\
\mathcal{J}_{i}=\left\\{i,\ldots,n\right\\},$ $\displaystyle
u(\boldsymbol{s})$
$\displaystyle\sim\mathcal{GP}\left\\{0,\text{M}_{\nu}(\cdot;\sigma,\rho)\right\\},\
\text{Cov}[u(\boldsymbol{s}+\boldsymbol{h}),u(\boldsymbol{s})]=\text{M}_{\nu}(\left\lVert\boldsymbol{h}\right\rVert;\sigma,\rho),\boldsymbol{s}\in\mathcal{M},\boldsymbol{h}\in\mathbb{R}^{2}.$

Because $\mathcal{J}_{1}=[n]$ we have $k,l\in\mathcal{J}_{1}$ for every
$k,l\in[n]$ and hence
$\boldsymbol{C}_{\boldsymbol{\eta}}(\boldsymbol{\eta},\boldsymbol{\theta}_{1})$
is a fully dense matrix (§3).

For the spatially-varying Gaussian process $u(\cdot)$ we use a Matern
covariance function as in §6.1, however we do not use a piecewise constant
approximation to $u(\cdot)$ with the sparse approximation to its precision
matrix, and instead evaluate the dense covariance matrix at the observed point
locations. The unknown spatial parameter is therefore
$\boldsymbol{u}=(u(\boldsymbol{s}_{i}))_{i\in[n]}$ with dimension $d=n$. The
model for $\boldsymbol{u}$ is:

$\displaystyle\boldsymbol{u}|\boldsymbol{\theta}_{2}$
$\displaystyle\sim\text{N}\left\\{\boldsymbol{0},\boldsymbol{\Sigma}(\boldsymbol{\theta}_{2})\right\\},\
\boldsymbol{\Sigma}_{jk}=\text{M}_{\nu}(\left\lVert\boldsymbol{s}_{j}-\boldsymbol{s}_{k}\right\rVert;\sigma,\rho),j,k\in[n],$

where, as in §6.1, $\text{M}_{\nu}(\cdot;\sigma,\rho)$ is a Matern covariance
function with fixed shape $\nu=1$, and standard deviation $\sigma$ and range
$\rho$ in the parametrization of Brown (2015). The nonlinear parameters are
$\boldsymbol{\theta}=(\log\sigma,\log\rho)$. Again as in §6.1 we follow
Fuglstad et al. (2019) and choose independent Exponential priors on $\sigma$
and $\rho^{-1}$ which satisfy
$\mathbb{P}(\sigma>1)=\mathbb{P}(\rho<20\text{km})=0.5$.

To approximate the joint predictive distribution
$\boldsymbol{u}^{*}|\boldsymbol{y}$ of
$\boldsymbol{u}^{*}=(u(\boldsymbol{s}^{*}_{l}))_{l\in[L]}$ for any new
locations
$(\boldsymbol{s}^{*}_{l})_{l\in[L]}\subset\mathcal{M},L\in\mathbb{N}$, we make
use of the full joint posterior approximation by first drawing
$\boldsymbol{w}=(\boldsymbol{u},\boldsymbol{\beta})$ from
$\widetilde{\pi}(\boldsymbol{w}|\boldsymbol{y})$ according to Algorithm 1. We
then draw from $\boldsymbol{u}^{*}|\boldsymbol{u}$ using methods implemented
in the RandomFields package (Schlather et al., 2015).

Table 4 shows the posterior medians and credible regions for the coefficients
of age, sex, Cell Count (WBC) and Deprivation (TPI), along with those for the
the nonlinear parameters. The use of partial likelihood means that these are
obtained without assumptions on the form of the baseline hazard. Figure 4
shows the approximate posterior distributions of the two nonlinear parameters
in the Matern covariance function as well as the approximate posterior mean
and $20\%$ exceedence probabilities obtained from $100$ draws from
$\boldsymbol{u}^{*}|\boldsymbol{y}$ for a fine grid of points
$\boldsymbol{s}^{*}_{1},\ldots,\boldsymbol{s}^{*}_{L},L=58,215$ obtained by
laying a $400\times 200$ square over the study area and taking the
intersection of this square with the map. The use of exact point locations and
simulation on a fine grid leads to a higher-resolution estimate than that of
Henderson et al. (2002) and Lindgren et al. (2011), though the overall pattern
appears to match that reported by these previous analyses.

(a) $\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\sigma|\boldsymbol{y})$

(b) $\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\rho|\boldsymbol{y})$

(c) $\mathbb{E}[\exp\left\\{u(\boldsymbol{s})\right\\}|\boldsymbol{y}]$

(d) $\mathbb{P}[\exp\left\\{u(\boldsymbol{s})\right\\}>1.2|\boldsymbol{y}]$

Figure 4: (Top) Priors (- - -) and approximate posterior distributions (—) for (a) the spatial standard deviation and (b) range (KM, density $\times 10^{5}$). (Bottom) maps showing (c) point locations of deaths $(\cdot)$ and approximate posterior mean excess spatial variation and (d) approximate posterior probability of exceeding $20\%$ excess spatial variation for the spatial survival analysis of §6.2 Table 4:  Posterior moments and quantiles for fixed effects and nonlinear parameters for the Leukaemia data of Example 6.2. Variable | Mean | SD | $2.5\%$ | $97.5\%$
---|---|---|---|---
Age | 0.029 | 0.002 | 0.025 | 0.033
Male | 0.029 | 0.070 | -0.107 | 0.171
Cell Count | 0.005 | 0.001 | 0.003 | 0.004
Deprivation | 0.027 | 0.010 | 0.009 | 0.046
SD, $\sigma$ | 0.257 | 0.095 | 0.115 | 0.499
Range, $\rho$ (KM) | 44.2 | 32.0 | 9.77 | 134

###  6.3 Estimating the mass profile of the Milky Way with measurement error

In a series of papers, Eadie and Harris (2016), Eadie et al. (2017), Eadie et
al. (2018) and Eadie and Juric (2019) develop and test a model for estimating
the mass of the Milky Way Galaxy using position and velocity measurements of
star clusters in its orbit. They use a physical model which implies a
probability distribution for the position and velocity of star clusters, which
depends on parameters that determine the mass of the Galaxy at any radial
distance from its centre. This distribution is used in a hierarchical model
which incorporates measurement uncertainties for the multivariate response,
and missing data. Inference is based on posterior distributions for the
parameters of interest, upon which estimates and uncertainty quantification
for the mass of the Galaxy are based.

Throughout these papers, posteriors are computed using MCMC which requires
extensive user tuning (Eadie et al., 2017). Here we present similar results to
Eadie et al. (2017) using the ELGM procedure with a running time faster than
that reported by Eadie et al. (2017) and no user tuning, although we simplify
some aspects of inference for purposes of illustration.

We describe only the relevant statistical details of this example; the
physical details are well beyond the scope of this paper. Let
$\boldsymbol{y}_{i}=(y_{i,1},y_{i,2},y_{i,3},y_{i,4})$ denote the four
kinematic measurements taken for the $i^{th}$ star cluster: position, line-of-
sight velocity, proper motion in right-ascention corrected for declination,
and proper motion in declination, and let
$\boldsymbol{y}=(\boldsymbol{y}_{i})_{i\in[n]}$. These heliocentric
measurements are taken with respect to the position and motion of the sun, and
are converted to galactocentric position and motion relative to the centre of
the Galaxy using a deterministic, nonlinear transformation $\Omega(\cdot)$.
The probability distribution relating the measurements to the mass of the
Galaxy is defined in this galactocentric frame of reference. Prior to
transformation, each heliocentric measurement is subject to a random
measurement error $\boldsymbol{u}_{i}$ with mean zero and fixed standard
deviations $\boldsymbol{\delta}_{i}$ which are reported as part of the
measurement process.

For simplicity we consider only clusters with complete kinematic measurements
and include measurement errors on the position and line-of-sight velocity
only, defining
$\boldsymbol{u}_{i}=(u_{i1},u_{i2}),\boldsymbol{w}=(\boldsymbol{u}_{i})_{i\in[n]}\in\mathbb{R}^{2n}$,
and letting
$\boldsymbol{y}_{i}^{*}=(y_{i,1}+u_{i,1},y_{i,2}+u_{i,2},y_{i,3},y_{i,4})$ for
$i\in[n]$. There are $n=70$ star clusters with complete data. A probability
distribution
$\pi\left\\{\Omega(\boldsymbol{y}_{i})|\boldsymbol{\theta}_{1}\right\\}$ over
the galactocentric measurements $\Omega(\boldsymbol{y}_{i}),i\in[n]$ depends
on parameters $\boldsymbol{\theta}_{1}=(\Psi_{0},\gamma,\alpha,\beta)$, and
the mass of the Galaxy at radial distance $r$ kiloparsecs (kpc) from its
centre, $M_{\boldsymbol{\theta}_{1}}(r)=\gamma\Psi_{0}r^{1-\gamma}$, is the
object of inferential interest. Parameter transformations and strongly-
informative priors on $\boldsymbol{\theta}_{1}$ are set according to the
extensive discussion in Eadie and Harris (2016). For further detail, refer to
Eadie and Harris (2016) and Eadie et al. (2017).

The model is as follows:

$\displaystyle\pi(\boldsymbol{y}|\boldsymbol{w},\boldsymbol{\theta}_{1})$
$\displaystyle=\prod_{i=1}^{n}\pi\left\\{\Omega(\boldsymbol{y}_{i}^{*})|\boldsymbol{\theta}_{1}\right\\},\
\mathcal{J}_{i}=\left\\{\left\\{i,1\right\\},\left\\{i,2\right\\}\right\\},\
\boldsymbol{u}_{i}\overset{ind}{\sim}\text{N}\left\\{\boldsymbol{0},\text{diag}(\boldsymbol{\delta}_{i})\right\\},i\in[n],$
(10)

where $i\in[n]$. The measurement standard deviations
$\boldsymbol{\delta}_{i}=(\delta_{i1},\delta_{i2}),i\in[n]$ are reported for
each cluster and taken as fixed and known. The strongly-informative priors are
$\Psi_{0}\sim\text{Unif}(1,200),\gamma\sim\text{Unif}(0.3,0.7),\alpha-3\sim\text{Gamma}(1,4.6)$
and $\beta\sim\text{Unif}(-0.5,1)$ (Eadie and Harris, 2016). The index sets
$\mathcal{J}_{i},i\in[n]$ have $\left|\mathcal{J}_{i}\right|=2>1$ for each
$i\in[n]$, because the observations are multivariate. This is therefore an
example of an ELGM compatible with our methodology, but not an LGM.

(a) $\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\Psi_{0}|\boldsymbol{Y})$

(b) $\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\gamma|\boldsymbol{Y})$

(c) $\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\alpha|\boldsymbol{Y})$

(d) $\widetilde{\pi}_{\text{\tiny{\tt LA}}}(\beta|\boldsymbol{Y})$

Figure 5: (a) – (d) prior (- - -) and approximate posterior (—) distributions
for the four parameters for the astronomy data of Example 6.3. Figure 6:
Posterior mean mass of the Milky Way Galaxy relative to the mass of the sun as
a function of radial distance from Galaxy centre in kiloparsecs (kps)
($M_{\boldsymbol{\theta}_{1}}(r)$, —), with one- (dark) and two- (light)
standard deviation bands for the astronomy data of Example 6.3.

We fit this model using to Algorithm 1, obtaining results that are broadly
comparable to those reported by Eadie and Harris (2016) and Eadie et al.
(2017) accounting for the priors and data inclusion rules used. Figure 5 shows
the posterior distributions for $\Psi_{0},\gamma,\alpha,\beta$ and Figure 6
shows the estimated cumulative mass profile $M_{\boldsymbol{\theta}_{1}}(r)$
of the Galaxy for chosen values of $r$. While the posterior mean and standard
deviation of $M_{\boldsymbol{\theta}_{1}}(r)$ are straightforward to obtain
using quadrature, an approximation to the posterior density of a nonlinear,
trans-dimensional transformation of $\boldsymbol{\theta}_{1}$ is in general
difficult to obtain using our procedure as no algorithm is available to draw
samples from $\widetilde{\pi}_{\text{\tiny{\tt
LA}}}(\boldsymbol{\theta}_{1}|\boldsymbol{Y})$, so we instead report the
posterior mean along with pointwise one and two standard deviation bands for
each value of $r=1,2,\ldots,150$.

There are several computational challenges in this example which we leave to
future work. We do not incorporate missing data into the hierarchical model as
done by Eadie et al. (2017) and Eadie and Juric (2019). Further, the
parameters $\boldsymbol{w}$ and $\boldsymbol{\theta}_{1}$ are subject to
complicated nonlinear constraints implied by the underlying physics, and these
constraints depend both on $\boldsymbol{\theta}_{1}$ and
$\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}_{1}}$. Advanced methods for
optimization in the presence of such nonlinear constraints are readily
available (see the IPOPT software of Wachter and Biegler (2006) and
corresponding ipoptr R package). However, in the presence of measurement
errors this would require the derivatives of the constraints which involves
differentiating through $\widehat{\boldsymbol{w}}_{\boldsymbol{\theta}_{1}}$
with respect to $\boldsymbol{\theta}_{1}$. This is a challenging task which
has recently been investigated in a related context by Margossian et al.
(2020), and we leave its implementation here to future work.

##  7 Discussion

We have defined a novel class of Extended Latent Gaussian Models and developed
approximate Bayesian inference methodology for this class. The method relies
less on sparse matrix algorithms but also depends on matrices of smaller size
than previous approaches, and we have provided numerical evidence of our model
providing faster run times than R-INLA and MCMC, as well as scaling to large
sample sizes (§5). We have proved that the error in our approximation
converges to zero as the sample size increases, as long as the assumptions
ensuring convergence of the error in each component are satisfied (§4).
Further, we demonstrated three challenging examples of ELGMs and fit them with
our procedure: inference for a continuous spatial field using aggregated point
process data, including a comparison to MCMC (§6.1); a Cox Proportional
hazards model with partial likelihood for mapping the spatial variation in
Leukaemia survival times, for which we are unaware of any other method for
making Bayesian inferences (§6.2); and an astrophysical model for estimating
the mass of the Milky Way galaxy accounting for multivariate measurement
uncertainties (§6.3). The core method (Algorithm 1) is implemented in the open
source aghq package in the R language.

There are a number of compelling avenues for future research. While we are
unaware of any theoretical results pertaining to the finite-sample performance
of Gaussian and Laplace approximations to posterior distributions, several
authors (Rue and Martino, 2007; Rue et al., 2009) have suggested that such
approximations may be inaccurate in finite samples, and this merits further
attention. Rue et al. (2009) and Wood (2020) introduce approximations which
are empirically more accurate than the Gaussian, however they work for
_marginal_ posterior distributions only. In contrast, our procedure depends on
access to fast independent sampling from the approximate _joint_ posterior
(Algorithm 1), which we use to compute complicated posterior summaries in a
straightforward manner (§6.1, §6.2). As noted in the discussion to Rue et al.
(2009), this is a compelling practical advantage of MCMC algorithms, and we
believe it is important to retain this advantage when developing methods for
approximate Bayesian inference.

Another potential method for comparison would be Variational Inference (VI)
methods based on Gaussian approximations (Blei et al., 2017). Gaussian VI
methods find the Gaussian approximation which minimizes a lower bound on the
Kullback-Leibler divergence to the true posterior—a complicated, non-convex
optimization—in contrast to our approach which matches the mode and curvature
of the true posterior. Computationally, our approach differs from VI in that
we require only convex optimizations of a tractable objective function, which
are therefore computationally feasible (§3.2), while VI requires a challenging
non-convex optimization of an objective function whose value and gradient
cannot be calculated analytically, leading to difficulty in implementation. We
provide code for fitting the model of §5 with VI using rstan, finding that the
model either crashes or provides unreliable answers as determined by the
Pareto $\widehat{k}$ diagnostic (Yao et al., 2018; Vehtari et al., 2021;
Vehtari et al., 2017), for a variety of sample sizes and priors. While it is
possible that an expert user could tune the implementation to produce more
reliable answers, our approach does not require such user tuning, again owing
to its reliance on stable, tractable convex optimization. We leave further
comparisons to VI for future work.

To our knowledge, Theorem 1 is the first formal theoretical result pertaining
to a nested approximation, of the type made popular in the INLA method of Rue
et al. (2009). In our proof, we explicitly handle the non-standard application
of three different types of approximations, including a misspecified Gaussian,
a marginal Laplace, and the re-use of adapted quadrature points and weights
for approximating two different integrals. It is possible that with further
work, this technique could be expanded to apply to other nested approximation
algorithms, like the ones proposed by Rue et al. (2009) and Wood (2020).

## Acknowledgements

The authors are grateful for helpful comments provided by two referees and an
associate editor, as well as Blair Bilodeau, Yanbo Tang, and Gwen Eadie. Parts
of this manuscript appeared in the first author’s PhD thesis for which he is
grateful for comments from Gwen Eadie, Linbo Wang, Nancy Reid and Anthony
Davison. This work was supported by the Natural Sciences and Engineering
Research Council of Canada’s Discovery Grant and Postgraduate Scholarships
programs.

## Conflict of Interest

The authors report there are no competing interests to declare.

## References

  * Baghishania and Mohammadzadeh (2012) Baghishania, H. and M. Mohammadzadeh (2012).  Asymptotic normality of posterior distributions for generalized linear mixed models.  Journal of Multivariate Analysis 111, 66–77.
  * Bates and Maechler (2019) Bates, D. and M. Maechler (2019).  Matrix: sparse and dense matrix classes and methods.  R package version 1.2-18.
  * Bilodeau et al. (2021) Bilodeau, B., A. Stringer, and Y. Tang (2021).  Stochastic convergence rates and applications of adaptive quadrature in Bayesian inference.
  * Blei et al. (2017) Blei, D. M., A. Kucukelbirb, and J. D. McAuliffe (2017).  Variational inference: a review for statisticians.  Journal of the American Statistical Association 112(518), 859–877.
  * Braun (2014) Braun, M. (2014).  trustOptim: An R Package for trust region optimization with sparse Hessians.  Journal of Statistical Software 60(4), 1–16.
  * Breslow and Clayton (1993) Breslow, N. E. and D. G. Clayton (1993).  Approximate inference in generalized linear mixed models.  Journal of the American Statistical Association 88(421), 9–25.
  * Brown (2015) Brown, P. (2015).  Model-based Geostatistics the easy way.  Journal of Statistical Software 63(12).
  * Brown and Stafford (2021) Brown, P. and J. Stafford (2021).  A root Gaussian Cox process for spatially-aggregated disease incidence data.  In prepartion.
  * Carpenter et al. (2017) Carpenter, B., A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell (2017).  Stan: a probabilistic programming language.  Journal of statistical software 76(1).
  * Diggle et al. (2013) Diggle, P., P. Moraga, B. Rowlingson, and B. Taylor (2013).  Spatial and spatio-temporal log-Gaussian Cox processes: extending the geostatistical paradigm.  Statistical Science 28(4), 542–563.
  * Diggle and Giorgi (2016) Diggle, P. J. and E. Giorgi (2016).  Model-based Geostatistics for prevalence mapping in low-resource settings.  Journal of the American Statistical Association $\bf{111}$, 1096–1120.
  * Diggle et al. (2010) Diggle, P. J., I. Kaimi, and R. Abellana (2010).  Partial-Likelihood Analysis of Spatio-Temporal Point-Process Data.  Biometrics 66, 347–354.
  * Ding et al. (2018) Ding, H., Y. Lee, I. Sato, and M. Sugiyama (2018).  Variational inference for Gaussian process with panel count data.  arXiv:1803.04232 [stat.ML].
  * Eadie and Harris (2016) Eadie, G. M. and W. E. Harris (2016).  Bayesian mass estimates of the Milky Way: the dark and light sides of parameter assumptions.  The Astrophysical Journal 829(108).
  * Eadie and Juric (2019) Eadie, G. M. and M. Juric (2019).  The cumulative mass profile of the milky way as determined by globular cluster kinematics from Gaia DR2.  The Astrophysical Journal 875(159).
  * Eadie et al. (2018) Eadie, G. M., B. Keller, and W. E. Harris (2018).  Estimating the Milky Way’s mass via hierarchical Bayes: a blind test on MUGS2 simulated galaxies.  The Astrophysical Journal 865(72).
  * Eadie et al. (2017) Eadie, G. M., A. Springford, and W. E. Harris (2017).  Bayesian mass estimates of the Milky Way: including measurement uncertainties with hierarchical Bayes.  The Astrophysical Journal 835(167).
  * Fing et al. (2010) Fing, Y., H. Rue, and J. Wakefield (2010).  Bayesian inference for generalized linear mixed models.  Biostatistics 11(3), 397–412.
  * Fuglstad et al. (2019) Fuglstad, G.-A., D. Simpson, F. Lindgren, and H. Rue (2019).  Constructing priors that penalize the complexity of Gaussian random fields.  Journal of the American Statistical Association 114(525), 445–452.
  * Geyer (2020) Geyer, C. J. (2020).  trust: Trust Region Optimization.  R package version 0.1-8.
  * Hall (2000) Hall, D. B. (2000).  Zero-Inflated Poisson and Binomial Regression with Random Effects: A Case Study.  Biometrics 54(4), 1030 – 1039.
  * Heiss and Winschel (2008) Heiss, F. and V. Winschel (2008).  Likelihood approximation by numerical integration on sparse grids.  Journal of Econometrics 144, 62–80.
  * Henderson et al. (2002) Henderson, R., S. Shimakura, and D. Gorst (2002).  Modelling spatial variation in Leukaemia survival.  Journal of the American Statistical Association 97(460), 965–972.
  * Hoffman and Gelman (2014) Hoffman, M. D. and A. Gelman (2014).  The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo.  Journal of Machine Learning Research 15, 1351–1381.
  * ICPSR (2015) ICPSR (2015).  Treatment episode data set – discharges (teds-d) – concatenated, 2006 to 2011.
  * Jin and Andersson (2020) Jin, S. and B. Andersson (2020).  A note on the accuracy of adaptive Gauss-Hermite quadrature.  Biometrika 00(0), 1–8.
  * Kass et al. (1990) Kass, R. E., L. Tierney, and J. B. Kadane (1990).  The validity of posterior expansions based on Laplace’s method.  Bayesian and Likelihood Methods in Statistics and Econometrics, 473–488.
  * Kleijn and van der Vaart (2012) Kleijn, B. and A. van der Vaart (2012).  The Bernstein von-Mises theorem under misspecification.  Electronic Journal of Statistics 6, 354–381.
  * Knorr-Held and Rue (2002) Knorr-Held, L. and H. Rue (2002).  On block updating in Markov random field models for disease mapping.  Scandinavian Journal of Statistics $\bf{29}$, 597–614.
  * Kristensen et al. (2016) Kristensen, K., A. Nielson, C. W. Berg, H. Skaug, and B. M. Bell (2016).  TMB: automatic differentiation and Laplace approximation.  Journal of statistical software 70(5).
  * Law et al. (2018) Law, H. C. L., D. Sejdinovic, E. Cameron, T. C. Lucas, S. Flaxman, K. Battle, and K. Fukumizu (2018).  Variational learning on aggregate outputs with Gaussian processes.  arXiv:1805.08463 [stat.ML].
  * Lee et al. (2017) Lee, J. S., P. Nguyen, P. E. Brown, J. Stafford, and N. Saint-Jacques (2017).  A local-EM algorithm for spatio-temporal disease mapping with aggregated data.  Spatial Statistics 21, 75–95.
  * Li et al. (2012) Li, Y., P. Brown, D. C. Gesink, and H. Rue (2012).  Log Gaussian Cox processes and spatially aggregated disease incidence data.  Statistical Methods in Medical Research 21(5), 479 – 507\.
  * Lindgren and Rue (2015) Lindgren, F. and H. Rue (2015).  Bayesian spatial modelling with R-INLA.  Journal of Statistical Software 63(19), 1–25.
  * Lindgren et al. (2011) Lindgren, F., H. Rue, and J. Lindstroom (2011).  An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society, Series B (Statistical Methodology) 73(4), 423–498.
  * Margossian et al. (2020) Margossian, C. C., A. Vehtari, D. Simpson, and R. Agrawal (2020).  Hamiltonian Monte Carlo using an adjoint-differentiated Laplace approximation.  arXiv:2004.12550v3 [stat.CO].
  * Martino et al. (2011) Martino, S., R. Akerkar, and H. Rue (2011).  Approximate bayesian inference for survival models.  Scandinavian Journal of Statistics 38(3), 514 – 528.
  * Martins et al. (2013) Martins, T. G., D. Simpson, F. Lindgren, and H. Rue (2013).  Bayesian computing with INLA: new features.  Computational Statistics and Data Analysis 67, 68–83.
  * Monnahan and Kristensen (2018) Monnahan, C. and K. Kristensen (2018).  No-U-turn sampling for fast Bayesian inference in ADMB and TMB: introducing the adnuts and tmbstan R packages.  PloS one 13(5).
  * Murray et al. (2010) Murray, I., R. P. Adams, and D. J. C. MacKay (2010).  Elliptical Slice Sampling.  Proceedings of the 13th International Conference on Artificial Intelligence and Statistics, 541–548.
  * Nandi et al. (2020) Nandi, A. K., T. C. D. Lucas, R. Arambepola, P. Gething, and D. Weiss (2020).  disaggregation: an R package for Bayesian spatial disaggregation modelling.  arxiv.
  * Ogden (2017) Ogden, H. (2017).  On asymptotic validity of naive inference with an approximate likelihood.  Biometrika 104.
  * Ogden (2021) Ogden, H. (2021).  On the error in laplace approximations of high-dimensional integrals.  Stat 10.
  * Petras (2003) Petras, K. (2003).  Smolyak cubature of given polynomial degree with few nodes for increasing dimension.  Numerische Mathematik 93, 729–753.
  * Pinheiro and Bates (1995) Pinheiro, J. C. and D. M. Bates (1995).  Approximations to the log-likelihood function in the nonlinear mixed effects models.  Journal of Computational and Graphical Statistics 4(1), 12–35.
  * Rue (2001) Rue, H. (2001).  Fast sampling of Gaussian Markov random fields.  Journal of the Royal Statistical Society, Series B: Statistical Methodology 63(2), 325–338.
  * Rue and Martino (2007) Rue, H. and S. Martino (2007).  Approximate Bayesian inference for hierarchical Gaussian Markov random field models.  Journal of Statistical Planning and Inference 137, 3177 – 3192.
  * Rue et al. (2009) Rue, H., S. Martino, and N. Chopin (2009).  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  Journal of the Royal Statistical Society. Series B (Statistical Methodology) 71(2), 319 – 392.
  * Rue et al. (2017) Rue, H., A. Riebler, S. H. Sørbye, J. B. Illian, D. P. Simpson, and F. K. Lindgren (2017).  Bayesian computing with INLA: a review.  Annu. Rev. Stat. Appl. 4, 395–421.
  * Schlather et al. (2015) Schlather, M., A. Malinowski, P. J. Menck, M. Oesting, and K. Strokorb (2015).  Analysis, simulation and prediction of multivariate random fields with package RandomFields.  Journal of Statistical Software 63(8).
  * Shun and McCullagh (1995) Shun, Z. and P. McCullagh (1995).  Laplace approximation of high dimensional integrals.  Journal of the Royal Statistical Society, Series B: Statistical Methodology 57(4), 749–760.
  * Simpson et al. (2016) Simpson, D., J. B. Illian, F. Lindgren, S. H. Sorbye, and H. Rue (2016).  Going off grid: computatioanlly efficient inference for log-Gaussian Cox processes.  Biometrika 103(1), 49–70.
  * Simpson et al. (2017) Simpson, D., H. Rue, T. G. Martins, A. Riebler, and S. H. Sørbye (2017).  Penalising model component complexity: A principled, practical approach to constructing priors.  Statistical Science 32(1).
  * Stringer et al. (2021) Stringer, A., P. Brown, and J. Stafford (2021).  Approximate Bayesian inference for case crossover models.  Biometrics 77(3), 785–795.
  * Tierney and Kadane (1986) Tierney, L. and J. B. Kadane (1986).  Accurate approximations to posterior moments and marginal densities.  Journal of the American Statistical Association 81(393).
  * van der Vaart (1998) van der Vaart, A. (1998).  Asymptotic Statistics.  Cambridge University Press.
  * van Niekerk et al. (2021) van Niekerk, J., H. Bakka, and H. Rue (2021).  Competing risks joint models using r-inla.  Statistical Modelling 21(1-2), 56–71.
  * Vehtari et al. (2017) Vehtari, A., A. Gelman, and J. Gabry (2017).  Practical bayesian model evaluation using leave-one-out cross-validation and waic.  Statistics and Computing 27, 1413–1432.
  * Vehtari et al. (2021) Vehtari, A., D. Simpson, A. Gelman, Y. Yao, and J. Gabry (2021).  Pareto smoothed importance sampling.  arXiv:1507.02646v7 [stat.CO].
  * Wachter and Biegler (2006) Wachter, A. and L. T. Biegler (2006).  On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming.  Mathematical Programming 106(1), 25 – 57.
  * Walker (1969) Walker, A. M. (1969).  On the asymptotic behaviour of posterior distributions.  Journal of the Royal Statistical Society, Series B (Statistical Methodology) 31(1), 80–88.
  * Weiser (2016) Weiser, C. (2016).  mvQuad: Methods for Multivariate Quadrature. (R package version 1.0-6).
  * Wood (2011) Wood, S. (2011).  Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models.  Journal of the Royal Statistical Society, Series B (Statistical Methodology) 73(1), 3 – 36.
  * Wood (2016) Wood, S. (2016).  Smoothing parameter and model selection for general smooth models.  Journal of the American Statistical Association 111, 1548 – 1575.
  * Wood (2020) Wood, S. (2020).  Simplified integrated nested Laplace approximation.  Biometrika 107(1), 223–230.
  * Yao et al. (2018) Yao, Y., A. Vehtari, D. Simpson, and A. Gelman (2018).  Yes, but did it work?: evaluating variational inference.  Proceedings of the 35th International Conference on Machine Learning.
  * Zhang (2004) Zhang, H. (2004).  Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics.  Journal of the American Statistical Association 99(465), 250–261.