Markov Chain Monte Carlo with the Integrated Nested Laplace Approximation
Virgilio Gómez-Rubio Håvard Rue


Human-Robot Mutual Adaptation in Shared Autonomy 


    
      Stefanos Nikolaidis
 
          Carnegie Mellon University 
            
      snikolai@cmu.edu
      Yu Xiang Zhu
    
          Carnegie Mellon University 
            
      yuxiangz@cmu.edu
  David Hsu
                       
    National University of Singapore 
            
    dyhsu@comp.nus.edu.sg  
      Siddhartha Srinivasa
                       
    Carnegie Mellon University 
            
    siddh@cmu.edu

    December 30, 2023
==================================================================================================================================================================================================================================================================================================================================================================================================================================================





    The Integrated Nested Laplace Approximation (INLA) has established
    itself as a widely used method for approximate inference on
    Bayesian hierarchical models which can be represented as a latent
    Gaussian model (LGM). INLA is based on producing an accurate
    approximation to the posterior marginal distributions of the
    parameters in the model and some other quantities of interest by
    using repeated approximations to intermediate distributions and
    integrals that appear in the computation of the posterior
    marginals.

    INLA focuses on models whose latent effects are a Gaussian Markov
    random field (GMRF). For this reason, we have explored alternative
    ways of expanding the number of possible models that can be fitted
    using the INLA methodology. In this paper, we present a novel
    approach that combines INLA and Markov chain Monte Carlo (MCMC).
    The aim is to consider a wider range of models that cannot be
    fitted with INLA unless some of the parameters of the model have
    been fixed. Hence, conditioning on these parameters the model
    could be fitted with the R-INLA package. We show how new
    values of these parameters can be drawn from their posterior by
    using conditional models fitted with INLA and standard MCMC
    algorithms, such as Metropolis-Hastings. Hence, this will extend
    the use of INLA to fit models that can be expressed as a
    conditional LGM. Also, this new approach can be used to build
    simpler MCMC samplers for complex models as it allows sampling
    only on a limited number parameters in the model.

    We will demonstrate how our approach can extend the class of
    models that could benefit from INLA, and how the R-INLA
    package will ease its implementation. We will go through simple
    examples of this new approach before we discuss more advanced
    problems with datasets taken from relevant literature.


Keywords: Bayesian Lasso, INLA, MCMC, Missing Values, Spatial Models





§ INTRODUCTION


Bayesian inference for complex hierarchical models has almost entirely
relied upon computational methods, such as Markov chain Monte Carlo
<cit.>. <cit.> propose a
new paradigm for Bayesian inference on hierarchical models that can be
represented as latent Gaussian models (LGMs), that focuses on
approximating marginal distributions for the parameters in the model.
This new approach, the Integrated Nested Laplace Approximation (INLA,
henceforth), uses several approximations to the conditional
distributions that appear in the integrals needed to obtain the
marginal distributions. See Section <ref> for details.

INLA is implemented as an R package, called R-INLA, that allows
us to fit complex models often in a matter of seconds. Hence, this is
often much faster than fitting the same model using MCMC methods.
Fitting models using INLA is restricted, in practice, to the classes
of models implemented in the R-INLA package. Several authors
have provided ways of fitting other models with INLA by fixing some of
the parameters in the model so that conditional models are fitted with
R-INLA. We have included a brief summary below.


<cit.> provide an early application of the idea of
fitting conditional models on some of the model parameters with
R-INLA. They developed this idea for a very specific example on
spatiotemporal models in which some of the models parameters are fixed
at their maximum likelihood estimates, which are then plugged-in the
overall model, thus ignoring the uncertainty about these parameters
but greatly reducing the dimensionality of the model. However, they do
not tackle the problem of fitting the complete model to make inference
on all the parameters in the model.


<cit.> propose an approach to extend
the type of models that can be fitted with R-INLA and apply
their ideas to fit some spatial models. They note how some models can
be fitted after conditioning on one or several parameters in the
model. For each of these conditional models R-INLA reports the
marginal likelihood, which can be combined with a set of priors for
the parameters to obtain their posterior distribution. For the
remainder of the parameters, their posterior marginal distribution can
be obtained by Bayesian model averaging <cit.> the
family of models obtained with R-INLA.

Although <cit.> focus on some spatial
models, their ideas can be applied in many other examples. They apply
this to estimate the posterior marginal of the spatial autocorrelation
parameter in some models, and this parameter is known to be bounded,
so that computation of its marginal distribution is easy because the
support of the distribution is a bounded interval.

For the case of unbounded parameters, the previous approach can be
applied, but a previous search may be required. For example, the
(conditional) maximum log-likelihood plus the log-prior could be
maximised to obtain the mode of the posterior marginal. This will mark
the centre of an interval where the values for the parameter are taken
from and where the posterior marginal can be evaluated.

