Mixture model fitting using conditional models and modal Gibbs sampling
Virgilio Gomez-Rubio
stat.CO stat.ME


#1 1


1
  
  
  Mixture model fitting using conditional models and modal Gibbs sampling
    Virgilio Gómez-Rubio
  This work has been supported by grant PPIC-2014-001-P, 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.

    Department of Mathematics, Universidad de Castilla-La Mancha (Spain)
    Received <date>; Accepted <date>
=================================================================================================================================================================================================================================================================================================================

 

0
  

    Mixture model fitting using conditional models and modal Gibbs sampling


   













Mixture models are a convenient way of modeling data using a convex combination
of different parametric distributions.  In this paper, we present a novel
approach to fitting mixture models based on estimating first the posterior
distribution of the auxiliary variables that assign each observation to a group
in the mixture.  The posterior distributions of the remainder of the parameters
in the mixture is obtained by averaging over their conditional posterior
marginals on the auxiliary variables using Bayesian model averaging.

A new algorithm based on Gibbs sampling is used to approximate the posterior
distribution of the auxiliary variables without sampling any other parameter in
the model. In particular, the modes of the full conditionals of the parameters
of the densities in the mixture are computed and these are plugged-in to the
full conditional of the auxiliary variables to draw samples. This approximation,
that we have called 'modal' Gibbs sampling,
reduces the computational burden in the Gibbs sampling algorithm and still
provides very good estimates of the posterior distribution of the auxiliary
variables.  Conditional models on the auxiliary variables are fitted using the
Integrated Nested Laplace Approximation (INLA) to obtain the conditional
posterior distributions, including modes, of the distributional parameters in the mixtures. 

This approach is general enough to consider mixture models with discrete or
continuous outcomes from a wide range of distributions and latent models as
conditional model fitting is done with INLA. This presents several other
advantages, such as fast fitting of the conditional models, not being
restricted to the use of conjugate priors on the model parameters and being
less prone to label switching.  Within this framework, computing the marginal
likelihood of the mixture model when the number of groups in the mixture is
known is easy and it can be used to tackle selection of the number of
components.  Finally, a simulation study has been carried out to assess the
validity of 'modal' Gibbs sampling and two examples on well-known datasets are
discussed using a mixture of Gaussian and Poisson distributions, respectively.




Keywords:  Bayesian model averaging; INLA; Marginal likelihood; Mixture models; Model selection


1.45 



§ INTRODUCTION


Mixture models are a convenient way of describing data when these are thought to
come from different groups or components which are represented by
different distributions.  Some well known examples of this type of data include
the velocities of galaxies <cit.>, waiting and
eruption times of the Old Faithful geyser in the Yellowstone National Park
<cit.> and the number of fetal movements in lambs 
<cit.>.  In all these cases, the data
generating process is made of different components with different parametric
densities and the observed data does not include information about the
component to which observation belongs. For this reason, mixture models are
often represented using latent auxiliary variables to indicate to which
component each observation belongs and the parametric distribution of each
component. 


In this paper a new approach to fitting mixture models is developed by focusing
on estimating the posterior distribution of the auxiliary variables first. Once
this posterior distribution has been obtained, the posterior marginal of the
remainder of the parameters in the model is obtained by Bayesian model
averaging over conditional  posterior distributions given the values of the
auxiliary variables. This is possible because, given the value of the
auxiliary variables, a mixture model becomes a model with different
independent components and the posterior marginals of their parameters
can easily be obtained.

Hence, model fitting will be split in two steps. First of all, the posterior
distribution of the latent auxiliary variables will be estimated.  This will be
done by means of what we have called 'modal' Gibbs sampling because some of the
parameters will not be sampled but the modes of their full conditional
distributions will be used instead.
Secondly, conditional models on these auxiliary variables will be fitted by
considering values of the auxiliary variables with non-negligible posterior
probability.  These conditional models will be fitted with INLA during 'modal'
Gibbs sampling and the resulting conditional marginals will be combined using
Bayesian model averaging <cit.>.
As conditional models are fitted with INLA, it is possible
to consider a wide range of distributions for the mixture components and
the use of non-conjugate priors for their parameters.
Furthermore, this approach seems to be less likely to suffer from label
switching and it allows for the computation of measures for the
selection of the optimal number of components, such as the marginal likelihood.

The analysis of mixture models using Bayesian inference has been considered by
several authors. Including a set of auxiliary variables to describe to which
component each observations belongs simplifies model fitting with Markov chain
Monte Carlo (MCMC) methods, as described in Section <ref>.  Gibbs
sampling <cit.> and Metropolis-Hastings
<cit.> algorithms are popular MCMC methods
to fit mixture models. See, for example, <cit.> for a
discussion.  More recently, <cit.> describe how to combine
the Integrated Nested Laplace Approximation (INLA) within the
Metropolis-Hastings algorithm to fit models that can be expressed as a
conditional Gaussian Markov random field (GMRF).  Conditional on the values of
the auxiliary variables that assign each observation to a group, a mixture
model becomes a model with several likelihoods that can be easily fitted within
the INLA framework.

