# Flexible Bayesian modelling of concomitant covariate effects in mixture
models

Marco Berrettini1, Giuliano Galimberti1, Saverio Ranciati1 and
Thomas Brendan Murphy2

1Department of Statistical Sciences, University of Bologna
2School of Mathematics and Statistics, University College Dublin

###### Abstract

Mixture models provide a useful tool to account for unobserved heterogeneity
and are at the basis of many model-based clustering methods. In order to gain
additional flexibility, some model parameters can be expressed as functions of
concomitant covariates. In particular, component weights of the mixture can be
linked to the covariates through a multinomial logistic regression model,
where each component weight is a function of the linear predictor involving
one or more covariates. The proposed contribution extends this approach by
replacing the linear predictor, used for the component weights, with an
additive structure, where each term is a smooth function of the covariates
considered. An estimation procedure within the Bayesian paradigm is suggested.
In particular, a data augmentation scheme based on differenced random utility
models is exploited, and smoothness of the covariate effects is controlled by
suitable choices for the prior distributions of the spline coefficients. The
performance of the proposed methodology is investigated via simulation
experiments, and an application to an original dataset about UK parliamentary
votes on Brexit is discussed.

Keywords: P-splines, data agumentation, additive models, MCMC

##  1 Introduction

Mixture models are the basis of many model-based clustering methods, where
uncertainty can be accounted for in a probabilistic framework. Fields in which
mixture models have been successfully applied include agriculture, astronomy,
bioinformatics, biology, economics, engineering, genetics, imaging, marketing,
medicine, neuroscience, physics psychiatry, psychology and social sciences.
Extensive reviews of mixture models and their application are given in Everitt
and Hand (1981), Titterington et al. (1985), McLachlan and Basford (1988),
Lindsay (1995), Böhning (1999), McLachlan and Peel (2004), Frühwirth-Schnatter
(2006), Mengersen et al. (2011), and McNicholas (2016).

Mixtures of experts (MoE) models provide a way to extend mixture models,
allowing the parameters to depend on concomitant covariate information. This
nomenclature arises in Jacobs et al. (1991), whose authors model the prior
probabilities of latent group membership (i.e. the component weights) as a
logistic function of the covariates. However, some models belonging to the
general class of MoE were already present in the statistical literature under
the name of switching regression models (Quandt, 1972), concomitant variable
latent-class models (Dayton and Macready, 1988), clusterwise regression models
(DeSarbo and Cron, 1988) and mixed models (Wang et al., 1996). Any mixture
model which incorporates covariates or concomitant variables falls within the
family of MoE models, which ranges from the “full” MoE model, in which all
model parameters are functions of covariates, to the special cases where some
of the model parameters do not depend on covariates. A graphical model
representation of the suite of MoE models is provided by Murphy and Murphy
(2020). The MoE framework allows for a wide range of applications: among
others, rank data (Gormley and Murphy, 2008), network data (Gormley and
Murphy, 2010), time series data (Frühwirth-Schnatter et al., 2012), non-normal
data (Chamroukhi, 2015; Nguyen and McLachlan, 2016) and longitudinal data
(Tang and Qu, 2016) have been developed. Recent reviews of MoE models is
provided by Nguyen and Chamroukhi (2018) and Gormley and Frühwirth-Schnatter
(2019).

Jacobs et al. (1991) and Jordan and Jacobs (1994) derive maximum likelihood
(ML) estimates for MoE models via the expectation-maximization (EM) algorithm.
However, the EM algorithm for fitting MoE models is often difficult in
practice, usually due to a complex component density and/or component weights
model or a large parameter set. Inference for MoE models can also be carried
out within the Bayesian paradigm, using a Markov chain Monte Carlo (MCMC)
algorithm. Bayesian approaches for mixture models are summarized by Frühwirth-
Schnatter (2006). Both the Gibbs sampler (Geman and Geman, 1984) and the
Metropolis-Hastings (MH) (Metropolis et al., 1953) algorithm are typically
required, though sampling the component weight parameters through a MH-
algorithm brings issues, such as choosing suitable proposal distributions and
tuning parameters. Alternatively, Frühwirth-Schnatter et al. (2012) introduced
data augmentation of the multinomial logit regression model, based on the
difference random utility model (dRUM) representation, in the context of MoE
models. Following Frühwirth-Schnatter and Frühwirth (2010), the authors
approximate the logistic distribution by a finite scale mixture of Gaussian
distributions in order to avoid any MH step.

This paper focuses on the gating network MoE model (Murphy and Murphy, 2020),
where only the component weights depend on concomitant covariates. In
particular, a more flexible specification of these covariate effects is
considered, in terms of a sum of smooth functions as in a generalized additive
model (Hastie and Tibshirani, 1990). Several proposals are available for
modeling and estimating the smooth functions, see, e.g., Hastie et al. (2009)
and Wood (2017) for an overview. In order to achieve a parsimonious
representation of these smooth functions Bayesian P-splines are used, as
suggested by Lang and Brezger (2004). Following Frühwirth-Schnatter et al.
(2012), data augmentation is exploited based on the dRUM scheme, thus
introducing a set of auxiliary variables.

The remainder of the paper is organized as follows. The data used in the
empirical analysis are introduced in Section 2; specification of the model is
provided in Section 3, and in Section 4 the associated Bayesian inference
procedure is elicited; results from simulation studies are presented in
Section 5, while the ones from real data are reported in Section 6. Finally,
Section 7 is devoted to discussion and conclusions.

##  2 Data

The proposed empirical analysis is based on data about the Parliamentary votes
on Brexit, sometimes referred to as “meaningful votes”, that are the
parliamentary votes under the terms of Section 13 of the United Kingdom’s
European Union (Withdrawal) Act 2018 (Parliament, 2018). This Act requires the
government of the United Kingdom (UK) to bring forward an amendable
parliamentary motion at the end of the Article 50 negotiations between the
government and the European Union (EU) to ratify the Brexit withdrawal
agreement.

The sample considered consists of $n=638$ members of the UK Parliament (MPs)
who voted 16 divisions during the period 25/03/2019 – 01/04/2019. The voting
data were collated using the “hansard” R package (Odell, 2017).

In the United Kingdom, a member of Parliament (MP) is an individual elected to
serve in the House of Commons, the lower house of the Parliament of the United
Kingdom. The Commons is an elected body consisting of 650 members representing
constituencies (electoral areas). Among the MPs not taken into account in this
study there are 4 Speakers (one speaker and 3 deputies) who neither took part
nor voted, 7 members of the political party Sinn Féin which followed a policy
of abstentionism (refusing to attend Parliament or vote on bills) and one MP
who passed away on February 17th, 2019 and was replaced after April 4th, 2019.

At the end of March 2019, the government had not won any of the meaningful
votes. This led to a series of non-binding “indicative votes” on potential
options for Brexit, and also to a delay of the departure date. The amendment
tabled by Conservative MP Sir Oliver Letwin on March 25th, to take power to
control business in the Commons away from the government on March 27th, to
allow MPs to put forward business motions relating to Brexit, passed 329 –
302. Instead, the one tabled by Labour MP Dame Margaret Beckett was defeated
311 – 314. This amendment would have required Parliament to vote favourably
for a “no deal” Brexit, or request an extension to Article 50 if the
government was without a deal within seven days of leaving the European Union.
The amended main motion (Letwin but not Beckett) passed 327 – 300. As a result
of the Letwin amendment’s success, indicative votes on Parliament’s preferred
Brexit options were held on March 27th. Eight propositions were voted upon, of
which all eight failed:

  * •

(B) No deal – Conservative MP Mr John Baron’s proposal to immediately leave
the EU without any deal (For: 160, Against: 400);

  * •

(D) Common market 2.0 – Conservative MP Mr Nicholas Boles’ proposal to join
the Single Market and a customs union (For: 189 – against: 283);

  * •

(E) EFTA and EEA – Conservative MP Mr George Eustice’s proposal to remain in
the Single Market outside of a customs union (For: 64 – Against: 377);

  * •

(J) Customs union – Conservative MP Mr Kenneth Clarke’s proposal for a
permanent customs union (For: 265 – Against: 271);

  * •

(K) Labour’s plan – Labour’s alternative position proposed by MP Mr Jeremy
Corbyn, including a comprehensive customs union with the EU, close alignment
with the Single Market, dynamic alignment on rights and protections,
commitments on participation in EU agencies and funding programmes, and clear
agreements on the detail of future security arrangements (For: 237 – Against:
307);

  * •

(L) Revocation to avoid no deal – Scottish National Party MP Ms Joanna
Cherry’s proposal’s to revoke Article 50 (For: 184 – Against: 293);

  * •

(M) Confirmatory public vote – Labour MP Dame Margaret Beckett’s proposal for
a public vote on any withdrawal bill (For: 268 – Against: 295);

  * •

(O) Managed no deal – Conservative MP Mr Marcus Fysh’s proposal to immediately
leave the EU seeking a tariff-free trade agreement (For: 139 – Against: 422).

