text:
\section{Introduction}\label{introduction}}

Bayesian inferences are based upon the posterior distribution. Computing
this distribution and summaries based upon it requires evaluating a
difficult integral, called the normalizing constant or marginal
likelihood, a computation which is infeasible to do exactly in all but
the simplest of problems. Modern approaches for inference in a wide
range of models are based on approximations to the posterior. These
approximations, in turn, rely on approximations to the normalizing
constant.

The normalizing constant is an integral, and in principle any quadrature
technique may be used to approximate it. In practice, because the
location and shape of the posterior are random variables which vary with
the data, techniques which adapt to these properties are desirable. One
such adaptive technique is Adaptive Gauss-Hermite Quadrature
(\texttt{AGHQ}{}). Since \texttt{AGHQ}{} was introduced in the
statistical literature
\citep{nayloradaptive,laplace,adaptive_GH_1994,adaptive_GH_2020} it has
been used in a variety of applications including latent variable models
\citep{aghqmle} and as the default option for approximating the likelihood
in Generalized Linear Mixed Models in the \pkg{glmer} function in the
popular \pkg{lme4} package \citep{lme4}. Both the theoretical and
practical advantages of \texttt{AGHQ}{} as a tool for Bayesian inference
in complicated models have recently been studied by \citet{aghqus}.

Like all quadrature techniques, \texttt{AGHQ}{} in principle only
applies to models in which the dimension of the parameter space is very
small. However, \texttt{AGHQ}{} has the special property that when run
with only a single quadrature point, it becomes a Laplace approximation
and hence may be used when the dimension of the parameter space is
large. Further, while it has been demonstrated that there are difficult
problems for which \texttt{AGHQ}{} can be used on its own
\citep{aghqus}, it also forms an important part of modern approximate
Bayesian inference techniques for more general classes of models with
high dimensional parameter spaces \citep{noeps}. There is a wide range
of Bayesian inference problems for which \texttt{AGHQ}{} is a useful
tool.

Despite its popularity, no general purpose \pkg{R} package exists that
is tailored to the use of \texttt{AGHQ}{} for Bayesian inference. While
quadrature is conceptually straightforward to implement and general
purpose packages exist in \texttt{R} for this purpose---such as the
\pkg{mvQuad} package \citep{mvquad}---posterior summaries such as
moments, quantiles, and marginal densities are nontrivial to obtain once
the approximate normalizing constant has been computed. Marginal
densities are especially nontrival, see \citet{nayloradaptive} or
\citet{aghqus}, Appendix E. A package which provides an interface for
normalizing the posterior using \texttt{AGHQ}{} and then using this
normalized posterior for inference is desirable for scientists for whom
\texttt{AGHQ}{} may be faster than existing tools for their specific
problems, and for statisticians who wish to use \texttt{AGHQ}{} to use
approximate Bayesian inference to fit new models. I provide examples of both types of
use in this manuscript.

There are several existing packages which implement \texttt{AGHQ}{} to
some extent, however they are all either too specific or too general,
applying only to specific models or lack functions for
using the results to compute posterior summaries or plots. For example,
the \pkg{fastGHQuad} package \citep{fastghquad} implements
\texttt{AGHQ}{} for general scalar functions only, requires the user to
supply their own optimization results, and does not include any
functions for post-processing the results. The \pkg{GLMMadaptive}
package \citep{glmmadaptive} handles optimization and has a formula
interface and an API for post-processing model results, but only fits
Generalized Linear Mixed Models with independent random effects. The
\pkg{LaplacesDemon} package provides an option to use \texttt{AGHQ}{} to
improve moment estimation in models fit either with \texttt{MCMC}{} or a
Laplace Approximation, but does not offer further use of the technique.
A package for \texttt{AGHQ}{} which sits in the middle, providing a
flexible interface for quadrature and Bayesian inference in a wide range
of models, is desirable.

In this paper I describe how to use the \pkg{aghq} package to
approximate the normalizing constant using \texttt{AGHQ}{}, and make
Bayesian inferences based on the result. I demonstrate that the package
is simple enough to cover many problems, yet flexible enough to be used
as a building block within more complicated frameworks for Bayesian
inference. For simpler problems, the user simply provides a function
which computes the un-normalized posterior, and the package handles the
optimization; for more complicated problems, the user may supply their
own optimization results if desired. Standard \texttt{print},
\texttt{summary}, and \texttt{plot} methods are provided for the output,
as are several helper functions which act on \texttt{aghq} objects to
compute moments, quantiles, and marginal densities and distribution
functions. The use of the basic features of the package for inference is demonstrated through the detailed reconstruction of the two low-dimensional examples from \citet{aghqus}.

I further demonstrate how the more advanced features of the package are used to
implement approximate Bayesian inference for complicated high-dimensional models. The popular \pkg{TMB} package \citep{tmb}
offers an interface for automatic differentiation and unnormalized marginal Laplace
approximations in \pkg{R}, and the \pkg{tmbstan} package \citep{tmbstan} offers an
interface for Bayesian inference via MCMC for objects created using \pkg{TMB}. The
\pkg{aghq} package is designed to work with the output of \pkg{TMB} in a similar fashion to \pkg{tmbstan}. However, instead of MCMC samplng, \pkg{aghq} offers
the ability to take a log-posterior defined using \pkg{TMB}, normalize
it using AGHQ, and make approximate Bayesian inferences. I demonstrate this
functionality through several examples including non-Gaussian geostatistical modelling
with and without zero-inflation, and a zero-inflated, over-dispersed Poisson regression.
The latter example makes use of the \pkg{glmmTMB} package to define the \pkg{TMB}
function template, providing a convenient user interface that can be adapted in other problems as well.

The remainder of this paper is organized as follows. In \S\ref{sec:basic-use} I describe the basic use of the package in
detail using the simple example from the simulations in \citet{aghqus}.
In \S\ref{sec:background} I give the necessary background on
\texttt{AGHQ}{} including explicit formulas and a brief statement of
relevant mathematical properties. In \S\ref{sec:examples} I implement
the examples from \citet{aghqus} using the \pkg{aghq} package, exposing
the computational details of their implementation and hence providing
users with a detailed template describing the use of the package for
difficult Bayesian inference problems. In \S\ref{sec:elgms} I implement
the Generalized Linear Geostatistical Model (GLGM) from \citet{geostatsp} from
scratch using \pkg{aghq}, and then extend it to fit the much more challenging
zero-inflated binomial
geostatistical model described by \citet{geostatlowresource}. All of the examples
of \S\ref{sec:examples} and \S\ref{sec:elgms} make use of \pkg{TMB} for defining the models. In
\S\ref{sec:glmmTMB} I take the zero-inflated, overdispersed Poisson
model implemented in the \pkg{glmmTMB} package \citep{glmmTMB} and use
the \pkg{aghq} package to make fully Bayesian inferences for it,
leveraging their existing code. This provides a template for doing
similar extensions using other existing models. I conclude in
\S\ref{sec:discussion} with a discussion of impact and future work.

\hypertarget{sec:basic-use}{%
\section{Basic Use}\label{sec:basic-use}}

In this section I describe the basic use of the \pkg{aghq} package by
way of analyzing a simple example in which the answer is known. I focus
on the use and features of the package, with detailed discussion of
Bayesian inference and quadrature in \texttt{R} deferred to
\S\ref{sec:background}. All notation and terms are defined in detail in
that section.

The following conjugate model is used by \citet{aghqus}:
\begin{equation}\begin{aligned}
Y_i \ \vert \ \lambda &\overset{ind}{\sim} \text{Poisson}(\lambda), i\in[n], \\
\lambda &\sim \text{Exponential}(1), \lambda > 0, \\
\implies \lambda \ \vert \ \boldsymbol{Y}&\sim \text{Gamma}\left(1 + \sum_{i=1}^{n}y_{i},n+1\right), \\
\pi(\lambda,\boldsymbol{Y}) &= \frac{\lambda^{\sum_{i=1}^{n}y_{i}}e^{-n\lambda}}{\prod_{i=1}^{n}y_{i}!}\times e^{-\lambda}, \\
\pi(\boldsymbol{Y}) &= \frac{\Gamma\left(1 + \sum_{i=1}^{n}y_{i}\right)}{(n+1)^{\left(1 + \sum_{i=1}^{n}y_{i}\right)}\prod_{i=1}^{n}y_{i}!}
\end{aligned}\end{equation} Although \texttt{AGHQ}{} may be applied to
models with constrained parameter spaces, in practice it
is desirable to remove simple box constraints like \(\lambda > 0\) using
transformations. Further, \citet{nayloradaptive} argue that any
transformation which brings the posterior closer to spherical symmetry
will improve either the empirical accuracy of the quadrature, the
stability and speed of the optimization, or both. Here we set:
\begin{equation}\begin{aligned}
\eta &= \log\lambda, \\
\implies \log\pi(\eta,\boldsymbol{Y}) &= \eta\sum_{i=1}^{n}y_{i} - (n+1)e^{\eta} - \sum_{i=1}^{n}\log(y_{i}!) + \eta, \\
\end{aligned}\end{equation} which, because the jacobian
\(\left\lvert \mathrm{d}\lambda/\mathrm{d}\eta \right\rvert\) is
included, does not change the value of \(\pi(\boldsymbol{Y})\).

To approximate \(\pi(\boldsymbol{Y})\) and then make Bayesian inferences
based on \(\pi(\lambda|\boldsymbol{Y})\), first install and load the
\pkg{aghq} package:

\begin{CodeChunk}
\begin{CodeInput}
R> # Install stable version from CRAN:
R> # install.packages("aghq")
R> # or development version from github:
R> # devtools::install_github("awstringer1/aghq")
R> library(aghq)
\end{CodeInput}
\end{CodeChunk}

We will simulate some data from the model:

\begin{CodeChunk}
\begin{CodeInput}
R> set.seed(84343124)
R> y <- rpois(10,5) # True lambda = 5, n = 10
\end{CodeInput}
\end{CodeChunk}

The main function is the \texttt{aghq::aghq} function. The user supplies
a list containing the log-posterior and its first two derivatives:

\begin{CodeChunk}
\begin{CodeInput}
R> # Define the log posterior, log(pi(eta,y)) here
R> logpietay <- function(eta,y) {
+   sum(y) * eta - (length(y) + 1) * exp(eta) - sum(lgamma(y+1)) + eta
+ }
R> objfunc <- function(x) logpietay(x,y)
R> objfuncgrad <- function(x) numDeriv::grad(objfunc,x)
R> objfunchess <- function(x) numDeriv::hessian(objfunc,x)
R> # Now create the list to pass to aghq()
R> funlist <- list(
+   fn = objfunc,
+   gr = objfuncgrad,
+   he = objfunchess
+ )
\end{CodeInput}
\end{CodeChunk}

Now perform the quadrature. The only other required input is the number
of quadrature points, \(k\), and a starting value for the optimization:

\begin{CodeChunk}
\begin{CodeInput}
R> # AGHQ with k = 3
R> # Use eta = 0 as a starting value
R> thequadrature <- aghq::aghq(ff = funlist,k = 3,startingvalue = 0)
\end{CodeInput}
\end{CodeChunk}

The object \texttt{thequadrature} has class \texttt{aghq}, with
\texttt{summary} and \texttt{plot} methods:

\begin{CodeChunk}
\begin{CodeInput}
R> summary(thequadrature)
\end{CodeInput}
\begin{CodeOutput}
AGHQ on a 1 dimensional posterior with  3 quadrature points

The posterior mode is: 1.49

The log of the normalizing constant/marginal likelihood is: -23.3

The posterior Hessian at the mode is:
     [,1]
[1,]   49

The covariance matrix used for the quadrature is...
       [,1]
[1,] 0.0204

...and its Cholesky is:
      [,1]
[1,] 0.143

Here are some moments and quantiles for theta:

       mean median mode    sd 2.
theta1 1.48   1.48 1.49 0.142  1.2  1.76
\end{CodeOutput}
\begin{CodeInput}
R> plot(thequadrature)
\end{CodeInput}

\begin{center}\includegraphics{01-aghq-softwarepaper_files/figure-latex/doaghq2-1} \end{center}

\begin{center}\includegraphics{01-aghq-softwarepaper_files/figure-latex/doaghq2-2} \end{center}

\end{CodeChunk}

The actual object is a list with elements (see \S\ref{sec:background}
for notation):

\begin{itemize}
\item
  \texttt{normalized\_posterior}: the output of the
  \texttt{aghq::normalize\_posterior} function, which is itself a list
  with elements:

  \begin{itemize}
  \item
    \texttt{nodesandweights}: a dataframe containing the nodes and
    weights from the quadrature and the normalized and un-normalized log
    posterior at the nodes,
  \item
    \texttt{thegrid}: a \texttt{NIGrid} object from the \texttt{mvQuad}
    package, see \texttt{?mvQuad::createNIGrid},
  \item
    \texttt{lognormconst}: the numeric value of
    \(\log\widetilde{\pi}_{\texttt{AGHQ}}(\boldsymbol{Y})\).
  \end{itemize}
\item
  \texttt{marginals}: a list of the same length as
  \texttt{startingvalue} of which element \texttt{j} is the result of
  calling \texttt{aghq::marginal\_posterior} with that \texttt{j}. This
  is a \texttt{tbl\_df/tbl/data.frame} containing the normalized
  marginal posterior for the \(j^{th}\) element of
  \(\boldsymbol{\theta}\). In this one-dimensional example, this is
  exactly the same information that is present in
  \texttt{thequadrature\$normalized\_logpost\$nodesandweights}, however
  this object will become important in the multidimensional examples in
  \S\ref{sec:examples}, \S\ref{sec:elgms}, and \S\ref{sec:glmmTMB},
\item
  \texttt{optresults}: information and results from the optimization of
  \(\log\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y})\). This is itself a
  list with elements:

  \begin{itemize}
  \item
    \texttt{ff}: the function list passed to \texttt{aghq},
  \item
    \texttt{mode}: \(\widetilde{\boldsymbol{\theta}}_n\),
  \item
    \texttt{hessian}: \(\boldsymbol{H}\),
  \item
    \texttt{convergence}: the convergence code, specific to the
    optimizer used.
  \end{itemize}
\end{itemize}

We can observe the structure of the object and compare the results to
the truth for this simple example:

\begin{CodeChunk}
\begin{CodeInput}
R> # The normalized posterior:
R> thequadrature$normalized_posterior$nodesandweights
\end{CodeInput}
\begin{CodeOutput}
  theta1 weights logpost logpost_normalized
1   1.25   0.267   -23.7             -0.357
2   1.49   0.239   -22.3              1.027
3   1.74   0.267   -23.9             -0.605
\end{CodeOutput}
\begin{CodeInput}
R> # The log normalization constant:
R> options(digits = 6)
R> thequadrature$normalized_posterior$lognormconst
\end{CodeInput}
\begin{CodeOutput}
[1] -23.3212
\end{CodeOutput}
\begin{CodeInput}
R> # Compare to the truth:
R> lgamma(1 + sum(y)) - (1 + sum(y)) * log(length(y) + 1) - sum(lgamma(y+1))
\end{CodeInput}
\begin{CodeOutput}
[1] -23.3195
\end{CodeOutput}
\begin{CodeInput}
R> # Quite accurate with only n = 10 and k = 3; this example is very simple.
R> # The mode found by the optimization:
R> thequadrature$optresults$mode
\end{CodeInput}
\begin{CodeOutput}
[1] 1.49393
\end{CodeOutput}
\begin{CodeInput}
R> # The true mode:
R> log((sum(y) + 1)/(length(y) + 1))
\end{CodeInput}
\begin{CodeOutput}
[1] 1.49393
\end{CodeOutput}
\begin{CodeInput}
R> options(digits = 3)
\end{CodeInput}
\end{CodeChunk}