In this paper, we will propose a different approach based on Markov
chain Monte Carlo techniques. Instead of trying to obtain the
posterior marginal of the parameters we condition on, we show how to
draw samples from their posterior distribution by combining MCMC
techniques and conditioned models fitted with R-INLA. This
provides several advantages, as described below.


This will increase the number of models that can be fitted using INLA
and its associated R package R-INLA. In particular,
models that can be expressed as a conditional LGM could be fitted. The
implementation of MCMC algorithms will also be simplified as only the
important parameters will be sampled, while the remaining parameters
are integrated out with INLA and R-INLA.
<cit.> have also effectively combined MCMC and
INLA for efficient variable selection and model choice.

The paper is structured as follows. The Integrated Nested Laplace
Approximation is described in Section <ref>. Markov chain
Monte Carlo methods are summarised in Section <ref>. Our
proposed combination of MCMC and INLA is detailed in Section
<ref>. Some simple examples are developed in Section
<ref> and some real applications are provided in Section
<ref>. Finally, a discussion and some final remarks
are provided in Section <ref>.




§ INTEGRATED NESTED LAPLACE APPROXIMATION



We will now describe the types of models that we will be considering
and how the Integrated Nested Laplace Approximation method works. We
will assume that our vector of n observed data
𝐲 = (y_1,…,y_n) are observations from a distribution
in the exponential family, with mean μ_i. We will also assume that
a linear predictor on some covariates plus, possibly, other effects
can be related to mean μ_i by using an appropriate link function.
Note that this linear predictor η_i may be made of linear terms
on some covariates plus other types of terms, such as non-linear
functions on the covariates, random effects, spatial random effects,
etc. All these terms will define some latent effects 𝐱.

The distribution of 𝐲 will depend on a vector of
hyperparameters θ_1. Because of the approximation that INLA
will use, we will also assume that the vector of latent effects
𝐱 will have a distribution that will depend on a vector of
hyperparameters θ_2. Altogether, the hyperparameters can be
represented using a single vector θ=(θ_1, θ_2).

From the previous formulation, it is clear that observations are
independent given the values of the latent effects 𝐱 and
the hyperparameters θ. That is, the likelihood of our model
can be written down as

    π(𝐲|𝐱,θ) =
        ∏_i∈ℐπ(y_i|x_i,θ)

Here, i is indexed over a set of indices
ℐ⊆{1,…,n} that indicates observed
responses. Hence, if the value of y_i is missing then
i ∉ℐ (but the predictive distribution y_i could
be computed once the model is fitted).

Under a Bayesian framework, the aim is to compute the posterior
distribution of the model parameters and hypermeters using Bayes'
rule. This can be stated as

    π(𝐱, θ|𝐲) ∝π(𝐲|𝐱,θ) π(𝐱,θ)

Here, π(𝐱,θ) is the prior distribution of the latent
effects and the vector of hyperparameters. As the latent effects
𝐱 have a distribution that depends on θ_2, it is
convenient to write this prior distribution as
π(𝐱,θ) = π(𝐱|θ) π(θ).

Altogether, the posterior distribution of the latent effects and
hyperparameters can be expressed as

    π(𝐱, θ|𝐲) ∝π(𝐲|𝐱,θ) π(𝐱,θ) =
      π(𝐲|𝐱,θ) π(𝐱|θ)
      π(θ) = 
    π(𝐱|θ) π(θ) ∏_i∈ℐπ(y_i|x_i,θ)

The joint posterior, as presented in Equation (<ref>), is
seldom available in a closed form. For this reason, several estimation
methods and approximations have been developed over the years.

Recently, <cit.> have provided approximations
based on the Laplace approximation to estimate the marginals of all
parameters and hyperparameters in the model. They develop this
approximation for the family of latent Gaussian Markov random fields
models. In this case, the vector of latent effects is a Gaussian
Markov random field (GMRF). This GMRF will have zero mean (for
simplicity) and precision matrix Q(θ_2).

Assuming that the latent effects are a GMRF will let us develop
Equation (<ref>) further. In particular, the posterior
distribution of the latent effects 𝐱 and the vector of
hyperparameters θ can be written as

    π(𝐱, θ|𝐲) ∝π(θ) |𝐐(θ)|^n/2exp{-1/2𝐱^T 𝐐(θ)
        𝐱+∑_i∈ℐlog(π(y_i|x_i, θ)) }.

With INLA, the aim is not the joint posterior distribution
π(𝐱, θ|𝐲) but the marginal
distributions of latent effects and hyperparameters. That is,
π(x_j|𝐲) and π(θ_k|𝐲), where indices
j and k will take different ranges of values depending on the
number of latent effects and hyperparameters.