In this paper we follow a different approach from <cit.> to
estimate the posterior distribution of the auxiliary variables.  A variation of
the Gibbs sampling algorithm is developed to use the modes of the full
conditional distributions for some of the parameters of the model and these are
plugged-in to the full conditional distribution of the auxiliary variables to
draw samples.


This
provides a good approximation to their posterior distribution, as confirmed by
a simulation study and the examples developed later.

The paper is organized as follows. In Section <ref> an introduction
to mixture models and Bayesian inference is provided. Next, INLA is described
in Section <ref>. Fitting mixture models with INLA is detailed in
Section <ref>. The computation of the marginal likelihood of a
model for the selection of the number of components in the mixture is tackled
in Section <ref>.  A simulation study is carried out in
Section <ref>.  Several examples are developed in
Section <ref>. Finally, the developments and main results of this
paper are discussed in Section <ref>.



§ BAYESIAN INFERENCE ON MIXTURE MODELS





 §.§ Mixture models


Given a set of n observations 𝐲 = (y_1,…,y_n), a mixture model with
K components or groups is usually defined as follows:


    y_i ∼∑_j=1^K w_j f_j(y_i|θ_j), i=1,…,n

Here, f_j(·|θ_j) is the distribution of group or component j in
the mixture, which is defined by a set of parameters θ_j. Parameters
w_j indicate the weight of each component in the mixture and these are taken to sum up to 1,
i.e., ∑_j=1^K w_j = 1. Weights depend on the number of observations
that belong to each component and the vector of weights {w_j}_j=1^K 
will be denoted by 𝐰.
The ensemble of parameters {θ_j}_j=1^K will be denoted by vector
θ and a single element of it by θ_∙∈θ.

Typically, the form of distributions f_j(·|θ_j), and their
associated parameters θ_j, will depend on the actual problem. For
example, for continuous observations, f_j(·|θ_j) can be a Normal
distribution with parameters θ_j = (μ_j, τ_j), where μ_j is the
mean and τ_j the precision. When the observations are count data,
distributions f_j(·|θ_j) could well be Poisson with different means,
i.e., θ_j = λ_j.



An alternative formulation of the model, very useful for Bayesian inference
with MCMC, is to express the mixture by including auxiliary variables
{z_i}_i=1^n to indicate to which component observation i belongs.
These auxiliary variables are defined so that the probability of observation
i belonging to component j, conditional on 𝐰, is:


    π(z_i = j|𝐰) =  w_j, j=1,…,K

This means that, given z_i = j, the distribution of y_i is 
f_z_i(·|θ_z_i). The ensemble of auxiliary variables
(z_1, … , z_n) will be denoted by 𝐳.







 §.§ Bayesian inference




We have followed <cit.> in order to provide a brief summary of
Bayesian inference for mixture models below and the description of the Gibbs
sampling algorithm for mixture models described in Section <ref>. First of all, the graphical
representation of a mixture model in Figure <ref> shows the
relationship between the different parameters in the mixture model.



Observations y_i are independent given z_i and
θ_z_i, hence:


    π(𝐲|𝐳, θ) =
      ∏_i=1^n f_z_i(y_i|θ_z_i)



Note also that the prior distribution on 𝐳 is π(𝐳|𝐰), which is the distribution of 𝐳
given weights 𝐰 and it is defined as:


    π(𝐳|𝐰) = ∏_j=1^K w_j^n_j

where n_j is the number of observations assigned to component j.

Finally, in order to complete our Bayesian formulation of the model, priors on
θ and 𝐰 must be defined.  For convenience, conjugate
priors are often used as this will make inference using MCMC easier (see,
Section <ref> below for details).  However, as we will see in
Section <ref>, sampling on θ and 𝐰 will
not be done in practice.  The prior on θ will depend on the
distributions in the mixture. For example, for a mixture of Normal
distributions, means {μ_j}_j=1^K can be assigned a vague Normal prior
(e.g., centered at zero with a small precision)
whilst precisions {τ_j}_i=1^K can be assigned vague Gamma priors.
For a mixture of Poisson distributions, means {μ_j}_j=1^K can be
assigned a vague Gamma prior. However, given that model fitting will be
done with INLA there is a wider choice of priors available.

Weights 𝐰 can be assigned a Dirichlet prior with parameters
(α_1,…, α_K).  A convenient vague prior on 𝐰 is a
Dirichlet with parameters α_1 = … = α_K = 2, which also
ensures that the mode of this prior distribution exists.

It is worth noting that the prior on θ can be used to include some
identifiability constraints in the model <cit.>.
For example, in a mixture model with two components defined by Normal
distributions with means μ_1 and μ_2, respectively, the prior may be
informative so that μ_1 < μ_2, or include a constraint to make sure that
μ_1 < μ_2. This will make the two components fully identifiable. Hence,
the prior in this case could be proportional to π(μ_1)π(μ_2)I(μ_1 <
μ_2), where I(·) is the indicator function.





 §.§ Gibbs sampling



Gibbs sampling has been a convenient approach to obtain samples from the
posterior distribution of the model parameters in mixture models
<cit.>. <cit.> have derived the full conditionals of
the mixture model described in Section <ref>.