Of course, in practical problems, the true mode and normalizing constant
won't be known.

The \pkg{aghq} package provides further routines for computing moments,
quantiles, and distributions/densities. These are especially useful when
working with transformations like we are in the example, since interest
here is in the original parameter \(\lambda = e^{\eta}\). The main
functions are:

\begin{itemize}
\item
  \texttt{compute\_pdf\_and\_cdf}: compute the density and cumulative
  distribution function for a marginal posterior distribution of a
  variable and (optionally) a smooth monotone transformation of that
  variable,
\item
  \texttt{compute\_moment}: compute the posterior moment of any
  function,
\item
  \texttt{compute\_quantiles}: compute posterior quantiles for a
  marginal posterior distribution.
\end{itemize}

To compute the approximate marginal
posterior for \(\lambda\),
\(\widetilde{\pi}_{\texttt{AGHQ}}(\lambda|\boldsymbol{Y})\), first
compute the marginal posterior for \(\eta\) and then transform:
\begin{equation}
\widetilde{\pi}_{\texttt{AGHQ}}(\lambda|\boldsymbol{Y}) = \widetilde{\pi}_{\texttt{AGHQ}}(\eta|\boldsymbol{Y})\left\lvert \frac{\mathrm{d}\eta}{\mathrm{d}\lambda} \right\rvert
\end{equation}

The \texttt{aghq::compute\_pdf\_and\_cdf} function has an option,
\texttt{transformation}, which allows the user to specify a parameter
transformation that they would like the marginal density of. The user
specifies functions \texttt{fromtheta} and \texttt{totheta} which
convert from and to the transformed scale (on which quadrature was
done), and the function returns the marginal density on both scales,
making use of a numerically differentiated jacobian. This is all done as
follows:

\begin{CodeChunk}
\begin{CodeInput}
R> # Compute the pdf for eta
R> pdfwithlambda <- compute_pdf_and_cdf(
+   thequadrature$marginals[[1]],
+   transformation = list(
+     totheta = log,
+     fromtheta = exp
+   )
+ )
R> glimpse(pdfwithlambda)
\end{CodeInput}
\begin{CodeOutput}
Rows: 1,000
Columns: 5
$ theta          <dbl> 0.999, 1.000, 1.001, 1.002, 1.003, 1.004, 1.005, 1.0...
$ pdf            <dbl> 0.00860, 0.00881, 0.00902, 0.00923, 0.00945, 0.00968...
$ cdf            <dbl> 0.00000000, 0.00000873, 0.00001766, 0.00002681, 0.00...
$ transparam     <dbl> 2.72, 2.72, 2.72, 2.72, 2.73, 2.73, 2.73, 2.73, 2.74...
$ pdf_transparam <dbl> 0.00317, 0.00324, 0.00331, 0.00339, 0.00347, 0.00355...
\end{CodeOutput}
\begin{CodeInput}
R> # Plot along with the true posterior
R> lambdapostplot <- pdfwithlambda
+   ggplot(aes(x = transparam,y = pdf_transparam)) +
+   theme_classic() +
+   geom_line() +
+   stat_function(fun = dgamma,
+                 args = list(shape = 1+sum(y),rate = 1 + length(y)),
+                 linetype = 'dashed') +
+   labs(x = expression(lambda),y = "Density")
R> lambdapostplot
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{01-aghq-softwarepaper_files/figure-latex/postforlambda1-1}

}

\caption[Approximate (---) and true (- - -) posterior for $\lambda$]{Approximate (---) and true (- - -) posterior for $\lambda$}\label{fig:postforlambda1}
\end{figure}
\end{CodeChunk}

We may compute the posterior mean of \(\lambda = e^{\eta}\),
\(\mathop{\mathrm{\mathbb{E}}}(\lambda|\boldsymbol{Y})\) by using the
\texttt{compute\_moment} function. This can also be used to confirm that
the posterior is, in fact, properly normalized:

\begin{CodeChunk}
\begin{CodeInput}
R> options(digits = 6)
R> # Check if the posterior integrates to 1, by computing the "moment" of "1":
R> compute_moment(thequadrature$normalized_posterior,
+                ff = function(x) 1)
\end{CodeInput}
\begin{CodeOutput}
[1] 1
\end{CodeOutput}
\begin{CodeInput}
R> # Posterior mean for eta:
R> compute_moment(thequadrature$normalized_posterior,
+                ff = function(x) x)
\end{CodeInput}
\begin{CodeOutput}
[1] 1.48374
\end{CodeOutput}
\begin{CodeInput}
R> # Posterior mean for lambda = exp(eta)
R> compute_moment(thequadrature$normalized_posterior,
+                ff = function(x) exp(x))
\end{CodeInput}
\begin{CodeOutput}
[1] 4.45441
\end{CodeOutput}
\begin{CodeInput}
R> # Compare to the truth:
R> (sum(y) + 1)/(length(y) + 1)
\end{CodeInput}
\begin{CodeOutput}
[1] 4.45455
\end{CodeOutput}
\begin{CodeInput}
R> options(digits = 3)
\end{CodeInput}
\end{CodeChunk}

Quantiles are computed using the \texttt{compute\_quantiles} function.
For example, to get the first and third percentiles along with median
and first and third quartiles:

\begin{CodeChunk}
\begin{CodeInput}
R> # Get the quantiles for eta:
R> etaquants <- compute_quantiles(
+   thequadrature$marginals[[1]],
+   q = c(.01,.25,.50,.75,.99)
+ )
R> # The quantiles for lambda = exp(eta) are obtained directly:
R> exp(etaquants)
\end{CodeInput}
\begin{CodeOutput}

3.17 4.00 4.40 4.85 6.15
\end{CodeOutput}
\begin{CodeInput}
R> # and compared with the truth:
R> qgamma(c(.01,.25,.50,.75,.99),shape = 1+sum(y),rate = 1+length(y))
\end{CodeInput}
\begin{CodeOutput}
[1] 3.11 4.01 4.42 4.87 6.07
\end{CodeOutput}
\end{CodeChunk}

The estimation of quantiles, especially extreme quantiles, is much more
difficult than that of moments, and this is reflected in the small
differences in accuracy observed between the two in this example.

\hypertarget{sec:background}{%
\section{Background and Details}\label{sec:background}}

In this section I expand upon the underlying methodology of
\texttt{AGHQ}{}, including a brief review of Bayesian inference in
general, a description of quadrature methods in \texttt{R}, and a
discussion of the calculation of posterior summary statistics. The notation follows \citet{aghqus} who describe the methodological background in much greater detail. I focus here on the computational details.

\hypertarget{bayesian-inference}{%
\subsection{Bayesian Inference}\label{bayesian-inference}}

Let \(\boldsymbol{\theta}\in\Theta \subseteq \mathbb{R}^{d}\) denote a
parameter, let \(\pi(\boldsymbol{\theta})\) denote its prior
density, and let \(\boldsymbol{Y}= \left\{ Y_{i}:i\in[n] \right\}\)
denote a sample with likelihood
\(\pi(\boldsymbol{Y}|\boldsymbol{\theta})\).
Bayesian inference is based upon the posterior distribution:
\begin{equation}\label{eqn:posterior}
\pi(\boldsymbol{\theta}|\boldsymbol{Y}) = \frac{\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y})}{\int_{\Theta}\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y})\mathrm{d}\boldsymbol{\theta}},
\end{equation} where
\(\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y}) = \pi(\boldsymbol{Y}|\boldsymbol{\theta})\pi(\boldsymbol{\theta})\).
The denominator of (\ref{eqn:posterior}) is called the normalizing
constant, denoted by \(\pi(\boldsymbol{Y})\): \begin{equation}
\pi(\boldsymbol{Y}) = \int_{\Theta}\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y})\mathrm{d}\boldsymbol{\theta},
\end{equation} and this is infeasible to compute directly unless the model is very simple. The \pkg{aghq} package provides a framework for
Bayesian inference in which \(\pi(\boldsymbol{Y})\) is approximated
using \texttt{AGHQ}{}.

\hypertarget{subsec:quadrature}{%
\subsection{Quadrature}\label{subsec:quadrature}}

For any function \(f:\mathbb{R}^{d}\to\mathbb{R}\), a quadrature
approximation to its integral with \(k\) nodes is given by
\begin{equation}
\int_{A}f(\boldsymbol{x})\mathrm{d}\boldsymbol{x} \approx \sum_{j=1}^{k}f(\boldsymbol{x}_{j})\omega_{j},
\end{equation} where \(A\subseteq \mathbb{R}^{d}\),
\(\left\{ \boldsymbol{x}_{j} \right\}_{j=1}^{k}\) are \emph{nodes}, and
\(\left\{ \omega_{j} \right\}_{j=1}^{k}\) are \emph{weights}. In
\pkg{R}, quadrature is performed using the flexible interface provided
by the \pkg{mvQuad} package:

\begin{CodeChunk}
\begin{CodeInput}
R> # Integrate exp(-x^2) on (-Inf,Inf)
R> f <- function(x) exp(-(x^2)/2)
R> # Define a grid for Gauss-Hermite quadrature
R> # with k = 3 points
R> gg <- mvQuad::createNIGrid(1,"GHe",3)
R> # Perform the quadrature
R> mvQuad::quadrature(f,gg)
\end{CodeInput}
\begin{CodeOutput}
[1] 2.51
\end{CodeOutput}
\begin{CodeInput}
R> # True value is sqrt(2pi)
R> sqrt(2*pi)
\end{CodeInput}
\begin{CodeOutput}
[1] 2.51
\end{CodeOutput}
\end{CodeChunk}

The \pkg{mvQuad} package computes the nodes and weights for a given
quadrature rule and order \(k\). For the above rule, these are

\begin{CodeChunk}
\begin{CodeInput}
R> mvQuad::getNodes(gg)
\end{CodeInput}
\begin{CodeOutput}
          [,1]
[1,] -1.73e+00
[2,] -1.39e-16
[3,]  1.73e+00
\end{CodeOutput}
\begin{CodeInput}
R> mvQuad::getWeights(gg)
\end{CodeInput}
\begin{CodeOutput}
     [,1]
[1,] 1.87
[2,] 1.67
[3,] 1.87
\end{CodeOutput}
\end{CodeChunk}

In multiple dimensions, the nodes are repeated in each dimension, and
the weights are multiplied:

\begin{CodeChunk}
\begin{CodeInput}
R> # Integrate  exp(-(x+y)^2 / 2) on R^2
R> f <- function(x) apply(x,1,function(y) prod(exp(-(y^2)/2)))
R> # Define a grid for Gauss-Hermite quadrature
R> # with k = 3 points
R> gg2 <- mvQuad::createNIGrid(2,"GHe",3)
R> # Nodes, and manual construction:
R> mvQuad::getNodes(gg2)
\end{CodeInput}
\begin{CodeOutput}
           [,1]      [,2]
 [1,] -1.73e+00 -1.73e+00
 [2,] -1.39e-16 -1.73e+00
 [3,]  1.73e+00 -1.73e+00
 [4,] -1.73e+00 -1.39e-16
 [5,] -1.39e-16 -1.39e-16
 [6,]  1.73e+00 -1.39e-16
 [7,] -1.73e+00  1.73e+00
 [8,] -1.39e-16  1.73e+00
 [9,]  1.73e+00  1.73e+00
\end{CodeOutput}
\begin{CodeInput}
R> # Equivalent:
R> expand.grid(mvQuad::getNodes(gg),mvQuad::getNodes(gg))
\end{CodeInput}
\begin{CodeOutput}
       Var1      Var2
1 -1.73e+00 -1.73e+00
2 -1.39e-16 -1.73e+00
3  1.73e+00 -1.73e+00
4 -1.73e+00 -1.39e-16
5 -1.39e-16 -1.39e-16
6  1.73e+00 -1.39e-16
7 -1.73e+00  1.73e+00
8 -1.39e-16  1.73e+00
9  1.73e+00  1.73e+00
\end{CodeOutput}
\begin{CodeInput}
R> # Weights:
R> mvQuad::getWeights(gg2)
\end{CodeInput}
\begin{CodeOutput}
      [,1]
 [1,] 3.51
 [2,] 3.13
 [3,] 3.51
 [4,] 3.13
 [5,] 2.79
 [6,] 3.13
 [7,] 3.51
 [8,] 3.13
 [9,] 3.51
\end{CodeOutput}
\begin{CodeInput}
R> # Equivalent:
R> cbind(apply(expand.grid(mvQuad::getWeights(gg),mvQuad::getWeights(gg)),1,prod))
\end{CodeInput}
\begin{CodeOutput}
      [,1]
 [1,] 3.51
 [2,] 3.13
 [3,] 3.51
 [4,] 3.13
 [5,] 2.79
 [6,] 3.13
 [7,] 3.51
 [8,] 3.13
 [9,] 3.51
\end{CodeOutput}
\begin{CodeInput}
R> # Perform the quadrature
R> mvQuad::quadrature(f,gg2)
\end{CodeInput}
\begin{CodeOutput}
[1] 6.28
\end{CodeOutput}
\begin{CodeInput}
R> # True value is 2pi
R> 2 * pi
\end{CodeInput}
\begin{CodeOutput}
[1] 6.28
\end{CodeOutput}
\end{CodeChunk}

Applying a quadrature rule with \(k\) points per dimension to a
\(d\)-dimensional integrand therefore requires \(k^{d}\) function
evaluations. Methods based on sparse grids \citep{sparsegrids} or
employing more challenging computations \citep{designedquadrature} may
alleviate this challenge somewhat, however these are not yet shown to be
appropriate for the specific task of Bayesian inference and are not
implemented in the \pkg{aghq} package.

\hypertarget{gauss-hermite-quadrature}{%
\subsection{Gauss-Hermite Quadrature}\label{gauss-hermite-quadrature}}

Gauss-Hermite Quadrature (\texttt{GHQ}{}, \citet{numint}) is a useful
quadrature rule based on the theory of polynomial interpolation. Define
the \emph{Hermite polynomials}
\[H_{k}(x) = (-1)^k e^{x^{2}/2}\frac{d^{k}}{dx^{k}}e^{-x^{2}/2}\] where
\(k\in\mathbb{N}\). The \texttt{GHQ}{} nodes and weights with \(k\)
quadrature points are: \begin{equation}\begin{aligned}
\left\{ x_{j} \right\}_{j=1}^{k} &= \left\{ x\in\mathbb{R}: H_{k}(x) = 0 \right\} \\
\omega_{j} &= \frac{k!}{H_{k+1}(x_{j})^{2} \times \phi(x_{j})}, j\in[k]
\end{aligned}\end{equation} where \(\phi(\cdot)\) is a standard normal
density. These may be computed in \pkg{R}:

\begin{CodeChunk}
\begin{CodeInput}
R> # Gauss Hermite, 1 dimension, 3 points:
R> gg <- mvQuad::createNIGrid(1,"GHe",3)
R> nn <- as.numeric(mvQuad::getNodes(gg))
R> ww <- as.numeric(mvQuad::getWeights(gg))
R> round(nn,14) # Print zero properly
\end{CodeInput}
\begin{CodeOutput}
[1] -1.73  0.00  1.73
\end{CodeOutput}
\begin{CodeInput}
R> ww
\end{CodeInput}
\begin{CodeOutput}
[1] 1.87 1.67 1.87
\end{CodeOutput}
\begin{CodeInput}
R> # Check zeroes of Hermite polynomial:
R> Hk <- as.function(mpoly::hermite(3,kind = "he")) # Hermite polynomial defined above
R> round(Hk(nn),14) # All zero
\end{CodeInput}
\begin{CodeOutput}
[1] 0 0 0
\end{CodeOutput}
\end{CodeChunk}