Before computing these marginal distributions, INLA will obtain an
approximation to π(θ|y),
π̃(θ|𝐲). This approximation will later
be used to compute an approximation to marginals
π(x_j|𝐲). Given that the marginal can be written down as

    π(x_j|𝐲) = ∫π(x_j|θ, 𝐲)
        π(θ| 𝐲) dθ,

the approximation is as follows:

    π̃(x_j|𝐲)=
        ∑_g π̃(x_j|θ_g, 𝐲)×π̃(θ_g|𝐲)×Δ_g.

Here, π̃(x_j|θ_g, 𝐲) is an approximation to
π (x_j|θ_g, 𝐲), which can be obtained using
different methods <cit.>.
θ_g refers to an ensemble of hyperparameters, that
take values on a grid (for example), with weights Δ_g.

INLA is a general approximation that can be applied to a large number
of models. An implementation for the R programming language
is available in the R-INLA package at www.r-inla.org,
which provides simple access to model fitting. This includes a simple
interface to choose the likelihood, latent effects and priors. The
implementation provided by R-INLA includes the computation of
other quantities of interest. The marginal likelihood
π(𝐲) is approximated, and it can be used for model
choice. As described in <cit.>, the approximation
to the marginal likelihood provided by INLA is computed as

    π̃(𝐲) = ∫π(θ, 𝐱,
        𝐲)/π̃_G(𝐱|θ,𝐲)|_𝐱=𝐱^*(θ)d θ.

Here,
π(θ, 𝐱, 𝐲) =
π(𝐲|𝐱, θ)
π(𝐱|θ) π(θ),
π̃_G(𝐱|θ,𝐲) is a
Gaussian approximation to π(𝐱|θ,𝐲)
and 𝐱^*(θ) is the posterior mode of 𝐱 for a
given value of θ. This approximation is reliable when
the posterior of θ is unimodal, as it is often the
case for latent Gaussian models. Furthermore,
<cit.> demonstrate that this approximation is
accurate for a wide range of models.


Other options for model choice and assessment include the Deviance
Information Criterion <cit.> and the
Conditional Predictive Ordinate <cit.>. Other
features in the R-INLA package include the use of different
likelihoods in the same model, the computation of the posterior
marginal of a certain linear combination of the latent effects and
others <cit.>.



§ MARKOV CHAIN MONTE CARLO



In the previous Section we have reviewed how INLA computes an
approximation of the marginal distributions of the model parameters
and hyperparameters. Instead of focusing on an approximation to the
marginals, Markov chain Monte Carlo methods could be used to
obtain a sample from the joint posterior marginal
π( x, θ| y). To simplify the notation, we will
denote the vector of latent effects and hyperparameters by
z = ( x, θ). Hence, the aim now is to estimate
π( z| y) or, if we are only interested on the posterior
marginals, π(z_i| y).

Several methods to estimate or approximate the posterior distribution
have been developed over the years <cit.>. In the case
of MCMC, the interest is in obtaining a
Markov chain whose limiting distribution is π( z |  y). We
will not provide a summary of MCMC methods here, and the reader is
referred to <cit.> for a detailed description.

The values generated using MCMC are (correlated) draws from
π( z |  y) and, hence, can be used to estimate quantities of
interest. For example, if we are interested in marginal inference on
z_i, the posterior mean from the N sampled values
{z_i^(j)}_j=1^N can be estimated using the empirical mean
of {z^(j)_i}_j=1^N. Similarly, estimates of the posterior
expected value of any function on the parameters f(𝐳) can be found
using that 


    E[f(𝐳)| y] ≃1/N∑_j = 1 ^N f(𝐳^(j))


Multivariate estimates inference can be made by using the multivariate
nature of vector z^(j). For example, the posterior covariance
between parameters z_k and z_l could be computed by considering
samples {(z_k^(j)), z_l^(j))}_j=1^N.



 §.§ The Metropolis-Hastings algorithm


This algorithm was firstly proposed by <cit.> and
<cit.>. The Markov chain is generated by proposing new
moves according to a proposal distribution q(·|·). The new
point is accepted with probability

    α = min{
          1,
          π( z^(j+1)| y) q( z^(j) | z^(j+1))/π( z^(j)| y) q( z^(j+1)| z^(j))}

In the previous acceptance probability, the posterior probabilities of
the current point and the proposed new point appear as
π(z^(j)| y) and π(z^(j+1)| y), respectively. These
two probabilities are unknown, in principle, but using Bayes' rule
they can be rewritten as

    π( z| y) = π( y | z)π( z)/π( y)

Hence, the acceptance probability α can be rewritten as

    α = min{
          1,
          π( y | z^(j+1))π( z^(j+1)) q( z^(j)
              |
               z^(j+1))/π( y | z^(j))π( z^(j))
              q( z^(j+1)| z^(j))}