It is worth noting that only the first and the last propositions express clear
pro-leave positions, while the other 6 aim at mitigating the effects of
Brexit, or even stopping it.

As Parliament had agreed to an extension of Article 50 to June 30th, the
possibility of a third meaningful vote was raised and took place on March
29th, 2019. Mrs Theresa May promised she would resign as Prime Minister if the
Withdrawal Agreement passed. In the end, Mrs Theresa May’s deal was voted down
again (For: 286 – Against: 344), albeit by a smaller margin than in the
previous two votes that took place on January, 15th, and March, 12th,
respectively. Further indicative votes were held on April 1st on four
propositions chosen by the Speaker, all of which failed:

  * •

(C) Customs union by Conservative MP Mr Kenneth Clarke (For: 273 – Against:
276);

  * •

(D) Common market 2.0 by Conservative MP Mr Nicholas Boles (For: 261 –
Against: 282);

  * •

(E) Confirmatory public vote by Labour MPs Mr Peter Kyle and Mr Phil Wilson
(For: 280 – Against: 292);

  * •

(G) Parliamentary supremacy by Scottish National Party MP Ms Joanna Cherry
(For: 191 – Against: 292).

Note that all the four proposals were modified versions of those put to the
vote on March 27th, although only the proposers of the confirmatory public
vote and the title of Ms Cherry’s one changed (previously, “revocation to
avoid no deal”). The proposal of a third round of indicative votes, to be held
on April, 8th, was then rejected. Due to huge opposition to the fourth
withdrawal agreement, on May, 24th, Mrs Theresa May announced she would resign
as Conservative Party leader and Prime Minister on June, 7th.

The main purpose of the analysis described in this paper is to identify groups
of MPs whose opinions about Brexit, in terms of votes for the aforementioned
divisions, are similar. Furthermore, the influence of some concomitant
information related to the MPs themselves, or the constituencies they
represent, on group membership is considered. In particular, $J=3$ concomitant
covariates are included in the analysis: (i) age of the MP; (ii) share of
Leave votes at the Brexit referendum in parliamentary constituencies; (iii)
“safeness” of the seat of each MP. Regarding the second covariate, it is worth
noting the Brexit referendum vote was not counted by constituency, except in
Northern Ireland. Some local councils (districts) republished local results by
electoral ward or constituency. Some constituencies are coterminous with
(overlap) their local government district. For the others, Hanretty (2017)
estimated through a demographic model for the Leave and Remain vote.

About the third covariate, Apostolova et al. (2017) analyzed and made
available the results of the 2017 UK general election and, in particular, the
number of votes taken by each party for each of the 650 constituencies. In
this dataset, 12 main parties are considered, while all the others are
gathered together, unless one of these won the seat: the votes taken by the
winning party are counted separately and placed in another category, for a
total of $P=13$ categories. To quantify how much an MP, or the party they
represent, was supported in the constituency they were elected into, a measure
of the degree of heterogeneity of votes among parties in that constituency is
considered. In particular, by denoting with $\omega_{pi}$ the share of votes
taken by party $k\,(k=1,\dots,P)$, the entropy (Shannon, 1948) of the votes in
the $i$-th constituency can be computed as:

$\mbox{EN}(\bm{\omega}_{i})=-\sum_{k=1}^{P}\omega_{ki}\log(\omega_{ki}).$

In this case, the entropy $\mbox{EN}(\bm{\omega}_{i})$ quantifies the
uncertainty in predicting the number of votes taken by a party drawn at random
in a given constituency. It ranges from a minimum of 0, which corresponds to a
situation of no heterogeneity (i.e. all votes are taken by a single party), to
a maximum of $\log(P)$ when there is equidistribution of votes between the
parties. Thus, by considering $\exp[\mbox{EN}(\bm{\omega}_{i})]$, a quantity
that can be interpreted as the effective number of competing political parties
(or candidates) in a given constituency is obtained.

##  3 Model specification

Consider an independent and identically distributed sample of outcome
observations $\\{\mathbf{y}_{i}\\}$, with $i=1,\dots,n$, from a population
clustered by a $G$ components finite mixture model. Each component $g$, for
$g=1,\dots,G$, is modelled by the probability (density) function
$f(\mathbf{y}_{i}|\bm{\theta}_{g})$ with parameters $\bm{\theta}_{g}$ and
weight $p_{g}$ , such that $0<p_{g}<1$ and $\sum_{g=1}^{G}p_{g}=1$. All
component densities $f(\mathbf{y}|\bm{\theta}_{g})$, for $g=1,\dots,G$, are
assumed to arise from the same parametric distribution family, depending on
the specific nature of the outcome variable. For the discrete case, one writes
$y_{iqc}=1,c=1,\dots,C_{q},$ if the $i$-th unit presents the $c$-th category
out of $C_{q}$ possibilities for the $q$-th variable, $q=1,\dots,Q$. For each
group $g$, conditional independence is assumed between these $Q$ manifest
variables, leading to the following multinomial distribution:

$f(\mathbf{y}_{i}|\bm{\theta}_{g})=\prod_{q=1}^{Q}\prod_{c=1}^{C_{q}}(\theta_{gqc})^{y_{iqc}}.$

In the Brexit data application, vector
$\bm{\theta}_{gq}=(\theta_{gq1},\theta_{gq2},\theta_{gq3})$ contains the
probabilities of all the possible outcomes (“absent”, “aye” and “no”
respectively) for each of the $Q=16$ divisions and each group $g=1,\dots,G$.

Observation $i$ has $J$ associated covariates $\mathbf{x}_{i}=$ $(1,$
$x_{i1},$ $\dots,$ $x_{iJ^{*}-1},$ $x_{iJ^{*}},$ $x_{iJ^{*}+1},$ $\dots,$
$x_{iJ})$, of which the last $J-J^{*}$ are metrical, with $J^{*}\in[1,J]$. The
gating network MoE model extends the finite mixture model by allowing the
distribution of the latent variable to depend on the concomitant variables:

$f(\mathbf{y}_{i}|\mathbf{x}_{i})=\sum_{g=1}^{G}p_{g}(\mathbf{x}_{i})f(\mathbf{y}_{i}|\bm{\theta}_{g}).$
(1)

Jacobs et al. (1991) model the components’ weights using a multinomial logit
regression model. Arbitrarily selecting a “reference” class - for example, the
$G$-th one, one can assume that the log-odds of the latent class (prior)
membership $p_{g}$, with respect to that class $G$, are linear functions of
the covariates:

$\log\,\frac{p_{g}(\mathbf{x}_{i})}{p_{G}(\mathbf{x}_{i})}=\mathbf{x}_{i}\bm{\gamma}_{g},\quad
g=1,\dots,G-1,$

where $\bm{\gamma}_{g},\,g=1,\dots,G-1,$ denotes the vector of coefficients
corresponding to the $g$-th latent class. Following some simple algebra, this
produces the general result that:

$p_{g}(\mathbf{x}_{i})=\frac{\exp(\mathbf{x}_{i}\bm{\gamma}_{g})}{1+\sum_{g=1}^{G-1}\exp(\mathbf{x}_{i}\bm{\gamma}_{g})}.$

The main goal of this paper is to introduce extra-flexibility in the effects
of the covariates on $p_{g}(\mathbf{x})$ by extending the linear predictor
according to the additive paradigm (see, for example, Green and Silverman
(1993) or Hastie et al. (2009)). More precisely, the linear predictor is
defined as:

$\ln\frac{p_{g}(\mathbf{x}_{i})}{p_{G}(\mathbf{x}_{i})}=\eta_{gi}=\sum_{j=0}^{J^{*}}\gamma_{gj}x_{ij}+\sum_{j=J^{*}+1}^{J}s_{gj}(x_{ij}).$

Here $s_{g,J^{*}+1},\dots,s_{gJ}$ are unknown smooth functions of the metrical
covariates. These smooth functions can be expressed as:

$s_{gj}(x_{ij})=\sum_{\rho=1}^{m}\beta_{gj\rho}B_{j\rho}(x_{ij}),$

were $\bm{\beta}_{g}$ denotes the vector containing the parameters associated
to the nonlinear part of the predictor for the $g$-th component
($g=1,\dots,G-1$), and $B_{j\rho}(\cdot)$ is a B-spline basis function for a
cubic spline ($j=J^{*}+1,\dots,J$ and $\rho=1,\dots,m$) . To ensure
identifiability of the additive predictors, each function
$s_{gj}(\mathbf{x}_{j})$ is constrained to have zero means, that is,

$\frac{1}{\mbox{range}(\;{x}_{j})}\int s_{gj}(x_{j})\,\mbox{d}x_{j}=0,\quad
j=J^{*}+1,\dots,J.$

This can be incorporated into estimation via iterative algorithms by centering
the functions $s_{gj}(\mathbf{x}_{j})$ about their means in every iteration
and thus adding the means to the intercept $\mathbf{\gamma}_{g0}$. See Section
4.1 for further details about Bayesian constrained sampling.

