colorlinks=true, citecolor=blue, urlcolor=blue, linkcolor=blue,
filecolor=blue, backref=page

# Gaussian Parsimonious Clustering Models
with Covariates

Keefe Murphy1,2 Thomas Brendan Murphy1,2
mailto:me@keefe.murphy<EMAIL_ADDRESS>
mailto:me@brendan.murphy<EMAIL_ADDRESS>
###### Abstract

We consider model-based clustering methods for continuous, correlated data
that account for external information available in the presence of mixed-type
fixed covariates by proposing the MoEClust suite of models. These allow
covariates influence the component weights and/or component densities by
modelling the parameters of the mixture as functions of the covariates. A
familiar range of constrained eigen-decomposition parameterisations of the
component covariance matrices are also accommodated. This paper thus addresses
the equivalent aims of including covariates in Gaussian Parsimonious
Clustering Models and incorporating parsimonious covariance structures into
the Gaussian mixture of experts framework. The MoEClust models demonstrate
significant improvement from both perspectives in applications to univariate
and multivariate data sets.

Keywords: Model-based clustering, mixtures of experts, EM algorithm,
parsimony, multivariate response, covariates

##  1 Introduction

In many analyses using the standard mixture model framework, a clustering
method is typically implemented on the outcome variables only. Reference is
not made to the associated covariates until the structure of the produced
clustering is investigated in light of the information present in the
covariates. Therefore, interpretations of the values of the model parameters
within each component are guided by covariates that were not actually used in
the construction of the clusters. It is desirable to have covariates
incorporated into the clustering process and not only into the interpretation
of the clustering structure and model parameters, thereby making them
endogenous rather than exogenous to the clustering model. This both exploits
clustering capabilities and provides richer insight into the type of
observation which characterises each cluster.

When each observation consists of a response variable $y_{i}$ on which the
clustering is based and covariates $x_{i}$, there are, broadly speaking, two
main approaches in the literature to having covariates guide construction of
the clusters, neatly summarised by Lamont et al. (2016) and compared in
Ingrassia et al. (2012). Letting $z_{i}$ denote the latent cluster membership
indicator vector, where $z_{ig}=1$ if observation $i$ belongs to cluster $g$
and $z_{ig}=0$ otherwise, the first approach assumes that $z_{i}$ affects the
distribution of $x_{i}$. In probabilistic terms, this means to replace the
actual group-specific conditional distribution
$f\left(y_{i}\,|\,x_{i},z_{ig}=1\right)\mathrm{P}\left(z_{ig}=1\right)~{}\mbox{with}~{}f\left(y_{i}|x_{i},z_{ig}=1\right)f\left(x_{i}\,|\,z_{ig}=1\right)\mathrm{P}\left(z_{ig}=1\right)$.
The name ‘cluster-weighted model’ (CWM) is frequently given to this approach,
e.g. Dang et al. (2017) and Ingrassia et al. (2015); the latter provides a
recent extension allowing for mixed-type covariates. The use of the
alternative term ‘mixtures of regressions with _random_ covariates’ to
describe CWMs (e.g. Hennig (2000)), provides opportunity to clarify that the
remainder of this paper focuses on the second approach, with _fixed_
covariates affecting cluster membership via
$f\left(y_{i}\,|\,x_{i},z_{ig}=1\right)\mathrm{P}\left(z_{ig}=1\,|\,x_{i}\right)$.

This is achieved using the mixture of experts (MoE) paradigm (Jacobs et al.,
1991), in which the parameters of the mixture are modelled as functions of
fixed, potentially mixed-type covariates. However, for multivariate,
continuous, correlated responses, a unifying framework combining the special
cases of the Gaussian MoE model with the flexibility afforded by the various
parsimonious covariance parameterisations within the Gaussian Parsimonious
Clustering Models (GPCM) family of finite mixture models (Banfield & Raftery,
1993; Celeux & Govaert, 1995) has to date been lacking. Indeed, the main
contribution of this paper is in addressing the aim of incorporating
covariates into the GPCM family and the equivalent aim of incorporating GPCM
covariance constraints into the Gaussian MoE framework, by proposing the
MoEClust family of models: the name comes from the interest in employing MoE
models chiefly for clustering purposes. From both perspectives, MoEClust
models demonstrate significant improvement in applications to univariate and
multivariate response data.

A software implementation for the full suite of MoEClust models is provided by
the associated R package MoEClust (Murphy & Murphy, 2018), which is available
from www.r-project.orgwww.r-project.org (R Core Team, 2018), with which all
results were obtained. The syntax of the popular mclust package (Scrucca et
al., 2016) is closely mimicked, with formula interfaces for specifying
covariates in the gating and/or expert networks.

The structure of the paper is as follows; for both Gaussian mixtures of
experts and MoEClust models, the modelling frameworks and inferential
procedures are described, respectively, in Section 2 and Section 3. Section 4
discusses some practical issues affecting performance, namely the
initialisation of the EM algorithm used to fit the models (Section 4.1), and
issues around model selection (Section 4.2). The performance of the proposed
models is illustrated in Section 5 with applications to CO2 emissions data
(Section 5.1) and data from the Australian Institute of Sports (Section 5.2).
Finally, the paper concludes with a brief discussion in Section 6, with some
additional results deferred to the Appendices.

##  2 Modelling

This section builds up the MoEClust models by first describing the mixture of
experts (MoE) modelling framework in Section 2.1 – elaborating on the four
special cases of the MoE model in Section 2.1.1 – and then extending to the
family of MoEClust models comprising Gaussian mixture of experts models with
parsimonious covariance structures from the GPCM family in Sections 2.2 and
2.2.1. Finally, a brief review of existing models and software is given in
Section 2.3.

###  2.1 Mixtures of Experts

The mixture of experts model (Jacobs et al., 1991) extends the mixture model
used to cluster response data $y_{i}$ by allowing the parameters of the model
for observation $i$ to depend on covariates $x_{i}$, where
$y_{1},\ldots,y_{n}$ is an independent sample of response/outcome variables of
dimension $p$, modelled by a $G$ component finite mixture model where the
model parameters depend on associated covariate inputs $x_{1},\ldots,x_{n}$ of
dimension $d$. The MoE model is often referred to as a conditional mixture
model (Bishop, 2006) because, given the set of covariates $x_{i}$, the
distribution of the response variable $y_{i}$ is a finite mixture model:

$f\left(y_{i}\,|\,x_{i}\right)=\sum_{g=1}^{G}\tau_{g}\left(x_{i}\right)f\left(y_{i}\,|\,\theta_{g}\left(x_{i}\right)\right)$

where each component is modelled by a probability density function
$f\left(y_{i}\,|\,\theta_{g}\left(x_{i}\right)\right)$, with component-
specific parameters $\theta_{g}\left(x_{i}\right)$, and mixing proportions
$\tau_{g}\left(x_{i}\right)$. As usual, $\tau_{g}\left(x_{i}\right)>0$ and
$\sum_{g=1}^{G}\tau_{g}\left(x_{i}\right)=1$.

This MoE framework facilitates flexible modelling. While the response variable
$y_{i}$ is modelled via a finite mixture, model parameters are modelled as
functions of other, related, covariates $x_{i}$ from the context under study.
Thus, both the mixing proportions and the parameters of component densities
can depend on the covariates $x_{i}$. It is worth noting that the terminology
used to describe MoE models in the machine learning literature often refers to
the component densities $f\left(y_{i}\,|\,\theta_{g}\left(x_{i}\right)\right)$
as ‘experts’ or the ‘expert network’, and to the mixing proportions
$\tau_{g}\left(x_{i}\right)$ as ‘gates’ or the ‘gating network’, hence the
nomenclature _mixture of experts_.

In the original formulation of the MoE model (Jacobs et al., 1991), the mixing
proportions (gating network) are modeled using multinomial logistic regression
(MLR), though this need not strictly be the case; Geweke & Keane (2007) impose
a multinomial probit structure here instead. The mixture components (expert
networks) are generalised linear models (GLM) (McCullagh & Nelder, 1983).
Thus:

$\hat{\tau}_{g}\left(x_{i}\right)=\frac{\exp\left(\hat{\beta}_{g}^{\top}\tilde{x}_{i}\right)}{\sum_{h=1}^{G}\exp\left(\hat{\beta}_{h}^{\top}\tilde{x}_{i}\right)}\qquad\mbox{and}\qquad\hat{\theta}_{g}\left(x_{i}\right)=\left\\{\psi\left(\hat{\gamma}_{g}^{\top}\tilde{x}_{i}\right),\hat{\Sigma}_{g}\right\\}\vspace{-0.25em}$
(1)