Some interesting properties of \texttt{GHQ}{} are as follows
\citep{numint}:

\begin{itemize}
\item $x = 0$ is a node if and only if $k$ is odd,
\item $-x$ is a node if and only $x$ is,
\item $\omega_{j} > 0$ for each $j\in[k]$,
\item for a polynomial $P_{\alpha}(x)$ of degree $\alpha \leq 2k-1$, \texttt{GHQ}{}
satisfies $$\int_{\mathbb{R}}P_{\alpha}(x)\phi(x)\mathrm{d}x = \sum_{j=1}^{k}P_{\alpha}(x_{j})\omega_{j}\phi(x_{j}),$$
\item for any suitable function $f(\cdot)$, there exists some $\xi\in\mathbb{R}$ such that
$$\int_{\mathbb{R}}f(x)\phi(x)\mathrm{d}x \equiv \sum_{j=1}^{k}f(x_{j})\phi(x_{j})\omega_{j} + f^{(2k)}(\xi)E_{k},$$
where $E_{k}$ is an error term depending on $k$ with $E_{k}\downarrow 0$ as $k\uparrow\infty$.
\end{itemize}

\hypertarget{adaptive-gauss-hermite-quadrature}{%
\subsection{Adaptive Gauss-Hermite
Quadrature}\label{adaptive-gauss-hermite-quadrature}}

In Bayesian inference, \(\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y})\)
has a location and curvature which vary stochastically with
\(\boldsymbol{Y}\), and no fixed quadrature rule will perform well for
all datasets and models. A rule which adapts to the changing shape and
location of \(\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y})\) is
desirable.

Adaptive Gauss-Hermite Quadrature (\texttt{AGHQ}{}; \citet{nayloradaptive,adaptive_GH_1994,adaptive_GH_2020,aghqus}) shifts and scales
the \texttt{GHQ}{} nodes and weights according to the mode and curvature
of \(\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y})\). Define
\begin{equation}\begin{aligned}
\widetilde{\boldsymbol{\theta}}_n&= \text{argmax}_{\boldsymbol{\theta}}\log\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y}), \\
\boldsymbol{H}&= -\partial^{2}_{\boldsymbol{\theta}}\log\pi^{*}(\widetilde{\boldsymbol{\theta}}_n,\boldsymbol{Y}), \\
\boldsymbol{H}^{-1} &= \boldsymbol{L}\boldsymbol{L}^{\tiny{\texttt{T}}}
\end{aligned}\end{equation} where \(\boldsymbol{L}\) is the (unique)
lower Cholesky triangle. The \AGHQ{} nodes and weights are defined by shifting the \GHQ{}
nodes by $\widetilde{\boldsymbol{\theta}}_n$ and rescaling the nodes and weights by $\mb{L}$. This can be implemented making use of the
\texttt{mvQuad::rescale} function, as in the following example with a
fake \(\widetilde{\boldsymbol{\theta}}_n\) and \(\boldsymbol{H}\):

\begin{CodeChunk}
\begin{CodeInput}
R> # Mode
R> theta_hat <- c(2,3)
R> # Hessian
R> thehess <- matrix(c(3,1,1,5),nrow=2,ncol=2)
R> thehessinverse <- solve(thehess)
R> # Int grid for rescaling using software
R> k <- 3
R> intgridsoftware <- mvQuad::createNIGrid(2,"GHe",k)
R> mvQuad::rescale(intgridsoftware,m = theta_hat,C = thehessinverse,dec.type=2)
R> # Int grid for rescaling manually
R> intgridmanual <- mvQuad::createNIGrid(2,"GHe",k)
R>
R> # Nodes
R> # Note that chol() gives L^T, upper triangle
R> nnrescaled <- mvQuad::getNodes(intgridmanual)
R> nnrescaled[ ,1] <- nnrescaled[ ,1] + theta_hat[1]
R> nnrescaled[ ,2] <- nnrescaled[ ,2] + theta_hat[2]
R>
R> mvQuad::getNodes(intgridsoftware)
\end{CodeInput}
\begin{CodeOutput}
       [,1] [,2]
 [1,] 0.965 2.43
 [2,] 2.000 2.23
 [3,] 3.035 2.02
 [4,] 0.965 3.21
 [5,] 2.000 3.00
 [6,] 3.035 2.79
 [7,] 0.965 3.98
 [8,] 2.000 3.77
 [9,] 3.035 3.57
\end{CodeOutput}
\begin{CodeInput}
R> nnrescaled
\end{CodeInput}
\begin{CodeOutput}
       [,1] [,2]
 [1,] 0.965 2.43
 [2,] 2.000 2.23
 [3,] 3.035 2.02
 [4,] 0.965 3.21
 [5,] 2.000 3.00
 [6,] 3.035 2.79
 [7,] 0.965 3.98
 [8,] 2.000 3.77
 [9,] 3.035 3.57
\end{CodeOutput}
\begin{CodeInput}
R> # Weights
R> wwrescaled <- mvQuad::getWeights(intgridmanual) * det(chol(thehessinverse))
R> mvQuad::getWeights(intgridsoftware)
\end{CodeInput}
\begin{CodeOutput}
       [,1]
 [1,] 0.937
 [2,] 0.836
 [3,] 0.937
 [4,] 0.836
 [5,] 0.746
 [6,] 0.836
 [7,] 0.937
 [8,] 0.836
 [9,] 0.937
\end{CodeOutput}
\begin{CodeInput}
R> wwrescaled
\end{CodeInput}
\begin{CodeOutput}
       [,1]
 [1,] 0.937
 [2,] 0.836
 [3,] 0.937
 [4,] 0.836
 [5,] 0.746
 [6,] 0.836
 [7,] 0.937
 [8,] 0.836
 [9,] 0.937
\end{CodeOutput}
\end{CodeChunk}

Some interesting properties of \texttt{AGHQ}{} are as follows
\citep{aghqus}:

\begin{itemize}
\item
  The mode \(\widetilde{\boldsymbol{\theta}}_n\) is a quadrature point
  if and only if \(k\) is odd,
\item
  When \(k = 1\), \texttt{AGHQ}{} is a Laplace approximation,
\item
  When \(k > 1\), \texttt{AGHQ}{} achieves higher-order asymptotic
  accuracy: \citet{aghqus} show that if
  \(\widetilde{\pi}_{\texttt{AGHQ}}(\boldsymbol{Y})\) is an
  \texttt{AGHQ}{} approximation to \(\pi(\boldsymbol{Y})\) using \(k\)
  points, then: \begin{equation}
  \widetilde{\pi}_{\texttt{AGHQ}}(\boldsymbol{Y}) = \pi(\boldsymbol{Y})\left[1 + \mathcal{O}_{p}\left( n^{-\lfloor (k+2)/3\rfloor}\right) \right],
  \end{equation} where the convergence is in probability measured with
  respect to the distribution of \(\boldsymbol{Y}\). Table
  \ref{tab:ratetable} shows the asymptotic accuracy achieved for
  selected \(k\). Note that the well-known \(n^{-1}\) rate for the
  Laplace approximation is recovered when \(k = 1\).
\end{itemize}

\begin{table}
\centering
\begin{tabular}{|l|r|}
\hline
$k$ & Rate \\
\hline
1 & $\mathcal{O}_{p}\left(n^{-1}\right)$ \\
3 & $\mathcal{O}_{p}\left(n^{-1}\right)$ \\
5 & $\mathcal{O}_{p}\left(n^{-2}\right)$ \\
7 & $\mathcal{O}_{p}\left(n^{-3}\right)$ \\
9 & $\mathcal{O}_{p}\left(n^{-3}\right)$ \\
11 & $\mathcal{O}_{p}\left(n^{-4}\right)$ \\
13 & $\mathcal{O}_{p}\left(n^{-5}\right)$ \\
15 & $\mathcal{O}_{p}\left(n^{-5}\right)$ \\
\hline
\end{tabular}
\caption{Asymptotic rates for \texttt{AGHQ}{}}
\label{tab:ratetable}
\end{table}

\hypertarget{posterior-summaries}{%
\subsection{Posterior summaries}\label{posterior-summaries}}

With the \texttt{AGHQ}{} approximation
\(\widetilde{\pi}_{\texttt{AGHQ}}(\boldsymbol{Y})\) to
\(\pi(\boldsymbol{Y})\), the approximate posterior density is:
\begin{equation}
\widetilde{\pi}_{\texttt{AGHQ}}(\boldsymbol{\theta}|\boldsymbol{Y}) = \frac{\pi^{*}(\boldsymbol{\theta},\boldsymbol{Y})}{\widetilde{\pi}_{\texttt{AGHQ}}(\boldsymbol{Y})},
\end{equation} and, by definition, satisfies: \begin{equation}
\sum_{j=1}^{k}\omega_{j}\widetilde{\pi}_{\texttt{AGHQ}}(\boldsymbol{\theta}_{j}|\boldsymbol{Y}) = 1.
\end{equation} For any function \(f:\mathbb{R}^{d}\to\mathbb{R}\), the
posterior moment
\(\mathop{\mathrm{\mathbb{E}}}(f(\boldsymbol{\theta})|\boldsymbol{Y})\)
may be approximated by \begin{equation}
\mathop{\mathrm{\mathbb{E}}}(f(\boldsymbol{\theta})|\boldsymbol{Y})\approx \sum_{j=1}^{k}f(\boldsymbol{\theta}_{j})\omega_{j}\widetilde{\pi}_{\texttt{AGHQ}}(\boldsymbol{\theta}_{j}|\boldsymbol{Y}).
\end{equation} \citet{aghqus} prove that the density and moment
approximations attain the same relative rate of asymptotic convergence
as the approximation to the normalizing constant, when computed in this
manner.

Marginal posteriors and quantiles are more difficult to compute, as described in
detail by \citet{nayloradaptive} and \citet{aghqus}, Appendix E. The challenge is that the quadrature
rule is only defined at specific points
\(\boldsymbol{\boldsymbol{\theta}}_{j}\), and these points are not on a
regular grid, due to having been multiplied by \(\boldsymbol{L}\). It is
therefore only straightforward to evaluate marginal posteriors at \(k\)
distinct points, and this only works for the first element of
\(\boldsymbol{\theta}\). To mitigate these challenges, the parameter
vector has to be re-ordered and the joint normalization re-computed for
each desired marginal distribution, and a polynomial interpolant based
on these points is used for plotting and computing quantiles. The
\pkg{aghq} package handles these difficulties for the user, as I now
describe.

For any \(t\in[d]\), let \(\boldsymbol{\theta}^{-t}\) be the parameter
vector with the \(t^{th}\) index removed. The approximate marginal
posterior
\(\widetilde{\pi}_{\texttt{AGHQ}}^{t}(\theta_{t}|\boldsymbol{Y})\) is
approximated as follows: \begin{equation}
\widetilde{\pi}_{\texttt{AGHQ}}^{t}(\theta_{t}|\boldsymbol{Y}) \approx \sum_{j=1}^{k}\omega_{j}\widetilde{\pi}_{\texttt{AGHQ}}\left[(\theta_{t},\boldsymbol{\theta}^{-t}_{j})|\boldsymbol{Y}\right],
\end{equation} where the parameter vector has been re-ordered so
\(\theta_{t}\) is the first component. This only provides \(k\) unique evaluations of
\(\widetilde{\pi}_{\texttt{AGHQ}}^{t}(\theta_{t}|\boldsymbol{Y})\) and
is insufficient for plotting and for computation of the marginal CDF and
quantiles. To obtain evaluations at any value of \(\theta_{t}\), a
\(k^{th}\)-degree Lagrange polynomial interpolant is fit to
\(\log\widetilde{\pi}_{\texttt{AGHQ}}^{t}(\theta_{t}|\boldsymbol{Y})\)
using the \texttt{poly.calc} function in the \pkg{polynom} package
\citep{polynom}. This is also useful in the case where each each
evaluation of
\(\widetilde{\pi}_{\texttt{AGHQ}}(\boldsymbol{\theta}|\boldsymbol{Y})\)
is computationally expensive, such as in \S\ref{sec:elgms}.

To obtain the marginal cumulative distribution function
\(F_{t}(x) = \mathop{\mathrm{\mathbb{P}}}(\theta_{t} \leq x|\boldsymbol{Y})\),
a fine grid \(\left\{ x_{l} \right\}_{l=1}^{L}\) and simple quadrature
are used: \begin{equation}
F_{t}(x) \approx \sum_{x_{l}\leq x}\widetilde{\pi}_{\texttt{AGHQ}}^{t}(x_{l}|\boldsymbol{Y})(x_{l} - x_{l-1})
\end{equation} The default grid has \(L = 1000\) and covers the range of
the quadrature points plus a buffer, and the user may input their own
grid if this gives undesirable results. This can easily be verified in
any specific application by checking that \(F_{t}(x_{1})\approx 0\),
\(F_{t}(x_{L})\approx 1\), and \(F_{t}(x_{l}) \leq F_{t}(x_{l+1})\) for
each \(l\in[L-1]\), properties which are immediately revealed by a plot
of \(F_{t}(\cdot)\).

With the marginal CDF available, computing quantiles is straightforward.
Let \(q_{t}^{\alpha}\) be the \(\alpha\) quantile of
\(\widetilde{\pi}_{\texttt{AGHQ}}^{t}(\theta_{t}|\boldsymbol{Y})\), such
that
\(\int_{-\infty}^{q_{t}^{\alpha}}\widetilde{\pi}_{\texttt{AGHQ}}^{t}(\theta_{t}|\boldsymbol{Y})\mathrm{d}\theta_{t} = \alpha\).
Then we compute \begin{equation}
q_{t}^{\alpha} \approx \text{min}\left\{ x_{l}: F_{t}(x_{l}) \geq \alpha, l\in[L] \right\}
\end{equation} If
\(\widetilde{\alpha} = \int_{-\infty}^{\infty}\cdots\int_{-\infty}^{q_{t}^{\alpha}}\pi(\boldsymbol{\theta}|\boldsymbol{Y})\mathrm{d}\theta_{t}\mathrm{d}\boldsymbol{\theta}^{-t}\)
is the true coverage probability of \(q_{t}^{\alpha}\) and
\(\widetilde{q}_{t}^{\alpha}\) is the corresponding \(\alpha\)-quantile
of the true marginal posterior \(\pi(\theta_{t}|\boldsymbol{Y})\), then
for any \(0 \leq \alpha \leq 1\), \citet{aghqus} guarantee that
\(\left\lvert \alpha - \widetilde{\alpha} \right\rvert\) and
\(\left\lvert q_{t}^{\alpha} - \widetilde{q}_{t}^{\alpha} \right\rvert\)
both converge to zero in probability at the same rate as the normalizing
constant. These theorems do not account for the further approximation
introduced to compute the quantiles, however.

The computation of the approximate marginal CDF, PDF and quantiles is
done using the \texttt{compute\_pdf\_and\_cdf} and
\texttt{compute\_quantiles} functions as described in
\S\ref{sec:basic-use}. Internally, the computations proceed as follows:

\begin{CodeChunk}
\begin{CodeInput}
R> # Get the marginal posterior
R> (margpost <- thequadrature$marginals[[1]])
\end{CodeInput}
\begin{CodeOutput}
# A tibble: 3 x 3
  theta1 logmargpost     w
   <dbl>       <dbl> <dbl>
1   1.25      -0.357 0.267
2   1.49       1.03  0.239
3   1.74      -0.605 0.267
\end{CodeOutput}
\begin{CodeInput}
R> # aghq::interpolate_marginal_posterior() constructs the polynomial
R> # interpolant, returning a function.
R> margpostinterp <- interpolate_marginal_posterior(margpost)
R> # Construct a fine grid covering the quadrature points, plus a bit.
R> rn <- range(margpost$theta1)
R> rnl <- diff(rn)
R> thetarange <- c(min(rn) - rnl/2,max(rn) + rnl/2)
R> finegrid <- c(seq(thetarange[1],thetarange[2],length.out=1000))
R> # Approximate the CDF (and PDF)
R> approxdens <- dplyr::tibble(theta = finegrid,
+                pdf = exp(margpostinterp(.data$theta)),
+                cdf = cumsum(.data$pdf * c(0,diff(.data$theta))))
R> # Check it goes from 0 to 1...
R> range(approxdens$cdf)
\end{CodeInput}
\begin{CodeOutput}
[1] 0.000 0.999
\end{CodeOutput}
\begin{CodeInput}
R> # ...and that it is (numerically) nondecreasing:
R> all(approxdens$cdf[2:1000] > approxdens$cdf[1:999])
\end{CodeInput}
\begin{CodeOutput}
[1] TRUE
\end{CodeOutput}
\begin{CodeInput}
R> # Approximate the quantiles
R> q <- c(.01,.25,.50,.75,.99)
R> quants <- numeric(length(q))
R> names(quants) <- paste0(as.character(100 * q),
R>
R> for (i in 1:length(q)) {
+   quants[i] <- approxdens$theta[max(which(approxdens$cdf < q[i]))]
+ }
R> quants
\end{CodeInput}
\begin{CodeOutput}

1.15 1.39 1.48 1.58 1.82
\end{CodeOutput}
\begin{CodeInput}
R> # For comparison to before:
R> exp(quants)
\end{CodeInput}
\begin{CodeOutput}

3.17 4.00 4.40 4.85 6.15
\end{CodeOutput}
\end{CodeChunk}

It is an open question whether computing the quantiles in this way
induces non-negligible approximation error,
and whether there is a better method for computing them. As such methods
are developed, they will be implemented in the \pkg{aghq} package.

\hypertarget{sec:examples}{%
\section{Examples}\label{sec:examples}}

In this section I give the compuational details for the two difficult examples considered by \citet{aghqus} of models for which
\texttt{AGHQ}{} provides a complete methodology for approximate Bayesian inference.
In both cases the dominant method in the literature is \texttt{MCMC}{},
and in the first example I provide an explicit comparison between the
two where I demonstrate the computational advantages of \texttt{AGHQ}{}.
In the second example the focus is on how to approach a difficult
constrained optimization and compute interesting posterior summary
statistics. A more detailed discussion of both these examples appears in \citet{aghqus}, and the focus of this section is to provide readers complete examples of complicated models
that they can immediately implement themselves using the \pkg{aghq} package and the provided code.

\hypertarget{sec:inf-disease}{%
\subsection{Infectious disease models}\label{sec:inf-disease}}

\citet{aghqus} implement \texttt{AGHQ}{} for a
\emph{Susceptible, Infectious, Removed} (SIR) model for the spread of
infectious disease and compared to \texttt{MCMC}{} for the same model as
implemented in the \pkg{EpiILMCT} \texttt{R} package \citep{epi}. Here
we illustrate how to use the \pkg{aghq} package to fit this model. The
run time is less than half a second and no user tuning is required. In
contrast, \texttt{MCMC}{} does require tuning and although the run time
is of course dependent on the number of iterations (which the user also
has to somehow choose), tends to run far slower than the run time of \AGHQ{}.

We consider their example of an outbreak of the Tomato Spotted Wilt Virus in
\(n = 520\) plants, of which \(n_{0} = 327\) were infected. For observed
infection times \(I_{1} \leq \cdots \leq I_{n_{0}}\) and
\(I_{i} = \infty, i = n_{0}+1,\ldots,n\), and removal times
\(R_{i},i\in[n]\), the likelihood is: \begin{equation}
\pi(\boldsymbol{I},\boldsymbol{R} | \alpha,\beta) = \prod_{j=2}^{n_{0}}\left( \sum_{i:I_{i} < I_{j} \leq R_{i}}\lambda_{ij}\right)\exp\left\{ -\sum_{i=1}^{n_{0}}\sum_{j=1}^{n}[\text{min}(R_{i},I_{j}) - \text{min}(I_{i},I_{j})]\lambda_{ij} \right\},
\end{equation} where \(\lambda_{ij} = \alpha d_{ij}^{-\beta}\) and
\(d_{ij}\) is the Euclidean distance between plants \(i\) and \(j\).
Independent \(\text{Exponential}(0.01)\) priors are placed on the
parameters of interest \((\alpha,\beta)\), and quadrature is performed
on the transformed scale with
\(\boldsymbol{\theta}= (\log\alpha,\log\beta)\). See \citet{aghqus} and
\citet{epi} for further information.

To implement \texttt{AGHQ}{} for this problem using the \pkg{aghq}
package, we

\begin{enumerate}
\def\arabic{enumi}.{\arabic{enumi}.}
\item
  Load the \texttt{tswv} data from the \pkg{EpiILMCT} package,
\item
  Define a list of functions \texttt{ff} which compute the log-posterior
  and its first two derivatives, and
\item
  Call the \texttt{aghq} function.
\end{enumerate}

To obtain \texttt{ff} we utilize the automatic differentiation
capabilities of the \texttt{TMB} package \citep{tmb}. The user writes a
\texttt{c++} template which computes the log-posterior and passes it
along with the data to the \texttt{TMB::MakeADFun} function, which
outputs the list \texttt{ff} in the correct format for use within
\texttt{aghq::aghq}. This is illustrated here, with the \texttt{TMB}
function template deferred to Appendix \ref{app:inf-dis-tmb}. The
following illustration focusses on how to run the code; a more detailed
methodological comparison between \texttt{AGHQ}{} and \texttt{MCMC}{}
for this example is given by \citet{aghqus}.

\begin{CodeChunk}
\begin{CodeInput}
R> library(TMB)
R> # Get the Tomato disease data
R> data("tswv", package = "EpiILMCT")
R> # Use TMB
R> # Put the objective function code in a subdirectory code/ of your
R> # working directory
R> precompile()
\end{CodeInput}
\begin{CodeOutput}
Removing: libTMB.cpp libTMBdbg.cpp libTMBomp.cpp
\end{CodeOutput}
\begin{CodeInput}
R> # Compile: only need to do if source has changed
R> # compile("code/02_disease.cpp")
R> dyn.load(dynlib("code/02_disease"))
R> # Create the functions
R> dat <- tswv$tswvsir
R> dat$epidat <- dat$epidat[order(dat$epidat[ ,4]), ]
R>
R> I <- dat$epidat[ ,4]
R> R <- dat$epidat[ ,2]
R> infected <- !is.infinite(I)
R>
R> datlist <- list(
+   D = as.matrix(dist(dat$location[dat$epidat[ ,1], ])),
+   I = I,
+   R = R,
+   infected = as.numeric(infected[infected])
+ )
R>
R> ff <- MakeADFun(data = datlist,
+                 parameters = list(theta1 = 0,theta2 = 0),
+                 DLL = "02_disease",
+                 ADreport = FALSE,
+                 silent = TRUE)
R> # Do the quadrature
R> tm <- Sys.time()
R> tmp <- capture.output(# Suppress optimization printing
+   ilmquadrature <- aghq::aghq(ff,7,c(0,0))
+ )
R> runtime <- difftime(Sys.time(),tm,units = "secs")
R> cat("Runtime of AGHQ with k = 7:",runtime,"seconds\n")
\end{CodeInput}
\begin{CodeOutput}
Runtime of AGHQ with k = 7: 0.485 seconds
\end{CodeOutput}
\begin{CodeInput}
R> summary(ilmquadrature)
\end{CodeInput}
\begin{CodeOutput}
AGHQ on a 2 dimensional posterior with  7 7 quadrature points

The posterior mode is: -4.39 0.291

The log of the normalizing constant/marginal likelihood is: -1088

The posterior Hessian at the mode is:
     [,1] [,2]
[1,]  327 -532
[2,] -532  948

The covariance matrix used for the quadrature is...
       [,1]   [,2]
[1,] 0.0349 0.0196
[2,] 0.0196 0.0120

...and its Cholesky is:
      [,1]   [,2]
[1,] 0.187 0.1048
[2,] 0.000 0.0325

Here are some moments and quantiles for theta:

         mean median   mode    sd    2.
theta1 -4.439 -4.426 -4.386 0.200 -4.8808 -4.096
theta2  0.258  0.268  0.291 0.121 -0.0175  0.459
\end{CodeOutput}
\end{CodeChunk}

I have loaded the \texttt{MCMC}{} results from disk (code in Appendix
\ref{app:inf-dis}) and saved posterior samples in a \texttt{data.frame}
called \texttt{postsamps}. Now, compare moments and the density plots:

\begin{CodeChunk}
\begin{CodeInput}
R> # Posterior mean and quantiles for alpha = exp(theta1) and beta = exp(theta2)
R> # AGHQ:
R> ilm_means <- compute_moment(ilmquadrature$normalized_posterior,function(x) exp(x))
R> names(ilm_means) <- c('alpha','beta')
R> ilm_means
\end{CodeInput}
\begin{CodeOutput}
alpha  beta
0.012 1.304
\end{CodeOutput}
\begin{CodeInput}
R> # MCMC:
R> postsamps
\end{CodeInput}
\begin{CodeOutput}
# A tibble: 1 x 2
   alpha  beta
   <dbl> <dbl>
1 0.0120  1.30
\end{CodeOutput}
\begin{CodeInput}
R> # Quantiles, AGHQ, alpha:
R> ilm_quants_alpha <- exp(compute_quantiles(ilmquadrature$marginals[[1]]))
R> ilm_quants_beta <- exp(compute_quantiles(ilmquadrature$marginals[[2]]))
R> ilm_quants_alpha
\end{CodeInput}
\begin{CodeOutput}
   2.
0.00759 0.01665
\end{CodeOutput}
\begin{CodeInput}
R> # MCMC:
R> quantile(postsamps$alpha,c(.025,.975))
\end{CodeInput}
\begin{CodeOutput}
   2.
0.00742 0.01682
\end{CodeOutput}
\begin{CodeInput}
R> # and for beta:
R> ilm_quants_beta
\end{CodeInput}
\begin{CodeOutput}
 2.
0.983 1.582
\end{CodeOutput}
\begin{CodeInput}
R> quantile(postsamps$beta,c(.025,.975))
\end{CodeInput}
\begin{CodeOutput}
 2.
0.974 1.593
\end{CodeOutput}
\begin{CodeInput}
R> # Density plots
R> translist <- list(totheta = function(x) log(x),fromtheta = function(x) exp(x))
R> ilmpdfandcdf_alpha <- compute_pdf_and_cdf(ilmquadrature$marginals[[1]],translist)
R> ilmpdfandcdf_beta <- compute_pdf_and_cdf(ilmquadrature$marginals[[2]],translist)
R>
R> mcmcplotalpha <- postsamps
+   ggplot(aes(x = alpha)) +
+   theme_classic() +
+   geom_histogram(aes(y = ..density..),
+                  bins = 100,
+                  colour = "transparent",
+                  fill = "grey",
+                  alpha = .5) +
+     geom_line(
+       data = ilmpdfandcdf_alpha,
+       mapping = aes(x = transparam,y = pdf_transparam)
+     ) +
+     labs(title = "",x = expression(alpha),y = "Density") +
+     coord_cartesian(xlim = c(0.005,0.0255),ylim = c(0,200))
R>
R> mcmcplotbeta <- postsamps
+   ggplot(aes(x = beta)) +
+   theme_classic() +
+   geom_histogram(aes(y = ..density..),
+                  bins = 100,
+                  colour = "transparent",
+                  fill = "grey",
+                  alpha = .5) +
+     geom_line(
+       data = ilmpdfandcdf_beta,
+       mapping = aes(x = transparam,y = pdf_transparam)
+     ) +
+     labs(title = "",x = expression(beta),y = "Density") +
+     coord_cartesian(xlim = c(0.8,2.05),ylim = c(0,3))
R>
R> mcmcplotalpha | mcmcplotbeta
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{01-aghq-softwarepaper_files/figure-latex/ilmsummary1-1}

}

\caption[MCMC ($\textcolor{lightgray}{\blacksquare}$) and AGHQ (---) results for infectious disease model]{MCMC ($\textcolor{lightgray}{\blacksquare}$) and AGHQ (---) results for infectious disease model.}\label{fig:ilmsummary1}
\end{figure}
\end{CodeChunk}

This code forms the basis for the extended comparison between
\texttt{MCMC}{} and \texttt{AGHQ}{} for this example given by
\citet{aghqus}.

Joint posterior moments may also be approximated: for example, for two
plants which are \(d_{ij} = 2\) units apart, the infectivity rate is
\(\lambda_{ij} = \alpha\cdot2^{-\beta}\). We can approximate
\(\mathop{\mathrm{\mathbb{E}}}(\lambda_{ij}|\boldsymbol{Y})\) as:

\begin{CodeChunk}
\begin{CodeInput}
R> compute_moment(
+   ilmquadrature$normalized_posterior,
+   function(x) exp(x)[1] * 2^(-exp(x)[2])
+ )
\end{CodeInput}
\begin{CodeOutput}
[1] 0.00481
\end{CodeOutput}
\begin{CodeInput}
R> # Compare to MCMC
R> mean(postsamps$alpha * 2^(-postsamps$beta))
\end{CodeInput}
\begin{CodeOutput}
[1] 0.00481
\end{CodeOutput}
\end{CodeChunk}

\hypertarget{subsec:astro}{%
\subsection{Estimating the mass of the Milky Way
Galaxy}\label{subsec:astro}}

Astrophysicists are interested in estimating the mass of the Milky Way
galaxy. \citet{gwen2} describe a statistical model which relates
multivariate position/velocity measurements of star clusters orbiting
the Galaxy to parameters which in turn determine its mass according to a
physical model. Their specific interest in this work is in assessing the
impact of various prior choices and data inclusion rules on inference,
and they fit and assess many models. Despite this, they rely on
\texttt{MCMC}{} for inference, with long run times and requiring user
tuning for each new fit. In contrast, \citet{aghqus} show that
\texttt{AGHQ}{} is fast and stable for this problem and much better
suited to the type of rapid prototyping done by \citet{gwen2}.

Implementing \texttt{AGHQ}{} for this problem requires care, however.
The parameters are subject to various box and nonlinear constraints and
in certain cases have a sharply peaked mode lying close to the boundary
of the parameter space. Here we demonstrate how to do a constrained
optimization using \pkg{TMB} and the \pkg{ipoptr} package which
interfaces the \texttt{IPOPT} software \citep{ipopt} and then pass the
results to \pkg{aghq} for inference. This is an example where the model
is too complex for \texttt{aghq::aghq()} to handle it on its own, but
the \pkg{aghq} package still provides a flexible interface for Bayesian
inference.