By defining the $n\times m$ design matrices $\mathbf{B}_{j}$, where the
element in row $i$ and column $\rho$ is given by $B_{j\rho}(x_{ij})$, the
predictor can be rewritten in matrix notation:

$\bm{\eta}_{g}=\mathbf{X}\bm{\gamma}_{g}+\sum_{j=J^{*}+1}^{J}\mathbf{B}_{j}\bm{\beta}_{gj},$

where
$\bm{\beta}_{gj}=(\beta_{gj1},\dots,\beta_{gj1})^{\prime},j=J^{*}+1,\dots,J$,
correspond to the vectors of unknown regression coeffcients. The matrix
$\mathbf{X}$ is the design matrix of fixed effects.

##  4 Bayesian Inference

Lang and Brezger (2004) propose to set a high number of knots to ensure enough
flexibility, and define the priors for the regression parameters $\beta_{gj}$
in terms of a random walk:

$\beta_{gj\rho}=\beta_{gj,\rho-1}+u_{gj\rho},\quad
u_{gj\rho}\sim\mbox{N}(0,\tau_{gj}^{2}).$ (2)

The amount of smoothness is controlled by the additional variance parameters
$\tau_{gj}^{2}$, which correspond to the inverse smoothing parameters in the
frequentist approach; this parameter $\tau_{gj}^{2}$ protects against
possibile overfitting if a large number of knots is chosen. The priors can be
equivalently written in the form of global smoothness priors

$\bm{\beta}_{gj}|\tau_{gj}^{2}\propto\exp\left(-\frac{1}{2\tau_{gj}^{2}}\bm{\beta}_{gj}^{\prime}\mathbf{K}_{j}\bm{\beta}_{gj}\right),$

where the penalty matrix $\mathbf{K}_{j}$ is given by
$\mathbf{K}_{j}=\bm{\Delta}_{1}^{\prime}\bm{\Delta}_{1}$, with
$\bm{\Delta}_{1}$ being the first order difference matrix (Rue and Held, 2005,
Chapter 2, p. 52). Because $\mathbf{K}_{j}$ is rank deficient with
$\mbox{rank}(\mathbf{K}_{j})=m-1$ for a first-order random walk, the prior is
improper. It is worth mentioning that these kind of models are usually
referred to in the literature as Intrinsic Gaussian Markov Random Fields (Rue
and Held, 2005). Lang and Brezger (2004) suggest a number of knots between 20
and 40. Hyperpriors are assigned to the variances $\tau_{gj}^{2}$, selecting
Inverse Gamma distributions $\tau_{gj}^{2}\sim\mbox{IG}(a_{gj},b_{gj})$, with
$a_{gj}=1$ and a small value for $b_{gj}$, for example, $b_{gj}=0.005$ leading
to almost diffuse priors for $\tau_{gj}^{2}$. Regarding the coefficients
involved in $\eta_{gi}$, Gaussian priors are assumed for the fixed effects
parameters:
$\bm{\gamma}_{g}\sim\mbox{N}(\mathbf{0}_{J^{*}+1},v\,\mathbf{I}_{J^{*}+1})$,
with variance hyperparameter $v$ set high (e.g. equal to 100), enough to make
the prior non informative.

Following Frühwirth-Schnatter et al. (2012), any multinomial logit model can
be written as a binary one in the partial dRUM representation:

$z_{gi}=\eta_{gi}-\log\left(\sum_{l\neq
g}\lambda_{li}\right)+\epsilon_{gi},\quad D_{gi}=\mathds{1}(z_{gi}>0),$ (3)

where $\mathds{1}$ denotes the indicator function, $z_{gi}$ is a latent
variable, $\lambda_{gi}=\exp(\eta_{gi})$ and $\epsilon_{gi}$ are i.i.d. errors
following a logistic distribution.

Given $\lambda_{1i},\dots,\lambda_{Gi}$ and the latent indicator variables
$D_{1i},\dots,D_{Gi}$, the latent variables $z_{1i},\dots,z_{Gi}$ are
distributed according to an Exponential distribution and can be easily sampled
in a data augmented implementation. To avoid any Metropolis-Hastings step,
Frühwirth-Schnatter and Frühwirth (2010) approximate, for each
$\epsilon_{gi}$, the logistic distribution by a finite scale mixture of H
normal distributions with zero means and variances $\omega_{1},\dots,$ drawn
with fixed probabilities $w_{1},\dots,w_{H}$. The same authors obtained their
finite scale mixture approximation by minimizing the Kullback-Leibler
divergence between the densities, and recommend choosing $H=3$ in larger
applications, where computing time matters, and to work with $H=6$ whenever
possible. In a second step of data augmentation, the component indicator
$r_{gi}$, taking value $h=1,\dots,H$, is introduced as yet another latent
variable. Conditional on the latent variables $\mathbf{z}_{g}$ and the
indicators $\mathbf{r}_{g}=(r_{g1},\dots,r_{gn})$, the binary logit regression
model (3)reduces to a Gaussian regression model.

Regarding the parameters of each component, appropriate full conditionals can
be exploited in order to obtain samples from the posterior distribution. For
the data application, a Dirichlet distribution is assumed for the
probabilities of each division’s possible outcome $\bm{\theta}_{gq}$ (for
$g=1,\dots,G;\,q=1,\dots,Q$), with hyperparameters
$\delta_{1},\dots,\delta_{C_{q}}$ all set equal to 1.

###  4.1 MCMC Algorithm

Based on the representation of Section 4, a new MCMC algorithm is implemented
for fixed $G$ by integrating the scheme proposed by Frühwirth-Schnatter et al.
(2012), with the Bayesian P-spline approach by Lang and Brezger (2004). A
sketch of the algorithm is comprised of the following steps:

  1. 1.

sample the regression coefficients $\bm{\beta}_{g}$ conditional on
$\mathbf{z}_{g}$ and $\mathbf{r}_{g}$, for $g=1,\dots,G-1$. Using the prior in
(2), the conditional posterior of $\bm{\beta}_{gj}$ is given by a multivariate
normal density. Straightforward calculations (Brezger and Lang, 2006) show
that the precision matrix $\mathbf{P}_{gj}$ and the mean $\mathbf{m}_{gj}$ of
$\bm{\beta}_{gj}|\cdot$ are given by

$\begin{split}&\mathbf{P}_{gj}=\mathbf{B}_{j}^{\prime}\mathbf{W}_{g}^{-1}\mathbf{B}_{j}+\frac{1}{\tau_{gj}^{2}}\mathbf{K}_{j},\\\
&\mathbf{m}_{gj}=\mathbf{P}_{gj}^{-1}\mathbf{B}_{j}^{\prime}\mathbf{W}_{g}^{-1}\left(\mathbf{z}_{g}-\bm{\tilde{\eta}}_{g,-j}+\log\sum_{l\neq
g}\bm{\lambda}_{l}\right),\end{split}$

where $\tilde{\bm{\eta}}_{g,-j}$ is the part of the predictor associated with
all but the $j$-th effects in the model, and $\mathbf{W}_{g}$ is a $n\times n$
diagonal matrix with nonzero elements equal to the randomly drawn variances
$\omega_{1g}=s_{r_{g1}}^{2},\dots,\omega_{ng}=s_{r_{gn}}^{2}$ for the $g$-th
group;

  2. 2.

center each smooth function $s_{gj}(\mathbf{x}_{j})$: imposing the constraint
in (3) is equivalent to sampling
$\bm{\beta}_{gj}|(\mathbf{1}_{n}^{\prime}\mathbf{B}_{j}\bm{\beta}_{gj}=\mathbf{0})$
for each $j=J^{*}+1,\dots,J$. Following Algorithm 2.6 in Rue and Held (2005),
this can be done by trasforming each vector of coefficients $\bm{\beta}_{gj}$
as follows:

$\bm{\tilde{\beta}}_{gj}=\bm{\beta}_{gj}-\mathbf{P}_{gj}^{-1}\mathbf{B}_{j}^{\prime}\mathbf{1}_{n}\left(\mathbf{1}_{n}^{\prime}\mathbf{B}_{j}\mathbf{P}_{gj}^{-1}\mathbf{B}_{j}^{\prime}\mathbf{1}_{n}\right)^{-1}\mathbf{1}_{n}^{\prime}\mathbf{B}_{j}\bm{\beta}_{gj};$

  3. 3.

sample the fixed effects parameters $\bm{\gamma}_{g}$ from a multivariate
normal density with precision matrix $\mathbf{P}_{\bm{\gamma}_{g}}$ and the
mean vector $\mathbf{m}_{\bm{\gamma}_{g}}$:

$\begin{split}&\mathbf{P}_{\bm{\gamma}_{g}}=\mathbf{X}^{\prime}\mathbf{W}_{g}^{-1}\mathbf{X}+v^{-1}\mathbf{I}_{J^{*}+1},\\\
&\mathbf{m}_{\bm{\gamma}_{g}}=\mathbf{P}_{\bm{\gamma}_{g}}^{-1}\mathbf{X}^{\prime}\mathbf{W}_{g}^{-1}\left(\mathbf{z}_{g}-\bm{\tilde{\eta}_{g,-\gamma}}+\log\sum_{l\neq
g}\bm{\lambda}_{l}\right),\end{split}$