for some link function $\psi\left(\cdot\right)$, a collection of parameters
$\big{\\{}\hat{\gamma}_{g},\hat{\Sigma}_{g}\big{\\}}$ in the component
densities, a collection of parameters $\hat{\beta}_{g}$ in the gates, and
$\tilde{x}_{i}=\left(1,x_{i}\right)^{\top}$. Note that expert network
covariates influence only the component means, and not the component
covariance matrices $\Sigma_{g}$. Henceforth, we restrict our attention to
continuous outcome variables as per the GPCM family. Therefore, component
densities are assumed to be the $p$-dimensional multivariate Gaussian
$\phi\left(\cdot\right)$, and the link function $\psi\left(\cdot\right)$ is
simply the identity, such that covariates are linearly related to the response
variables, i.e.:

$\displaystyle f\left(y_{i}\,|\,x_{i}\right)$
$\displaystyle=\sum_{g=1}^{G}\tau_{g}\left(x_{i}\right)\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}\right)=\left\\{\gamma_{g}^{\top}\tilde{x}_{i},\Sigma_{g}\right\\}\right)$
(2)

####  2.1.1 The MoE Family of Models

It’s possible that some, none, or all model parameters depend on the
covariates. This leads to the four special cases of the Gaussian MoE model
shown in Figure 1, with the following interpretations, due to Gormley & Murphy
(2011):

  1. (a)

in the _mixture model_ the distribution of $y_{i}$ depends on the latent
cluster membership variable $z_{i}$, the distribution of $z_{i}$ is
independent of the covariates $x_{i}$, and $y_{i}$ is independent of the
covariates $x_{i}$ conditional on $z_{i}$:

$f\left(y_{i}\right)=\sum_{g=1}^{G}\tau_{g}\phi\left(y_{i}\bigm{|}\theta_{g}=\left\\{\mu_{g},\Sigma_{g}\right\\}\right)$

  2. (b)

in the _expert network mixture of experts model_ the distribution of $y_{i}$
depends on both the covariates $x_{i}$ and the latent cluster membership
variable $z_{i}$ and the distribution of the latent variable is independent of
the covariates $x_{i}$:

$f\left(y_{i}\,|\,x_{i}\right)=\sum_{g=1}^{G}\tau_{g}\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}\right)=\left\\{\gamma_{g}^{\top}\tilde{x}_{i},\Sigma_{g}\right\\}\right)$

  3. (c)

in the _gating network mixture of experts model_ the distribution of $y_{i}$
depends on the latent cluster membership variable $z_{i}$, the distribution of
the latent variable depends on the covariates $x_{i}$, and the distribution of
the outcome variable is independent of $x_{i}$ conditional on $z_{i}$:

$f\left(y_{i}\right)=\sum_{g=1}^{G}\tau_{g}\left(x_{i}\right)\phi\left(y_{i}\bigm{|}\theta_{g}=\left\\{\mu_{g},\Sigma_{g}\right\\}\right)$

  4. (d)

in the _full mixture of experts model_ , given by (2), the distribution of
$y_{i}$ depends on both the covariates $x_{i}$ and on the latent cluster
membership variable $z_{i}$, and the distribution of the latent variable
$z_{i}$ depends on the covariates $x_{i}$.

For models (c) and (d), $z_{i}$ has a multinomial distribution with a single
trial and probabilities equal to $\tau_{g}\left(x_{i}\right)$. The full MoE
model thus has the following latent variable representation:
$\left(y_{i}\,|\,x_{i},z_{ig}=1\right)\sim\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}\right)=\left\\{\gamma_{g}^{\top}\tilde{x}_{i},\Sigma_{g}\right\\}\right),~{}\mathrm{P}\left(z_{ig}=1\,|\,x_{i}\right)=\tau_{g}\left(x_{i}\right)$.

$x$$z$$y$$\tau$$\theta=\left\\{\mu,\Sigma\right\\}$ (a) Mixture model.

$x$$z$$y$$\tau$$\theta=\left\\{\gamma,\Sigma\right\\}$ (b) Expert network MoE
model.

$x$$z$$y$$\beta$$\theta=\left\\{\mu,\Sigma\right\\}$ (c) Gating network MoE
model.

$x$$z$$y$$\beta$$\theta=\left\\{\gamma,\Sigma\right\\}$ (d) Full MoE model.

Figure 1: The graphical model representation of the mixture of experts models.
The differences between the special cases are due to the presence or absence
of edges between the covariates $x$ and the latent variable $z$ and/or
response variable $y$.

At this point, it’s worth reflecting back on CWMs by noting that they most
fundamentally differ from MoE models in their handling of the mixing
proportions $\tau_{g}$ and in how the joint density
$f\left(x_{i},z_{ig}=1\right)$ is treated, either as
$\mathrm{P}\left(z_{ig}=1\,|\,x_{i}\right)=\tau_{g}\left(x_{i}\right)$ (MoE)
or $f\left(x_{i}\,|\,z_{ig}=1\right)\mathrm{P}\left(z_{ig}=1\right)$ (CWM). In
other words, the direction of the edge between $x$ and $z$ in the full MoE
model in Figure 1(1(d)) is reversed under CWMs (Ingrassia et al., 2012). By
virtue of modelling the distribution of the covariates, CWMs are also
inherently less parsimonious. In principle, the same covariate(s) can enter
both parts of full MoE models. Such a model can provide a useful estimation of
the conditional density of the outcome given the covariates, but the
interpretation of the clustering model and the effect of the covariates
becomes more difficult in this case. Conversely, allowing different covariates
enter different parts of the model further differentiates MoE models from
CWMs. It is common to distinguish among the overall set of covariates between
_concomitant_ gating network variables and _explanatory_ expert network
variables. Thus, for clarity, $x_{i}^{\left(G\right)}$ and
$x_{i}^{\left(E\right)}$ will henceforth respectively refer to the possibly
overlapping subsets of gating and expert network covariates. Higher order
terms and interaction effects between covariates are also allowed in both
parts of the model.

###  2.2 Parsimonious Model-Based Clustering

Parsimonious covariance matrix parameterisations are obtained in GPCMs – of
which mclust is perhaps the most well known implementation – by means of an
eigen-decomposition of the form
$\Sigma_{g}=\lambda_{g}\mathcal{D}_{g}\mathcal{A}_{g}\mathcal{D}_{g}^{\top}$,
where $\lambda_{g}$ is a scalar controlling the volume, $\mathcal{A}_{g}$ is a
diagonal matrix specifying the shape of the density contours, with
$\det(\mathcal{A}_{g})=1$, and $\mathcal{D}_{g}$ is an orthogonal matrix
governing the corresponding ellipsoid’s orientation (Banfield & Raftery, 1993;
Celeux & Govaert, 1995). Imposing constraints reduces the number of free
covariance parameters from $Gp\left(p+1\right)/2$ in the VVV model. This is
desirable when $p$ is even moderately large. Table 1 summarises these GPCM
constraints and their geometric characteristics, which are then shown in
Figure 2.

Table 1: Parameterisations of the within-group covariance matrix $\Sigma_{g}$ available under GPCMs. ‘$\boldsymbol{\dagger}$’ indicates availability in the Gaussian mixture of experts context before and after the introduction of the MoEClust family; ‘$\bullet$’ indicates other models available in the MoEClust family. The horizontal line divides models for univariate and multivarate responses. While all models are possible when $G=1$, they are all equivalent to one of the highlighted available models, otherwise missing entries correspond to models which are never available. The other central columns refer to $G>1$ settings. Name | Model | $G=1$ | $n>p$ | $n\leq p$ | Distribution | Volume | Shape | Orientation
---|---|---|---|---|---|---|---|---
E | $\sigma$ | $\boldsymbol{\dagger}$ | $\bullet$ |  | (univariate) | equal |  |
V | $\sigma_{g}$ |  | $\boldsymbol{\dagger}$ |  | (univariate) | variable |  |
EII | $\lambda\mathcal{I}$ | $\boldsymbol{\dagger}$ | $\bullet$ | $\bullet$ | spherical | equal | equal | —
VII | $\lambda_{g}\mathcal{I}$ |  | $\bullet$ | $\bullet$ | spherical | variable | equal | —
EEI | $\lambda\mathcal{A}$ | $\bullet$ | $\bullet$ | $\bullet$ | diagonal | equal | equal | coordinate axes
VEI | $\lambda_{g}\mathcal{A}$ |  | $\bullet$ | $\bullet$ | diagonal | variable | equal | coordinate axes
EVI | $\lambda\mathcal{A}_{g}$ |  | $\bullet$ | $\bullet$ | diagonal | equal | variable | coordinate axes
VVI | $\lambda_{g}\mathcal{A}_{g}$ |  | $\boldsymbol{\dagger}$ | $\boldsymbol{\dagger}$ | diagonal | variable | variable | coordinate axes
EEE | $\lambda\mathcal{D}\mathcal{A}\mathcal{D}^{\top}$ | $\bullet$ | $\bullet$ |  | ellipsoidal | equal | equal | equal
EVE | $\lambda\mathcal{D}\mathcal{A}_{g}\mathcal{D}^{\top}$ |  | $\bullet$ |  | ellipsoidal | equal | variable | equal
VEE | $\lambda_{g}\mathcal{D}\mathcal{A}\mathcal{D}^{\top}$ |  | $\bullet$ |  | ellipsoidal | variable | equal | equal
EEV | $\lambda\mathcal{D}_{g}\mathcal{A}\mathcal{D}_{g}^{\top}$ |  | $\bullet$ |  | ellipsoidal | equal | equal | variable
VEV | $\lambda_{g}\mathcal{D}_{g}\mathcal{A}\mathcal{D}_{g}^{\top}$ |  | $\bullet$ |  | ellipsoidal | variable | equal | variable
EVV | $\lambda\mathcal{D}_{g}\mathcal{A}_{g}\mathcal{D}_{g}^{\top}$ |  | $\bullet$ |  | ellipsoidal | equal | variable | variable
VVE | $\lambda_{g}\mathcal{D}\mathcal{A}_{g}\mathcal{D}^{\top}$ |  | $\bullet$ |  | ellipsoidal | variable | variable | equal
VVV | $\lambda_{g}\mathcal{D}_{g}\mathcal{A}_{g}\mathcal{D}_{g}^{\top}$ |  | $\boldsymbol{\dagger}$ |  | ellipsoidal | variable | variable | variable
Figure 2: Ellipses of isodensity for each of the 14 parsimonious eigen-
decomposition covariance parameterisations for multivariate data in GPCMs,
with three components in two dimensions.