Let \(\boldsymbol{Z}_{i} = (z_{i1},\ldots,z_{i4})\) denote the position
and three-component velocity measured relative to the sun
(\emph{heliocentric} measurements) for the \(i^{th}\) star cluster,
\(i\in[n], n = 70\). These are what is actually measured. These
measurements are then transformed into position and radial and
tangential velocity measurements relative to the centre of the Galaxy
(\emph{galactocentric} measurements)
\(\boldsymbol{Y}_{i} = (y_{i1},y_{i2},y_{i3})\). The probability density
for \(\boldsymbol{Y}_{i}\) is \begin{equation}
f(\boldsymbol{Y}_{i};\Psi_{0},\gamma,\alpha,\beta) = \frac{L_{i}^{-2\beta}\mathcal{E}_{i}^{\frac{\beta(\gamma-2)}{\gamma} + \frac{\alpha}{\gamma} - \frac{3}{2}}\Gamma\left( \frac{\alpha}{\gamma} - \frac{2\beta}{\gamma} + 1\right)} {\sqrt{8\pi^{3}2^{-2\beta}} \Psi_{0}^{-\frac{2\beta}{\gamma} + \frac{\alpha}{\gamma}}\Gamma\left( \frac{\beta(\gamma - 2)}{\gamma} + \frac{\alpha}{\gamma} - \frac{1}{2}\right)},
\end{equation} where \(L_{i} = y_{i1}y_{i3}\) and
\(\mathcal{E}_{i} = \Psi_{0}y_{i1}^{1-\gamma} - (y_{i2}^{2} + y_{i3}^{2})/2\).
The parameters \(\boldsymbol{\Xi} = (\Psi_{0},\gamma,\alpha,\beta)\)
determine the mass of the Galaxy at radial distance \(r\) kiloparsecs
(kpc) as: \begin{equation}
M_{\boldsymbol{\Xi}}(r) = \Psi_{0}\gamma r^{1-\gamma},
\end{equation} according to a physical model, and inference is done for
\(M_{\boldsymbol{\Xi}}(r)\) at various values of \(r\).

The parameters are subject to nonlinear and box constraints:
\begin{equation}\begin{aligned}
\alpha &>\gamma, \\
\alpha &> \beta(2 - \gamma) + \frac{\gamma}{2}, \\
\mathcal{E}_{i} &> 0, i\in[n], \\
\alpha &\geq 3, \\
1 \leq \Psi &\leq 200, \\
0.3 \leq \gamma &\leq 0.7, \\
-0.5 \leq \beta &\leq 1.
\end{aligned}\end{equation} Uniform priors on these ranges are used for
\(\Psi,\gamma,\beta\), and a \(\text{Gamma}(1,4.6)\) prior is used for
\(\alpha - 3\). We perform quadrature on the transformed scale
\(\boldsymbol{\theta}= (\theta_{1},\theta_{2},\theta_{3},\theta_{4})\)
with \begin{equation}\begin{aligned}
\theta_{1} &= \log\left(-\log\left[ \frac{\Psi - 1}{200 - 1}\right]\right), \\
\theta_{2} &= \log\left(-\log\left[ \frac{\gamma - 0.3}{0.7 - 0.3}\right]\right), \\
\theta_{3} &= \log(\alpha - 3), \\
\theta_{4} &= \log\left(-\log\left[ \frac{\beta + 0.5}{1 + 0.5}\right]\right). \\
\end{aligned}\end{equation} This transformation enhances the spherical
symmetry of the posterior and we find it to improve the stability of the
optimization and the empirical quality of the \texttt{AGHQ}{}
approximation, as suggested by \citet{nayloradaptive}. The nonlinear
constraints remain, however. We perform constrained quasi-newton
optimization using the \pkg{ipoptr} package. This requires the objective
function along with its gradient and Hessian, as well as the Jacobian
matrix of the nonlinear constraints. These can all be obtained from a
single \texttt{TMB} template as follows:

\begin{enumerate}
\def\arabic{enumi}.{\arabic{enumi}.}
\item
  The \texttt{TMB} objective function template is created as usual,
\item
  Inside this template, the user defines a vector which contains the
  values of each nonlinear constraint, and calls \texttt{ADREPORT()} on
  this vector,
\item
  \texttt{TMB::MakeADFun} is called twice: once with
  \texttt{ADreport=FALSE} and once with \texttt{ADreport=TRUE}. The
  former creates the \texttt{ff} function list and the latter returns a
  function list which evaluates the nonlinear constraints and their
  Jacobian matrix.
\end{enumerate}

The data are originally tabulated by \citet{gwen2} and are included in
the \pkg{aghq} package in both the raw format (variables
\(\boldsymbol{Z}_{i}\), \texttt{gcdata} object) and also after the
extensive, problem-specific pre-processing required for this example has
been completed (variables \(\boldsymbol{Y}_{i}\), \texttt{gcdatalist}
object). Appendix \ref{app:astrotmb} shows the \texttt{TMB} template for
this example. We then implement optimization, deferring code to Appendix
\ref{app:astro}.

In contrast to the other examples in this manuscript, the optimization
is extraordinarily difficult: it is sensitive to starting values, and
requires box constraints to ensure stability. In contrast to
\texttt{MCMC}{}, however, once these issues are solved using all of the
above transformations and constraints, they do not have to be revisited
each time the model is run, which facilitates rapid prototyping of
models. The run time is very fast:

\begin{CodeChunk}
\begin{CodeOutput}
Run time for mass model optimization: 0.462 seconds.
\end{CodeOutput}
\end{CodeChunk}

The \pkg{aghq} package does not try to support optimizations of this
difficulty; rather, it supports user-supplied optimization results as
follows. The \texttt{ff} object is the return value of
\texttt{TMB::MakeADFun}, as in Example \ref{sec:inf-disease}. The
\texttt{ipopt\_result} object containing the mode of the optimization
and other information.

\begin{CodeChunk}
\begin{CodeInput}
R> # Create the optimization template
R> useropt <- list(
+   ff = ff,
+   mode = ipopt_result$solution,
+   hessian = -1 * ff$he(ipopt_result$solution)
+ )
R> # Do the quadrature
R> tm <- Sys.time()
R> astroquad <- aghq::aghq(ff,3,thetastart,optresults = useropt)
R> runtime <- difftime(Sys.time(),tm,units = 'secs')
R> cat("Run time for mass model quadrature:",runtime,"seconds.\n")
\end{CodeInput}
\begin{CodeOutput}
Run time for mass model quadrature: 0.168 seconds.
\end{CodeOutput}
\begin{CodeInput}
R> summary(astroquad)
\end{CodeInput}
\begin{CodeOutput}
AGHQ on a 4 dimensional posterior with  3 3 3 3 quadrature points

The posterior mode is: 0.613 1.16 -4.75 0.0282

The log of the normalizing constant/marginal likelihood is: -706

The posterior Hessian at the mode is:
        [,1]    [,2]     [,3]     [,4]
[1,] 455.277  8.4792  0.63972 -0.56376
[2,]   8.479  8.7865  0.01546 -1.21210
[3,]   0.640  0.0155  1.00373 -0.00709
[4,]  -0.564 -1.2121 -0.00709 23.50307

The covariance matrix used for the quadrature is...
           [,1]      [,2]      [,3]       [,4]
[1,]  0.0022387 -0.002166 -0.001394 -0.0000584
[2,] -0.0021660  0.116725 -0.000376  0.0059677
[3,] -0.0013939 -0.000376  0.997183  0.0002478
[4,] -0.0000584  0.005968  0.000248  0.0428541

...and its Cholesky is:
       [,1]    [,2]     [,3]      [,4]
[1,] 0.0473 -0.0458 -0.02946 -0.001235
[2,] 0.0000  0.3386 -0.00509  0.017459
[3,] 0.0000  0.0000  0.99814  0.000301
[4,] 0.0000  0.0000  0.00000  0.206271

Here are some moments and quantiles for theta:

          mean median    mode     sd   2.
theta1  0.6057  0.606  0.6126 0.0463  0.515  0.696
theta2  1.1946  1.191  1.1585 0.3373  0.529  1.856
theta3 -5.1940 -5.208 -4.7475 0.9012 -6.939 -3.497
theta4  0.0401  0.039  0.0282 0.2098 -0.372  0.450
\end{CodeOutput}
\begin{CodeInput}
R> # plot(astroquad) # Produces 8 plots, too large for paper
R> # Manually plot
R> pdfandcdf <- vector(mode = 'list',length = 4)
R> par(mfrow = c(2,2))
R> for (j in 1:4) {
+   with(compute_pdf_and_cdf(astroquad$marginals[[j]]),
+        plot(pdf~theta,
+             type='l',
+             main = paste0('Theta ',j),
+             xlab = expression(theta))
+        )
+ }
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{01-aghq-softwarepaper_files/figure-latex/astroquad11-1}

}

\caption[Marginal posteriors on transformed scale for mass estimation example]{Marginal posteriors on transformed scale for mass estimation example.}\label{fig:astroquad11}
\end{figure}
\end{CodeChunk}

While the results are not meaningful on the transformed scale, it is
good to check that the marginal posteriors are visually ``nice'' and
log-concave. We may back-transform to get the results on the original
scale as follows:

\begin{CodeChunk}
\begin{CodeInput}
R> ## Define the parameter transformation functions ##
R> parambounds <- list(
+   Psi0 = c(1,200),
+   gamma = c(.3,.7),
+   alpha = c(3.0,3.7),
+   beta = c(-.5,1)
+ )
R>
R> get_psi0 <- function(theta)
+   (parambounds$Psi0[2] - parambounds$Psi0[1]) *
+     exp(-exp(theta)) + parambounds$Psi0[1]
R> get_theta1 <- function(Psi0)
+   log(-log( (Psi0 - parambounds$Psi0[1]) /
+               (parambounds$Psi0[2] - parambounds$Psi0[1]) ))
R>
R> get_gamma <- function(theta)
+   (parambounds$gamma[2] - parambounds$gamma[1]) *
+   exp(-exp(theta)) + parambounds$gamma[1]
R> # Add a little buffer, for stability
R> get_theta2 <- function(gamma)
+   log(-log( (gamma - parambounds$gamma[1] + 1e-03) /
+               (parambounds$gamma[2] - parambounds$gamma[1] + 1e-03) ))
R>
R> get_alpha <- function(theta)
+   exp(theta) + parambounds$alpha[1]
R> # Add a little buffer, for stability
R> get_theta3 <- function(alpha)
+   log(alpha - parambounds$alpha[1] + 1e-03)
R>
R> get_beta <- function(theta)
+   (parambounds$beta[2] - parambounds$beta[1]) *
+     exp(-exp(theta)) + parambounds$beta[1]
R> get_theta4 <- function(beta)
+   log(-log( (beta - parambounds$beta[1]) /
+               (parambounds$beta[2] - parambounds$beta[1]) ))
R> ## Compute the transformed pdfs ##
R> translist1 <- list(totheta = get_theta1,fromtheta = get_psi0)
R> translist2 <- list(totheta = get_theta2,fromtheta = get_gamma)
R> translist3 <- list(totheta = get_theta3,fromtheta = get_alpha)
R> translist4 <- list(totheta = get_theta4,fromtheta = get_beta)
R>
R> psi0pdf <- compute_pdf_and_cdf(astroquad$marginals[[1]],translist1)
R> gammapdf <- compute_pdf_and_cdf(astroquad$marginals[[2]],translist2)
R> alphapdf <- compute_pdf_and_cdf(astroquad$marginals[[3]],translist3)
R> betapdf <- compute_pdf_and_cdf(astroquad$marginals[[4]],translist4)
R>
R> ## Plots
R> psi0postplot <- psi0pdf
+   ggplot(aes(x = transparam,y = pdf_transparam)) +
+   theme_classic() +
+   geom_line() +
+   labs(x = expression(Psi),y = "Density")
R> gammapostplot <- gammapdf
+   ggplot(aes(x = transparam,y = pdf_transparam)) +
+   theme_classic() +
+   geom_line() +
+   labs(x = expression(gamma),y = "Density")
R> alphapostplot <- alphapdf
+   ggplot(aes(x = transparam,y = pdf_transparam)) +
+   theme_classic() +
+   geom_line() +
+   labs(x = expression(alpha),y = "Density")
R> betapostplot <- betapdf
+   ggplot(aes(x = transparam,y = pdf_transparam)) +
+   theme_classic() +
+   geom_line() +
+   labs(x = expression(beta),y = "Density")
R>
R> (psi0postplot | gammapostplot) / (alphapostplot | betapostplot)
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{01-aghq-softwarepaper_files/figure-latex/astrogamma1-1}

}

\caption[Marginal posteriors for original parameters in mass estimation example]{Marginal posteriors for original parameters in mass estimation example.}\label{fig:astrogamma1}
\end{figure}
\end{CodeChunk}

Comparing Figures \ref{fig:astrogamma1} and \ref{fig:astroquad11},
especially the marginal posteriors for \(\gamma\) and \(\alpha\), helps
to illustrate why it is helpful to do quadrature on a transformed scale.

The primary object of interest for inference in this example is
\(M_{\boldsymbol{\Xi}}(r)\), the mass of the Galaxy, as defined above.
It is a function of two variables. We can compute moments of it for any
\(r\) as follows:

\begin{CodeChunk}
\begin{CodeInput}
R> Mr <- function(r,theta) {
+   p = get_psi0(theta[1])
+   g = get_gamma(theta[2])
+
+   # Manual unit conversion into "mass of one trillion suns" (so awesome)
+   g*p*r^(1-g) * 2.325e09 * 1e-12
+ }
R>
R> rtodo <- 1:150
R> Mrout <- numeric(length(rtodo))
R> Mrsdout <- numeric(length(rtodo))
R> for (rr in 1:length(rtodo)) {
+   r <- rtodo[rr]
+
+   Mrout[rr] <- compute_moment(
+     astroquad$normalized_posterior,
+     function(x) Mr(r,x)
+   )
+   Mrsdout[rr] <- sqrt(compute_moment(
+     astroquad$normalized_posterior,
+     function(x) (Mr(r,x) - Mrout[rr])^2
+   ))
+ }
R> tibble(
+   r = rtodo,
+   Mr = Mrout,
+   Mrsd = Mrsdout
+ )
+   ggplot(aes(x = r)) +
+   theme_classic() +
+   geom_line(aes(y = Mr),size = 1) +
+   geom_ribbon(aes(ymin = Mrout - Mrsdout,
+                   ymax = Mrout + Mrsdout),
+               alpha = .5) +
+   geom_ribbon(aes(ymin = Mrout - 2*Mrsdout,
+                   ymax = Mrout + 2*Mrsdout),
+               alpha = .2) +
+   labs(title = "",
+        x = "r (kpc)",
+        y = bquote('M(r) ('~10^12~M[sun]~')')) +
+   scale_x_continuous(breaks = seq(0,150,by = 25)) +
+   scale_y_continuous(breaks = seq(0,1,by=.1))
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{01-aghq-softwarepaper_files/figure-latex/computemass1-1}

}

\caption[Mass estimate (---) and one ($\textcolor{gray}{\blacksquare}$) and two ($\textcolor{lightgray}{\blacksquare}$) standard deviation bands]{Mass estimate (---) and one ($\textcolor{gray}{\blacksquare}$) and two ($\textcolor{lightgray}{\blacksquare}$) standard deviation bands.}\label{fig:computemass1}
\end{figure}
\end{CodeChunk}