where $\bm{\tilde{\eta}_{g,-\gamma}}$ represents the nonlinear part of the
$g$-th predictor;

  4. 4.

sample the variance parameters $\tau_{gj}^{2}$ conditional on
$\bm{\tilde{\beta}}_{gj}$:

$\tau^{2}_{gj}|\bm{\beta}_{gj}\sim\mbox{IG}\left(a_{gj}+\frac{\mbox{rank}(\mathbf{K}_{j})}{2},b_{gj}+\frac{1}{2}\bm{\tilde{\beta}}_{gj}^{\prime}\mathbf{K}_{j}\bm{\tilde{\beta}}_{gj}\right);$

  5. 5.

sample all utilities $z_{1i},\dots,z_{G-1,i}$ simultaneously for each $i$
from:

$z_{gi}=\log\left(\frac{\lambda_{gi}}{\log\sum_{l\neq
g}\lambda_{li}}U_{gi}+D_{gi}\right)-\log\left(1-U_{gi}+\frac{\lambda_{gi}}{\log\sum_{l\neq
g}\lambda_{li}}D_{gi}\right),$

with $U_{gi}\sim\mbox{Unif}(0,1);$

  6. 6.

sample the component indicators $r_{gi}$ conditional on $z_{gi}$ from:

$\mbox{Pr}(r_{gi}=h|z_{gi},\bm{\beta}_{g},\bm{\gamma}_{g})\propto\frac{w_{h}}{s_{h}}\exp\left[-\frac{1}{2}\left(\frac{z_{gi}-\eta_{gi}+\log\sum_{l\neq
g}\lambda_{li}}{s_{h}}\right)^{2}\right];$

  7. 7.

sample the cluster-specific conditional probabilities
$\bm{\theta}_{g1},\dots,\bm{\theta}_{gQ}$ given the allocation vector
$\mathbf{D}_{g}=(D_{g1},\dots,D_{gn})^{\prime}$, from a total of $G\cdot Q$
Dirichlet distributions:

$\bm{\theta}_{gq}|\mathbf{D}_{g},\mathbf{y}\sim\mbox{Dir}\left(\delta_{1}+\sum_{i=1}^{n}D_{gi}\cdot
y_{iq1},\dots,\delta_{C_{q}}+\sum_{i=1}^{n}D_{gi}\cdot y_{iqC_{q}}\right);$

  8. 8.

Bayes’ classification for each individual $i$: draw the allocation vector
$\mathbf{D}_{i}=(D_{1i},\dots,D_{Gi})^{\prime}$ from the following discrete
probability distribution which combines the likelihood
$f(\mathbf{y}_{i}|\bm{\theta}_{g})$ and the prior $p_{g}(\mathbf{x}_{i})$:

$\mbox{Pr}(D_{gi}=1|\mathbf{y}_{i},\mathbf{x}_{i},\bm{\beta},\bm{\gamma},\bm{\theta})\propto\frac{\exp(\eta_{gi})}{1+\sum_{g=1}^{G-1}\exp(\eta_{gi})}\prod_{q=1}^{Q}\prod_{c=1}^{C_{q}}(\theta_{gqc})^{y_{iqc}}.$

###  4.2 Posterior inference and model selection

Once the MCMC algorithm has completed the prefixed number $T$ of iterations,
posterior inference is carried out by estimating each parameter’s posterior
mean over the last $T-T_{0}$ draws of the chains, with $T_{0}$ defining the
burn-in phase. Posterior quantities can be computed for the smooth functions
by considering them as linear combinations of spline bases and the
corresponding regression coefficients’ estimates. The uncertainty associated
to the smooth functions is quantified via their pointwise percentiles (usually
2.5 - 97.5 or 5 - 95), for each function over the last $T-T_{0}$ posterior
draws.

Observations can be allocated into the $G$ components using the maximum-a-
posteriori (MAP) rule. In particular, each unit $i=1,\dots,n$ is assigned to
the component $\hat{c}_{i}$ such that

$\hat{c}_{i}=\arg\max_{g}\left(\sum_{t=T_{0}}^{T}D_{1i}^{(t)},\dots,\sum_{t=T_{0}}^{T}D_{Gi}^{(t)}\right).$
(4)

where $\mathbf{D}_{i}^{(t)}=(D_{1i}^{(t)},\dots,D_{Gi}^{(t)})$ represents the
allocation vector for unit $i$ at iteration $t$. Sometimes, using the MAP
rule, one or more components could have no units assigned to them: thus, it
might be worth distinguishing between the number of components $G$ and the
number of non-empty components, denoted as

$\tilde{G}=\sum_{g=1}^{G}\mathds{1}\left(\sum_{i=1}^{n}\mathds{1}(\hat{c}_{i}=g)>0\right).$
(5)

Choosing the number of components in a mixture model is an important problem,
which originated many efforts in the statistical literature. In most
approaches, selecting $G$ is related to the number of free parameters, which
is not clear for the proposed model, due to the presence of regulariziation
induced by the prior distribution on the regression coefficients. One simple
solution in the Bayesian framework is given by the AICM (Raftery et al.,
2007), whose formula depends only on the log-likelihoods from the posterior
simulation:

$\mbox{AICM}=2(\bar{l}-s_{l}^{2}),$ (6)

where $\bar{l}$ and $s_{l}^{2}$ are the sample mean and variance of the
sequence of log-likelihoods
$f(\mathbf{y}_{i}|\bm{\theta}_{\mathbf{D}_{i}}^{(t)})$, for each iteration
$t=T_{0},\dots,T$, after the burn-in. AICM has already been applied
successfully in the mixture modelling context, for example by Erosheva et al.
(2007), Gormley and Murphy (2010), Gormley and Murphy (2011), and Mollica and
Tardella (2017).

###  4.3 Label switching

As for any finite mixture model, label switching may occur during MCMC
sampling; see Frühwirth-Schnatter (2006, Section 3.5) for a review. To
identify a mixture of experts model, Frühwirth-Schnatter et al. (2012) suggest
to focus on a subset of a group-specific parameter and apply $k$-means
clustering (with $G$ clusters) to the posterior draws. MCMC draws belonging to
the same group are assigned to the same cluster by $k$-means clustering, and
the resulting classification sequences
$\zeta_{t}=(S_{1}^{(t)},\dots,S_{G}^{(t)})$ – where each $S_{g}^{(t)}$,
$g=1,\dots,G$, is a classification index taking values in {1,…,G} – show how
to re-arrange the group-specific parameters for each iteration $t=1,\dots,T$,
even if label switching occurred during sampling. In particular, if the
mixture is not overfitting the number $G$ of groups, $\zeta_{t}$ is a
permutation of $\\{1,\dots,G\\}$, and a unique labeling is achieved by
reordering the draws in the following way:

  * •

relabel the hidden allocations $\mathbf{D}_{1},\dots,\mathbf{D}_{n}$ through
the inverse $\zeta_{t}^{-1}$: re-arrange $\mathbf{D}_{1},\dots,\mathbf{D}_{G}$
by $\mathbf{D}_{\zeta_{t}^{-1}(1)},\dots,\mathbf{D}_{\zeta_{t}^{-1}(G)}$,
respectively;

  * •

relabel the group-specific parameters through $\zeta_{t}^{-1}(1)$, $\dots$,
$\zeta_{t}^{-1}(G)$: re-arrange $\bm{\theta}_{1},\dots,\bm{\theta}_{G}$ by
$\bm{\theta}_{\zeta_{t}^{-1}(1)},\dots,\bm{\theta}_{\zeta_{t}^{-1}(G)}$;

  * •

relabel the regression coefficients $\bm{\gamma}$ and $\bm{\beta})$,
corresponding to the linear and nonlinear part of the predictor in the
multinomial logistic regression model: substitute
$\bm{\gamma}_{1},\dots,\bm{\gamma}_{G}$ by
$\bm{\gamma}_{\zeta_{t}^{-1}(1)}-\bm{\gamma}_{\zeta_{t}^{-1}(G)}$,
$\bm{\gamma}_{\zeta_{t}^{-1}(2)}-\bm{\gamma}_{\zeta_{t}^{-1}(G)}$, $\dots$,
$\bm{\gamma}_{\zeta_{t}^{-1}(G)}-\bm{\gamma}_{\zeta_{t}^{-1}(G)}=\mathbf{0}$
and, thanks to the fact that the additive predictors are still linear in the
parameters, substitute $\bm{\beta}_{1},\dots,\bm{\beta}_{G}$ by
$\bm{\beta}_{\zeta_{t}^{-1}(1)}-\bm{\beta}_{\zeta_{t}^{-1}(G)},\bm{\beta}_{\zeta_{t}^{-1}(2)}-\bm{\beta}_{\zeta_{t}^{-1}(G)},\dots,\bm{\beta}_{\zeta_{t}^{-1}(G)}-\bm{\beta}_{\zeta_{t}^{-1}(G)}=\mathbf{0}$.
Subtracting $\bm{\gamma}_{\zeta_{t}^{-1}(G)}$ and
$\bm{\beta}_{\zeta_{t}^{-1}(G)}$, respectively, from all draws ensures that
the regression coefficients of the baseline are all equal to 0 in the
identified model.