####  2.2.1 The MoEClust Family of Models

Interest lies in bringing parsimonious covariance structures to Gaussian MoE
models:

$f\left(y_{i}\,|\,x_{i}\right)=\sum_{g=1}^{G}\tau_{g}\left(x_{i}^{\left(G\right)}\right)\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}^{\left(E\right)}\right)=\left\\{\gamma_{g}^{\top}\tilde{x}_{i}^{\left(E\right)},\Sigma_{g}\right\\}\right)$

where $\Sigma_{g}$ can follow any of the GPCM constraints outlined in Table 1.
It is equivalent to say that interest lies in incorporating covariate
information into the GPCM model family. Using these covariance
parameterisations, combined with the four special cases of the MoE model in
Figure 1, yields the MoEClust family of models, which are capable of dealing
with correlated responses and offering additional parsimony in the component
densities compared to current implementations of Gaussian MoE models.

###  2.3 Existing Models and Software

A number of tools for fitting MoE models are available in the R programming
environment (R Core Team, 2018). These include flexmix (Grün & Leisch, 2007,
2008), mixtools (Benaglia et al., 2009), and others.

The flexmix package (Grün & Leisch, 2007, 2008) can accommodate the full MoE
model outlined in Section 2.1.1 in the case where $y_{i}$ is univariate,
though only models with unequal variance can be fitted. The user can specify
the form of the GLM and covariates (if any) to be used in the gating and
expert networks, for which the package has a similar interface to the (glm)
functions within R. In the case of multivariate continuous outcome variables,
there is functionality for multivariate Gaussian component distributions.
However, only the VVI and VVV parameterisations are facilitated.

The mixtools package (Benaglia et al., 2009) can also accommodate the full MoE
model. The package allows for nonparametric estimation of the functional form
for the mixing proportions (gating networks) and the component densities
(expert networks), so it offers further flexibility beyond flexmix. However,
the multivariate models in mixtools use the local independence assumption, so
it doesn’t directly offer the facility to model multivariate Gaussian
component densities with non-diagonal covariance. Again, for univariate data,
only models with unequal variance can be fitted.

The mclust package (Scrucca et al., 2016) can of course accommodate the full
range of covariance constraints in Table 1, and is thus an example of existing
software which can fit GPCMs, but only under the standard finite mixture model
framework (model (1(a)) in Figure 1); it doesn’t facilitate dependency on
covariates in any way.

Finally, it should be noted that eigen-decomposition parsimony has been
introduced to the alternative CWM framework, in which all covariates enter the
same part of the model, by Dang et al. (2017), for the multivariate Gaussian
distributions of both the response variables and the covariates, assuming only
continuous covariates. The flexCWM package (Mazza et al., 2018) allows GPCM
covariance structures in the distribution of the continuous covariates only,
though only univariate responses are accommodated. It also allows,
simultaneously or otherwise, covariates of other types, as well as omitting
the distribution for the covariates entirely, leading to non-parsimonious
mixtures of regressions, with or without concomitant variables.

##  3 Inference via EM

To conduct inference for MoEClust models, we focus on maximum likelihood
estimation using the EM algorithm (Dempster et al., 1977). This is outlined
first for MoE models in Section 3.1 and then extended to MoEClust models in
Section 3.2. A simple trick involving expert network residuals assists fitting
when using GPCM constraints.

###  3.1 MoE Inference

For the full mixture of experts model, the likelihood is of the form:

$\mathcal{L}\left(\beta,\gamma,\Sigma\right)=\prod_{i=1}^{n}\sum_{g=1}^{G}\tau_{g}\left(x_{i}^{\left(G\right)}\right)\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}^{\left(E\right)}\right)\right)$

where
$\tau_{g}\big{(}x_{i}^{\left(G\right)}\big{)}~{}\mbox{and}~{}\theta_{g}\big{(}x_{i}^{\left(E\right)}\big{)}$
are as defined by (1). The data are augmented by imputing the latent cluster
membership indicator $z_{i}=\left\\{z_{i1},\ldots,z_{iG}\right\\}$. Thus, the
conditional distribution of $\left(y_{i},z_{i}\,|\,x_{i}\right)$ is of the
form:

$f\left(y_{i},z_{i}\,|\,x_{i}\right)=\prod_{g=1}^{G}\left[\tau_{g}\left(x_{i}^{\left(G\right)}\right)\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}^{\left(E\right)}\right)\right)\right]^{z_{ig}}$

Hence, the complete data likelihood is of the form:

$\mathcal{L}_{c}\left(\beta,\gamma,\Sigma\right)=\prod_{i=1}^{n}\prod_{g=1}^{G}\left[\tau_{g}\left(x_{i}^{\left(G\right)}\right)\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}^{\left(E\right)}\right)\right)\right]^{z_{ig}}$

and the complete data log-likelihood has the form:

$\displaystyle\ell_{c}\left(\beta,\gamma,\Sigma\right)$
$\displaystyle=\sum_{i=1}^{n}\sum_{g=1}^{G}z_{ig}\left[\log\tau_{g}\left(x_{i}^{\left(G\right)}\right)+\log\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}^{\left(E\right)}\right)\right)\right]$
(3)
$\displaystyle=\sum_{i=1}^{n}\sum_{g=1}^{G}z_{ig}\log\tau_{g}\left(x_{i}^{\left(G\right)}\right)+\sum_{i=1}^{n}\sum_{g=1}^{G}z_{ig}\log\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}^{\left(E\right)}\right)\right)$

The iterative EM algorithm for MoE models follows in a similar manner to that
for mixture models. It consists of an E-step (expectation) which replaces for
each observation the missing data $z_{i}$ with their expected values
$\hat{z}_{i}$, followed by a M-step (maximisation), which maximises the
complete data log-likelihood, computed with the estimates
$\hat{Z}=\left(\hat{z}_{1},\ldots,\hat{z}_{n}\right)$, to provide estimates of
the component weight parameters $\hat{\tau}_{g}$ and the component parameters
$\hat{\theta}_{g}$. Aitken’s acceleration criteria is used to assess
convergence (Böhning et al., 1994). Parameter estimates produced on
convergence achieve at least a local maximum of the likelihood function. Upon
convergence, cluster memberships are estimated via the maximum _a posteriori_
(MAP) classification.

The E-step involves computing:

$\hat{z}_{ig}^{\left(t+1\right)}=\mathbb{E}\left(z_{ig}\bigm{|}y_{i},x_{i},\hat{\beta}^{\left(t\right)},\hat{\gamma}^{\left(t\right)},\hat{\Sigma}^{\left(t\right)}\right)=\frac{\hat{\tau}_{g}^{\left(t\right)}\left(x_{i}^{\left(G\right)}\right)\phi\left(y_{i}\bigm{|}\hat{\theta}_{g}^{\left(t\right)}\left(x_{i}^{\left(E\right)}\right)\right)}{\sum_{h=1}^{G}\hat{\tau}_{h}^{\left(t\right)}\left(x_{i}^{\left(G\right)}\right)\phi\left(y_{i}\bigm{|}\hat{\theta}_{h}^{\left(t\right)}\left(x_{i}^{\left(E\right)}\right)\right)}$