This is easier to compute as the acceptance probability depends on
known quantities, such as the likelihood π( y | z), the
prior on the parameters π( z) and the probabilities of the
proposal distribution. Note that the term π( y) that appears in
Equation (<ref>) is unknown but that it cancels out as it
appears both in the numerator and denominator.

In Equation (<ref>) we have described the move to sample from
the joint ensemble of model parameters. However, this can be applied
to individual paramaters one at a time, so that acceptance
probabilities will be

    α = min{
          1,
          π( y |z_i^(j+1))π(z_i^(j+1)) q(z_i^(j) |
              z_i^(j+1))/π( y |z_i^(j))π(z_i^(j))
              q(z_i^(j+1)| z_i^(j))}




§ INLA WITHIN MCMC



In this Section, we will describe how INLA and MCMC can be combined to
fit complex Bayesian hierarchical models. In principle, we will assume
that the model cannot be fitted with R-INLA unless some of the
parameters or hyperparameters in the model are fixed. This set of
parameters is denoted by z_c so that the full ensemble of
parameters and hyperparameters is z = ( z_c,  z_-c).
Here z_-c is used to denote all the parameters in z that
are not in z_c. Our assumptions are that the posterior
distribution of z can be split as

    π( z | y) ∝π( y| z_-c) π( z_-c| z_c) π( z_c)

and that π( y| z_-c) π( z_-c| z_c) is a latent
Gaussian model suitable for INLA. This means that conditional models
(on z_c) can still be fitted with R-INLA, i.e., we can
obtain marginals of the parameters in z_-c given z_c.
The conditional posterior marginals for the k-th element in vector
z_-c will be denoted by π(z_-c,k | z_c,  y). Also, the
conditional marginal likelihood π( y |  z_c) can be easily
computed with R-INLA.



 §.§ Metropolis-Hastings with INLA


We will now discuss how to implement the Metropolis-Hastings algorithm
to estimate the posterior marginal of z_c. Note that this is a
multivariate distribution and that we will use block updating in the
Metropolis-Hastings algorithm. Say that we start from an initial point
z_c^(0) then we can use the Metropolis-Hastings algorithm to
obtain a sample from the posterior of z_c.

We will draw a new proposal value for z_c, z_c^(1),
using the proposal distribution q(·|·). The acceptance
probability, shown in Equation (<ref>), becomes now:

    α = min{
          1,
          π( y | z_c^(j+1))π( z_c^(j+1))
              q( z_c^(j) | z_c^(j+1))/π( y |
               z_c^(j))π( z_c^(j))
              q( z_c^(j+1)| z_c^(j))}

Note that π( y | z_c^(j)) and
π( y | z_c^(j+1)) are the conditional marginal likelihoods
on z_c^(j) and z_c^(j+1), respectively. All these
quantities can be obtained by fitting a model with R-INLA with
the values of z_c set to z_c^(j) and z_c^(j+1).
Hence, at each step of the Metropolis-Hastings algorithm only a model
conditional on the proposal needs to be fitted.

π( z_c^(j)) and π( z_c^(j+1)) are the priors of
z_c evaluated at z_c^(j) and z_c^(j+1),
respectively, and they can be easily computed as the priors are known
in the model. Values q( z_c^(j) | z_c^(j+1)) and
q( z_c^(j+1)| z_c^(j)) can also be computed as the
proposal distribution is known. Hence, the Metropolis-Hastings
algorithm can be implemented to obtain a sample from the posterior
distribution of z_c. The marginal distribuions of the elements
of z_c can be easily obtained as well.

Regarding the marginals of z_-c,k, it is worth noting that at step
j of the Metropolis-Hastings algorithm a conditional marginal
distribution on z_c^(j) (and the data y) is obtained:
π(z_-c,k|  z_c^(j),  y). The posterior marginal can be
approximated by integrating over z_c as follows:

    π(z_-c,k|  y) = ∫π(z_-c,k|
         z_c,  y) π( z_c| y)
        d  z_c≃1/N∑_j=1^N π(z_-c,k|  z_c^(j),  y),

where N is the number of samples of the posterior distribution of
z_c. That is, the posterior marginal of z_-c,k can be
obtained by averaging the conditional marginals obtained at each
iteration of the Metropolis-Hastings algorithm.




 §.§ Effect of approximating the marginal likelihood


So far, we have ignored the fact that the conditional marginal
likelihood π( y |  z_c) used in the acceptance probability
α is actually an approximation. In this section, we will
discuss how this approximation will impact the validity of the
inference.

The situations where a Metropolis-Hastings algorithm has inexact
acceptance probabilities are often called pseudo-marginal MCMC
algorithms and were first introduced in <cit.> in the
context of statistical genetics where the likelihood in the acceptance
probability is approximated using importance sampling.
<cit.> provided a more general justification of
the pseudo-marginal MCMC algorithm, whose properties are further
studied in <cit.> and
<cit.>. These results show that if the (random)
acceptance probability is unbiased then the Markov chain will still
have as stationary distribution the posterior distribution of the
model parameters.