These results match closely to those reported by \citet{gwen2},
accounting for the priors and subset of data we used here.

\hypertarget{sec:elgms}{%
\section{Extended example: implementing Extended Latent Gaussian
Models}\label{sec:elgms}}

While I have demonstrated that there are some difficult models in which
\texttt{AGHQ}{} can provide the sole framework for Bayesian inference,
the real utility of the \pkg{aghq} package is its use within a more
complicated framework for approximate Bayesian inference in certain
high-dimensional models. \citet{aghqus} and \citet{noeps} describe methods which use \AGHQ{} in this manner with several supporting examples. I describe the computational details
of two such examples here, and how the \pkg{aghq} package is used to fit them. I first
implement a geostatistical binomial regression model---which is already
implemented in several existing packages---using the \pkg{aghq}
package. I then show how the code for this model is easily extended to
the more complicated case of a \emph{zero-inflated} geostatistical
binomial regression model defined by \citet{geostatlowresource}, which is not easily fit using existing
packages, but was fit by \citet{aghqus} using \pkg{aghq}.

\hypertarget{sec:loaloa}{%
\subsection{Geostatistical Binomial regression: the Loa Loa parasitic
roundworm in Nigeria and Cameroon}\label{sec:loaloa}}

\citet{geostatsp} considers counts of subjects who tested positive for a tropical
disease caused by the \emph{Loa loa} parasitic roundworm in \(n = 190\)
villages in Cameroon and Nigeria. We wish to model the association of
infection probability with elevation, vegitation index, and land type of
each village, as well as quantify excess spatial variation in it. These
data are available as the \texttt{loaloa} data in the \pkg{geostatsp}
package \citep{geostatsp}, and originally analyzed using \texttt{MCMC}{}
\citep{loaloa}. Consider the following geostatistical binomial
regression model \citep{geostatsp}: \begin{equation}\begin{aligned}
Y_{i} | p_{i} &\overset{ind}{\sim}\text{Binomial}(N_{i},p_{i}), \\
\log\left(\frac{p_{i}}{1-p_{i}}\right) &= \boldsymbol{x}_{i}^{T}\boldsymbol{\beta} + u(\boldsymbol{s}_{i}), \\
u(\cdot) | \boldsymbol{\theta} &\sim \mathcal{GP}(0,\text{C}_{\boldsymbol{\Xi}}), \\
\end{aligned}\end{equation} where \(Y_{i},N_{i}\) are the counts of
diseased and total individuals, \(\boldsymbol{x}_{i}\) is a vector
containing the covariates, \(\boldsymbol{s}_{i}\) is the spatial
location of village \(i\in[n]\), and
\(p_{i} \equiv p(\boldsymbol{s}_{i})\) is the infection probability at
location \(\boldsymbol{s}_{i}\). The unknown function \(u(\cdot)\)
governs excess spatial variation in logit infection risk, and is
modelled as a Gaussian Process \citep{gaussianprocesses} with Matern
covariance function \(\text{C}_{\boldsymbol{\Xi}}\) \citep{spde}. This
depends on parameters \(\boldsymbol{\Xi} = (\sigma,\rho)\) where
\(\sigma\) is the marginal standard deviation and \(\rho\) the range in
the parameterization of \citet{geostatsp}. Quadrature is done on the
transformed scale using the transformations introduced by
\citet{pcpriormatern}: \begin{equation}\begin{aligned}
\boldsymbol{\theta}&= (\log\kappa,\log\tau), \\
\kappa &= \sqrt{8\nu}/\rho, \\
\tau &= \sigma\kappa^{\nu} \sqrt{\frac{\Gamma\left(\nu + d/2\right)(4\pi)^{d/2}}{\Gamma(\nu)}} \\
\end{aligned}\end{equation} where \(\nu = d = 2\) in this example.
Priors are chosen according to \citet{pcpriormatern} and satisfy
\(P(\sigma < 4) = P(\rho < 50\text{km}) = 97.5\%\).

Define
\(\boldsymbol{U} = \left\{ u(\boldsymbol{s}_{i}):i\in[n] \right\}\) and
\(\boldsymbol{W} = (\boldsymbol{U},\boldsymbol{\beta})\). The posterior
distributions of interest are: \begin{equation}\begin{aligned}
\pi(\boldsymbol{W}|\boldsymbol{Y}) &= \int\pi(\boldsymbol{W}|\boldsymbol{\theta},\boldsymbol{Y})\pi(\boldsymbol{\theta}|\boldsymbol{Y})d\boldsymbol{\theta}, \\
\pi(\boldsymbol{\theta}|\boldsymbol{Y}) &= \frac{\int\pi(\boldsymbol{W},\boldsymbol{\theta},\boldsymbol{Y})d\boldsymbol{W}}{\int\int\pi(\boldsymbol{W},\boldsymbol{\theta},\boldsymbol{Y})d\boldsymbol{W}d\boldsymbol{\theta}}.
\end{aligned}\end{equation} The method for inference using \AGHQ{} described by \citet{aghqus} and \citet{noeps}
is to:

\begin{enumerate}
\item Approximate $\int\pi(\boldsymbol{W},\boldsymbol{\theta},\boldsymbol{Y})d\boldsymbol{W} \approx \widetilde{\pi}^{*}_{\texttt{LA}}(\boldsymbol{\theta},\boldsymbol{Y})$ using \texttt{AGHQ}{} with $k=1$
(a Laplace approximation), saving the mode $\widehat{\boldsymbol{W}}_{\boldsymbol{\theta}} = \text{argmax}_{\boldsymbol{W}}\log\pi(\boldsymbol{W},\boldsymbol{\theta},\boldsymbol{Y})$
and Hessian $\boldsymbol{H}(\boldsymbol{\theta}) = -\partial^{2}_{\boldsymbol{W}}\log\pi(\widehat{\boldsymbol{W}}_{\boldsymbol{\theta}},\boldsymbol{\theta},\boldsymbol{Y})$ each time this
approximation is evaluated,
\item Normalize $\widetilde{\pi}^{*}_{\texttt{LA}}(\boldsymbol{\theta},\boldsymbol{Y})$ using \texttt{AGHQ}{}
with $k > 1$, obtaining $\widetilde{\pi}_{\texttt{AGHQ}{}}(\boldsymbol{\theta}|\boldsymbol{Y})$,
\item Approximate:
\begin{equation}\begin{aligned}
\pi(\boldsymbol{W}|\boldsymbol{Y}) \approx \widetilde{\pi}(\mb{W}|\mb{Y}) &\approx
\int\widetilde{\pi}_{\texttt{G}}(\boldsymbol{W}|\boldsymbol{\theta},\boldsymbol{Y})
\widetilde{\pi}_{\texttt{AGHQ}{}}(\boldsymbol{\theta}|\boldsymbol{Y})d\mb{\theta} \\
&\approx \sum_{j=1}^{k}\widetilde{\pi}_{\texttt{G}}(\boldsymbol{W}|\boldsymbol{\theta}^{j},\boldsymbol{Y})
\widetilde{\pi}_{\texttt{AGHQ}{}}(\boldsymbol{\theta}^{j}|\boldsymbol{Y})\omega_{j}
\end{aligned}\end{equation}
where $\widetilde{\pi}_{\texttt{G}}(\boldsymbol{W}|\boldsymbol{\theta},\boldsymbol{Y})$ is a Gaussian approximation
to $\pi(\boldsymbol{W}|\boldsymbol{\theta},\boldsymbol{Y})$ for fixed $\boldsymbol{\theta}$ having mean $\widehat{\boldsymbol{W}}_{\boldsymbol{\theta}}$
and precision matrix $\boldsymbol{H}(\boldsymbol{\theta})$. This integral has the form of a moment with respect to $\widetilde{\pi}_{\texttt{AGHQ}{}}(\boldsymbol{\theta}|\boldsymbol{Y})$ and is further approximated using the same \AGHQ{} points and weights, in the manner of \S\ref{sec:basic-use} and \S\ref{sec:background}.
\end{enumerate}

With these approximations available, inference any set of points in the
study region \(\boldsymbol{U}^{*}\) is approached using the posterior predictive
distribution
\(\pi(\boldsymbol{U}^{*}|\boldsymbol{Y})\).
A sample is drawn from this distribution as follows:

\begin{enumerate}
\item Draw $t\sim\text{Multinomial}(\lambda_{1},\ldots,\lambda_{k})$ where $\lambda_{j} = \widetilde{\pi}_{\texttt{AGHQ}{}}(\boldsymbol{\theta}^{j}|\boldsymbol{Y})\omega_{j}, j\in[k]$,
\item Draw $\boldsymbol{W} = (\boldsymbol{U},\boldsymbol{\beta})\sim\widetilde{\pi}_{\texttt{G}}(\boldsymbol{W}|\boldsymbol{\theta}^{t},\boldsymbol{Y})$, and retain $\boldsymbol{U}$,
\item Draw $\boldsymbol{U}^{*}$ from $\pi(\boldsymbol{U}^{*}|\boldsymbol{U})$,
\end{enumerate}

which makes further use of the \texttt{AGHQ}{} approximation. Steps 1.
and 2. together correspond to drawing a sample from
\(\widetilde{\pi}(\boldsymbol{W}|\boldsymbol{Y})\) and are implemented
in \texttt{aghq::sample\_marginal}. Step 3 is specific to this problem
and unrelated to \texttt{AGHQ}{} as a method, and I use methods from the
\texttt{RandomFields} package \citep{randomfields} for this. A similar
approach is taken to sample from the posterior of the spatial risk.

The \pkg{aghq} package is used to fit this model. A function computing
\(\log\pi(\boldsymbol{W},\boldsymbol{\theta},\boldsymbol{Y})\) is
defined:

\begin{CodeChunk}
\begin{CodeInput}
R> ## Log posterior ##
R> # See appendix for full code
R> log_posterior_joint <- function(W,theta) {
+   Q <- Q_matrix(theta) # Prior precision matrix
+   log_prior_theta(theta) + log_prior_W(W,theta,Q) + log_likelihood(W)
+ }
R>
R> ff <- list(
+   fn = log_posterior_joint,
+   gr = grad_log_posterior_W,
+   he = function(W,theta) -1 * H_matrix(W,theta)
+ )
R>
R> # Starting values taken from Brown (2011)
R> startingsig <- .988
R> startingrho <- 4.22*1e04
R> startingtheta <- c(
+   log(get_kappa(startingsig,startingrho)),
+   log(get_tau(startingsig,startingrho))
+ )
\end{CodeInput}
\end{CodeChunk}

This is then passed to the \texttt{aghq::marginal\_laplace} function:

\begin{CodeChunk}
\begin{CodeInput}
R> loaloaquad <- aghq::marginal_laplace(ff,3,list(W = rep(0,196),theta = startingtheta))
\end{CodeInput}
\end{CodeChunk}

The \texttt{marginal\_laplace} function returns an object of class
\texttt{aghq} as usual, with one additional list item
\texttt{modesandhessians} containing the mode
\(\widehat{\boldsymbol{W}}_{\boldsymbol{\theta}}\) and Hessian
\(\boldsymbol{H}(\boldsymbol{\theta})\) for each quadrature point
\(\boldsymbol{\theta}\). Inference for \(\boldsymbol{\theta}\) can be
done using the usual \texttt{summary} and \texttt{plot} methods, however
interest in this example is for \(\boldsymbol{W}\). The \texttt{aghq}
package contains a function \texttt{sample\_marginal} which takes the
output of \texttt{marginal\_laplace} and efficienty draws a specified
number of samples from
\(\widetilde{\pi}(\boldsymbol{W}|\boldsymbol{Y})\):

\begin{CodeChunk}
\begin{CodeInput}
R> ## Sample the full spatial field ##
R> # Sampling is very fast, but the next step with the maps is very slow
R> # so only do 10 here.
R> fieldsamps <- sample_marginal(loaloaquad,10)
R> # fieldsamps$samps[1:190, ] = W, fieldsamps$samps[191:196, ] = beta
R> UU <- fieldsamps$samps[1:190, ] # Excess spatial variation in logit-prob of infection
R> EE <- apply(fieldsamps$samps,2,make_eta) # logit prob of infection
\end{CodeInput}
\end{CodeChunk}

Each column of \texttt{fieldsamps} is one sample from the posterior of
\(\boldsymbol{W} = (\boldsymbol{U},\boldsymbol{\beta})\), of dimension
\(196 = 190 + 6\).

The final step required to draw the maps is unique to this specific
problem, and we use a custom helper function to draw
\(\boldsymbol{U}^{*}\) using \texttt{R}'s spatial data structures.
Appendix \ref{app:elgm} shows full code for defining all the helper
functions.

\begin{CodeChunk}
\begin{CodeInput}
R> # Get the samples for U
R> # Problem-specific helper function, Appendix A.3
R> fieldbrickU <- simulate_spatial_fields(
+     UU,
+     fieldsamps$theta,
+     loaloa,
+     resolution = list(nrow = 200,ncol = 400)
+   )
R> # Now for eta
R> fieldbricketa <- simulate_spatial_fields(
+     EE,
+     fieldsamps$theta,
+     loaloa,
+     resolution = list(nrow = 200,ncol = 400)
+   )
\end{CodeInput}
\end{CodeChunk}

With posterior samples in hand, we may estimate any quantity of
interest:

\begin{CodeChunk}
\begin{CodeInput}
R> # Average over the samples
R> simfieldsmeanU <- mean(fieldbrickU)
R> simfieldsmeanlambda <- calc(fieldbricketa,function(x) mean(1 / (1 + exp(-x))))
\end{CodeInput}
\end{CodeChunk}

Figure \ref{fig:makemaps1} shows the resulting maps of
\texttt{simfieldsmeanU} and \texttt{simfieldsmeanlambda}.

\begin{CodeChunk}
\begin{figure}

{\centering \includegraphics{01-aghq-softwarepaper_files/figure-latex/makemaps1-1}

}

\caption[Posterior mean for $u(\cdot)$ (top) and $p(\cdot)$ (bottom) for loaloa example]{Posterior mean for $u(\cdot)$ (top) and $p(\cdot)$ (bottom) for loaloa example.}\label{fig:makemaps1}
\end{figure}
\end{CodeChunk}

\hypertarget{sec:loaloazip}{%
\subsection{Zero-inflated Geostatistical Binomial regression: extending
the Loa loa model}\label{sec:loaloazip}}

\citet{geostatlowresource} describe an extension of the geostatistical
binomial regression model which accounts for the notion that certain
spatial locations may be unsuitable for disease transmission, and fit
this model using frequentist inference. \citet{aghqus} fit this model using approximate
Bayesian inference with \AGHQ{} in a similar manner to that described in \S\ref{sec:loaloa}. The \emph{zero-inflated} geostatistical binomial regression model
is: \begin{equation}\begin{aligned}
\mathbb{P}\left(Y_{i} = y_{i} | p_{i},\phi_{i}\right) &= \left(1 - \phi_{i}\right)\text{I}\left(y_{i} = 0\right) + \phi_{i}\text{Binomial}(y_{i};N_{i},p_{i}), \\
\log\left(\frac{p_{i}}{1-p_{i}}\right) &= \boldsymbol{x}_{i}^{T}\boldsymbol{\beta}_{\texttt{risk}} + u(\boldsymbol{s}_{i}), \\
\log\left(\frac{\phi_{i}}{1-\phi_{i}}\right) &= \boldsymbol{x}_{i}^{T}\boldsymbol{\beta}_{\texttt{zi}} + v(\boldsymbol{s}_{i}), \\
u(\cdot) | \boldsymbol{\theta} &\sim \mathcal{GP}(0,\text{C}_{\boldsymbol{\Xi},u}), \\
v(\cdot) | \boldsymbol{\theta} &\sim \mathcal{GP}(0,\text{C}_{\boldsymbol{\Xi},v}), \\
\end{aligned}\end{equation} where
\(\phi_{i} \equiv \phi(\boldsymbol{s}_{i})\) is the probability that
location \(\boldsymbol{s}_{i}\) is suitable for disease transmission,
and \(u(\cdot)\) and \(v(\cdot)\) are independent spatial processes.
\citet{loaloazero} consider a related model for these data specifically.
Inference for this model is very challenging, and not easily done using
any existing packages. For illustration we set the standard deviation
and range parameters \((\sigma,\rho)\) for these two processes to be
equal, and ignore covariates, setting
\(\boldsymbol{\beta}_{\texttt{risk}} = \boldsymbol{\beta}_{\texttt{zi}} = 0\).