where
$\big{\\{}\hat{\beta}^{\left(t\right)},\hat{\gamma}^{\left(t\right)},\hat{\Sigma}^{\left(t\right)}\big{\\}}$
are the estimates of the parameters in the gating and expert networks on the
$t$-th iteration of the EM algorithm.

For the M-step, we notice that the complete data log-likelihood in (3) can be
considered as a separation into the portion due to the gating network and the
portion due to the expert network, which can be maximised separately under the
EM framework:

$\mathbb{E}\left[\ell_{c}\left(\hat{\beta}^{\left(t\right)},\hat{\gamma}^{\left(t\right)},\hat{\Sigma}^{\left(t\right)}\right)\right]=\sum_{i=1}^{n}\sum_{g=1}^{G}\hat{z}^{\left(t+1\right)}_{ig}\log\tau_{g}\left(x_{i}^{\left(G\right)}\right)+\sum_{i=1}^{n}\sum_{g=1}^{G}\hat{z}^{\left(t+1\right)}_{ig}\log\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}^{\left(E\right)}\right)\right)$
(4)

The first term is of the same form as MLR, here written with component 1 as
the baseline reference level, for identifiability reasons:

$\log\frac{\tau_{g}(x_{i}^{\left(G\right)})}{\tau_{1}(x_{i}^{\left(G\right)})}=\log\frac{\mathbb{P}\big{(}z_{ig}=1\big{)}}{\mathbb{P}\big{(}z_{i1}=1\big{)}}=\beta_{g}^{\top}\tilde{x}_{i}^{\left(G\right)}\quad\forall~{}g=\left\\{2,\ldots,G\right\\},~{}\mbox{where}~{}\beta_{1}=\left\\{0,\ldots,0\right\\}^{\top}$

and thus methods for fitting such models can be used to maximise this term and
find estimates of the parameters in the gating network. The second term is of
the same form as fitting $G$ separate weighted multivariate GLMs, and thus
methods for fitting such models can be used to estimate the expert network
parameters. Note that these are multivariate in the sense of a multivariate
outcome $y_{i}$, with associated covariates having dimension $d$ rendering
these regressions possibly also multivariate in terms of the explanatory
variables. Thus, fitting MoE models is straight forward in principle.

###  3.2 MoEClust Inference

Maximising the second term in (4), corresponding to the expert network, gives
rise to the following expression:

$-\frac{1}{2}\left(p\log
2\pi+\sum_{i=1}^{n}\sum_{g=1}^{G}\hat{z}_{ig}^{\left(t+1\right)}\log\lvert\Sigma_{g}\rvert+\sum_{i=1}^{n}\sum_{g=1}^{G}\hat{z}_{ig}^{\left(t+1\right)}\left(y_{i}-\gamma_{g}^{\top}\tilde{x}_{i}^{\left(E\right)}\right)^{\top}\Sigma_{g}^{-1}\left(y_{i}-\gamma_{g}^{\top}\tilde{x}_{i}^{\left(E\right)}\right)\right)$
(5)

from which it can be shown that $\hat{\gamma}_{g}$ does not depend on
$\Sigma_{g}$, much like a Seemingly Unrelated Regression model (SUR) (Zellner,
1962). We first estimate $\hat{\gamma}_{g}$ and then $\hat{\Sigma}_{g}$.
Fitting $G$ separate multivariate regressions (weighted by $\hat{z}_{ig}$),
yields $G$ sets of $n\times p$ SUR residuals
$\hat{r}_{ig}=y_{i}-\hat{\gamma}_{g}^{\top}\tilde{x}_{i}^{\left(E\right)}$.
Crucially, these residuals satisfy $\sum_{i=1}^{n}\hat{z}_{ig}\hat{r}_{ig}=0$.
Thus, maximising (5) is equivalent to minimising the expression below:

$\sum_{i=1}^{n}\sum_{g=1}^{G}\hat{z}_{ig}^{\left(t+1\right)}\log\lvert\Sigma_{g}\rvert+\sum_{i=1}^{n}\sum_{g=1}^{G}\hat{z}_{ig}^{\left(t+1\right)}\hat{r}_{ig}^{\top}\Sigma_{g}^{-1}\hat{r}_{ig}$
(6)

which is of the same form as the criterion used in the M-step of a standard
Gaussian finite mixture model with component covariance matrices $\Sigma_{g}$,
component means equal to zero, and new augmented data set $\hat{R}$. Thus,
when estimating the component covariance matrices via (6), the same M-step
function as used within mclust can be applied to augmented data, constructed
so that each observation is represented as follows:

  1. 1.5em  1. Stack the $G$ sets of SUR residuals into the $\left(n\times G\right)\times p$ matrix $\hat{R}$: $\hat{R}={\normalsize\begin{bmatrix}\hat{r}_{111}&\hat{r}_{121}&\dots&\hat{r}_{1p1}\\\ \hat{r}_{211}&\hat{r}_{221}&\dots&\hat{r}_{2p1}\\\ \vdots&\vdots&\ddots&\vdots\\\ \hat{r}_{n11}&\hat{r}_{n21}&\dots&\hat{r}_{np1}\\\ \hline\cr\hat{r}_{112}&\hat{r}_{122}&\dots&\hat{r}_{1p2}\\\ \hat{r}_{212}&\hat{r}_{222}&\dots&\hat{r}_{2p2}\\\ \vdots&\vdots&\ddots&\vdots\\\ \hat{r}_{n12}&\hat{r}_{n22}&\dots&\hat{r}_{np2}\\\ \hline\cr\vdots&\vdots&\ddots&\vdots\\\ \hline\cr\hat{r}_{11G}&\hat{r}_{12G}&\dots&\hat{r}_{1pG}\\\ \hat{r}_{21G}&\hat{r}_{22G}&\dots&\hat{r}_{2pG}\\\ \vdots&\vdots&\ddots&\vdots\\\ \hat{r}_{n1G}&\hat{r}_{n2G}&\dots&\hat{r}_{npG}\\\ \end{bmatrix}}$ 2. Create the $\left(n\,\times\,G\right)\,\times\,G$ block-diagonal matrix $\hat{\zeta}$ from the columns of $\hat{Z}$: $\hat{\zeta}={\normalsize\begin{bmatrix}\hat{z}_{11}&0&\dots&0\\\ \hat{z}_{21}&0&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ \hat{z}_{n1}&0&\dots&0\\\ \hline\cr 0&\hat{z}_{12}&\dots&0\\\ 0&\hat{z}_{22}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ 0&\hat{z}_{n2}&\dots&0\\\ \hline\cr\vdots&\vdots&\ddots&\vdots\\\ \hline\cr 0&0&\dots&\hat{z}_{1G}\\\ 0&0&\dots&\hat{z}_{2G}\\\ \vdots&\vdots&\ddots&\vdots\\\ 0&0&\dots&\hat{z}_{nG}\\\ \end{bmatrix}}$

Structuring the model in this manner allows GPCM covariance structures to be
imposed on Gaussian MoE models with gating and/or expert network covariates,
hence the nomenclature _MoEClust_. In the end, the M-step involves three sub-
steps, each using the current estimate of $\hat{Z}$: i) estimating the gating
network parameters $\hat{\beta}_{g}$ and hence the component weights
$\hat{\tau}_{g}\big{(}x_{i}^{\left(G\right)}\big{)}$ via MLR, ii) estimating
the expert network parameters $\hat{\gamma}_{g}$ and hence the component-
specific means via weighted multivariate multiple linear regression, and iii)
estimating the component covariance matrices using the augmented data set
comprised of SUR residuals, as outlined above.

In the absence of covariates in the gating and/or expert networks, under the
special cases outlined in Section 2.1.1, their respective contribution to the
expected complete log-likelihood in (4) is maximised as per the corresponding
term in a standard GPCM. In other words, the gating and expert networks,
without covariates, can be seen as regressions with only an intercept term.
Thus, for instance, the formula for estimating the gates in the absence of
concomitant variables is $\hat{\tau}_{g}=n^{-1}\sum_{i=1}^{n}\hat{z}_{ig}$,
rather than (1). It is also possible and sometimes useful to constrain
$\tau_{g}$ to be equal across components, when gating concomitants are
omitted, thereby estimating fewer parameters. Similarly, removing the
corresponding regression intercept(s) from the part(s) of the model where
covariates enter can yield further parsimony in appropriate settings.

The MoEClust R package also facilitates modelling with an additional noise
component to encompass clusters with non-Gaussian distributions, distributed
as a homogeneous spatial Poisson process, via either of the following models:

$\displaystyle f\left(y_{i}\,|\,x_{i}\right)$
$\displaystyle=\sum_{g=1}^{G}\tau_{g}\left(x_{i}^{\left(G\right)}\right)\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}^{\left(E\right)}\right)=\left\\{\gamma_{g}^{\top}\tilde{x}_{i}^{\left(E\right)},\Sigma_{g}\right\\}\right)+\frac{\tau_{0}\left(x_{i}^{\left(G\right)}\right)}{V}$
(7) $\displaystyle f\left(y_{i}\,|\,x_{i}\right)$
$\displaystyle=\sum_{g=1}^{G}\tau_{g}\left(x_{i}^{\left(G\right)}\right)\phi\left(y_{i}\bigm{|}\theta_{g}\left(x_{i}^{\left(E\right)}\right)=\left\\{\gamma_{g}^{\top}\tilde{x}_{i}^{\left(E\right)},\Sigma_{g}\right\\}\right)+\frac{\tau_{0}}{V}$
(8)

as per Banfield & Raftery (1993), with $V$ estimated by the hypervolume of the
convex hull, ellipsoid hull, or smallest hypercube enclosing the data, and the
extension that concomitant variables are allowed to affect (7) or not affect
(8) the noise component’s gates, but such mixtures with a noise component are
not considered further here.

##  4 Practical Issues

In this section, factors affecting the performance of MoEClust models are
discussed; namely, the necessity of a good initial partition to prevent the EM
algorithm from converging to a suboptimal local maximum (Section 4.1), and the
issue of model selection with regard to where and what covariates enter the
model, if any (Section 4.2).

###  4.1 EM Initialisation

With regards to initialisation of the EM algorithm for MoEClust models,
agglomerative hierarchical clustering and quantile-based clustering have been
found to be suitable for multivariate and univariate data, respectively.
Random initialisation of the allocations have also been found to perform well,
despite the obvious computational drawback of the need to run the EM algorithm
from multiple random starting points. However, when explanatory variables
$x_{i}^{\left(E\right)}$ enter the expert network, it is useful to use the
explanatory variables to augment the initialisation strategy with an extra
step, which takes the initial partition of the data (whether obtained by
hierarchical clustering, random initialisation, or some other method) and
iteratively reallocates observations in such a way that each subset can be
well-modelled by a single expert. Algorithm 1 outlines the proposed
initialisation strategy, similar to that of Ning et al. (2008).

1 Start from some non-overlapping partition
$\Omega_{1},\Omega_{2},\ldots,\Omega_{G}$.

2 Estimate the expert network regressor
$\eta_{g}\left(\gamma_{g},\cdot\right)$ on every subset
$\left\\{\Omega_{g}\right\\}_{g=1}^{G}$.

3 For every observation, compute the Mahalanobis distance between $y_{i}$ and
the fitted values
$\hat{y}_{ig}=\eta_{g}\left(\gamma_{g},x_{i}^{\left(E\right)}\right)$:
$M_{ig}=d\left(y_{i},\hat{y}_{ig}\right)=\sqrt{\left(y_{i}-\hat{y}_{ig}\right)^{\top}\hat{\Sigma}_{g}^{-1}\left(y_{i}-\hat{y}_{ig}\right)}$.

4 Let $k_{i}=\arg\,\min_{g}M_{ig}$.

5 Reassign observation $i$ to subset $\Omega_{k_{i}}$.

Repeat Steps 2-5 until convergence is achieved.

Algorithm 1 Iterative reallocation initialisation with expert network
covariates

If at any stage a level is dropped from a categorical variable in subset
$\Omega_{g}$ the variable itself is dropped from the corresponding regressor
for the observations with missing levels. Convergence of the algorithm is
guaranteed and the additional computational burden incurred is negligible. By
using the Mahalanobis distance metric (Mahalanobis, 1936), each observation is
assigned to the cluster corresponding to the Gaussian ellipse to which it is
closest. This has the added advantage of speeding up the running of the EM
algorithm. The estimates of $\eta_{g}\left(\gamma_{g},\cdot\right)$ at
convergence are used as starting values for the expert network. The gating
network is initialised by considering the partition itself at convergence as a
discrete approximation of the gates.

When using a deterministic approach to obtain starting partitions
$\Omega_{1},\ldots,\Omega_{G}$ for Algorithm 1, initialisation can be further
improved by considering information in the expert network covariates to find a
good clustering of the joint distribution of
$\big{(}y_{i},x_{i}^{\left(E\right)}\big{)}$. When $x_{i}^{\left(E\right)}$
includes categorical or ordinal covariates, the model-based approach to
clustering data of mixed type of McParland & Gormley (2016) could be employed
at this stage; however, the present implementation only considers the
continuous covariates within $x_{i}^{\left(E\right)}$ for this step. Figure 3
illustrates the necessity of this procedure using a toy data set, with a
single continuous covariate and a univariate response clearly arising from a
mixture of two linear regressions, which otherwise would not be discerned
without including the covariate in the initialisation routine via Algorithm 1.

Our initialisation strategy has the same limitation that the result may
represent a suboptimal local maximum. However, the problem is transferred from
the difficult task of initialising the EM algorithm to initialising Algorithm
1. Thus, it is feasible to repeat the algorithm with many different partitions
and choose the best result to initialise the EM algorithm, since it converges
very fast, requires much less computational effort than the EM itself, and
generally reduces the number of required EM iterations.

(a) Hierarchical clustering of the
response variables only.

(b) Random allocation.

(c) Partition obtained via Algorithm 1.

Figure 3: Initial 2-component partitions on univariate data clearly arising
from a mixture of two linear regressions, obtained using (3(a)) agglomerative
hierarchical clustering, (3(b)) random allocation, and (3(c)) Algorithm 1
applied to the initialisation in (3(b)), demonstrating the improvement
achieved by incorporating expert network covariates into the initialisation
strategy. Fitted lines are also shown.

###  4.2 Model Selection

Whether a variable should be considered as a covariate or part of the response
is usually clear from the context of the data being clustered. However, within
the suite of MoE models outlined above, it is natural to then question which
covariates, if any, are to be included, and if so in which of the
gating/expert networks. Unless the manner in which covariates enter the MoE
model is guided by the question of interest in the application under study,
this is a challenging problem as the space of MoE models is potentially very
large, once variable selection for the covariates entering the gating and
expert networks is considered. Thus in practice only models where covariates
enter all mixture components or all component weights in a linear manner are
typically considered in order to restrict the size of the model search space.
However, even within this reduced model space, there are $2^{d}+(G-1)\times
2^{2d}$ models to consider. If the number of components $G$ is unknown, the
model search space increases further.

Model comparison for the MoEClust class of models is even more challenging, as
for multivariate data there are potentially 14 different GPCM covariance
parameterisations to consider for models with $G\geq 2$, and 3 otherwise.
However, model selection can still be implemented in a similar manner to other
model-based clustering methods: the Bayesian Information Criterion (BIC)
(Schwarz, 1978) and Integrated Completed Likelihood (ICL) (Biernacki et al.,
2000) have been shown to give suitable model selection criteria, both for the
number of component densities (and thus clusters) required and for selecting
variables to include in the model. Furthermore, Gormley & Murphy (2010)
demonstrates how such criteria can be employed to select the appropriate model
structure, to choose the appropriate number of components, and to guide the
inclusion of covariates across the four MoE special cases shown in Figure 1.

For MoEClust models involving mixtures of GLMs, standard variable selection
approaches can be used to find the optimal covariates for inclusion in either
the multinomial logistic regression (the gating network) or the multivariate
weighted least squares regression (the expert network). However, the thorny
variable selection issue is complicated further when one considers that the
selected variables may only be optimal for the given $G$ and the given set of
GPCM covariance matrix constraints. Even stepwise model selection approaches
are therefore difficult, because any step could involve adding or removing a
single covariate – in either the gating or expert networks – adding or
removing a component, or altering the covariance parameterisation, with
backwards stepwise selection being particularly cumbersome (see Appendix A).

Thus for MoEClust models, the recommended approach closely follows Gormley &
Murphy (2010): fix the covariates for inclusion in the component weights and
component densities and then find the $G$ value and GPCM covariance
parameterisation which together optimise some criterion. Different fits with
different combinations of covariates can then be compared according to the
same criterion. The number of free parameters in the penalty term for the BIC
and ICL criteria of course depends on the included gating and expert network
covariates and the GPCM constraints employed.

##  5 Results

The clustering performance of the MoEClust family of models is illustrated
with application to two well known benchmark data sets: univariate CO2 data
(Section 5.1), and multivariate data from the Australian Institute of Sports
(Section 5.2). Henceforth, any mention of the initialisation of the
allocations, where covariates enter the expert network, refers to the method
used to find an initial partition for Algorithm 1. All results were obtained
using the R package MoEClust (Murphy & Murphy, 2018). BIC was used to
determine the optimal number of components, choose the covariance type, and
select the best subset of covariates, as well as where to put them, in all
cases.