In our case, the error in the acceptiance rate is coming from a
deterministic estimate of the conditional marginal likelihood, hence
the framework of pseudo-marginal MCMC does not apply. However, since
it is deterministic, our MCMC chain will converge to a
stationary distribution. This limiting distribution will be

    π̃( z_c |  y) ∝π( z_c)
        π̃( y |  z_c)

where the “∼” indicates an approximation. R-INLA returns
an approximation to the conditional marginal likelihood term, which
implies an approximation to π( z_c |  y). This leaves the
question, about how good this approximation is, for which we have to
rely on asymptotic results, heuristics and numerical experience.

The conditional marginal likelihood estimate returned from
R-INLA is based on numerical integration and uses a sequence of
Laplace approximations <cit.>. This
estimate is more accurate than the classical estimate using one
Laplace approximation. This approximation has, with classical
assumptions, relative error 𝒪(n^-1) <cit.>,
where n is the number of replications in the observations. For our
purpose, this error estimate is sufficient, as it demonstrates that

    π̃( z_c |  y)/π( z_c |  y)∝π̃( y |  z_c)/π( y|  z_c) =
        1 + 𝒪(n^-1)

for plausible values of z_c. However, as discussed by
<cit.>, the classical assumptions are
rarely met in practice due to “random effects”, smoothing etc.
Precise error estimates under realistic assumptions are difficult to
obtain; see <cit.> for a more detailed discussion of this
issue.

About numerical experience with the conditional marginal likelihood
estimate, <cit.> have studied empirically its
properties and accuracy for a wide range of latent Gaussian models.
They have compared the estimates with those obtained using MCMC, and
in all their cases the approximates of the marginal likelihood
provided by INLA were very accurate. For this reason, we believe that
the approximate stationary distribution π̃( z_c |  y)
should be close to the true one, without being able to quantify this
error in more details.

Although the error in Equation (<ref>) is pointwise, we do expect
the error would be smooth in z_c. This is particularly
important, as in most cases we are interested in the univariate
marginals of π̃( z_c |  y). These marginals will
typically have less error as the influence of the approximation error
will be averaged out integrating out all the other components. A final
renormalization would also remove constant offset in the error.

Additionally, we will validate the approximation error in a simulation
study in Section <ref> where we fit various models using
INLA, MCMC and INLA within MCMC and very similar posterior
distributions are obtained. Furthermore, the real applications in
Section <ref> also support that the approximations to
the marginal likelihood are accurate.



 §.§ Some remarks


Common sense is still not out of fashion, hence there is an implicit
assumption that our INLA within MCMC approach should be only for models for
which it is reasonable to use the INLA-approach to do the inference
for the conditional model. The procedure that we have just shown will
allow INLA to be used together with the Metropolis-Hastings algorithm
(and, possibly, other MCMC methods) to obtain the posterior
distribution (and marginals) of z_c and the posterior marginals
of the elements in z_-c. Hence, this will allow INLA to be
used to fit models not implemented in the R-INLA package as well
as providing other options for model fitting, that we summarise here.

The Metropolis-Hastings algorithm will allow any choice of the priors
on the set of parameters z_c. This is an advantage (as shown in
the example in Section <ref>) of combining MCMC and INLA
because priors that are not implemented in R-INLA can be used in
the model. In particular, improper flat priors, multivariate priors
and objective priors can be used.

The framework of conditional LGMs that we now can fit using our new
approach is quite rich. It includes models with missing covariates
that are imputed at each step of the Metroplis-Hastings algorithm (see
example in Section <ref>), models with complex
non-linear effects in the linear predictor (see example in
Section <ref>) or models that have a mixture of effects in
the linear predictor <cit.>.



§ SIMULATION STUDY



In this section we develop simple examples to illustrate the method
proposed in the previous sections, and we investigate how this new
approach works in practice.



 §.§ Bivariate linear regression



The first example is based on a linear regression with two covariates.
Our aim is to use our proposed method to obtain the posterior
distribution of the coefficients of the two covariates and then
compare the estimated marginals to the results obtained when the full
model is fitted with MCMC and INLA.

The simulated dataset contains 100 observations of a response variable
y and covariates u_1 and u_2. The model used to
generate the data is a typical linear regression, i.e.,


    y_i = α + β_1 u_1i + β_2 u_2i +
        ε_i; i = 1, …, 100

Here, ε_i is a Gaussian error term with zero mean and
precision τ. The dataset has been simulated using α = 3,
β_1 = 2, β_2 = -2 and τ = 1. Covariates u_1i and
u_2i have also been simulated using a uniform distribution between
0 and 1 in both cases.