We modify the helper functions from \S\ref{sec:loaloa} to compute the
new log-likelihood, again with the help of the \pkg{TMB} package for
computing derivatives. With the likelihood changed, the code for fitting
the model, sampling from the posteriors, and drawing the maps of the
suitability \(\phi(\boldsymbol{s})\) and infection
\(\phi(\boldsymbol{s})p(\boldsymbol{s})\) probabilities is nearly
identical to the simpler model of \S\ref{sec:loaloa}. The only practical
differences owe to the difficulty of the model, and are unrelated to the
\texttt{AGHQ}{} methodology and software: there are few zeroes in these
data, making it challenging to identify \(\phi(\boldsymbol{s})\); and
the posterior for \(\boldsymbol{\Xi}\) is very flat, making it very
challenging to infer separate parameters for the two spatial processes.
Even when fixing them to be equal, the surface is so flat that the
default \texttt{BFGS} optimization routine for
\(\widetilde{\pi}_{\texttt{LA}}(\boldsymbol{\theta}|\boldsymbol{Y})\)
diverges, and we switch to the slower but highly stable \texttt{trust}
trust optimization routine.

\begin{CodeChunk}
\begin{CodeInput}
R> ## Quadrature ##
R>
R> log_posterior_joint <- function(W,theta) {
+   log_prior_theta(theta) +  log_posterior_W(W,theta)
+ }
R>
R> ff <- list(
+   fn = log_posterior_joint,
+   gr = grad_log_posterior_W,
+   he = function(W,theta) -1 * H_matrix(W,theta)
+ )
R>
R> startingsig <- .988
R> startingrho <- 4.22*1e04
R> startingtheta <- c(
+   log(get_kappa(startingsig,startingrho)),
+   log(get_tau(startingsig,startingrho))
+ )
R> # Get good starting values for W by optimizing once for the starting theta,
R> # which speeds up evaluation of log pi(theta|Y)
R> ffW <- list(
+   fn = function(W) ff$fn(W,startingtheta),
+   gr = function(W) ff$gr(W,startingtheta),
+   he = function(W) ff$he(W,startingtheta)
+ )
R> tmp <- capture.output({
+   Wopt <- optimize_theta(ffW,
+                          rep(0,model_data$Wd),
+                          control = list(method = 'sparse_trust'))
+ })
\end{CodeInput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> loaloazipquad <- aghq::marginal_laplace(
+   ff,
+   3,
+   list(W = Wopt$mode,theta = startingtheta),
+   control = list(method = 'trust',
+                  inner_method = 'sparse_trust')
+ )
\end{CodeInput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> # Post samples
R> loazippostsamples <- sample_marginal(loaloazipquad,10)
R> nn <- nrow(loazippostsamples$samps)
R> UUzip <- loazippostsamples$samps[1:(nn/2), ]
R> UUrisk <- loazippostsamples$samps[(1+nn/2):nn, ]
\end{CodeInput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> fieldbrickzip <- simulate_spatial_fields(
+   UUzip,
+   loazippostsamples$theta,
+   loaloa,
+   resolution = list(nrow = 200,ncol = 400)
+ )
R>
R> fieldbrickrisk <- simulate_spatial_fields(
+   UUrisk,
+   loazippostsamples$theta,
+   loaloa,
+   resolution = list(nrow = 200,ncol = 400)
+ )
\end{CodeInput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> simfieldsmeanzip <- calc(fieldbrickzip,function(x) mean(1 / (1 + exp(-x))))
R> simfieldsmeanrisk <- calc(
    (1/(1+exp(-fieldbrickzip))) *
    (1 /(1+exp(-fieldbrickrisk))),
  mean
)
\end{CodeInput}
\end{CodeChunk}

Figure \ref{fig:loaloazipmaps1} shows the resulting maps for the
suitability probabilities \(\phi(\cdot)\) as well as the predicted
incidence probabilities \(\phi(\cdot)\times p(\cdot)\).

\begin{CodeChunk}
\begin{figure}

{\centering \includegraphics{01-aghq-softwarepaper_files/figure-latex/loaloazipmaps1-1}

}

\caption[Posterior mean for $\phi(\cdot)$ (top) and $\phi(\cdot)\times p(\cdot)$ (bottom) for loaloa zero-inflated example]{Posterior mean for $\phi(\cdot)$ (top) and $\phi(\cdot)\times p(\cdot)$ (bottom) for loaloa zero-inflated example in the manner of \citet{aghqus}.}\label{fig:loaloazipmaps1}
\end{figure}
\end{CodeChunk}

\hypertarget{sec:glmmTMB}{%
\section{Bayesian Inference for Existing Frequentist
Models}\label{sec:glmmTMB}}

The \pkg{aghq} package can also be used to extend existing frequentist
methods to be Bayesian, allowing the users to specify priors over and
quantify uncertainty in parameters which are otherwise difficult to
infer. We consider an examples of a zero-inflated, overdispersed Poisson
regression with random effects fit using the \pkg{glmmTMB} package, and
illustrate how Bayesian inferences can be made for this model using
either \texttt{aghq} or \texttt{marginal\_laplace}.

\hypertarget{zero-inflated-overdispersed-poisson-regression}{%
\subsection{Zero-inflated, overdispersed Poisson
regression}\label{zero-inflated-overdispersed-poisson-regression}}

Consider the following example from the \pkg{glmmTMB} package
\citep{glmmTMB}. Counts \(Y_{i}, i = 1,\ldots,n\) of salamanders in 23
streams were measured 4 times. Some streams are affected by mining, and
this is thought to affect the mean number of salamanders observed.
Further, the mining status of a site may render the stream inhabitable
which may result in many \emph{structural zero} counts. Some zero counts
would be expected under a standard Poisson model, however
(\emph{non-structural zeros}), and the probability of an observed count
being a structural zero is modeled to depend on whether the site was
mined. Finally, the variance in the counts is thought to be higher than
expected under the Poisson model, and this \emph{overdispersion} is
allowed to depend on the day of the year. The full hierarchical model
is: \begin{equation}\begin{aligned}
Y_{i} | \lambda_{i},d_{i}=0,\sigma^{2}_{i} &\overset{ind}{\sim} \text{quasi-Poisson}(\lambda_{i},\sigma^{2}_{i}), \\
\log\lambda_{i} &= \beta_{0}^{\texttt{(mean)}} + \beta_{\texttt{mined}}\cdot x_{\texttt{mined},i} + u_{i}, \\
u_{i} &\overset{iid}{\sim}\text{Normal}(0,\sigma^{2}_{u}), \\
\sigma^{2}_{i} &= \lambda_{i}\left(1 + \frac{\lambda_{i}}{\theta_{i}}\right), \\
\log\left(\frac{p_{i}}{1-p_{i}}\right) &= \beta_{0}^{\texttt{(zi)}} + \beta_{\texttt{mined}}^{\texttt{(zi)}}\cdot x_{\texttt{mined},i}, \\
\log\theta_{i} &= \beta_{0}^{\texttt{(disp)}} + \beta_{\texttt{DOY}}\cdot x_{\texttt{DOY}}
\end{aligned}\end{equation} where
\(\lambda_{i} = \mathbb{E}(Y_{i}|u_{i},d_{i} = 0)\), \(d_{i} = 1\) if
\(Y_{i}\) is a structural zero such that
\(\mathbb{P}(Y_{i} = 0|d_{i} = 1) = 1\),
\(p_{i} = \mathbb{P}(d_{i} = 1)\),
\(\sigma^{2}_{i} = \text{Var}(Y_{i}|\lambda_{i},\theta_{i})\),
\(\theta_{i}\) is the quasi-Poisson/negative Binomial overdispersion
parameter, \(x_{\texttt{mined},i} = 1\) if site \(i\) was mined and
\(0\) otherwise, \(u_{i}\) is a random intercept for site \(i\) with
fixed, unknown variance \(\sigma^{2}_{u}\), and \(x_{\texttt{DOY}}\) is
the (centred and scaled) day of the year on which measurement \(i\) was
obtained.

Frquentist inference for this model is done using the \pkg{glmmTMB}
package. A Laplace approximation to the likelihood
\(L(\beta,\sigma^{2}_{u});\boldsymbol{Y} \equiv \int L(\boldsymbol{u},\beta,\sigma^{2}_{u};\boldsymbol{Y})d\boldsymbol{u}\)
is maximized yielding point estimates for \(\beta\) and
\(\sigma^{2}_{u}\). Uncertainty for \(\beta\) is approximated using the
inverse curvature of the log-Laplace approximated likelihood at the
maximum, and uncertainty for \(\sigma^{2}_{u}\) is not reported or
accounted for. The model may be easily fit using this package as
follows:

\begin{CodeChunk}
\begin{CodeInput}
R> library(glmmTMB)
R> ## ZIP example from glmmTMB ##
R> data("Salamanders",package = "glmmTMB")
R> zipmod <- glmmTMB(
+   count ~ mined + (1|site),
+   zi = ~mined,
+   disp = ~DOY,
+   data = Salamanders,
+   family = nbinom2
+ )
\end{CodeInput}
\begin{CodeInput}
R> summary(zipmod)
\end{CodeInput}
\begin{CodeOutput}
 Family: nbinom2  ( log )
Formula:          count ~ mined + (1 | site)
Zero inflation:         ~mined
Dispersion:             ~DOY
Data: Salamanders

     AIC      BIC   logLik deviance df.resid
    1735     1767     -861     1721      637

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 site   (Intercept) 0.134    0.366
Number of obs: 644, groups:  site, 23

Conditional model:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.540      0.376   -1.44     0.15
minedno        1.424      0.365    3.90 0.000098 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)    0.256      0.487    0.53   0.5988
minedno       -2.244      0.745   -3.01   0.0026 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dispersion model:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.0278     0.3177    0.09     0.93
DOY          -0.3947     0.1537   -2.57     0.01 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{CodeOutput}
\end{CodeChunk}

\hypertarget{bayesian-inference-adaptive-quadrature}{%
\subsection{Bayesian inference: adaptive
quadrature}\label{bayesian-inference-adaptive-quadrature}}

We extend inference this model to be fully Bayesian using the \pkg{aghq}
package, allowing in particular to assess uncertainty in
\(\sigma^{2}_{u}\) and account for this uncertainty when reporting
intervals for \(\beta\). We put independent weakly-informative
\(\text{Normal}(0,1000)\) priors for each component of \(\beta\) as well
as \(\sigma_{u} \sim \text{Exponential}(\phi)\) with \(\phi\) chosen so
that \(\mathbb{P}(\sigma_{u} > 1) = .5\), a penalized complexity prior
\citep{pcprior}. The \texttt{zipmod} object returned by
\texttt{glmmTMB::glmmTMB} contains a function which computes the
negative laplace-approximated log-likelihood and its gradient, using the
automatic differentation interface of the \pkg{TMB} package. We add the
priors to this and then pass the result to \pkg{aghq}:

\begin{CodeChunk}
\begin{CodeInput}
R> logprior <- function(x) {
+   sum(dnorm(x[1:6],0,log(1e03))) + dexp(exp(x[7]),log(2)) + x[7]
+ }
R>
R> logpost <- function(x) {
+   # x = beta1 ... beta6, log(sigma)
+   -1 * zipmod$obj$fn(x) + logprior(x)
+ }
R> logpostgr <- function(x) {
+   -1 * zipmod$obj$gr(x) + numDeriv::grad(logprior,x)
+ }
R> # TMB does not give AD hessians for laplace-approximate likelihoods
R> logposthe <- function(x) {
+   -1 * numDeriv::jacobian(zipmod$obj$gr,x) + numDeriv::hessian(logprior,x)
+ }
R>
R> ff <- list(
+   fn = function(x) logpost(x),
+   gr = function(x) logpostgr(x),
+   he = function(x) logposthe(x)
+ )
R>
R> # zipmod$fit$par contains the MLE; use as a starting value
R> zipmodquad <- aghq(ff,3,zipmod$fit$par)
\end{CodeInput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> summary(zipmodquad)
\end{CodeInput}
\begin{CodeOutput}
AGHQ on a 7 dimensional posterior with  3 3 3 3 3 3 3 quadrature points

The posterior mode is: -0.616 1.49 0.173 -2.22 0.00758 -0.393 -0.874

The log of the normalizing constant/marginal likelihood is: -864

The posterior Hessian at the mode is:
       [,1]   [,2]   [,3]    [,4]   [,5]   [,6]   [,7]
[1,]  75.18 51.117 -17.25  -5.043  -2.66 -2.603  6.637
[2,]  51.12 51.118  -5.04  -5.043  -2.15 -3.213  0.581
[3,] -17.25 -5.043  19.61   5.957 -19.06 -1.991  1.277
[4,]  -5.04 -5.043   5.96   5.958 -13.70 -0.767  0.379
[5,]  -2.66 -2.145 -19.06 -13.697  51.73 13.079 -2.510
[6,]  -2.60 -3.213  -1.99  -0.767  13.08 50.999  0.243
[7,]   6.64  0.581   1.28   0.379  -2.51  0.243 11.323

The covariance matrix used for the quadrature is...
        [,1]     [,2]    [,3]    [,4]    [,5]     [,6]     [,7]
[1,]  0.1380 -0.12194  0.1482 -0.0118  0.0516 -0.00790 -0.07932
[2,] -0.1219  0.13555 -0.1223  0.0772 -0.0229  0.00424  0.07057
[3,]  0.1482 -0.12230  0.2441 -0.0208  0.0859 -0.01253 -0.08811
[4,] -0.0118  0.07716 -0.0208  0.6324  0.1714 -0.03111  0.02282
[5,]  0.0516 -0.02292  0.0859  0.1714  0.1018 -0.01888 -0.02155
[6,] -0.0079  0.00424 -0.0125 -0.0311 -0.0189  0.02335  0.00218
[7,] -0.0793  0.07057 -0.0881  0.0228 -0.0215  0.00218  0.13555

...and its Cholesky is:
      [,1]   [,2]   [,3]    [,4]   [,5]    [,6]     [,7]
[1,] 0.371 -0.328 0.3989 -0.0319 0.1390 -0.0213 -0.21355
[2,] 0.000  0.167 0.0519  0.4002 0.1363 -0.0164  0.00277
[3,] 0.000  0.000 0.2870 -0.1005 0.0814 -0.0111 -0.01073
[4,] 0.000  0.000 0.0000  0.6791 0.1906 -0.0388  0.02036
[5,] 0.000  0.000 0.0000  0.0000 0.1447 -0.0373  0.03286
[6,] 0.000  0.000 0.0000  0.0000 0.0000  0.1400 -0.00300
[7,] 0.000  0.000 0.0000  0.0000 0.0000  0.0000  0.29718

Here are some moments and quantiles for theta:

          mean  median     mode    sd   2.
theta1 -0.7364 -0.7380 -0.61569 0.361 -1.441 -0.0297
theta2  1.5808  1.5760  1.49046 0.376  0.853  2.3011
theta3 -0.0330 -0.0587  0.17257 0.483 -0.922  0.7979
theta4 -2.5326 -2.6041 -2.22076 0.767 -4.093 -1.0984
theta5 -0.0679  0.0197  0.00758 0.238 -0.482  0.5861
theta6 -0.3903 -0.3942 -0.39264 0.144 -0.676 -0.1102
theta7 -0.7811 -0.7939 -0.87430 0.357 -1.537 -0.0483
\end{CodeOutput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> # Plot the posterior of sigma_u
R> # and compare confidence intervals for beta
R> sigmaplot <- compute_pdf_and_cdf(
+   zipmodquad$marginals[[7]],
+   list(
+     totheta = log,
+     fromtheta = exp
+   )
+ )
+   ggplot(aes(x = transparam,y = pdf_transparam)) +
+   theme_classic() +
+   geom_line() +
+   geom_vline(xintercept = as.numeric(sqrt(summary(zipmod)$varcor$cond$site)),
+              linetype = "dotdash") +
+   stat_function(fun = dexp,args = list(rate = log(2)),linetype = "dashed") +
+   labs(x = expression(sigma),y = "Density")
R> varnames <- c('beta0','beta1','beta0zi','beta1zi','beta0disp','beta1disp')
R> quadmodbetaconf <- summary(zipmodquad)$summarytable[1:6,c('2.
R> quadmodbetaconf$var <- varnames
R> quadmodbetaconf$Method <- 'AGHQ'
R> quadmodbetaconf <- as_tibble(quadmodbetaconf)
+   rename(lower = `2.
R> zipmodbetaconf <- summary(zipmod)$coefficients
+   reduce(rbind)
+   as_tibble()
+   dplyr::select(estimate = Estimate,
+          sd = `Std. Error`)
+   mutate(lower = estimate - 2*sd,
+          upper = estimate + 2*sd,
+          var = varnames,
+          Method = 'glmmTMB')
+   dplyr::select(-sd)
R>
R> betaplot <- bind_rows(quadmodbetaconf,zipmodbetaconf)
+   ggplot(aes(x = var,colour = Method)) +
+   theme_bw() +
+   geom_point(aes(y = estimate)) +
+   geom_errorbar(aes(ymin = lower,ymax = upper),width=.1) +
+   labs(colour = "Method",x = "Coefficient",y = "Estimate and 9
R>
R> sigmaplot / betaplot
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{01-aghq-softwarepaper_files/figure-latex/zeromod3aghq-1}

}

\caption[(Top) Posterior (---) and prior (- - -) for random effects standard deviation computed with AGHQ, with point estimate (- $\cdot$ -) from glmmTMB]{(Top) Posterior (---) and prior (- - -) for random effects standard deviation computed with AGHQ, with point estimate (- $\cdot$ -) from glmmTMB. (Bottom) Comparison of point and interval estimates for the regression coefficients estimated using glmmTMB and AGHQ.}\label{fig:zeromod3aghq}
\end{figure}
\end{CodeChunk}

Accounting for the uncertainty in \(\sigma_{u}\) has little effect on
the confidence (credible) intervals for \(\beta\), however it has a
marked effect on the inference for \(\sigma_{u}\) which we see is
somewhat imprecisely estimated. This information is not available from
the frequentist model used by \pkg{glmmTMB}.

We note that inference in this example is more similar to what is done
using \texttt{aghq::marginal\_laplace}, because the ``likelihood'' we
pass to \texttt{aghq::aghq} is already computed using a Laplace
approximation. Note also that this technique could, in principle, be
applied this easily to any model fit using the \pkg{glmmTMB} package.

\hypertarget{bayesian-inference-marginal-laplace-approximation}{%
\subsection{Bayesian inference: marginal Laplace
approximation}\label{bayesian-inference-marginal-laplace-approximation}}

Using \texttt{AGHQ}{} on the raw output of \pkg{glmmTMB} is feasible in
this example because there are a small number of fixed-effect regression
coefficients. In order to improve computational efficiency we include \(\beta\) as a
parameter to which the marginal Laplace approximation is applied. We may
do this by making explicit use of the \texttt{aghq::marginal\_laplace}
function.

We first construct an object which computes the likelihood using the
interface in \pkg{glmmTMB}:

\begin{CodeChunk}
\begin{CodeInput}
R> Salamanders2 <- Salamanders
R> Salamanders2$intercept <- rep(1,nrow(Salamanders2))
R> Salamanders2$site <- as.character(Salamanders2$site)
R> noremod <- glmmTMB(
+   count ~ -1 + site + intercept + mined,
+   zi = ~mined,
+   disp = ~DOY,
+   data = Salamanders2,
+   family = nbinom2
+ )
\end{CodeInput}
\end{CodeChunk}

As before, we create functions which compute the un-normalized joint
posterior and its derivatives with respect to \(\boldsymbol{W}\), which
contains the \texttt{site} random effects and the regression
coefficients. In contrast to the previous example we compute the prior
derivatives analytically and because the underlying \texttt{TMB} model
has no random effects, an automatically-differentiated log-likelihood
Hessian is available:

\begin{CodeChunk}
\begin{CodeInput}
R> nsite <- length(unique(Salamanders$site))
R> prior_prec_mat <- function(theta) {
+   # prior for W = (beta0,U,beta)
+   # theta = log(sigma), 1/sigma^2 = exp(-2*theta)
+   bdiag(
+     Diagonal(n = nsite,x = exp(-2*theta)),
+     Diagonal(n = 6,x = 1e-03)
+   )
+ }
R> log_posterior <- function(W,theta) {
+   # W = (U1,...,U23,beta1,betazi0,betazi1,betad0,betad1)
+   Q <- prior_prec_mat(theta)
+   as.numeric(
+   -1 * noremod$obj$fn(W) +
+     (1/2) * determinant(Q,logarithm = TRUE)$modulus +
+     (-1/2)*crossprod(W,crossprod(Q,W)) +
+     dexp(exp(theta),log(2),log = TRUE) + theta
+   )
+ }
R>
R> log_posterior_gr <- function(W,theta) {
+   Q <- prior_prec_mat(theta)
+   as.numeric(-1 * noremod$obj$gr(W)) +
+       -1 * as.numeric(crossprod(Q,W))
+ }
R>
R> log_posterior_he <- function(W,theta) {
+   Q <- prior_prec_mat(theta)
+   as(-1 * noremod$obj$he(W),'sparseMatrix') - Q
+ }
R>
R> ff <- list(
+   fn = log_posterior,
+   gr = log_posterior_gr,
+   he = log_posterior_he
+ )
\end{CodeInput}
\end{CodeChunk}

We now, as usual, pass the \texttt{ff} object to
\texttt{aghq::marginal\_laplace}:

\begin{CodeChunk}
\begin{CodeInput}
R> margquadmod <- marginal_laplace(
+   ff,
+   7,
+   startingvalue = list(W = noremod$fit$par,theta = -1)
+ )
\end{CodeInput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> summary(margquadmod)
\end{CodeInput}
\begin{CodeOutput}
AGHQ on a 1 dimensional posterior with  7 quadrature points

The posterior mode is: -0.705

The log of the normalizing constant/marginal likelihood is: -892

The posterior Hessian at the mode is:
     [,1]
[1,] 9.22

The covariance matrix used for the quadrature is...
      [,1]
[1,] 0.108

...and its Cholesky is:
      [,1]
[1,] 0.329

Here are some moments and quantiles for theta:

         mean median   mode    sd  2.
theta1 -0.806 -0.777 -0.705 0.382 -1.71 -0.163
\end{CodeOutput}
\begin{CodeInput}
R> # Plot sigma
R> sigmaplotdatmarg <- compute_pdf_and_cdf(
+   margquadmod$marginals[[1]],
+   transformation = list(
+     totheta = log,
+     fromtheta = exp
+   )
+ )
R> compute_pdf_and_cdf(
+   zipmodquad$marginals[[7]],
+   list(
+     totheta = log,
+     fromtheta = exp
+   )
+ )
+   ggplot(aes(x = transparam,y = pdf_transparam)) +
+   theme_classic() +
+   geom_line(linetype = 'dotdash') +
+   stat_function(fun = dexp,args = list(rate = log(2)),linetype = 'dashed') +
+   geom_line(data = sigmaplotdatmarg,
+             mapping = aes(x = transparam,y = pdf_transparam)) +
+   labs(x = expression(sigma),y = "Density") +
+   coord_cartesian(xlim = c(0,2))
\end{CodeInput}
\begin{figure}

{\centering \includegraphics{01-aghq-softwarepaper_files/figure-latex/doquadmarg2-1}

}

\caption[(Top) Posterior (---) and prior (- - -) for random effects standard deviation computed with the marginal Laplace approximation, compared to that previously obtained using full adaptive quadrature (- $\cdot$ -)]{(Top) Posterior (---) and prior (- - -) for random effects standard deviation computed with the marginal Laplace approximation, compared to that previously obtained using full adaptive quadrature (- $\cdot$ -).}\label{fig:doquadmarg2}
\end{figure}
\end{CodeChunk}

Because the regression coefficients and \texttt{site} random effects are
now being marginalized over, the quadrature output only contains
information about \(\theta = \log\sigma_{u}\), the random effects
standard deviation. We can use the \texttt{aghq::sample\_marginal}
function to generate samples from the approximate mixture posterior
\(\widetilde{\pi}(\boldsymbol{W}|\boldsymbol{Y})\) and then use these
for inference for \(\boldsymbol{W}\):

\begin{CodeChunk}
\begin{CodeInput}
R> # Draw 100,000 samples (very fast):
R> margsamps <- sample_marginal(margquadmod,1e05)
R> # Get the means, standard deviations, and 9
R> sm <- apply(margsamps$samps,1,mean)
R> ssd <- apply(margsamps$samps,1,sd)
R> cilower <- apply(margsamps$samps,1,quantile,probs = .025)
R> cihigher <- apply(margsamps$samps,1,quantile,probs = .975)
R> # Regression coefficients
R> betamargmean <- sm[24:29]
R> betamargsd <- ssd[24:29]
R> betamargcilower <- cilower[24:29]
R> betamargciupper <- cihigher[24:29]
R> glmmtmbest <- reduce(summary(zipmod)$coefficients,rbind)
R>
R> data.frame(
+   margmean = betamargmean,
+   freqest = glmmtmbest[ ,'Estimate'],
+   margsd = betamargsd,
+   freqsd = glmmtmbest[ ,'Std. Error'],
+   margci2.5 = betamargcilower,
+   freqci2.5 = glmmtmbest[ ,'Estimate'] - 2 * glmmtmbest[ ,'Std. Error'],
+   margci97.5 = betamargciupper,
+   freqci97.5 = glmmtmbest[ ,'Estimate'] + 2 * glmmtmbest[ ,'Std. Error']
+ )
\end{CodeInput}
\begin{CodeOutput}
  margmean freqest margsd freqsd margci2.5 freqci2.5 margci97.5 freqci97.5
1   -0.625 -0.5403  0.356  0.376    -1.347    -1.293     0.0411     0.2120
2    1.526  1.4236  0.370  0.365     0.849     0.693     2.2993     2.1545
3    0.173  0.2563  0.470  0.487    -0.789    -0.718     1.0471     1.2304
4   -1.948 -2.2442  0.623  0.745    -3.140    -3.733    -0.6945    -0.7551
5    0.171  0.0278  0.314  0.318    -0.443    -0.608     0.7836     0.6631
6   -0.435 -0.3947  0.165  0.154    -0.759    -0.702    -0.1098    -0.0874
\end{CodeOutput}
\begin{CodeInput}
R> # Random effects
R> # They have the global intercept added to them
R> # in the TMB model, and aren't in the same order...
R> Uest <- sm[1:23] + sm[24]
R> Uglmmtmb <- coef(zipmod)$cond$site[
+   match(
+     stringr::str_remove(
+       rownames(summary(noremod)$coefficients$cond),
+       "site"
+     )[1:23],
+     rownames(coef(zipmod)$cond$site)
+   ),1]
R>
R> data.frame(
+   frequentist = Uglmmtmb,
+   marginal_laplace = Uest
+ )
\end{CodeInput}
\begin{CodeOutput}
   frequentist marginal_laplace
1      -0.4730           -0.555
2      -0.6860           -0.780
3      -0.9436           -1.068
4      -0.6656           -0.748
5      -0.4341           -0.504
6      -0.5401           -0.605
7      -0.1561           -0.206
8      -0.4038           -0.468
9      -0.2616           -0.328
10     -0.7531           -0.850
11     -0.5099           -0.581
12     -0.7369           -0.812
13     -0.5689           -0.624
14     -0.8573           -1.087
15     -0.4634           -0.484
16      0.0826            0.168
17     -0.9176           -1.201
18     -0.6547           -0.753
19     -0.8070           -0.995
20     -0.8214           -1.012
21     -0.1966           -0.137
22     -0.1780           -0.122
23     -0.5671           -0.629
\end{CodeOutput}
\end{CodeChunk}

The posterior means of the random intercepts from the marginal Laplace
approximation are mostly larger in absolute value than the point
estimates from the frequentist model because their standard deviation
\(\sigma_{u}\) is estimated to be smaller in the latter compared to the
former, despite the former having a prior that shrinks it towards zero.
Relative to their standard errors, the regression coefficient point
estimates are similar to their posterior means, which is to be expected.

\hypertarget{sec:discussion}{%
\section{Discussion}\label{sec:discussion}}

The \pkg{aghq} package provides a flexible interface for Bayesian
inference where the posterior is normalized using adaptive quadrature.
Practitioners who use \texttt{MCMC}{} for their models and statisticians
wishing to use \texttt{AGHQ}{} in their research may use the package to
automate many of the difficult and tedious computations associated with
implementing the method for their problem. I have demonstrated the use
of the package both in low-dimensional examples where it provides a
comprehesive set of tools for Bayesian inference at far faster run times
and easier use than \texttt{MCMC}{}, and for a high-dimensional models
where it is used as a key part of a more complicated method for
inference.

A limitation of the implementation at the time of writing is the
inefficency of the marginal Laplace approximation which forms the core
of all of the high-dimensional examples. Numerical derivatives of
\(\log\widetilde{\pi}(\boldsymbol{\theta}|\boldsymbol{Y})\) result in
many repeated optimizations of
\(\log\pi(\boldsymbol{W}|\boldsymbol{\theta},\boldsymbol{Y})\), and this
drastically increases computation time. Current research is focused on
more efficient optimization routines and derivative calculations for
this method, and these will be implemented in the \pkg{aghq} package as
they are developed.

\hypertarget{acknowledgements}{%
\section*{Acknowledgements}\label{acknowledgements}}
\addcontentsline{toc}{section}{Acknowledgements}

I wish to thank my PhD student colleagues Yanbo Tang and Blair Bilodeau
for their helpful input, Professor Gwen Eadie for help with the
astronomy example, and Professors Patrick Brown and Emanuele Giorgi for
their helpful comments. This work was supported by a Postdoctoral
Fellowship from the Natural Sciences and Engineering Research Council of
Canada.

redpajama_set_name: RedPajamaArXiv