###  5.1 CO2 Data

As a univariate example of an application of MoEClust, we consider data on the
CO2 emissions and Gross National Product (GNP) for $n=28$ countries in the
year 1996 (Hurn et al., 2003). There is interest in studying the relationship
between CO2 and GNP. Thus, it appears that a suitable model for these data
would be a mixture of linear regression models, where the component density is
a Gaussian distribution.

Models with different numbers of components ($G\in\\{1,\ldots,9\\}$) with
either the equal variance (E) or unequal variance (V) models from Table 1 were
considered. Consideration was also given to inclusion/exclusion of GNP as a
covariate in the gating and/or expert networks. When GNP is present as a
gating concomitant, only models with $G\geq 2$ are fitted. Results are
reported based on the best of 20 random initial cluster allocations in each
case. The final selected model had $G=2$, unequal variance, and GNP entered
the expert network; as such, this is an _expert network mixture of experts
model_. This same model maximised both the BIC and ICL criteria.

A table of BIC and ICL values for the top model under each of the four special
cases of the MoE model are given in Table 2. The parameters of the optimal
model are detailed in Table 3 and its fit is exhibited in Figure 4. The first
row of Table 2 is equivalent to a standard finite Gaussian mixture. It’s worth
noting where GNP is included and the E model is optimal that this constraint
on the variance is a feature of MoEClust, which neither MoE nor GPCM models
can presently capture. While the V models can, however, also be fitted using
flexmix and mixtools, our initialisation strategy ultimately leads to higher
log-likelihood estimates in those instances. Table 3 and Figure 4 suggest the
relationship between CO2 and GNP is clustered around two different linear
regression lines; one small, compact cluster of six countries with a large
slope value and a larger, more diffuse one with a slope of almost zero.

(a) Fitted lines of the expert network GLMs, with dashed-line 95% confidence
intervals. Text label size is proportional to a country’s probability of
belonging to the cluster to which it was assigned. Clustering uncertainty is
also indicated by vertical bars relating to the second y-axis.

(b) Heat map of the conditional density of the outcome variable CO2,
accounting for the gating and expert networks, the latter of which includes
GNP as a covariate.

Figure 4: Scatter plots of GNP against CO2 emissions for a number of countries with two linear regression components from the optimal unequal variance MoEClust model. Table 2: The MoEClust BIC and ICL values of the top models under the four MoE special cases for the CO2 data. Each row is optimal with respect to $G$ and GPCM type, given the included covariates. Gating | Expert | G | GPCM | BIC | ICL
---|---|---|---|---|---
|  | 2 | E | -163.16 | -163.91
GNP |  | 2 | E | -166.05 | -166.68
| GNP | 2 | V | -157.20 | -160.04
GNP | GNP | 2 | V | -159.25 | -161.47
Table 3: Estimated parameters of the optimal MoEClust model fit to the CO2 data. Parameter | Component 1 | Component 2
---|---|---
Proportion | 0.75 | 0.25
(Intercept) | 8.68 | 1.41
GNP | -0.02 | 0.68
$\sigma_{g}^{2}$ | 4.20 | 0.66

###  5.2 Australian Institute of Sport (AIS) Data

Various physical and hematological (blood) measurements were made on 102 male
and 100 female athletes at the Australian Institute of Sport (AIS) (Cook &
Weisberg, 1994). The thirteen variables recorded in the study are detailed in
Table 4.

Table 4: Australian Institute of Sports data variables. The five in the first column are hematological response variables and the others are physical measurements for the athlete. Response | Description | Covariate | Description
---|---|---|---
RCC | red cell count | BMI | body mass index
WCC | white cell count | SSF | sum of skin folds
Hc | Hematocrit | Bfat | body fat percentage
Hg | Hemoglobin | LBM | lean body mass
Fe | plasma ferritin concentration | Ht | height
|  | Wt | weight
|  | sex | a factor with levels: female, male
|  | sport | a factor with levels: Basketball, Field,
|  |  | Gym, Netball, Rowing, Swimming, Tennis,
|  |  | Track 400m, Track Sprint, Water Polo

A MoEClust model can be used to investigate the clustering structure in the
hematological measurements of the athletes and investigate how covariates may
influence these measurements and the clusters. A range of MoEClust models were
fitted to this data, with $G\in\\{1,\ldots,9\\}$, with either, neither, or
both of the covariates ‘sex’ and ‘BMI’ allowed to enter either, neither, or
both of the gating and expert networks. Concomitant gating covariates were
only considered for models with $G\geq 2$. With 3 permissible covariance
parameterisations for the single component models, and 14 otherwise, and 16
different possible combinations of gating and/or expert network covariate
settings, this amounts to an exhaustive search over a total of 1,804 models
under consideration. However, not all models converged, particularly for
higher values of $G$. Cluster allocations were initialised using agglomerative
hierarchical clustering. Table 5 gives the BIC and ICL values of a selection
of these fitted models, representing the optimal models according to the BIC
criterion for certain special cases of interest. Similar results for a
stepwise model search are given in Appendix A.

Table 5: The MoEClust BIC and ICL values for a selection of MoEClust models fitted to the Australian Institute of Sports data. Rows 1 and 2 give the optimal model using only the covariance parameterisations implemented previously in flexmix and mixtools. Among the more general MoEClust family, the final row gives the optimal model according to the ICL criterion and the remaining rows give the top models according to the BIC criterion for each of the four mixture of experts special cases, which include covariates in neither, either, or both of the gating and expert networks. As such, row 3 corresponds to the optimal model according to mclust. Rank (BIC) | Gating | Expert | G | GPCM | BIC | ICL
---|---|---|---|---|---|---
46 | BMI | sex | 2 | VVV | -4054.86 | -4113.98
622 | BMI + sex | sex | 4 | VVI | -4302.82 | -4329.00
233 |  |  | 2 | EVE | -4146.16 | -4201.61
2 |  | sex | 2 | EVE | -4015.35 | -4059.55
14 | sex |  | 3 | EVE | -4037.32 | -4066.66
1 | BMI | sex | 2 | EVE | -4013.40 | -4074.10
15 | BMI + sex |  | 3 | EEE | -4038.75 | -4043.01

Clearly, the inclusion of covariates improves the fit compared to GPCM models.
Similarly, the use of GPCM covariance constraints clearly improves the fit
compared to standard Gaussian MoE models. The optimal model according to BIC
is a 2-component EVE _full mixture of experts_ model. It is interesting to
note that different covariates enter in different parts of this model: BMI
enters the gating network and sex enters the expert network. The gating
network had an intercept value of -5.02 and a slope coefficient of 0.21; thus
higher BMI values increase the probability of belonging to Cluster 2 (see
Figure 6). Though every observation has its own mean parameter, in the
presence of expert network covariates, given in this case by the fitted values
of Table 6(a), the means are summarised by the posterior mean of the fitted
values,
$\hat{\mu}_{g}=\frac{\sum_{i=1}^{n}\hat{z}_{ig}\hat{y}_{i}}{\sum_{i=1}^{n}\hat{z}_{ig}}=\frac{\sum_{i=1}^{n}\hat{z}_{ig}\left(\sum_{g=1}^{G}\hat{z}_{ig}\left(\hat{\gamma}^{\top}_{g}\tilde{x}_{i}\right)\right)}{\sum_{i=1}^{n}\hat{z}_{ig}}$
in Table 6(b), along with the average gates. Further details of the model
parameters are given in Table B.1 in Appendix B.

Table 6: Estimated parameters of the optimal MoEClust fit to the AIS data.
6(a) gives the coefficients of the expert network regression for the component
densities, with ‘female’ used as the reference category for the explanatory
variable ‘sex’. While every observation has its own mean parameter, given by
the fitted values of the expert network in 6(a), the means are summarised by
the posterior mean of the fitted values of the model in 6(b), along with the
average gates.

(a) (a) | RCC | WCC | Hc | Hg | Fe
---|---|---|---|---|---
Component 1
(Intercept) | 4.47 | 6.71 | 40.82 | 13.61 | 45.92
sexmale | 0.64 | 0.01 | 5.03 | 1.94 | 15.89
| Component 2
(Intercept) | 4.30 | 7.46 | 39.94 | 13.47 | 74.79
sexmale | 0.64 | 0.29 | 5.51 | 2.08 | 57.20

(b) (b) | Component 1 | Component 2
---|---|---
$\overline{\hat{\tau}}_{g}$ | 0.56 | 0.44
RCC | 4.74 | 4.70
WCC | 6.84 | 7.45
Hc | 43.00 | 43.21
Hg | 14.48 | 14.67
Fe | 59.79 | 98.77