##  5 Simulation study

The performance of the proposed approach is investigated in a simulated
environment. In particular, two experiments are carried out, differing for the
true number of components and the distribution of the manifest variables in
the respective simulated datasets. The quality of the estimates for the
covariates’ effects is evaluated through a graphical comparison between the
true effects and the estimated posterior effects, obtained using both the
proposed semiparametric method and by restricting the additive predictor to be
linear. To assess the performance of the estimators of the unknown regression
functions $s_{gj}(\mathbf{x}_{j})$, also the square root of the average
squared errors (RASE) has been considered:

$RASE_{s_{gj}}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(\hat{s}_{gj}(x_{ji}^{*})-s_{gj}(x_{ji}^{*})\right)^{2}},$
(7)

where $\\{x_{ji}^{*},i=1,\dots,n\\}$ are the grid points taken evenly in the
range of each covariate.

Regarding the clustering performance, a comparison between the true
allocations and the estimated ones is made in terms of both Adjusted Rand
Index (ARI) and soft ARI (Flynt et al., 2019). While the ARI measures the
similarity between two data clusterings (hard partitions), the sARI can
incorporate the posterior allocation distributions (soft partitions). Both
indexes take values ranging from a maximum of +1, if the two partitions are
identical (apart from labelling), to a negative minimum if the similarity
between the two partitions is lower than the expected similarity of all
pairwise comparisons between clusterings specified by a random model. For this
second part of the analysis, also Bayesian latent class (BLC) models are
considered, that ignore the effects of concomitant information on the
components’ weights. For each method and each value of $G$, 4000 MCMC draws
have been simulated after a burn-in of 1000 draws. The optimal choice for $G$
is based on the AICM. For the proposed method, the number of parameters for
the cubic splines is set equal to 23.

###  5.1 First simulation experiment: $G=2$

A total of 100 independent datasets have been generated, with $n=1000$ from a
2-components mixture distribution for 5 categorical manifest variables (Table
1). The component weights are assumed to depend on 2 uniformly distributed
covariates $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ in $[0,1]$ as

$\eta_{1}(\mathbf{x}_{1},\mathbf{x}_{2})=2\left[\sin(3\pi\mathbf{x}_{1})e^{(-\mathbf{x}_{1})}+(3\mathbf{x}_{2}-1.5)^{2}\right]-0.5.$

Table 1: Cluster-specific conditional probabilities.  $\theta_{gqc}$ | $g=1$ | $g=2$
---|---|---
$\bm{\theta}_{g1}$ | (0.7, 0.1, 0.2)’ | (0.2, 0.7, 0.1)’
$\bm{\theta}_{g2}$ | (0.2, 0.8)’ | (0.7, 0.3)’
$\bm{\theta}_{g3}$ | (0.3, 0.6, 0.1)’ | (0.1, 0.3, 0.6)’
$\bm{\theta}_{g4}$ | (0.1, 0.1, 0.5, 0.3)’ | (0.5, 0.3, 0.1, 0.1)’
$\bm{\theta}_{g5}$ | (0.1, 0.1, 0.8)’ | (0.1, 0.8, 0.1)’

Figure 1 and Figure 17 in Appendix A show the marginal effects
$s_{11}(\mathbf{x}_{1})$ and $s_{12}(\mathbf{x}_{2})$ along with the average
posterior estimated effects across the 100 simulations obtained on the
simulated dataset using both the semiparametric and the parametric approach
respectively. It appears that the underlying trend can be caught through the
proposed method. Conversely, the parametric competitor clearly cannot properly
approximate the nonlinear non-monotonic trend. These results are summarized in
terms of RASE in Table 2. Goodness of fit (versus complexity) and quality of
the allocations for fixed $G=2$ are analyzed in Table 3. It can be noticed
that the semiparametric approach outperforms its competitors in tems of AICM,
ARI and sARI.

Figure 1: Average posterior marginal effects (green solid line) and 2.5-97.5 pointwise percentiles (green dotted lines) of the concomitant covariates on the log-odds of mixture weights estimated with the semiparametric approach over the 100 simulations. The true effect is represented with a black dashed line. Table 2: Mean (and standard deviation) of the RASEs computed on the estimated log-odds of the mixture weights over 100 simulations.  $RASE_{s_{gj}}$ | $s_{11}(\mathbf{x}^{*})$ | $s_{12}(\mathbf{x}^{*})$
---|---|---
Semiparametric | 0.255 (0.057) | 0.314 (0.088)
Parametric | 0.860 (0.007) | 1.346 (0.002)
Table 3: Average AICM, ARI and sARI (counting how many times each model prevails) over 100 simulations, for fixed $G=2$.  | AICM | (best) | ARI | (best) | sARI | (best)
---|---|---|---|---|---|---
Semiparametric MoE | 8466.5 | (96) | 0.834 | (99) | 0.760 | (100)
Parametric MoE | 8533.7 | (4) | 0.794 | (0) | 0.704 | (0)
BLCA | 9063.4 | (0) | 0.794 | (0) | 0.704 | (0)

For each dataset, the three algorithms have been run setting the number of
components equal to 2, 3 and 4. Table 4 shows the number of non-empty
components for the optimal models selected according to AICM by each method
for all the simulated dataset. It is evident that, in this scenario,
restricting the additive predictor to be linear leads to a wrong choice of the
number of non-empty components $\tilde{G}$ most of times. If one compares the
statistics in Table 3 with the ones computed with reference to the best models
selected according to AICM, the results do not change much for the BLC model.
There is a negligible difference for the semiparametric MoE models due to the
fact that, in one sample, the best model selected according to AICM formally
presents 3 components, though 1 of them is empty. Obviously, for the
parametric MoE models there is an improvement in terms of average AICM
(8510.0), while both the quantities that measures the accuracy of the
allocations get sensibly worse (ARI=0.516, sARI=0.460).

Table 4: Optimal number of non-empty components selected by each method, according to AICM . Method | $g=2$ | $g=3$ | $g=4$
---|---|---|---
Semiparametric MoE | 100 | - | -
Parametric MoE | 26 | 69 | 5
BLCA | 99 | 1 | -

###  5.2 Second simulation experiment: $G=6$

A total of 50 independent datasets have been generated, with $n=1000$ from a
6-components mixture distribution for 12 categorical manifest variables (Table
5). As for the simulation study in Section 5.2, the component weights are
assumed to depend on 2 uniformly distributed covariates $\mathbf{x}_{1}$ and
$\mathbf{x}_{2}$ in $[0,1]$. In particular, the log-odds of the weights are
defined as follows:

$\begin{split}\eta_{1}(\mathbf{x}_{1},\mathbf{x}_{2})&=0.7\left[\sin(3\pi\mathbf{x}_{1})e^{-\mathbf{x}_{1}}+(3\mathbf{x}_{2}-1.5)^{2}\right]-0.5;\\\
\eta_{2}(\mathbf{x}_{1},\mathbf{x}_{2})&=0.5e^{-\mathbf{x}_{1}^{2}}-0.8;\\\
\eta_{3}(\mathbf{x}_{1},\mathbf{x}_{2})&=0.5\sin(6\mathbf{x}_{1}-1)+e^{-16(3\mathbf{x}_{1}-0.5)^{2}}+\frac{e^{-30(\mathbf{x}_{2}-0.3)}}{1+e^{-30(\mathbf{x}_{2}-0.3)}};\\\
\eta_{4}(\mathbf{x}_{1},\mathbf{x}_{2})&=0.6\Bigg{(}3.4827\mathbf{\tilde{x}}_{1}-4.7422\mathbf{\tilde{x}}_{1}^{2}+3.3035\mathbf{\tilde{x}}_{1}^{3}-1.2605\mathbf{\tilde{x}}_{1}^{4}+0.251\mathbf{\tilde{x}}_{1}^{5}\\\
&-0.0204\mathbf{\tilde{x}}_{1}^{6}+\frac{e^{-20(\mathbf{x}_{2}-0.4)}}{1+e^{-20(\mathbf{x}_{2}-0.4)}}\Bigg{)},\quad\mathbf{\tilde{x}}_{1}=2.5\mathbf{x}_{1}+0.5;\\\
\eta_{5}(\mathbf{x}_{1},\mathbf{x}_{2})&=0.5\left(\frac{e^{-10\mathbf{x}_{1}}}{1+e^{-10\mathbf{x}_{1}}}+\frac{e^{-50(\mathbf{x}_{2}-0.3)}}{1+e^{-50(\mathbf{x}_{2}-0.3)}}\right).\end{split}$