For weight parameters 𝐰 the full conditional is a Dirichlet
distribution with parameters (α_1 + n_1,…, α_K + n_K), where
n_j represents the number of observations allocated to component j, with
j=1,…, K. Note that this full conditional only depends on 𝐳
and π(𝐰).  It may be worth setting α_1 = … = α_K
to 2 (instead of 1) so that the mode of the full conditional exists when some
n_j is set to zero (i.e., if there are no observations assigned to group j
at some point during Gibbs sampling).

The full conditional distribution for grouping variable z_i is a discrete
distribution that will provide the probability of observation i being assigned
to component j. This can be stated as


    π(z_i = j|𝐰,θ,𝐲) ∝ w_j f_j(y_i|θ_j), j=1,…,K

Actual probabilities are computed by obtaining the terms in the right hand
side first and then re-scaling the probabilities to sum up to one:


    π(z_i = j|𝐰,θ,𝐲) = w_j f_j(y_i|θ_j)/∑_j=1^K w_j f_j(y_i|θ_j)


Finally, the full conditionals of θ_∙∈θ will depend on
the actual distribution being used but they are not difficult to obtain for
many typical distributions. 

For example, consider that f_j(·|θ_j) is an Gaussian distribution
with mean μ_j and precision τ_j. If the prior on τ_j is a Gamma
distribution with parameters a and b, and the prior on μ_j|τ_j is a
Gaussian distribution with mean μ and precision τ_j, then the
full conditionals are given by


    π(μ_j|τ_j,𝐳, 𝐲) ∝ N(μ+s_j/1+n_j, 1/τ_j(1 + n_j))



    π(τ_j|μ_j, 𝐳, 𝐲) ∝ Ga(a + n_j+1/2,
      b + 0.5(μ_j - μ)^2 + 0.5 ss_j)

In the previous expressions, n_j is the number of observations assigned to
component j, s_j the sum of their values and ss_j is the sum of squares
of the values around mean μ_j.

For a mixture of Poisson with f_j(·|θ_j) a Poisson with mean
μ_j and a Gamma prior (i.e, π(μ_j) ∼ Ga(a, b)), the full
conditional distribution is 


    π(μ_j|𝐳, 𝐲) ∝ Ga(a + s_j, b + n_j)

Here, n_j is the number of observations assigned to component j and s_j
the sum of their values.

The full conditionals of θ_∙ presented before for Gaussian and
Poisson mixtures can be derived in a closed form because of the conjugate
priors on the elements in θ_j. As we shall see in
Section <ref>, it is not always necessary to use conjugate
priors as the full conditionals can be obtained using numerical approximations
with INLA. Before describing how to fit mixture models with INLA, we describe
the INLA methodology in the next Section.  Furthermore, collapsed Gibbs
sampling <cit.> can be used if the remainder of the parameters in
θ, denoted by θ_-∙, are integrated out in the full
conditional of θ_∙.

Figure <ref> shows the steps to run the collapsed Gibbs sampling
algorithm using the full conditionals derived before.  After a suitable burn-in
time, Gibbs sampling will provide (correlated) draws from the joint posterior
distribution of (𝐰, 𝐳, θ). Posterior inference
will be based on these samples, which can be used to obtain estimates of the
posterior marginal distributions, summary statistics, etc.






§ INTEGRATED NESTED LAPLACE APPROXIMATION



<cit.> describe a novel approach for Bayesian model fitting that
focuses on obtaining good approximations to the posterior marginals of the
model parameters.  In addition, they only consider models that can be expressed
as a latent Gaussian Markov random field (GMRF) because of their important
computational properties <cit.>.  The
approximation to the posterior marginals is based on repeated use of the
Laplace approximation and, hence, it has been termed Integrated Nested Laplace
Approximation (INLA).

This method can be summarized as follows.  We will consider a vector of
observations 𝐲 =(y_1,…,y_n) whose elements have distributions
family which depend on a vector of parameters θ_1. The means μ_i
of these distributions are conveniently linked to a linear predictor that
depends on a number of latent effects 𝐱, which is a GMRF that
depends on hyperparameters θ_2.  The ensemble of
hyperparameters will be denoted by θ = (θ_1, θ_2).  Elements in 𝐱 include the linear predictor that is
linked to the mean μ_i, as well as the coefficients of the fixed effects
and other types of latent effects. 

INLA assumes that the observations are independent given the latent effects
and the hyperparameters:


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

In the previous equation, ℐ is a index of observed responses and
x_i are the elements of 𝐱 that represent the linear predictor.
The likelihood associated to observation i is defined by index l(i). This
allows for several observations to have the same likelihood.
Hence, the joint
posterior distribution can be written as follows:


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


Given that 𝐱 is a  GMRF, equation (<ref>) can
be rewritten as


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

Here, 𝐐(θ) is the precision matrix of 𝐱.

INLA starts by finding a good approximation to π(θ|𝐲),
π̃(θ|𝐲), that can be used to approximate
 π(x_j|𝐲) because it can be written down as:


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

The approximation provided by INLA is 

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