This model can be easily fitted using R-INLA. Given that we are
using a Gaussian model, inference is exact in this case (up to
integration error). For this reason, we can compare the marginal
distribution provided of β_1 and β_2 by INLA to the ones
obtained with our combined approach. Note that the Metropolis-Hastings
algorithm will provide the joint posterior distribution of
β = (β_1, β_2) that can be use to obtain the
posterior marginals of β_1 and β_2. Furthermore, we can
also compare the marginals of α and τ, that will be
estimated by averaging the different conditional marginals obtained in
the Metropolis-Hastings steps.

In order to implement the Metropolis-Hastings algorithm to obtain a
sample from π(β| y) we have chosen a starting point of
β^(0) = (0, 0). The transition kernel to obtain a
candidate β^(t+1) at iteration t has been a bivariate
Guassian kernel centered at β^(t) with diagonal
variance-covariance matrix with values 1/0.75^2 in the diagonal as
this provided a resonable acceptance rate. The prior distribution of
β has been the product of two Gaussian distributions with
zero mean and precission 0.001 because these are the default priors
for linear effects in R-INLA.

Figure <ref> shows a summary of the results. Given
that both covariates are independent, their coefficients should show
small correlation and this can clearly seen in the plot of the joint
posterior distribution of β. Also, it can be seen how the
marginals obtained with INLA within MCMC for β_1 and β_2
match those obtained with INLA and MCMC. In addition, we have included
the estimates of the posterior marginals of the intercept α and
the precission τ. When using INLA within MCMC these are obtained
by Bayesian model averaging over the fitted models at every step of
the Metropolis-Hastings algorithm, whilst when computed with R-INLA these
are obtained by using INLA alone. The three estimation methods provide very
similar posterior distributions of the posterior marginals of the intercept and
the precission, which again confirms the accuracy of INLA within MCMC.







 §.§ Missing covariates


In the next example, we will discuss the case of missing covariates.
In this example we will consider a linear regression with covariate
u_1 only and we will assume that a number of values of the
covariates are missing. The aim is to include the imputation of this
variables into the model, so that the output is a marginal
distribution of the missing values. We will not discuss here the
different frameworks under which the values have gone missing, but
this is something that should take into account in the model. In
particular, we have removed the values of 9 covariates, which is
almost 10% of our data and summary plots can be included in a 3x3
plot. Hence, in this case the missingness mechanism is of the type
missing completely at random <cit.>.

Now, we will treat the missing values as if they were covariates. We
will use a block updating scheme as we can have a large number of
missing covariates. The transition kernel will be a multivariate
Gaussian with diagonal variance-covariance. The mean and variance for
all values are the mean and variance of the observed covariates,
respectively. The prior distribution is also a multivariate Gaussian,
but now with zero mean and diagonal variance-covariance matrix with
entries four times the variance of a uniform random variable in the
unit interval (the one used to simulate the covariates). This is done
so that the prior information is small compared to the information
provided by the covariates.



Figure <ref> shows the posterior marginals
obtained from the samples. As it can be seen, most of them are
centered at the actual values removed from the model. 
Note that this time the model with missing covariates cannot be fitted
with R-INLA so that we can only compare the marginals to those
obtained with MCMC. In all cases the marginals obtained with INLA
within MCMC and full MCMC are very similar.







 §.§ Poisson regression


In this example we consider a Poisson regression with two covariates:

    y_i ∼ Po(μ_i); log(μ_i) = α + β_1 u_1i + β_2 u_2i;    i = 1, …, 100.

The values of the parameters used to simulate the dataset are
α = 0.5, β_1 = 2 and β_2 = -2.

As in Section <ref>, our purpose is to estimate the
joint posterior distribution of (β_1, β_2).
The prior distribution used now is the
same as in the previous example. Hence, the posterior marginal of
α is obtained by combining the different conditional marginals
obtained at the different steps of the Metropolis-Hastings algorithm.

Figure <ref> shows the estimates of the marginal
distributions of the three parameters in the model, together with the
joint posterior distribution of β_1 and β_2. In all cases,
there is a very good agreement between the estimates obtained with
INLA, MCMC and INLA within MCMC of the marginals of the parameters in
the model.






§ APPLICATIONS



In this section, we will focus on some real life applications that provide a
more reallistic test of this methodology.  In all the examples, we have run
INLA within MC and MCMC for a total of 100500 simulations and discarded the
first 500.  Then we applied a thinning to keep one in ten iterations, to obtain
a final chain of 10000 samples. This includes samples from the missing
observations and fitted models. To fit the model using MCMC alone, we have used
rjags <cit.> with the same number of iterations and thinning. 



 §.§ Bayesian Lasso