Table 5: Cluster-specific conditional probabilities.  $\theta_{gqp}$ | $\bm{\theta}_{g1}=\bm{\theta}_{g2}=\bm{\theta}_{g3}$ | $\bm{\theta}_{g4}=\bm{\theta}_{g5}=\bm{\theta}_{g6}$ | $\bm{\theta}_{g7}=\bm{\theta}_{g8}=\bm{\theta}_{g9}$ | $\bm{\theta}_{g10}=\bm{\theta}_{g11}=\bm{\theta}_{g12}$
---|---|---|---|---
g=1 | (0.7, 0.1, 0.2) | (0.7, 0.1, 0.2) | (0.7, 0.1, 0.2) | (0.7, 0.1, 0.2)
g=2 | (0.7, 0.2, 0.1) | (0.7, 0.2, 0.1) | (0.2, 0.1, 0.7) | (0.2, 0.1, 0.7)
g=3 | (0.1, 0.2, 0.7) | (0.1, 0.2, 0.7) | (0.2, 0.7, 0.1) | (0.2, 0.7, 0.1)
g=4 | (0.7, 0.1, 0.2) | (0.2, 0.1, 0.7) | (0.1, 0.2, 0.7) | (0.1, 0.7, 0.2)
g=5 | (0.1, 0.7, 0.2) | (0.1, 0.7, 0.2) | (0.1, 0.7, 0.2) | (0.1, 0.7, 0.2)
g=6 | (0.1, 0.7, 0.2) | (0.1, 0.2, 0.7) | (0.2, 0.1, 0.7) | (0.7, 0.1, 0.2)

Figures 2 to 6 and Figures 18 to 22 in Appendix B show the marginal effects
$s_{g1}(\mathbf{x}_{1})$ (first row) and $s_{g2}(\mathbf{x}_{2})$ (second
row), for $g=1,\dots,G-1$, along with the average posterior estimated effects
across the 100 simulations obtained on the simulated dataset using both the
semiparametric and the parametric approach respectively. The latter provides
some decent approximations, since most of the underlying trends are
(approximately) linear or monotonic. Moreover, oversmoothing can be easily
noticed in the estimates obtained by the semiparametric MoE model, though
fitting on average seems to be better than that obtained from the parametric
MoE model, and the bands fully contain most of the effects. RASE scores
reported in Table 6 confirm these remarks. Table 7 shows that the advantages
in terms of goodness of fit do not correspond to such differences in terms of
quality of the allocations, when the true number of components if fixed $G=6$.
Nevertheless, the semiparametric MoE appears to prevail most of times,
especially in terms of sARI.

Figure 2: Average posterior marginal effects (grey solid line) and 2.5-97.5 pointwise percentiles (grey dotted lines) of the concomitant covariates on the log-odds of mixture weights for group 1, estimated with the semiparametric approach over the 100 simulations. The true effect is represented with a black dashed line. Figure 3: Average posterior marginal effects (grey solid line) and 2.5-97.5 pointwise percentiles (grey dotted lines) of the concomitant covariates on the log-odds of mixture weights for group 2, estimated with the semiparametric approach over the 100 simulations. The true effect is represented with a black dashed line. Figure 4: Average posterior marginal effects (grey solid line) and 2.5-97.5 pointwise percentiles (grey dotted lines) of the concomitant covariates on the log-odds of mixture weights for group 3, estimated with the semiparametric approach over the 100 simulations. The true effect is represented with a black dashed line. Figure 5: Average posterior marginal effects (grey solid line) and 2.5-97.5 pointwise percentiles (grey dotted lines) of the concomitant covariates on the log-odds of mixture weights for group 4, estimated with the semiparametric approach over the 100 simulations. The true effect is represented with a black dashed line. Figure 6: Average posterior marginal effects (grey solid line) and 2.5-97.5 pointwise percentiles (grey dotted lines) of the concomitant covariates on the log-odds of mixture weights for group 5, estimated with the semiparametric approach over the 100 simulations. The true effect is represented with a black dashed line. Table 6: Mean (and standard deviation) of the RASEs computed on the estimated log-odds of the mixture weights over 100 simulations.  $RASE_{s_{gj}}$ | Semiparametric | Parametric
---|---|---
$s_{11}(\mathbf{x}^{*})$ | 0.251 (0.070) | 0.378 (0.077)
$s_{12}(\mathbf{x}^{*})$ | 0.242 (0.089) | 0.505 (0.035)
$s_{21}(\mathbf{x}^{*})$ | 0.159 (0.077) | 0.208 (0.139)
$s_{22}(\mathbf{x}^{*})$ | 0.122 (0.058) | 0.106 (0.079)
$s_{31}(\mathbf{x}^{*})$ | 0.326 (0.094) | 0.426 (0.121)
$s_{32}(\mathbf{x}^{*})$ | 0.210 (0.095) | 0.289 (0.201)
$s_{41}(\mathbf{x}^{*})$ | 0.143 (0.064) | 0.213 (0.132)
$s_{42}(\mathbf{x}^{*})$ | 0.164 (0.074) | 0.163 (0.077)
$s_{51}(\mathbf{x}^{*})$ | 0.139 (0.062) | 0.212 (0.122)
$s_{52}(\mathbf{x}^{*})$ | 0.167 (0.065) | 0.187 (0.064)
Table 7: Average AICM, ARI and sARI (counting how many times each model prevails) over 100 simulations, for fixed $G=6$.  | AICM | (best) | ARI | (best) | sARI | (best)
---|---|---|---|---|---|---
Semiparametric MoE | 20247.6 | (38) | 0.781 | (26) | 0.691 | (45)
Parametric MoE | 21359.6 | (12) | 0.777 | (17) | 0.684 | (5)
BLCA | 22360.5 | (0) | 0.776 | (7) | 0.682 | (0)

For each dataset, the three algorithms have been run setting the number of
components $G$ ranging from 2 to 8. All the competing methods detect the right
number of non-empty groups apart from 3 times: twice for the parametric MoE
model, once for the BLC model. Thus, the previous comparison based on
clustering performance is repeated for each competing model, the conclusions
do not change sensibly.

##  6 Application to Brexit Voting Data

To identify groups of MPs with similar opinion about Brexit, values of $G$
ranging from 1 to 15 are considered for the semiparametric MoE model described
in Section 3. For each value of $G$, 4000 MCMC draws are considered after a
burn-in of 1000 draws. The optimal number of components suggested by the AICM
is 11, although Figure 7 shows that the AICM curve is quite flat between 9 and
14. Clusters’ composition in terms of political party membership is shown in
Table 8. For space reasons, only the results relative to the most meaningful
clusters detected in the analysis are described.For interpretation reasons, it
is worth mentioning that the log-odds $\eta_{g}(\mathbf{x})$, $g=1,\dots,10$,
are expressed using Cluster 11 as the reference. This cluster is the most
numerous, with 217 MPs, and is characterized by an extreme pro-Leave (even
pro-no-deal) position, as shown in Figure 8. Looking at the composition, 211
out of these 217 MPs are conservatives, including Mr John Baron, proposer of
the no-deal, Mr Marcus Fysh, proposer of the managed no-deal, and Mr Boris
Johnson, who subsequently became prime minister, on July 24th, 2019.

Figure 7: AICM values corresponding to different number of components $G$ for the semiparametric gating network mixture of experts model Table 8: Posterior allocation and political party membership of the $n=638$ MPs. Adjusted Rand Index: 0.458. (C = Conservatives, DUP = Democratic Unionist Party, GP = Green Party, Ind = Independents, Lab = Labour, LD = Liberal Democrats, PC = Plaid Cymru, SNP = Scottish National Party).  | Cluster
---|---
Party | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11
C | 0 | 21 | 16 | 5 | 21 | 0 | 33 | 3 | 3 | 0 | 211
DUP | 0 | 0 | 0 | 10 | 0 | 0 | 0 | 0 | 0 | 0 | 0
GP | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0
Ind | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 2 | 12 | 0 | 1
Lab | 19 | 1 | 4 | 0 | 0 | 108 | 1 | 2 | 6 | 101 | 5
LD | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 9 | 1 | 0 | 0
PC | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 0 | 0 | 0
SNP | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 34 | 0 | 0 | 0
Figure 8: Vote estimates (posterior means) for Cluster 11.

Cluster 8 is the most heterogeneous in terms of political party membership. In
fact, 47 MPs out of 54 belonging to the Scottish National Party (including Ms
Joanna Cherry, who proposed the revocation of Article 50 to avoid no deal),
Plaid Cymru and the Liberal Democrats are in this group, together with 7 more
MPs belonging to 3 other different parties. Since there is no area where the
amount of votes for both the two leading parties (Conservatives and Labour) is
negligible, the constituencies represented in this group are usually
characterized by the presence of at least 3 political parties, i.e. the one
that actually won plus the two leading parties. In other words, the
probability for an MP to belong to this cluster is higher if the effective
number of competing candidates is greater than 3. This threshold (nonlinear)
effect is represented in the third plot of Figure 9. A mild nonlinearity
appears also in the second plot of Figure 9, representing the effect of the
opinion about Brexit in the constituency. In particular, the probability for
an MP to be assigned to Cluster 8 decreases as the fraction of Leave votes in
his constituency increases. In fact, most of the constituencies represented in
this group are known to be pro-Remain, and so are the aforementioned parties,
as confirmed by the posterior means of votes shown in Figure 10.