In the previous equation, π̃(x_j|θ_g, 𝐲) is an approximation to
π (x_j|θ_g, 𝐲). <cit.> describe
different ways in which this approximation can be obtained. θ_g
is a vector with values of θ over a grid, with associated weights
Δ_g.


INLA is implemented as a package for the R statistical software <cit.>
called R-INLA. This package implements a number of latent effects and allows for
an easy model fitting and visualization of the output. A recent review on
INLA and the R-INLA package can be found in <cit.>.


As stated before, INLA (and R-INLA) can handle different likelihoods.  We will
be using this feature to define and fit mixture models with INLA (see
Section <ref>). Furthermore, an approximation to the marginal
likelihood is also provided, so that model selection
<cit.> and Bayesian model averaging
<cit.> can be carried out. This approximation is based on 


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

<cit.> and <cit.> show that this
approximation is very accurate on a wide range of models.




§ FITTING MIXTURE MODELS WITH INLA




We will start by noting that, given 𝐳, INLA can be used to fit the
resulting model because all observations are assigned to a particular group in
the mixture. Hence, given 𝐳, the model can be expressed as a model
with several likelihoods. In this case, the approximations provided by INLA are
for (conditional) posterior marginals π(θ_∙|𝐲,
𝐳), with θ_∙∈θ. In addition, the
(conditional) marginal likelihood obtained is π(𝐲|𝐳).

Hence,
the posterior marginals of the
parameters in the mixture model can be obtained as follows: 


    π(θ_∙|𝐲) = ∑_z∈𝒵π(θ_∙|𝐲, 𝐳 = z)
    π(𝐳 = z|𝐲), θ_∙∈θ

Here, 𝒵 is the parameter space of the of the auxiliary variables,
which is the n-dimensional Cartesian product of {1,…, K}.

If an approximation π̂(𝐳|𝐲) to
π(𝐳|𝐲) is available, the posterior
marginal of θ_∙ can be approximated as


    π(θ_∙|𝐲) ≃∑_z∈𝒵π̃(θ_∙|𝐲, 𝐳 = z)
    π̂(𝐳 = z|𝐲)


Approximation π̂(𝐳|𝐲) can be obtained in a number of
ways. In the next Section, it will be described how Gibbs sampling can be used
together with INLA to obtain the posterior of 𝐳.






















 §.§ Gibbs sampling with INLA


Gibbs sampling can be implemented using the approximations to the conditional
marginals of θ_∙ provided by INLA. This algorithm is described
below and it will follow the lines described in Section <ref>,
with the difference that  some of the steps will be done with some of the
output provided by INLA. 

First of all, it is worth noting that instead of providing an 
approximation to full
conditional of θ_∙, π(θ_∙|𝐲,
𝐳, 𝐰, θ_-∙), INLA provides an
approximation to π(θ_∙|𝐲,
𝐳, 𝐰). However, using the approximation to the former
distribution can be regarded as using a 
collapsed Gibbs sampling <cit.> because π(θ_∙|𝐲,
𝐳, 𝐰) is obtained by integrating θ_-∙ out:


    π(θ_∙|𝐲, 𝐳, 𝐰) = 
    ∫π(θ|𝐲,
    𝐳, 𝐰) d θ_-∙











This simplifies model fitting as the approximation is computed by INLA,
which also allows for a wide range of distributions in the mixture components
and non-conjugate priors when defining the models.
This is summarized in Figure <ref>.










 §.§ 'Modal' Gibbs sampling with INLA



As stated previously, in order to fit a mixture model with INLA only the
posterior of 𝐳 is required.  We now propose a new sampling algorithm
which can be regarded as a collapsed Gibbs sampling <cit.> in the
sense that some of the parameters are integrated out in the full conditionals
and not all parameters in the model are sampled.  During Gibbs sampling
only new values of 𝐳 will be sampled given the (conditional) modes
of 𝐰 and θ, i.e., samples from 𝐰 and θ are replaced by the modes of their respective full conditional
distributions.

The full conditional distribution of 𝐰 is a Dirichlet distribution
with parameters (α_1 + n_1, …, α_K + n_K), which has a mode at ((α_1 +
n_1 - 1)/n^*, …, (α_K + n_K)/n^*), with n^* =n + ∑_i=1^K α_i
- K. For this reason, to make sure that the model is well defined it is better
  to take α_1=…=α_K equal to 2 (for example) instead of 1 in order to propose a vague
prior. The modes of the conditional distributions of θ used in the
Gibbs sampling algorithm are provided by numerical approximation of INLA and
these can be directly used.

Hence, the sampling process can be simplified as in Figure <ref>.






 §.§ Accuracy of 'model' Gibbs sampling



The critical point in this approach is the reliability of 'modal' Gibbs
sampling to obtain a good approximation to π(𝐳 | 𝐲).
This essentially means that all assignments with non-negligible probability are
explored  and that their probabilities are estimated with accuracy.

Let us denote by 𝒵^* ⊆𝒵 the subset of all
possible assignments explored by 'modal' Gibbs sampling. Because of the
structure of 'modal' Gibbs sampling in which different assignments are
simulated by conditioning on the modes of the parameters of the distributions
of the components, we believe that the parameter space 𝒵 is
conveniently explored as this mimics the data generating
process. All the assignments left out are likely to have a very small
probability of occurring under the data generating process <cit.>.