Though the plots in Figure 4 are suitable for univariate data with a single
continuous expert network covariate, visualising MoEClust results for
multivariate data with a mix of continuous and categorical covariates
constitutes a much greater challenge. For the optimal 2-component EVE model
fitted to the AIS data, the data and clustering results are shown using a
generalised pairs plot in Figure 5. The generalised pairs plot depicts the
pairwise relationships between the hematological response variables, the
included covariates, and the MAP classification, coloured according to the MAP
classification. The marginal distributions of each variable are given along
the diagonal. For the hematological response variables, ellipses with axes
corresponding to the within-cluster covariances are drawn. For the purposes of
visualising Figure 5, owing to the presence of a covariate in the expert
network of the model fitted to the AIS data, the MVN ellipses in panels
depicting two response variables are centered on the posterior mean of the
fitted values, as described above. The shape and size of these ellipses are
also modified for the same reason: they are derived by adding the extra
variability in the component means to $\hat{\Sigma}_{g}$.

Figure 5: Generalised pairs plot for the optimal MoEClust fit to the AIS data,
depicting the pairwise relationships between the hematological response
variables and the covariates BMI and sex. Colours and plotting symbols
correspond to the MAP classification. Mosaic plots are used to depict two
categorical variables, scatter plots are for two continuous variables, and a
mix of box-plots and jittered strip-plots are used for mixed pairs.

It’s clear from Figure 5 that the variable ‘plasma ferritin concentration’
(Fe) and the covariate ‘sex’ are driving much of the separation in the
clusters. Interestingly, the correspondence between the MAP cluster assignment
and the sex label is poor, with an Adjusted Rand Index (Hubert & Arabie, 1985)
of just 0.02; this is because BMI rather than sex influences the probability
of cluster membership (see Figure 6), and subjects are divided by sex within
each cluster, rather than the clusters capturing the athletes’ sex. Indeed,
Table 6(a) implies that males, on average, have elevated levels of all five
blood measurements in both clusters, though the magnitude of this effect is
more pronounced in Cluster 2, corresponding to athletes with higher average
BMI.

Figure 6: Gating network of the optimal fit to the AIS data, against the
concomitant variable BMI. Cluster 1 is shown by the downward sloping black
curve, Cluster 2 by the upward sloping red curve.

##  6 Discussion

The development of a suite of MoEClust models has been outlined, clearly
demonstrating the utility of mixture of experts models for parsimonious model-
based clustering where covariates are available. A novel means of visualising
such models has also been proposed. The ability of MoEClust models to jointly
model the response variable(s) and related covariates provides deeper and more
principled insight into the relations between such data in a mixture model
based analysis, and provides a principled method for both creating and
explaining the clustering, with reference to information contained in
covariates. Their demonstrated use to cluster observations and appropriately
capture heterogeneity in cross-sectional data provides only a glimpse of their
potential flexibility and utility in a wide range of settings. Indeed, given
that general MoE models have been used, under different names, in several
fields, including but not limited to statistics (Grün & Leisch, 2007, 2008),
biology (Wang et al., 1996), econometrics (Wang et al., 1998), marketing
(Wedel & Kamakura, 2012), and medicine (Thompson et al., ), MoEClust models
could prove useful in many domains.

Improvement over GPCM models has been introduced by accounting for external
information available in the presence of covariates. Similarly, improvement
over Gaussian mixture of experts models which incorporate fixed covariates has
been introduced by allowing GPCM family covariance structures in the component
densities. Due to sensitivity of the final solution obtained by the EM
algorithm to the initial values, an iterative reallocation procedure based on
the Mahalanobis distance has been proposed to mitigate against convergence to
a suboptimal local maximum. This initialisation algorithm converges quickly
and also speeds up convergence of the EM algorithm itself.

On a cautionary note, care must be exercised in choosing how and where
covariates enter when a MoEClust model is used as an analytic tool. The
interpretation of the analysis fundamentally depends on where covariates
enter, which of the suite of MoE models is invoked, and also on which GPCM
constraints are employed. Gating network MoEClust models may be of particular
interest to users of GPCMs; while concomitants influence the probability of
cluster membership, the correspondance thereafter between component densities
and clusters has the same interpretation as in standard GPCMs. When covariates
enter the component densities, we caution that observations with very
different response values can be clustered together, because they were being
modelled using the same GLM; similarly, regression distributions with distinct
parameters do not necessarily lead to well-separated clusters. Due to the
lower number of parameters in the parsimonious covariance matrices, MoEClust
models tend to favour including explanatory variables in the expert network
more than standard Gaussian MoE models would, as the number of parameters
added to the penalty term for the BIC and ICL by expert network covariates is
made up for by the reduction in the number of covariance parameters. Thus, in
cases where a mixture of experts model might elect to include a concomitant
variable in the gating network, a MoEClust model with fewer covariance
parameters may elect to include it as an explanatory expert network variable
instead; while this does lead to a better fit, it can complicate
interpretation.

Possible directions for future work in this area include exploring the use of
regularisation penalties in the gating and expert networks to help with
variable selection, investigating the utility of nonparametric estimation of
the gates, and introducing parsimonious factor-analytic covariance structures
for high-dimensional data (McNicholas & Murphy, 2008) into Gaussian mixture of
experts models using residuals in an equivalent fashion. Similarly, MoEClust
models could benefit from the heavier tails of the multivariate
$t$-distribution, and the robustness to outliers it affords, by considering
the associated $t$EIGEN family of covariance constraints (Andrews &
McNicholas, 2012). Finally, MoEClust models could also be extended to the
supervised or semi-supervised model-based classification settings, where some
or all of the observations are labelled.

## Acknowledgements

This research emanated from work conducted with the financial support of
Science Foundation Ireland under grant number SFI/12/RC/2289 in the Insight
Centre for Data Analytics in University College Dublin (UCD).