Figure 9: Estimated smooth effects for Cluster 8. Figure 10: Vote estimates
(posterior means) for Cluster 8.

Also Cluster 6 and 10 are characterized by a pro-Remain position towards
Brexit (see Figure 12 and Figure 14). Both are mostly made up by Labour and
represent constituencies where the competition is usually concentrated around
the two leading parties. Hence, in this case the effective number of competing
candidates has the opposite nonlinear effect on the probability of belonging
to these groups (see the third plot in Figure 11 and Figure 13), if compared
to the one observed for Cluster 9. The main difference between these two
clusters is due to the opinion towards Brexit of the respectively represented
constituencies. In particular, as it can be seen in the second plot of Figure
11 and Figure 13, the fraction of leave votes has a mild nonlinear effect on
the probability of belonging to Cluster 10 (similar to the one observed in
Figure 9), while it has no effect for Cluster 6. This seems to be reflected by
a more extreme position of the MPs of Cluster 10, especially towards the
revocation to avoid no deal (Motions L, G). The aforementioned Dame Margaret
Beckett, first proposer of the confirmatory public vote, as well as of the
homonymous amendment, belongs to this group, while the second co-proposers of
the confirmatory public vote, Mr Peter Kyle and Mr Phil Wilson are divided
between Clusters 10 and 6, respectively. Cluster 6 is also characterized by
the presence of Mr Jeremy Corbyn, proposer of the Labour’s alternative plan
for Brexit, and leader of the Labour Party until the 2019 United Kingdom
general election.

Figure 11: Estimated smooth effects for Cluster 6. Figure 12: Vote estimates
(posterior means) for Cluster 6. Figure 13: Estimated smooth effects for
Cluster 10. Figure 14: Vote estimates (posterior means) for Cluster 10.

There is a third Labour group, Cluster 1, which represents the pro-Leave
minority of the party. More precisely, Figure 15 shows that the MPs belonging
to this cluster distinguish themselves from party members assigned to
different clusters because they were elected in the few pro-Leave Labour
constituencies. This seems to affect their opinion in terms of divisions,
which sensibly departs from the party line. However, this cluster does not
look homogeneous, but rather divided into two fractions. In particular Figure
16 shows 5 divisions where the “no” fraction is close to 50%.

Figure 15: Estimated smooth effects for Cluster 1. Figure 16: Vote estimates
(posterior means) for Cluster 1.

##  7 Conclusions

In this article, a class of finite mixture models, which allows the mixing
proportions to be nonparametric functions of the covariates, has been
proposed. These models exploit an additive structure for the log-odds, and
resort to spline functions for approximating the smooth effect of the
concomitant variables. Simulated experiments proved that the proposed method
can be a useful tool for recollecting the underlying relation between
component weights and concomitant variables - especially when it is not linear
- and, consequently, for estimating models with a better goodness of fit and
leading to a more accurate allocation. The analysis of the determinants of
cluster membership of the United Kingdom MPs reported in Section 6 shows that
political party membership is not sufficient to explain their position towards
Brexit, because both the safeness of seat (in terms of effective number of
competing candidates at the previous general election) and the share of leave
vote in the constituency they represent have a strong - sometimes nonlinear -
effect.

Although the results shown in this paper seem encouraging, there are some
issues that might deserve further investigation. First of all, this paper
focused on models where manifest variables are categorical, but the proposed
methodology can be easily adapted to any other type of response variables, by
choosing an appropriate form for the component density
$f(\mathbf{y}|\bm{\theta}_{g})$.

One of the main advantage of the proposed MCMC algorithm is the absence of
Metropolis-Hastings steps. On the other hand, the use of mixture of Gaussians
to approximate the logistic distribution introduces an additional latent
variable that can increase the computational burden. Furthermore, the
implemented Gibbs sampler requires the number of components as input. If this
quantity is unknown, one needs to estimate it by running the algorithm many
times with different inputs, which might be time consuming, especially when
the “true” value is high. One solution could be incorporating the choice of
the number of components within the algorithm itself. As observed in the
simulation studies, it might happen that the proposed MCMC algorithm converges
to a solution that is characterised by empty components. This peculiar
behaviour could be exploited to devise a strategy similar to the one proposed
by Malsiner-Walli et al. (2016), which circumvents the issue of choosing the
optimal value for $G$ and focuses the attention on the posterior distribution
of the number of non-empty components, by combining a large value for $G$ with
appropriate prior distributions. Alternatively, a reversible jump MCMC
algorithm could be exploited (Richardson and Green, 1997), by designing
appropriate dimension-changing moves (such as split-and-merge moves and birth-
and-death moves).