INLA can be used to assess that parameter space 𝒵 has been
conveniently explored.
Assuming that only assignments in 𝒵^* have a non-negligible
probability, the posterior probability of π(𝐳|𝐲)
can be computed as


    π(𝐳|𝐲) ∝π(𝐲|𝐳)π(𝐳); 𝐳∈𝒵^*

Here, π(𝐲|𝐳) is the marginal likelihood, for which INLA
provides a good approximation. Hence, the posterior probabilities can be
computed as


    π(𝐳|𝐲) ≃π̃_I(𝐳|𝐲) =
    π̃(𝐲|𝐳)π(𝐳)/∑_z∈𝒵^*π̃(𝐲|𝐳 = z)π(𝐳 = z)


Note that this approximation is easy and fast to compute from the conditional
models fitted during 'modal' Gibbs sampling. Hence, this can be computed and
compared to the probabilities obtained with the samples from 'modal' Gibbs
sampling, π̂_G(𝐳|𝐲). If both estimates of the
posterior probabilities do not match then there is reason to suspect that
'modal' Gibbs sampling has not explored the parameter space 𝒵
conveniently. In this case, the posterior probabilities computed as in
equation (<ref>) can be used to approximate the posterior distribution
of 𝐳.

























All the previous algorithms will provide an approximation to the posterior of
𝐳, and a set of conditional posterior marginal distributions
π(θ_∙|𝐳, 𝐲) for the model parameters
θ_∙∈θ. The
actual posterior marginals of the model parameters can be approximated as


    π(θ_∙|𝐲) ≃∑_z∈𝒵π̃(θ_∙|𝐳 = z, 𝐲)
      π̂(𝐳 = z | 𝐲)

where π̂(𝐳 = z | 𝐲) is the approximation to
π(𝐳 = z | 𝐲) provided by Gibbs sampling, modal Gibbs
sampling or the one computed using the approximations to the conditional
marginal likelihoods provided by INLA, as in equation (<ref>) above.

The full conditional distribution of 𝐰 is the same as in
Section <ref>. i.e., a Dirichlet distribution. Similarly as before,
the posterior distribution of 𝐰 can be expressed as


    π(𝐰|𝐲) = ∑_z∈𝒵π(𝐰|𝐳 = z, 𝐲)
      π(𝐳 = z| 𝐲) ,

which can be approximated as


    π(𝐰|𝐲) ≃∑_z∈𝒵π(𝐰|𝐳 = z, 𝐲)
      π̂(𝐳 = z| 𝐲)

Note that the previous distribution is a mixture of Dirichlet distributions,
so it can be computed very efficiently.

<cit.> have proposed a similar algorithm to assign
documents to a mixture of topics. However, in their case they are able to
derive the full conditional of z_i on z_-i and 𝐲, so there is
no need to deal with parameters θ during Gibbs sampling.  Once they
have obtained the posterior distribution on the auxiliary variables
𝐳 they are able to compute the posterior distribution on the
remainder of the model parameters by averaging over 𝐳 using its
posterior distribution.

<cit.> follow a similar approach by they develop a fast
algorithm by arguing that, for a given document, posterior probabilities will
be concentrated on a number of topics. Hence, they first explore the set of 
topics fast and then refine the estimated posterior probabilities.

















Finally, the marginal likelihood of the model can also be expressed in a
similar way:


    π(𝐲) = ∑_z∈𝒵π(𝐲|𝐳 = z) π(𝐳 = z).

Given that π(𝐲|𝐳) is the conditional marginal
likelihood, which can be approximated with INLA, and π(𝐳) is
the prior distribution of 𝐳, so the marginal likelihood is easy to
compute. Hence, it can be approximated as


    π(𝐲) ≃π̃_I(𝐲) = ∑_z∈𝒵π̃(𝐲|𝐳 = z) π(𝐳 = z).



Note that  π(𝐳)  is <cit.>:


    π(𝐳) = ∫_𝒲π(𝐳|𝐰) π(𝐰) d𝐰 =
    Γ(∑_i=1^K α_i)/∏_i=1^K Γ(α_i)∏_i=1^K Γ(n_i + α_i)/Γ(n + ∑_i=1^K α_i)

Here, Γ(·) is the Gamma function and 𝒲 is the parametric
space of the weights vector 𝐰.




§ SELECTING THE NUMBER OF COMPONENTS



Finding the number of components is an important problem in mixture models. From
a Bayesian perspective, this implies computing the posterior probabilities of
the number of components  over a range of values. The problem of choosing
the number of components in the mixture can be 
regarded as a model selection problem as well, as discussed below. 

Given models ℳ_1, …,ℳ_p, we can identify each model
by its number of components (i.e., ℳ_1 has one component,
ℳ_2 has two components, and son on) and assign them prior
probabilities π(ℳ_k), k=1,…,p. Posterior probabilities
can be derived as:


    π(ℳ_k|𝐲) ∝π(𝐲|ℳ_k)π(ℳ_k)