## References

  * Andrews & McNicholas (2012) J. L. Andrews & P. D. McNicholas.  Model-based clustering, classification, and discriminant analysis via mixtures of multivariate $t$-distributions: the $t$EIGEN family.  _Statistics and Computing_ , 22, 2012.
  * Banfield & Raftery (1993) J. Banfield & A. E. Raftery.  Model-based Gaussian and non-Gaussian clustering.  _Biometrics_ , 49: 803–821, 1993.
  * Benaglia et al. (2009) T. Benaglia, D. Chauveau, D. R. Hunter, & D. Young.  mixtools: An R Package for Analyzing Finite Mixture Models.  _Journal of Statistical Software_ , 32(6): 1–29, 2009\.
  * Biernacki et al. (2000) C. Biernacki, G. Celeux, & G. Govaert.  Assessing a mixture model for clustering with the integrated completed likelihood.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 22(7): 719–725, 2000.
  * Bishop (2006) C. M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer-Verlag, New York, 2006.
  * Böhning et al. (1994) D. Böhning, E. Dietz, R. Schaub, P. Schlattmann, & B. G. Lindsay.  The distribution of the likelihood ratio for mixtures of densities from the one-parameter exponential family.  _Annals of the Institute of Statistical Mathematics_ , 46(2): 373–388, 1994.
  * Celeux & Govaert (1995) G. Celeux & G. Govaert.  Gaussian parsimonious clustering models.  _Pattern Recognition_ , 28: 781–793, 1995.
  * Cook & Weisberg (1994) R. D. Cook & S. Weisberg.  _An Introduction to Regression Graphics_.  John Wiley & Sons, New York, 1994.
  * Dang et al. (2017) U. J Dang, A. Punzo, P. D. McNicholas, S. Ingrassia, & R. P. Browne.  Multivariate response and parsimony for Gaussian cluster-weighted models.  _Journal of Classification_ , 34(1): 4–34, 2017.
  * Dempster et al. (1977) A. P. Dempster, N. M. Laird, & D. B. Rubin.  Maximum likelihood from incomplete data via the EM algorithm.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 39(1): 1–38, 1977.
  * Geweke & Keane (2007) J. Geweke & M. Keane.  Smoothly mixing regressions.  _Journal of Econometrics_ , 138(1): 252–290, 2007.
  * Gormley & Murphy (2010) I. C. Gormley & T. B. Murphy.  Clustering ranked preference data using sociodemographic covariates.  In S. Hess & A. Daly, editors, _Choice Modelling: The State-of-the-Art and the State-of-Practice_ , chapter 25, pages 543–569. Emerald, United Kingdom, 2010.
  * Gormley & Murphy (2011) I. C. Gormley & T. B. Murphy.  Mixture of Experts Modelling with Social Science Applications.  In K. Mengersen, C. Robert, & D. M. Titterington, editors, _Mixtures: Estimation and Applications_ , chapter 9, pages 101–121. John Wiley & Sons, Ltd., New York, 2011.
  * Grün & Leisch (2007) B. Grün & F. Leisch.  Fitting Finite Mixtures of Generalized Linear Regressions in R.  _Computational Statistics & Data Analysis_, 51(11): 5247–5252, 2007.
  * Grün & Leisch (2008) B. Grün & F. Leisch.  FlexMix Version 2: Finite Mixtures with Concomitant Variables and Varying and Constant Parameters.  _Journal of Statistical Software_ , 28(4): 1–35, 2008\.
  * Hennig (2000) C. Hennig.  Identifiability of models for clusterwise linear regression.  _Journal of Classification_ , 17(2): 273–296, 2000.
  * Hubert & Arabie (1985) L. Hubert & P. Arabie.  Comparing Partitions.  _Journal of Classification_ , 2(1): 193–218, 1985.
  * Hurn et al. (2003) M. Hurn, A. Justel, & C. P. Robert.  Estimating Mixtures of Regressions.  _Journal of Computational and Graphical Statistics_ , 12(1): 55–79, 2003.
  * Ingrassia et al. (2012) S. Ingrassia, S. C. Minotti, & G. Vittadini.  Local statistical modeling via the cluster-weighted approach with elliptical distributions.  _Journal of Classification_ , 29(3): 363–401, 2012.
  * Ingrassia et al. (2015) S. Ingrassia, A. Punzo, G. Vittadini, & S. C. Minotti.  The Generalized Linear Mixed Cluster-Weighted Model.  _Journal of Classification_ , 32(1): 85–113, 2015.
  * Jacobs et al. (1991) R. A. Jacobs, M. I. Jordan, S. J. Nowlan, & G. E. Hinton.  Adaptive mixtures of local experts.  _Neural Computation_ , 6: 79–87, 1991.
  * Lamont et al. (2016) A. E. Lamont, J. K. Vermunt, & M. L. Van Horn.  Regression mixture models: Does modeling the covariance between independent variables and latent classes improve the results?  _Multivariate Behavioural Research_ , 51(1): 35–52, 2016\.
  * Mahalanobis (1936) P. C. Mahalanobis.  On the generalised distance in statistics.  _Proceedings National Institute of Science, India_ , 2(1): 49–55, 1936.
  * Mazza et al. (2018) A. Mazza, A. Punzo, & S. Ingrassia.  flexCWM: A flexible framework for cluster-weighted models.  _Journal of Statistical Software_ , pages 1–27, 2018.
  * McCullagh & Nelder (1983) P. McCullagh & J. Nelder.  _Generalized Linear Models_.  Chapman and Hall, London, 1983.
  * McNicholas & Murphy (2008) P. D. McNicholas & T. B. Murphy.  Parsimonious Gaussian Mixture Models.  _Statistics and Computing_ , 18(3): 285–296, 2008.
  * McParland & Gormley (2016) D. McParland & I. C. Gormley.  Model Based Clustering for Mixed Data: clustmd.  _Advances in Data Analysis and Classification_ , 10(2): 155–169, 2016.
  * Murphy & Murphy (2018) K. Murphy & T. B. Murphy.  _MoEClust: Gaussian Parsimonious Clustering Models with Covariates_ , 2018.  URL https://cran.r-project.org/package=MoEClust.  R package version 1.2.1.
  * Ning et al. (2008) H. Ning, Y. Hu, & T. S. Huang.  Efficient initialization of Mixtures of Experts for human pose estimation.  In _Proceedings of the International Conference on Image Processing, ICIP 2008, October 12-15, 2008, San Diego, California, USA_ , pages 2164–2167, 2008.
  * R Core Team (2018) R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2018.
  * Schwarz (1978) G. Schwarz.  Estimating the Dimension of a Model.  _The Annals of Statistics_ , 6(2): 461–464, 1978.
  * Scrucca et al. (2016) L. Scrucca, M. Fop, T. B. Murphy, & A. E. Raftery.  mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models.  _The R Journal_ , 8(1): 289–317, 2016.
  * (33) T. J. Thompson, P. J. Smith, & J. P. Boyle.  Finite mixture models with concomitant information: assessing diagnostic criteria for diabetes.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 47: 393–404.
  * Wang et al. (1996) P. Wang, M. L. Puterman, I. Cockburn, & N. Le.  Mixed Poisson regression models with covariate dependent rates.  _Biometrics_ , 52: 381–400, 1996.
  * Wang et al. (1998) P. Wang, M. L. Puterman, & I. Cockburn.  Analysis of patent data – a mixed-Poisson regression-model approach.  _Journal of Business & Economic Statistics_, 16(1): 27–41, 1998.
  * Wedel & Kamakura (2012) M. Wedel & W. A. Kamakura.  _Market Segmentation: Conceptual and Methodological Foundations_.  International Series in Quantitative Marketing. Springer, US, 2012.
  * Zellner (1962) A. Zellner.  An efficient method of estimating seemingly unrelated regression equations and tests for aggregation bias.  _Journal of the American Statistical Association_ , 57: 348–368, 1962.

##  Appendix A AIS: Stepwise Model Search

The optimal model (according to the BIC criterion) for the AIS data,
highlighted in Table 5, is a $2$-component EVE model where BMI enters the
gating network and sex enters the expert network. This model was identified
after an exhaustive search over a range of $G$ values, the full range of GPCM
covariance parameterisations, and every possible combination of the BMI and
sex covariates in the gating and expert networks. Table A.1 shows the results
of a forward stepwise model selection algorithm, whereby each step can involve
any one of the following actions: adding or removing a component, adding or
removing a single covariate in either network, or changing the GPCM covariance
structure. The algorithm starts from a $1$-component model with no covariates,
with gating concomitants only considered when $G\geq 2$. The action which
yields the best improvement in terms of BIC is retained at each step, and the
algorithm terminates when no further improvement can be found. This forward
stepwise approach is much less computationally onerous than its equivalent in
the backwards direction. Note that the remaining covariates in Table 4 were
also considered for inclusion here, though still only BMI and sex remain in
the model at convergence. Indeed, the algorithm converged to the model chosen
via the BIC criterion under the exhaustive search.

Table A.1: Results of the forward stepwise model selection algorithm applied to the AIS data, detailing the number of steps and the optimal action, resulting model, and associated BIC value at each step. Step | Action | G | GPCM | Gating | Expert | BIC
---|---|---|---|---|---|---
1 | — | 1 | EEE | — |  | -4202.79
2 | Add explanatory variable (Expert) | 1 | EEE | — | sex | -4050.64
3 | Add component | 2 | EEE |  | sex | -4047.40
4 | Change GPCM type | 2 | EVE |  | sex | -4015.35
5 | Add concomitant variable (Gating) | 2 | EVE | BMI | sex | -4013.40
6 | Stop | 2 | EVE | BMI | sex | -4013.40

##  Appendix B AIS: Model Parameters

Given that there exists a binary variable, ‘sex’, in the expert network for
the optimal $G=2$ MoEClust model fit to the AIS data, there are effectively 4
mixture components. Components 1 and 2, corresponding to females and males in
Cluster 1, share the same covariance matrix but differ according to their
means. The same is true for females and males (Components 3 and 4) in Cluster
2. Table B.1 updates the parameters given in Table 6(b) to reflect this.

Table B.1: Estimated parameters of the optimal MoEClust fit to the AIS data, with further splitting of components due to the inclusion of the binary covariate ‘sex’ in the expert network, giving average gates and component mean values (for females and males) and the cluster-specific covariance matrices. |  | Female | Male |  $\hat{\Sigma}_{g}$ (EVE)
---|---|---|---|---
Cluster 1 | $\overline{\hat{\tau}}_{g}$ | 0.31 | 0.25 | RCC | WCC | Hc | Hg | Fe
RCC | 4.47 | 5.11 | 0.13 | 0.15 | 0.90 | 0.31 | 0.64
WCC | 6.71 | 6.71 |  | 1.78 | 1.21 | 0.37 | -1.98
Hc | 40.82 | 45.84 |  |  | 8.22 | 2.76 | 0.51
Hg | 13.61 | 15.56 |  |  |  | 1.16 | 0.62
Fe | 45.92 | 61.80 |  |  |  |  | 487.35
Cluster 2 | $\overline{\hat{\tau}}_{g}$ | 0.19 | 0.25 | RCC | WCC | Hc | Hg | Fe
RCC | 4.30 | 4.94 | 0.07 | 0.05 | 0.49 | 0.17 | 2.54
WCC | 7.46 | 7.74 |  | 4.46 | 0.14 | -0.09 | -7.82
Hc | 39.94 | 45.45 |  |  | 4.52 | 1.55 | 2.05
Hg | 13.47 | 15.55 |  |  |  | 0.62 | 2.45
Fe | 74.79 | 131.99 |  |  |  |  | 1927.93