## References

  * Apostolova et al. (2017) Apostolova, V., L. Audickas, C. Baker, A. Bate, R. Cracknell, N. Dempsey, O. Hawkins, R. McInnes, T. Rutherford, and E. Uberoi (2017).  General election 2017: results and analysis.  Briefing Paper no CBP 7979.
  * Böhning (1999) Böhning, D. (1999).  Computer-assisted analysis of mixtures and applications: meta-analysis, disease mapping and others, Volume 81.  CRC Press.
  * Brezger and Lang (2006) Brezger, A. and S. Lang (2006).  Generalized structured additive regression based on Bayesian P-splines.  Computational Statistics & Data Analysis 50(4), 967–991.
  * Chamroukhi (2015) Chamroukhi, F. (2015).  Non-normal mixtures of experts.  arXiv preprint arXiv:1506.06707.
  * Dayton and Macready (1988) Dayton, C. M. and G. B. Macready (1988).  Concomitant-variable latent-class models.  Journal of the American Statistical Association 83(401), 173–178.
  * DeSarbo and Cron (1988) DeSarbo, W. S. and W. L. Cron (1988).  A maximum likelihood methodology for clusterwise linear regression.  Journal of Classification 5(2), 249–282.
  * Erosheva et al. (2007) Erosheva, E. A., S. E. Fienberg, and C. Joutard (2007).  Describing disability through individual-level mixture models for multivariate binary data.  The Annals of Applied Statistics 1(2), 502–537.
  * Everitt and Hand (1981) Everitt, B. S. and D. J. Hand (1981).  Finite mixture distributions.  CRC Press.
  * Flynt et al. (2019) Flynt, A., N. Dean, and R. Nugent (2019).  sARI: A soft agreement measure for class partitions incorporating assignment probabilities.  Advances in Data Analysis and Classification 13(1), 303–323.
  * Frühwirth-Schnatter (2006) Frühwirth-Schnatter, S. (2006).  Finite mixture and Markov switching models.  Springer.
  * Frühwirth-Schnatter and Frühwirth (2010) Frühwirth-Schnatter, S. and R. Frühwirth (2010).  Data augmentation and MCMC for binary and multinomial logit models.  In Statistical modelling and regression structures, pp. 111–132. Springer.
  * Frühwirth-Schnatter et al. (2012) Frühwirth-Schnatter, S., C. Pamminger, A. Weber, and R. Winter-Ebmer (2012).  Labor market entry and earnings dynamics: Bayesian inference using mixtures-of-experts Markov chain clustering.  Journal of Applied Econometrics 27(7), 1116–1137.
  * Geman and Geman (1984) Geman, S. and D. Geman (1984).  Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.  IEEE Transactions on Pattern Analysis and Machine Intelligence 6(6), 721–741.
  * Gormley and Frühwirth-Schnatter (2019) Gormley, I. C. and S. Frühwirth-Schnatter (2019).  Mixture of experts models.  In Handbook of Mixture Analysis, pp. 271–307. CRC Press.
  * Gormley and Murphy (2008) Gormley, I. C. and T. B. Murphy (2008).  A mixture of experts model for rank data with applications in election studies.  The Annals of Applied Statistics 2(4), 1452–1477.
  * Gormley and Murphy (2010) Gormley, I. C. and T. B. Murphy (2010).  A mixture of experts latent position cluster model for social network data.  Statistical Methodology 7(3), 385–405.
  * Gormley and Murphy (2011) Gormley, I. C. and T. B. Murphy (2011).  Mixture of experts modelling with social science applications.  In Mixtures: Estimation and Applications, pp. 101–121. Wiley Online Library.
  * Green and Silverman (1993) Green, P. J. and B. W. Silverman (1993).  Nonparametric regression and generalized linear models: a roughness penalty approach.  Chapman and Hall/CRC.
  * Hanretty (2017) Hanretty, C. (2017).  Areal interpolation and the UK’s referendum on EU membership.  Journal of Elections, Public Opinion and Parties 27(4), 466–483.
  * Hastie et al. (2009) Hastie, T., R. Tibshirani, and J. Friedman (2009).  The elements of statistical learning: data mining, inference, and prediction.  Springer.
  * Hastie and Tibshirani (1990) Hastie, T. J. and R. J. Tibshirani (1990).  Generalized additive models, Volume 43.  CRC press.
  * Jacobs et al. (1991) Jacobs, R. A., M. I. Jordan, S. J. Nowlan, and G. E. Hinton (1991).  Adaptive mixtures of local experts.  Neural computation 3(1), 79–87.
  * Jordan and Jacobs (1994) Jordan, M. I. and R. A. Jacobs (1994).  Hierarchical mixtures of experts and the EM algorithm.  Neural computation 6(2), 181–214.
  * Lang and Brezger (2004) Lang, S. and A. Brezger (2004).  Bayesian P-splines.  Journal of Computational and Graphical Statistics 13(1), 183–212.
  * Lindsay (1995) Lindsay, B. G. (1995).  Mixture models: theory, geometry and applications.  IMS.
  * Malsiner-Walli et al. (2016) Malsiner-Walli, G., S. Frühwirth-Schnatter, and B. Grün (2016).  Model-based clustering based on sparse finite Gaussian mixtures.  Statistics and Computing 26, 303–324.
  * McLachlan and Basford (1988) McLachlan, G. J. and K. E. Basford (1988).  Mixture models: Inference and applications to clustering.  New York: Dekker.
  * McLachlan and Peel (2004) McLachlan, G. J. and D. Peel (2004).  Finite mixture models.  John Wiley & Sons.
  * McNicholas (2016) McNicholas, P. D. (2016).  Mixture model-based classification.  CRC press.
  * Mengersen et al. (2011) Mengersen, K. L., C. Robert, and M. Titterington (2011).  Mixtures: Estimation and Applications.  John Wiley & Sons.
  * Metropolis et al. (1953) Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller (1953).  Equation of state calculations by fast computing machines.  The Journal of Chemical Physics 21(6), 1087–1092.
  * Mollica and Tardella (2017) Mollica, C. and L. Tardella (2017).  Bayesian Plackett–Luce mixture models for partially ranked data.  Psychometrika 82(2), 442–458.
  * Murphy and Murphy (2020) Murphy, K. and T. B. Murphy (2020).  Gaussian parsimonious clustering models with covariates and a noise component.  Advances in Data Analysis and Classification 14, 293–325.
  * Nguyen and Chamroukhi (2018) Nguyen, H. D. and F. Chamroukhi (2018).  Practical and theoretical aspects of mixture-of-experts modeling: An overview.  Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery 8(4), e1246.
  * Nguyen and McLachlan (2016) Nguyen, H. D. and G. J. McLachlan (2016).  Laplace mixture of linear experts.  Computational Statistics & Data Analysis 93, 177–191.
  * Odell (2017) Odell, E. (2017).  hansard: Provides Easy Downloading Capabilities for the UK Parliament API.  R package version 0.8.0.
  * Parliament (2018) Parliament, U. (2018).  European union (withdrawal) act 2018.  Retrieved 2017/10/19.
  * Quandt (1972) Quandt, R. E. (1972).  A new approach to estimating switching regressions.  Journal of the American Statistical Association 67(338), 306–310.
  * Raftery et al. (2007) Raftery, A., M. Newton, J. Satagopan, and P. Krivitsky (2007).  Estimating the integrated likelihood via posterior simulation using the harmonic mean identity.  In Bayesian Statistics 8, pp. 371–416. Oxford University Press.
  * Richardson and Green (1997) Richardson, S. and P. J. Green (1997).  On Bayesian analysis of mixtures with an unknown number of components (with discussion).  Journal of the Royal Statistical Society: Series B (statistical methodology) 59(4), 731–792.
  * Rue and Held (2005) Rue, H. and L. Held (2005).  Gaussian Markov random fields: theory and applications.  Chapman and Hall/CRC.
  * Shannon (1948) Shannon, C. E. (1948).  A mathematical theory of communication.  The Bell System Technical Journal 27(3), 379–423.
  * Tang and Qu (2016) Tang, X. and A. Qu (2016).  Mixture modeling for longitudinal data.  Journal of Computational and Graphical Statistics 25(4), 1117–1137.
  * Titterington et al. (1985) Titterington, D. M., A. F. Smith, and U. E. Makov (1985).  Statistical analysis of finite mixture distributions.  Wiley,.
  * Wang et al. (1996) Wang, P., M. L. Puterman, I. Cockburn, and N. Le (1996).  Mixed Poisson regression models with covariate dependent rates.  Biometrics 52(2), 381–400.
  * Wood (2017) Wood, S. N. (2017).  Generalized additive models: an introduction with R.  CRC press.

##  Appendix A Simulation study $G=2$, further results

[ht!]

Figure 17: Average posterior marginal effects (green solid line) and 2.5-97.5
pointwise percentiles (green dotted lines) of the concomitant covariates on
the log-odds of mixture weights estimated with the parametric approach over
the 100 simulations. The true effect is represented with a black dashed line.

##  Appendix B Simulation study G=6, further results

Figure 18: Average posterior marginal effects (grey solid line) and 2.5-97.5
pointwise percentiles (grey dotted lines) of the concomitant covariates on the
log-odds of mixture weights for group 1, estimated with the semiparametric
approach over the 100 simulations. The true effect is represented with a black
dashed line. Figure 19: Average posterior marginal effects (grey solid line)
and 2.5-97.5 pointwise percentiles (grey dotted lines) of the concomitant
covariates on the log-odds of mixture weights for group 2, estimated with the
parametric approach over the 100 simulations. The true effect is represented
with a black dashed line. Figure 20: Average posterior marginal effects (grey
solid line) and 2.5-97.5 pointwise percentiles (grey dotted lines) of the
concomitant covariates on the log-odds of mixture weights for group 3,
estimated with the parametric approach over the 100 simulations. The true
effect is represented with a black dashed line. Figure 21: Average posterior
marginal effects (grey solid line) and 2.5-97.5 pointwise percentiles (grey
dotted lines) of the concomitant covariates on the log-odds of mixture weights
for group 4, estimated with the parametric approach over the 100 simulations.
The true effect is represented with a black dashed line. Figure 22: Average
posterior marginal effects (grey solid line) and 2.5-97.5 pointwise
percentiles (grey dotted lines) of the concomitant covariates on the log-odds
of mixture weights for group 5, estimated with the parametric approach over
the 100 simulations. The true effect is represented with a black dashed line.

##  Appendix C Brexit data further results

Cluster 4 is mostly made up by the whole Democratic Unionist Party. A
threshold effect of the number of competing candidates similar (slightly
smoother) to the one observed for Cluster 8 is present (Figure 23).

Figure 23: Estimated smooth effects for Cluster 4 Figure 24: Vote estimates
(posterior means) for Cluster 4

Cluster 9 is the last remaining group which is not mostly composed by MPs
belonging to the leading parties. In particular it includes a large portion of
the Independents group, made up by MPs who quitted the party they were elected
with. Despite this, they seem to reflect the pro-Remain position of the
constituencies where they were elected from (Figures 25 and 26).

Figure 25: Estimated smooth effects for Cluster 9 Figure 26: Vote estimates
(posterior means) for Cluster 9

Cluster 7 includes all of the MPs of the Cabinet (of the time), including Mrs
Theresa May. A high rate of abstensionism is present in this group, apart from
the Letwin-Beckett amendment and the third meaningful vote (Figure 28).
According to Figure 27, covariates seem to have no effect here.

Figure 27: Estimated smooth effects for Cluster 7 Figure 28: Vote estimates
(posterior means) for Cluster 7

Cluster 2 is (mostly) conservative as well, but the opinion of its MPs looks
more heterogenous with respect to Clusters 11 and 7. In particular, a clear
anti no-deal position emerges in Figure 30. This might be due to the
characteristics of the MPs themselves and the constituencies they represent.
In fact, although Conservatives, most of these MPs come from constituencies
that expressed a pro-Remain position and this mild tendency is reflected in
the second plot of Figure 29. Both Mr Nicholas Boles and Mr Kenneth Clarke,
proposers of the 2 couples of divisions for a customs union with the EU
(“Common market 2.0” and, indeed, “Customs union”), belong to this group.

Figure 29: Estimated smooth effects for Cluster 2 Figure 30: Vote estimates
(posterior means) for Cluster 2

The remaining conservative clusters, Cluster 3 and Cluster 5, mainly
distinguish themselves for their moderate positions expressed through the
divisions (Figures 32 and 34) rather than their concomitant information
(Figures 31 and 33). In particular, the latter is characterized by the
presence of Mr George Eustice, whose proposal to remain in the Single Market
outside of a customs union resulted in the division titled “(E) EFTA and EEA”,
which has been favorably voted only by this cluster and Cluster 2.

Figure 31: Estimated smooth effects for Cluster 3 Figure 32: Vote estimates
(posterior means) for Cluster 3 Figure 33: Estimated smooth effects for
Cluster 5 Figure 34: Vote estimates (posterior means) for Cluster 5