Here, π(𝐲|ℳ_k) is the marginal likelihood of model
ℳ_k and the actual posterior probabilities can be computed as
follows:


    π(ℳ_k|𝐲) = π(𝐲|ℳ_k)π(ℳ_k)/∑_k=1^p π(𝐲|ℳ_k)π(ℳ_k)


The marginal likelihood π(𝐲|ℳ_k) of a given model
ℳ_k can be approximated with INLA when the model is completely
fitted with it. However, in the case of mixture models, this is not possible
because INLA is combined with MCMC. For this reason, we need to resort to other
methods to compute the marginal likelihood.

<cit.> and <cit.> describe several approaches to
compute the marginal likelihood of a model using the MCMC output. They note
that the marginal likelihood is the scaling constant when writing the posterior
distribution using Bayes' rule. For mixture models, this can be written
by conditioning on 𝐳 as 


    π(𝐲) = π(𝐲 | 𝐳)π(𝐳)/π(𝐳 | 𝐲)

Note that this holds for any 𝐳 but, as pointed out by
<cit.> and <cit.>, values with a high posterior
probability (or density, if 𝐳 is a continuous random variable) are
preferred since π(𝐳 | 𝐲) will be away from zero. In this
case, posterior mode 𝐳^𝐦 will be used to compute the
log-marginal likelihood as:


    log(π(𝐲)) = log(π(𝐲 | 𝐳^𝐦)) + log(π(𝐳^𝐦))
      - log(π(𝐳^𝐦 | 𝐲))

Note that the first and second terms in the right hand side are the conditional
marginal likelihood (for which an accurate approximation is provided by INLA)
and the prior evaluated at the posterior mode of 𝐳, respectively.
The last term needs to be obtained from the MCMC output. In practice, the
following approximation will be used:


    log(π(𝐲)) ≃log(π̃(𝐲 | 𝐳^𝐦)) + log(π(𝐳^𝐦))
      - log(π̂(𝐳^𝐦 | 𝐲))


Note that the last term in equation (<ref>) can be approximated in
different ways, that will lead to different approximations to
π(𝐲). If it is the estimate provided by modal Gibbs sampling then
the approximation to the marginal likelihood will be denoted by
π̂_G(𝐲).  If the last term is computed using the conditional
marginal likelihoods provided by INLA, as in equation (<ref>), the
approximation will be denoted by π̃_M(𝐲).

In addition to these two ways of approximating the marginal likelihood,
equation (<ref>) provides another way to estimate the marginal
likelihood, denoted by π̃_I(𝐲). Computing these three
different estimates can be a way to perform an assessment of whether the
parametric space 𝒵 has been conveniently explored and the posterior
probabilities of its elements estimated with accuracy. In particular,
π̂_G(𝐲) will be very sensitive to this issue and large
differences with π̃_M(𝐲) and π̃_I(𝐲) will
indicate that the posterior probabilities provided by modal Gibbs sampling
may not be accurate.

Estimates of the posterior probability of a given model, as defined in
equation (<ref>), can be obtained with any of the three estimates of
the marginal likelihood and they will be denoted similarly.  For example,
π̂_G(ℳ_i | 𝐲), i=1,…,p are the estimates of
the posterior probabilities of the different models based on estimates
{π̂_G(𝐲|ℳ_k)}_k=1^p of their marginal likelihoods.




§ SIMULATION STUDY



In order to assess the performance of the methods presented to 
fit mixture models with INLA we have conducted a simulation study. The
aim of this study is twofold, as we are interested in estimating both the
actual number of components and the posterior distributions of the model
parameters. Furthermore, in this simulation study we have considered 
mixtures of Gaussian and Poisson distributions.




 §.§ Gaussian data


The first dataset comprises simulated observations from three Gaussian
distributions centered at 0, 5 and 10, respectively, with precision one.  We
have generated 50 observations from each Gaussian distribution to obtain the
final dataset. Note that this will produce a mixture with overlapping
distributions between groups and, hence, provides a good testing framework.
Mixture models with up to 5 components will be fitted to the data to assess
model choice using the marginal likelihood and the estimation of the model
parameters.  A density estimate of the simulated data can be found in
Figure <ref> (left plot).

In principle, we will not assume that all K components have the same
variances and the components will be defined by K Gaussian distributions with
different means and variances.  Furthermore, we will compare the results
obtained with  our current approach to those obtained with MCMC  with the JAGS
software <cit.>. For this, we have used a burn-in of 200 iterations
plus other 10000 iterations, of which we have only have kept one every ten, so
that a final 1000 samples have been used for inference.  The means of the
Gaussian distributions in the mixture have been ordered in increasing order at
every step of the sampling process with JAGS in order to reclass the assignment
to the mixture components to reduce label switching. 



In order to define the full Bayesian model, we have used a Gaussian prior
on the means with zero mean and precision 1/1000. For the precisions,
we have used a Gamma distribution with parameters 0.5 and 0.5. Also,
all models have been assigned the same probability a priori.