The Lasso <cit.> is a popular regression and variable
selection method for variable selection. It has the nice property of
providing coefficient estimates that are exactly zero and, hence, it
performs model fitting and variable selection at the same time. For a
linear model with a Gaussian likelihood, the Lasso is trying to
estimate the regression coeffcient by minimising

    ∑_i=1^n ( y_i - β_0 - ∑_j=1^p β_j
      x_ij)^2 + λ∑_j=1^p |β_j|

Here, y_i is the response variable and x_ij are associated
covariates. n is the number of observations and p the number of
covariates. λ is a non-negative penalty term to control how
the shrinkage of the coefficients is done. If λ = 0 then the
fitted coefficients are those obtained by maximum likelihood, whilst
higher values of λ will shrink the estimates towards zero.

The Lasso is closely related to Bayesian inference as it can be
regarded as a standard regression model with Laplace priors on the
variable coefficients. The Laplace distribution is defined as

    f(β) = 1/2σexp(-|β -
          μ|/σ), x∈ℝ

where μ and σ, a positive number, are parameters of
location and scale, respectively. The Laplace prior distribution is
not available for (parts of) the latent field in R-INLA.
However, conditioning on the values of the β-coefficients
the model can be easily fitted with R-INLA.

We will apply the methodology described in this paper to implement the
Bayesian Lasso by combining INLA and MCMC. We will be using the
Hitters dataset described in <cit.>. This
dataset records several statistics about players in the Major League
Baseball, including salary in 1987, number of times at bat in 1986 and
other variables. Our aim is to build a model to predict the player's
salary in 1987 on some of the other variables recorded in 1986 (the
previous season).

We will focus on a smaller model than the one described in
<cit.> and will consider predicting salary in 1987 on
only five variables measured from the 1986 season: number of times at
bat (AtBat), number of hits (Hits), the number of home runs (HmRun),
number of runs (Runs) and the number of runs batted (RBI).


For our implementation of the Bayesian Lasso, we
will be fitting models conditioning on the covariate coefficients
β = (β_1,…β_p). Also, we will assume that
β and the error term precision τ are independent a
priori, i.e., π(β, τ) = π(β) π(τ). This
will provide a simpler way to compare our results with the Lasso and
it will also make computations a bit simpler. However, note that
choosing it is also possible to choose a prior so that
π(β, τ) = π(β|τ)π(τ) <cit.> The posterior distribution of these
variables will be obtained using MCMC.

The summary of the Lasso estimates are available in Table
<ref> and the posterior distribution of the coefficients is
in Figure <ref>. In all cases, there is agreement between the
Lasso and Bayesian Lasso estimates. Also, the posterior distribution
of the model coefficients is the same for MCMC and combining INLA with
MCMC. For those coefficients with a zero estimate with the Lasso, the
posterior distribution obtained with the Bayesian Lasso is centered at
zero.







 §.§ Imputation of Missing Covariates



<cit.> describe the package mice that implements
several multiple imputation methods. We will be using the
nhanes dataset to illustrate how our approach can be used to
provide imputation of missing coariates in a real dataset. This
dataset contains data from <cit.> on age, body mass
index (bmi), hypertension status (hyp) and cholesterol
level (chl). Age is divided into three groups: 20-39,
40-59, 60+.


Our aim is to impute missing covariates in order to fit a model that
explains the cholesterol level on age and body mass index. Although
the values of age have been completely observed, there are missing
values in body mass index and cholesterol level. INLA can handle
missing values in the response (and will provide a predictive
distribution) but, as already stated, is not able to handle models
with missing values in the covariates.

We will consider a very simple imputation mechanism by assigning a
Gaussian prior to the missing values of body mass index. This Gaussian
is centred at the average value of the observed values (26.56) and
variance four times the variance of the observed values (71.07,
altogether). With this, we expect to provide some guidance on how the
imputed values should be but allowing for a wide range of variation.
More complex imputation mechanisms could be considered <cit.>. As in previous examples, we will fit
the same model using MCMC in order to compare both results. The model
that we will fit is:

    [                                                                     chl_i                                                                         = β_0 + β_1 bmi_i +
                               β_2 age2_i+ β_3 age3_i + ε_i;                                                                       β_0                                                                         ∝                                                                         1;                                                                       β_k                                                                         ∝                                                  N(0, 0.001); k = 1, 2, 3;                                                                       ε_i                                                                         ∼                                                                   N(0, τ);                                                                         τ                                                                         ∼                                                            Ga(1, 0.00005) ]



Figure <ref> shows the posterior
marginal distributions of the imputed values of the body mass index.
Both MCMC and our approach provide very similar point estimates. Table
<ref> summarises the model parameters obtained both with MCMC
and our approach and Figure <ref> displays the posterior
marginals of the model parameters obtained with our approach and MCMC.
In all cases, the marginals agree, and the point estimates look very
similar.











 §.§ Spatial econometrics models