Table <ref> shows the results of the simulation study for the
Gaussian data. The marginal likelihoods favor the model with 3 components,
which is the actual number of components in the mixture. All three methods to
approximate the marginal likelihood give very similar values and their
respective posterior probabilities agree to choose the model with 3 components.
Also, this model provides very good estimates of the means and precisions of
the Gaussian distributions in the mixture.  Figure <ref> shows the posterior
marginals of the model parameters obtained with MCMC and INLA (by Bayesian
model averaging on the conditional posterior marginals). In all cases the
agreement is quite high, which confirms that our approach provides similar
estimates to MCMC.





 §.§ Poisson data



A similar simulation study based on a mixture of Poisson distributions has been
conducted. A mixture of three Poisson distributions with means 1, 15 and 45,
respectively, has been used and 50 observations have been simulated from each
distribution.  A histogram of the simulated data can be found in
Figure <ref> (right plot).  Similarly as in the Gaussian case,
mixture models up to 5 components have been fitted to the data. Model fitting
has also been done using MCMC with the JAGS software, for which we have used
200 burn-in iterations, plus 1000 iterations for inference (after keeping one
in 10 iterations from the total 10000 iterations).



Table <ref> summarizes the results obtained when fitting the
different mixture models to the data. Similarly as in the previous example, all
three methods to compute the marginal likelihoods provide very similar values,
which lead to very similar posterior probabilities of the models.  The
preferred mixture model is clearly the one with three components.
Figure <ref> shows the posterior marginal distributions of the
means of the different components for the mixture with three components. Again,
there is a very good agreement between the estimates provided by INLA and MCMC.






 §.§ Final remarks


Although the simulation study only considers Gaussian and Poisson examples, it
shows that the proposed approach based on modal Gibbs sampling to fit mixture
models with INLA is feasible and works well in practice. The estimates of
the posterior marginals of the model parameters obtained with modal Gibbs
sampling match those obtained with MCMC. The selection of the number of
components in the mixture has also been tackled by providing a new way to
estimate the marginal likelihood of a given mixture model. In the simulation
study this criterion selects the model with the same number of components 
used to simulate the data.

However, fitting mixture models is hard in practice <cit.>. In
particular, if the number of components is not specified correctly most
algorithms will struggle.  For example, when the number of components is higher
than three in the simulation study it can be seen how the components with the
highest means are accurately estimated, whilst the components with lower means
are not properly estimated. This can happen because the algorithm converges to
a local optimum. When we fit a model with more than the actual number of
components (three, in our examples) some components have a mean close to zero.
We believe that this is because there is not enough data allocated to 
components with the lowest means and that the prior (centered at zero) has a
stronger effect on the posterior marginal.

When the fitted model has less components than the actual model, then some
groups will be merged together in one component. Usually, this will make
observations coming from 'neighbor' components being grouped together. Hence,
there is a potential identifiability problem because it is not clear which
components will be merged together.  When the number of components in the model
is larger than its actual number, then some of the components may not have
enough observations allocated as to obtain good estimates of the parameters of
that component. 

Identifiability problems are often observed as label switching during model
fitting. However, we have not observed this problem when fitting the models in
the simulation study. We also believe that using conditional modes for the
model parameters during modal Gibbs sampling protects from label switching
because the components are better identified. Although this sampling framework
may also lead to a loss of variability in the samples obtained of the auxiliary
variables, the results in the simulation study suggest that it is not the case
and the parametric space of the auxiliary variables is conveniently explored.



§ EXAMPLES





 §.§ Gaussian mixture models: galaxy data


<cit.> and <cit.> have studied the velocities of 82
galaxies in the Corona Borealis region using mixture models.  They have fitted
mixture models with different components to try to estimate the optimal number
of components in the mixture using different methods.  <cit.>
compute the posterior probability of each model using MCMC methods, whilst
<cit.> focuses on computing the marginal likelihood to derive Bayes
factors for model choice. We will fit mixture models with a different number of
components in order to obtain the parameter estimates and the marginal
likelihood of each model, so that model choice can be performed.

Figure <ref> shows a kernel density estimate of the data.
Three groups seem to appear in the data. Hence, it makes sense to fit a mixture
model to the data. In order to consider several models, 
mixture models with up to 4 components have been fitted to the data.





In this case, all models are a mixture of Gaussian distributions. We have
considered a vague Gaussian prior (with zero mean and precision 0.001) for the
means, and a Gamma (with parameters 0.5 and 0.5) for the precisions.
Furthermore, all precisions have been considered to be equal, i.e., all
components have the same precision <cit.>. For all
models, simulations included 200 burn-in iterations followed by other 1000
iterations for inference, obtained after thinning (by keeping one in ten iterations).





Table <ref> summarizes the results obtained after fitting the
different mixture models. The values of the marginal likelihood with all three
methods are very similar, as well as the posterior probabilities. In addition,
<cit.> report the marginal likelihoods of the model with two and
three components and these values are very similar to the ones shown in
Table <ref>.  The small discrepancies may be related to the choice
of the priors.  Furthermore,  <cit.> report point estimates of
the model parameters for the models with 3 and 4 components, which are very
similar to the ones in Table <ref> as well.  In this case, the
results favor the mixture model with three components. Hence, modal Gibbs
sampling seems to be a valid approach to fit Gaussian mixture models given that
similar results to other published bibliography have been obtained.




 §.§ Non-Gaussian mixture models: earthquake data


The next example will consider distributions of non-Gaussian mixtures. In
particular, the yearly number of major earthquakes (magnitude 7 or greater) from
1990 to 2006 has been considered.  Figure <ref> shows the
actual dataset and a histogram of the yearly counts.  This dataset has been
analyzed in <cit.> where it is suggested that this dataset
could be analyzed using a mixture model given its overdispersion.  Given that
now the response variable represents counts, it makes sense to use a mixture of
Poisson distributions.



For this reason, different mixture models with up to four components have been
fitted to the data.  The prior on the mean of each group is a vague Gaussian
distribution, with zero mean and precision equal to 0.001.  The results are
summarized in Figure <ref>. Modal Gibbs sampling
seems to provide a different estimate to the the marginal
likelihood to the other two methods. This may indicate that posterior
probabilities of the auxiliary variables may be unreliable in this case.
However, all estimates seem to favor the model with two components and
this is the one with the highest posterior probabilities in all cases.

<cit.> report results where the chosen model is the one with
three components. However, they use a parameterization of the model,
as well as different priors on the mean of the groups.
This may be the
reason why model choice differs. Point estimates of the means seem to be very
similar to those reported in <cit.> though.










§ DISCUSSION



This paper introduces a novel approach to fit mixture models based on obtaining
first the posterior distribution of the latent auxiliary variables using a
variation of Gibbs sampling. Then, the conditional (on the auxiliary variables)
posterior marginals of the remainder of the parameters in the mixture model are
obtained with the Integrated Nested Laplace Approximation. Finally, the
posterior marginals of these parameters are obtained by Bayesian model
averaging of their conditional posterior marginals marginals. Because INLA is
used to fit the conditional models, the distributions of the components
can be taken from a wide range of distributions, conjugate priors are not
required and label switching does not seem to be a problem.

This novel approach also provides a simple way to compute the posterior
probabilities of the number of components, by exploiting the marginal
likelihood provided by INLA and the output from sampling. The simulation study
carried out supports that this is a valid inferential framework for mixture
models.  The results in the simulation study and the comparison with MCMC
included in the case studies suggest that the estimates provided are very
accurate.

The Gibbs sampling algorithm employed in model fitting only draws samples from
the auxiliary variables. For the other parameters, instead of sampling from
their full conditional distributions their respective modes are used.  The
modes are available as part of the output provided by INLA at every step of the
Gibbs sampling algorithm when a conditional model on the current value of the
auxiliary variables is fitted. This simplifies the model fitting algorithm as
well as reduces computing time as actual sampling is reduced. Given the results
obtained in the simulation study, this simple approximation to the full
conditional of the model parameters appears to be enough to obtain a good
estimate of the posterior distribution of the indicator variables.  This
algorithm resembles collapsed Gibbs sampling <cit.> because some of
the model parameters are integrated out in the full conditionals.

Furthermore, if modal Gibbs sampling is thought to provide biased estimates of
the posterior probabilities of the auxiliary variables, then the posterior
probabilities can be recomputed with the estimates of the marginal likelihoods
provided by INLA. This is similar to what <cit.> have done
in the context of the classification of documents by topic. They develop a
fast Gibbs sampling algorithm for latent Dirichlet allocation that explores the
set of possible topics very quickly and then refine the posterior probabilities
because they argue that the probability mass will be concentrated in a small
number of topics.


The idea of obtaining estimates of the model parameter by averaging over the
indicator variables has been considered by several authors. See, for example,
<cit.>, for a general overview, and
<cit.>, for an early application of this idea on the
determination of topics covered in documents.  However, instead of averaging
over single values of the model parameters, this new approach performs Bayesian
model averaging on the (conditional) posterior distributions of the model
parameters obtained with INLA. Hence, no samples from the model parameters are
required, which simplifies model fitting. 

Finally, we believe that this approach can be applied to mixture models in a
variety of areas. For example, it can be used to fit cure rate models
<cit.> in survival analysis, which are expressed as a mixture
model with two components. One of the components is defined by a survival
model, which can be modeled with INLA as a typical parametric survival
distribution or a Cox model. Similarly, mixture models could be fitted to
classify areas which follow a shared spatio-temporal risk pattern or which
depart from this shared pattern in disease mapping <cit.>.
Here, spatial and temporal patterns in a given area can be modeled as the mixture of a
shared spatio-temporal pattern and an area-specific one. INLA provides several
latent effects to model spatial and temporal effects, and a mixture model
would fall within the current paradigm to fit mixture models with INLA.










§ ACKNOWLEDGMENTS


This work has been supported by grant PPIC-2014-001-P, 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. We are also grateful to
Luigi Spezia, assisted by Hector Perpiñán, for the excellent course on
Bayesian inference in hidden Markov models delivered at the University of
Valencia (Spain) in december 2016. Prof. Håvard Rue, Prof. Carmen Armero and
Elena Lázaro provided helpful comments on early versions of this paper.

chicago