<cit.> describe a novel approach to extend the
classes of models that can be fitted with R-INLA to fit some
spatial econometrics models. In particular, they fit several
conditional models by fixing the values of some of the parameters in
the model, and then they combine these models using a Bayesian model
averaging approach <cit.>. <cit.>
show a practical implementation with a spatial statistics model using
R package INLABMA. Some of these models have already been
included in R-INLA <cit.> but are still
considered as experimental.

In this example we will focus on one of the spatial econometrics
models described in <cit.> to illustrate how our new
approach to combine MCMC and R-INLA can be used to fit
unimplemented models. In particular, we will consider the spatial lag
model <cit.>:

    y = ρ W  y +  X β +  u;  u ∼ N(0,
    σ^2_u  I)

Here, y is a vector of observations at n areas, W is an
adjacency matrix, ρ a spatial autocorrelation parameter, X
a n× p matrix of covariates with associated coeffients
β = (β_1,…,β_p) and u = (u_1,…, u_n)
an error term. u_i, i=1,…,n, is Normally distributed with zero
mean and precision τ_u. This model can be rewritten as follows:

    y = ( I_n - ρ W)^-1 Xβ + ε;ε∼ N(0, 1/τ_u (( I_n-ρ W) (
    I_n-ρ W^'))^-1)

This model is difficult to fit with any standard software for
mixed-effects models because of parameter ρ. If the value of
ρ is fixed, then it is easy to fit the model with R-INLA as
it becomes a linear term on the covariates plus a random effects term
with a known structure. Hence, by conditioning on the value of ρ
we will be able to fit the model with R-INLA. In order to use
our new approach, we will be drawing values of ρ using MCMC and
conditioning on this parameter to fit the models with R-INLA.

Regarding prior distributions, ρ is assigned a uniform between
-1.5 and 1, β_i, i = 1,…, p a Gaussian prior with zero
mean and precision 0.001 (the default), and τ_u is assigned a
Gamma distribution with parameters 1 and 0.00005 (the default for
the precision of a 'generic0' latent class in R-INLA).

We have fitted this model to the Columbus dataset available in R
package spdep. This dataset contains information about 49
neighbourhoods in Columbus (Ohio) and we have considered a model with
crime rates as the response and household income and housing value as
covariates. We have also fitted the spatial lag model using a maximum
likelihood approach, the method proposed by <cit.>
and MCMC using an implementation of the model for the Jags software
included in package SEjags, which can be downloaded from Github.
The results are shown in Table <ref>. All Bayesian
approaches have very similar estimates, and these are also very
similar to the maximum likelihood estimates.






§ DISCUSSION



In this paper, we have developed a novel approach to extend the models
that can be fitted with INLA. For this, we have used INLA within the
Metropolis-Hastings algorithm, so that only a few number of parameters
are sampled.

We have shown three important applications. In the first one, we have
implemented a Bayesian Lasso for variable selection using Laplace
priors on the coefficient of the covariates. By following this
example, other priors could be easily used with INLA. This includes
not only univariate priors but multivariate priors, that are seldom
available in R-INLA.

In our second example we have tackled the problem of imputation of
missing covariates in model fitting. Here, we have included a very
simple imputation method for the missing values in the covariates, so
that model fitting and imputation were done at the same time. Compared
to fitting the same model with MCMC, we obtained the same posterior
estimates. In an ongoing work, <cit.> explore how
this can be extended to larger problems and how different imputation
models and missingness mechanisms can be properly addressed with INLA
and MCMC.

Finally, we have also shown how other models not included in the
R-INLA software can be fitted with INLA and MCMC. In particular,
we have fitted a spatial econometrics model by fitting conditional
models on the spatial autocorrelation parameter. This method can be
easily modified to suit any other models. In particular, Gibbs
sampling could be used if the full conditionals are available for a
subset of model parameters.

To sum up, we believe that this approach can be employed together with
INLA to fit more complex models and that it can also be combined with
other MCMC algorithms to develop simple samplers to fit complex
Bayesian hierarchical models. This method can work well when the
conditional models are hard to explore with current approaches for
which INLA provides a fast approximation, such as geostatistical
models. Furthermore, INLA could be embedded into a Reversible Jump
MCMC algorithm so that once the model dimension has been set, the
resulting model is approximated with INLA. See, for example,
<cit.> for a comprehensive list of MCMC algorithms that
could benefit from embedding INLA.



§ ACKNOWLEDGEMENTS


Virgilio Gómez-Rubio has been supported by grant PPIC-2014-001, funded
by Consejería de Educación, Cultura y Deportes (JCCM) and FEDER,
and grant MTM2016-77501-P, funded by Ministerio de Economía y
Competitividad.

Chicago 

