Overfitting Bayesian Mixtures of Factor Analyzers with an Unknown Number of Components

TL;DR

This paper introduces an overfitting Bayesian mixture of factor analyzers with an unknown number of components, demonstrating its effectiveness in clustering high-dimensional correlated data using MCMC and information criteria.

Contribution

It presents a novel Bayesian framework for overfitting mixtures of factor analyzers with fixed factors, addressing label switching and estimating the number of components and factors.

Findings

Effective in clustering high-dimensional data

Outperforms state-of-the-art software in simulations

Successfully applied to real datasets

Abstract

Recent advances on overfitting Bayesian mixture models provide a solid and straightforward approach for inferring the underlying number of clusters and model parameters in heterogeneous datasets. The applicability of such a framework in clustering correlated high dimensional data is demonstrated. For this purpose an overfitting mixture of factor analyzers is introduced, assuming that the number of factors is fixed. A Markov chain Monte Carlo (MCMC) sampler combined with a prior parallel tempering scheme is used to estimate the posterior distribution of model parameters. The optimal number of factors is estimated using information criteria. Identifiability issues related to the label switching problem are dealt by post-processing the simulated MCMC sample by relabelling algorithms. The method is benchmarked against state-of-the-art software for maximum likelihood estimation of mixtures of factor analyzers using an extensive simulation study. Finally, the applicability of the method is illustrated in publicly available data.

Figures

Tables

Table 1. Table 1: Prior parameter specification for the case of standardized data.

	$\alpha$	$\beta$	$\gamma$	$g$	$h$	$𝝃={({\xi }_1,\dots ,{\xi }_p)}^T$	$𝚿$
value	0.5	0.5	1	0.5	0.5	${(0,\dots ,0)}^T$	${𝐈}_p$

Table 2. Table 2: Adjusted Rand Index (ARI) between the estimated and true cluster assignments, estimate of the number of clusters ( K ^ ^ 𝐾 \widehat{K} ) and factors ( q ^ ^ 𝑞 \widehat{q} ) for the simulated datasets generated according to three scenarios with different levels of overlapping and covariance structure.

			fabMix			pgmm
	$K$	$q$	ARI	$\hat{K}$	$\hat{q}$	ARI	$\hat{K}$	$\hat{q}$
scenario 1	10	4	1	10	4	0.968	11	4
scenario 2	10	4	1	10	4	0.903	9	4
scenario 3	10	1	1	10	1	1	10	1

Table 3. Table 3: Confusion matrix between the estimated and true cluster assignments of the wave dataset. All methods select 1 1 1 factor and 3 3 3 clusters. The last two lines contain the adjusted and raw Rand index of each method with respect to the true cluster labels.

	fabMix			pgmm
	1	2	3	1	2	3
cluster 1	397	53	56	401	49	56
cluster 2	30	411	40	31	409	41
cluster 3	15	22	476	17	22	474
# factors (BIC)	1			1
ARI	0.616			0.615
RI	0.829			0.829

Table 4. Table 4: Confusion matrix between the estimated and true cluster assignments of the wine dataset. The last two lines contain the adjusted and raw Rand index of each method with respect to the true cluster labels.

	fabMix				pgmm
	1	2	3	4	1	2	3
cluster 1 (Barolo)	0	48	11	0	0	59	0
cluster 2 (Grignolino)	1	18	3	49	1	1	69
cluster 3 (Barbera)	48	0	0	0	48	0	0
# factors (BIC)	2				4
ARI	0.614				0.965
RI	0.830				0.984

Table 5. Table 5: Confusion matrix between the estimated and true cluster assignments of the coffee dataset. The last two lines contain the adjusted and raw Rand index of each method with respect to the true cluster labels.

	fabMix		pgmm
	1	2	1	2	3	4	5
cluster 1 (Arabica)	36	0	0	6	4	10	16
cluster 2 (Robusta)	0	7	7	0	0	0	0
# factors (BIC)	1		3
ARI	1		0.206
RI	1		0.508

Equations

$${\mathbf{x}}_i=\mathbf{\mu }+\mathbf{\Lambda }{\mathbf{y}}_i+{\mathbf{\epsilon }}_i.$$

$${\mathbf{y}}_i\sim {\mathcal{N}}_q(0,{\mathrm{I}}_q),$$

$${\mathbf{\epsilon }}_i\sim {\mathcal{N}}_p(0,\mathbf{\Sigma })$$

$${\mathbf{x}}_i|{\mathbf{y}}_i\sim {\mathcal{N}}_p(\mathbf{\mu }+\mathbf{\Lambda }{\mathbf{y}}_i,\mathbf{\Sigma }).$$

$${\mathbf{x}}_i\sim {\mathcal{N}}_p(\mathbf{\mu },\mathbf{\Lambda }{\mathbf{\Lambda }}^T+\mathbf{\Sigma }),$$

$$\frac{1}{2}p(p+1)-\left[p+pq-\frac{1}{2}q(q-1)\right]=\frac{1}{2}\left[(p-q{)}^2-(p+q)\right].$$

$${\zeta }_{rj}:={\lambda }_{rj}\frac{{\sum }_{i=1}^n{y}_{ij}}{n},$$

$$({\mathbf{x}}_i|z_i=k)={\mathbf{\mu }}_k+{\mathbf{\Lambda }}_k{\mathbf{y}}_i+{\mathbf{\epsilon }}_i.$$

$$\mathrm{P}(z_i=k)=w_k,k=1,\dots ,K,$$

$${\mathbf{\epsilon }}_i|z_i=k$$

$${\mathbf{x}}_i|z_i=k,{\mathbf{y}}_i$$

$${\mathbf{x}}_i\sim \sum _{k=1}^K{w}_k{\mathcal{N}}_p({\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k{\mathbf{\Lambda }}_k^T+{\mathbf{\Sigma }}_k),$$

$${\mathbf{\Sigma }}_k=\mathbf{\Sigma },k=1,\dots ,K.$$

$${\zeta }_{krj}:={\lambda }_{krj}\frac{{\sum }_{i=1}^{n}I(z_i=k)y_{ij}}{{\sum }_{i=1}^{n}I(z_i=k)},$$

$$\mathbf{w}$$

$${\mathbf{\mu }}_k$$

$${\mathbf{\Lambda }}_{kr\cdot }$$

$${\sigma }_{kr}^{-2}$$

$${\omega }_{\ell }^{-2}$$

$$f(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{w},\mathbf{\mu },\mathbf{\Sigma },\mathbf{\Lambda },\mathbf{\Omega })=f(\mathbf{x}|\mathbf{y},\mathbf{z},\mathbf{\mu },\mathbf{\Lambda },\mathbf{\Sigma })f(\mathbf{z}|\mathbf{w})f(\mathbf{\mu })f(\mathbf{\Sigma })$$

$$\times f(\mathbf{\Lambda }|\mathbf{\Omega })f(\mathbf{\Omega })f(\mathbf{w})f(\mathbf{y})$$

$$f_{K_0}(\mathbf{x})=\sum _{k=1}^{K_0}{w}_k{f}_k(\mathbf{x}|{\mathbf{\theta }}_k),$$

$$\mathbf{w}\sim \mathcal{D}\left(\frac{\gamma }{K},\dots ,\frac{\gamma }{K}\right)$$

$$\text{mbox}AIC(K_0,q)$$

$$\text{mbox}BIC(K_0,q)$$

$$\text{mbox}DIC(K_0,q)$$

$$A=\frac{f_i({\mathbf{\phi }}_j^{(t)}|\mathbf{x})f_j({\mathbf{\phi }}_i^{(t)}|\mathbf{x})}{f_i({\mathbf{\phi }}_i^{(t)}|\mathbf{x})f_j({\mathbf{\phi }}_j^{(t)}|\mathbf{x})}=\frac{f_i({\mathbf{\phi }}_j^{(t)})f_j({\mathbf{\phi }}_i^{(t)})}{f_i({\mathbf{\phi }}_i^{(t)})f_j({\mathbf{\phi }}_j^{(t)})}=\frac{{\tilde{f}}_i(w_j^{(t)}){\tilde{f}}_j(w_i^{(t)})}{{\tilde{f}}_i(w_i^{(t)}){\tilde{f}}_j(w_j^{(t)})},$$

$$\mathbf{w}\sim \mathcal{D}\left(\frac{{\gamma }_{(j)}}{K},\dots ,\frac{{\gamma }_{(j)}}{K}\right),$$

$$\frac{x_{ir}-{\bar{x}}_r}{\sqrt{s_r^2}},i=1,\dots ,n;r=1,\dots ,p$$

$${\gamma }_{(j)}=\gamma +\delta (j-1),j=1,\dots ,J,$$

$$({\mathbf{\Lambda }}_{kr\cdot }|\mathbf{\Omega },\mathbf{\Sigma },\mathbf{x},\mathbf{y},\mathbf{z})={\mathcal{N}}_q({\mathbf{\Gamma }}_{kr}{\mathbf{\Delta }}_{kr}^T,{\mathbf{\Gamma }}_{kr}),$$

$${\mathbf{\Gamma }}_{kr}$$

$${\mathbf{\Delta }}_{kr}$$

$${\dot{s}}_k(\mathbf{z})$$

$${\ddot{s}}_k(\mathbf{z},\mathbf{u})$$

Full text

url]panagiotis.papastamoulis@manchester.ac.uk

1 Introduction

Factor Analysis (FA) is a popular statistical model that aims to explain correlations in a high-dimensional space by dimension reduction. This is typically achieved by expressing the observed multivariate data as a linear combination of a smaller set of hypothetical and uncorrelated variables known as factors. The factors are not observed, so they are treated as missing data. The reader is referred to [1, 2] for an overview of factor analysis models, estimation techniques and applications.

However, when the observed data is not homogeneous, the typical FA model will not adequately fit the data. In such a case, a Mixture of Factor Analyzers (MFA) can be used in order to take into account the underlying heterogeneity. Thus, MFA models jointly treat two inferential tasks: model-based density estimation for high dimensional data as well as dimensionality reduction. Estimation of MFA models is straightforward by using the Expectation-Maximization (EM) algorithm [3, 4, 5, 6, 7]. The family of parsimonious Gaussian mixture models (PGMM) is introduced in [8, 9, 10], which is based on Gaussian mixture models with parsimonious factor analysis like covariance structures. Under a Bayesian setup, [11] estimate the number of mixture components and factors by simulating a continuous-time stochastic birth-death point process using a Birth-Death MCMC algorithm [12]. Their algorithm is shown to perform well in small to moderately scaled multivariate data.

Fully Bayesian approaches to estimate the number of components in a mixture model include the Reversible jump MCMC (RJMCMC) [13, 14, 15, 16], Birth-death MCMC (BDMCMC) [12] and allocation sampling [17] algorithms. In recent years there is a growing progress on the usage of overfitted mixture models in Bayesian analysis [18, 19]. An overfitting mixture model consists of a number of components which is much larger than its true (and unknown) value. Under a frequentist approach, overfitting mixture models is not a recommended practice. In this case, the true parameter lies on the boundary of the parameter space and identifiability of the model is violated due to the fact that some of the component weights can be equal to zero or some components may have equal parameters. Consequently, standard asymptotic Maximum Likelihood theory does not apply in this case [20]. Choosing informative prior distributions that bound the posterior away from unidentifiability sets can increase the stability of the MCMC sampler, however these informative priors tend to force too many distinct components and the possibility of reducing the overfitting mixture to the true model is lost (see Section 4.2.2 in [21]). Under suitable prior assumptions introduced by [18], it has been shown that asymptotically the redundant components will have zero posterior weight and force the posterior distribution to put all its mass in the sparsest way to approximate the true density. Therefore, the inference on the number of mixture components can be based on the posterior distribution of the “alive” components of the overfitted model, that is, the components which contain at least one allocated observation.

The simplicity of this approach is in stark contrast with the fully Bayesian approach of treating the number of clusters as a random variable. For example, in the RJMCMC algorithm the researcher has to design sophisticated move types that bridge models with different number of clusters. On the other hand, the allocation sampler is only applicable to cases where the model parameters can be analytically integrated out. Even in such cases though, the design of proper Metropolis-Hastings moves on the space of latent allocation variables of the mixture model is required to obtain a reasonable mixing of the simulated MCMC chain (see [17, 22]).

The contribution of this study is to utilise recent advances on overfitting mixture models [19] to the context of Bayesian MFA [11]. We use a Gibbs sampler [23, 24] which is embedded in a prior parallel tempering scheme in order to improve the mixing of the algorithm. In addition, we explore the usage of information criteria for estimating the number of factors. After estimating the number of clusters and factors, we perform inference on the chosen model by dealing with identifiability issues related to the label switching problem [25]. Our results indicate that overfitting Bayesian MFA models provide a simple and efficient approach to estimate the number of clusters in correlated high-dimensional data.

The rest of the paper is organized as follows. Section 2.1 reviews the basics of FA models. Finite mixtures of FA models are presented in Section 2.2 and a brief review of previous frequentist approaches is given in Section 2.3. The Bayesian formulation is presented in Section 2.4. The overfitting MFA model is introduced in Section 2.5. Section 2.6 deals with estimating the number of factors using information criteria. Section 2.7 presents the prior parallel tempering scheme which is incorporated into the MCMC sampler. Identifiability issues related to the label switching phenomenon are discussed in Section 2.8 and further details of the overall implementation are given in Section 2.9. Our method is illustrated and compared against the EM algorithm in Section 3 using a simulation study (Section 3.1) as well as three publicly available datasets (Sections 3.2, 3.3 and 3.4). The paper concludes in Section 4. Further technical details and simulation results are provided in the Appendix.

2 Methodology

At first we introduce some conventional guidelines that will be followed in our notation throughout this paper, unless explicitly stated otherwise. We will use bold face for vectors and matrices. The notation $\alpha_{k}$ will correspond to the $k$-th member of a vector $\boldsymbol{a}$. In addition, $\boldsymbol{\mathrm{A}}_{k}$ will denote the $k$-th member of a vector $\boldsymbol{\mathrm{A}}$ whose elements are matrices. The $(i,j)$ element of a matrix $\boldsymbol{\Sigma}$ will be denoted by the corresponding lower case letter, that is, $\sigma_{ij}$. The transpose matrix of $\boldsymbol{\Sigma}$ will be denoted as $\boldsymbol{\Sigma}^{T}$. We will not differentiate the notation between random variables and their specific realizations. We use $f(x|y)$ to denote the probability mass or density function of $x$ given $y$. For a discrete random variable $z$, the notation $\mbox{P}(z=k)$ will be also used to denote the probability of the event $\{z=k\}$. The $p\times p$ identity matrix is denoted as $\boldsymbol{\mathrm{I}}_{p}$, $p\in\mathbb{N}$.

2.1 Factor Analysis Model

Let $\boldsymbol{x}=(\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})$ denote a random sample of $p$ dimensional observations with $\boldsymbol{x}_{i}\in\mathbb{R}^{p}$; $i=1,\ldots,n$. We assume that $\boldsymbol{x}_{i}$ is expressed as a linear combination of a latent vector (factors) $\boldsymbol{x}_{i}\in\mathbb{R}^{q}$

[TABLE]

The unobserved random vector $\boldsymbol{y}_{i}$ lies on a lower dimensional space, that is, $q<p$ and it consists of uncorrelated features $y_{i1},\ldots,y_{iq}$. In particular, we assume that

[TABLE]

independent for $i=1,\ldots,n$ and $\boldsymbol{0}$ denotes a vector of zeros. The $p\times q$ dimensional matrix $\boldsymbol{\Lambda}=(\lambda_{rj})$ contains the factor loadings, while the $p$-dimensional vector $\boldsymbol{\mu}=(\mu_{1},\ldots,\mu_{p})$ contains the marginal mean of $\boldsymbol{x}_{i}$. For the error term $\boldsymbol{\varepsilon}_{i}$ assume that

[TABLE]

independent for $i=1,\ldots,n$, where $\boldsymbol{\Sigma}=\mbox{diag}(\sigma_{1}^{2},\ldots,\sigma_{p}^{2})$. Furthermore, $\boldsymbol{y}_{i}$ is assumed independent of $\boldsymbol{\varepsilon}_{i}$; $i=1,\ldots,n$. It follows that

[TABLE]

It can be easily derived that the marginal distribution of $\boldsymbol{x}_{i}$ is

[TABLE]

independent for $i=1,\ldots,n$.

According to the last expression, the covariance matrix of $\boldsymbol{x}_{i}$ is equal to $\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{T}+\boldsymbol{\Sigma}$. As shown in Equation (4), the knowledge of the missing data ($\boldsymbol{y}_{i}$) implies that the conditional distribution of $\boldsymbol{x}_{i}$ has a diagonal covariance matrix. This is the crucial characteristic of factor analysis models, where they aim to explain high-dimensional dependencies using a set of lower-dimensional uncorrelated factors.

Given that there are $q$ factors, the number of free parameters in the covariance matrix $\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{T}+\boldsymbol{\Sigma}$ is equal to $p+pq-\frac{1}{2}q(q-1)$ [26]. The number of free parameters in the uncostrained covariance matrix of $\boldsymbol{x}_{i}$ is equal to $\frac{1}{2}p(p+1)$. Hence, under the factor analysis model, the number of parameters in the covariance matrix is reduced by

[TABLE]

The last expression is positive if $q<\phi(p)$ where $\phi(p):=\frac{2p+1-\sqrt{8p+1}}{2}$, a quantity which is known as the Ledermann bound [27]. We assume that the number of latent factors does not exceed $\phi(p)$. When $q<\phi(p)$ it can be shown that $\boldsymbol{\Sigma}$ is almost surely unique [28]. Note however that this does not necessarily imply that the model is not identified if $q>\phi(p)$ [29].

Given identifiability of $\boldsymbol{\Sigma}$, a second source of identifiability problems is related to orthogonal transformations of the matrix of factor loadings. Indeed, consider a $p\times q$ matrix $\boldsymbol{\Gamma}$ with $\boldsymbol{\Gamma}^{T}\boldsymbol{\Gamma}=\boldsymbol{\mathrm{I}}_{q}$ and $\boldsymbol{\Gamma}\boldsymbol{\Gamma}^{T}=\boldsymbol{\mathrm{I}}_{p}$ (orthogonal matrix) and define $\tilde{\boldsymbol{y}}_{i}=\boldsymbol{\Gamma}\boldsymbol{y}_{i}$. It follows that the representation $\boldsymbol{x}_{i}=\boldsymbol{\mu}+\boldsymbol{\Lambda}\boldsymbol{\Gamma}^{T}\tilde{\boldsymbol{y}_{i}}+\boldsymbol{\epsilon}_{i}$ leads to the same marginal distribution of $\boldsymbol{x}_{i}$ as the one in Equation (5). Following [11], we preassign values to some entries of $\boldsymbol{\Lambda}$, in particular we set the entries of the upper diagonal of the first $(q-1)\times(q-1)$ block matrix of $\boldsymbol{\Lambda}$ equal to zero.

Another identifiability problem is related to the so-called “sign switching” phenomenon, see e.g. [30]. Observe that Equation (1) remains invariant when simultaneously switching the signs of a given row $r$ of $\boldsymbol{\Lambda}$; $r=1,\ldots,p$ and $\boldsymbol{y}_{i}$. Thus, when using MCMC samplers to explore the posterior distribution of model parameters, both $\boldsymbol{\Lambda}$ and $\boldsymbol{y}_{i}$; $i=1,\ldots,n$ are not marginally identifiable due to sign-switching across the MCMC trace. However, one could use the approximate maximum a posteriori estimate arising from the MCMC output in order to infer the mode of the posterior distribution of the specific parameters. These estimates correspond to the parameter values obtained at the iteration that maximizes the posterior distribution across the MCMC run. Another possibility is to restrict our attention to sign-invariant parameter functions. For example, notice that the covariance matrix $\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{T}+\boldsymbol{\Sigma}$ is invariant when switching the sign of each element in $\boldsymbol{\Lambda}$. Following [31], we also define the “regularised score” of variable $r$ to factor $j$ as

[TABLE]

for $r=1,\ldots,p$; $j=1,\ldots,q$. Notice that Equation (6) is invariant to simultaneously switching the signs of $\lambda_{rj}$ and $y_{ij}$’s, therefore, the estimation of $\mathbb{E}(\zeta_{rj}|\boldsymbol{x})$ is meaningful.

2.2 Finite Mixtures of Factor Analyzers

In this section the typical factor analysis model is generalized in order to take into account unobserved heterogeneity. Assume that there are $K$ underlying groups in the population, where the number of clusters $K\geqslant 1$ denotes a known integer. Each cluster is characterized by different structure, which is reflected in our model by assuming a distinct set of parameters $(\boldsymbol{\mu}_{k},\boldsymbol{\Lambda}_{k},\boldsymbol{\Sigma}_{k})$; $k=1,\ldots,K$. Consider the latent allocation parameters $z_{i}\in\{1\ldots,K\}$ which assign observation $\boldsymbol{x}_{i}$ to a cluster $k=1,\ldots,K$ for $i=1,\ldots,n$. Then, given the cluster allocations, each observation is expressed as

[TABLE]

Let $0\leqslant w_{k}\leqslant 1$; $k=1,\ldots,K$ and $\sum_{k=1}^{K}w_{k}=1$. A-priori each observation is generated from cluster $k$ with probability equal to $w_{k}$, that is,

[TABLE]

independent for $i=1,\ldots,n$. We will refer to $\boldsymbol{w}:=(w_{1},\ldots,w_{K})$ as the vector of mixing proportions of the model. Note that the allocation vector $\boldsymbol{z}:=(z_{1},\ldots,z_{n})$ is not observed, so it should be treated as missing data.

In general, the latent factor space can be different for each cluster, e.g. $\boldsymbol{y}_{i}|z_{i}=k\sim\mathcal{N}_{q_{k}}(\boldsymbol{0},\boldsymbol{\mathrm{I}}_{q_{k}})$, where $q_{k}$ represents the latent dimensionality of cluster $k$; $k=1,\ldots,K$. Following [11] we assume that $(\boldsymbol{y}_{i},z_{i})$ are independent and that the intrinsic latent dimension is the same for each cluster, so the marginal distribution of $\boldsymbol{y}_{i}$ is still described by Equation (2).

Now, we express Equations (3), (4) conditionally on the cluster membership

[TABLE]

independent for $i=1,\ldots,n$. Thus, Equation (5) becomes

[TABLE]

independent for $i=1,\ldots,n$. As shown in Equation (10), the marginal distribution of $\boldsymbol{x}_{i}$ is a finite mixture of distributions with $K$ components. Notice that when $q=0$, the parameters $\boldsymbol{\Lambda}_{k}$, $k=1,\ldots,K$ and latent data $\boldsymbol{y}_{i}$ are no longer present in the model. Thus, Equation (10) becomes $\boldsymbol{x}_{i}\sim\sum_{k=1}^{K}w_{k}\mathcal{N}_{p}(\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k})$, that is, a mixture model with diagonal variance structure per component.

A further assumption, also imposed by [11], is the restriction of the error variance, that is,

[TABLE]

In this case $(z_{i},\boldsymbol{\varepsilon}_{i})$ are independent and the marginal distribution of $\boldsymbol{\varepsilon}_{i}$ is given in Equation (3). Although we consider both cases, in the following sections we present the general model where the variance of errors is allowed to vary between clusters.

Finally, the generalization of Equation (6) to the case that there are $K$ clusters is straightforward. For cluster $k$, the regularized score of variable $r$ to factor $j$ is defined as

[TABLE]

for $k=1,\ldots,K$; $r=1,\ldots,p$ and $j=1,\ldots,q$, where $I(\cdot)$ denotes the indicator function. Note that Equation (12) is not defined in case where $\sum_{i=1}^{n}I(z_{i}=k)=0$, however this is not a problem in our implementation due to the fact that we only make inference on the “alive” clusters, that is, the subset of $\{1,\ldots,K\}$ defined as $\boldsymbol{K_{0}}=\{k=1,\ldots,K:\sum_{i=1}^{n}I(z_{i}=k)>0\}$.

2.3 EM-based approaches and available software

A compehensive perspective on the history and development of MFA models is given in Chapter 3 of the monograph by [32]. [4] applied the EM algorithm for estimating the MFA model (10), under the constraint (11). The general model where $\boldsymbol{\Sigma}_{k}$ is allowed to differ between components was estimated by [5]. [33] considered the case of isotropic error variance, that is, the covariance matrix of component $k$ is written as $\boldsymbol{\Lambda}_{k}\boldsymbol{\Lambda}_{k}^{T}+\sigma_{k}\boldsymbol{\mathrm{I}}_{p}$ for $\sigma_{k}>0$ (a model which is referred to as a mixture of probabilistic principal component analyzers).

[8] extended the covariance structure by considering the constraints: $\boldsymbol{\Lambda}_{k}=\boldsymbol{\Lambda}$, $\boldsymbol{\Sigma}_{k}=\boldsymbol{\Sigma}$ and $\boldsymbol{\Sigma}_{k}=\sigma_{k}\boldsymbol{\mathrm{I}}_{p}$. Furthermore, [10] introduced the parameterization $\boldsymbol{\Lambda}_{k}\boldsymbol{\Lambda}_{k}^{T}+\sigma_{k}\boldsymbol{\Delta}_{k}$, where $\sigma_{k}>0$ and $\boldsymbol{\Delta}_{k}=\mbox{diag}(\delta_{1},\ldots,\delta_{p})$ such that $|\boldsymbol{\Delta}_{k}|=1$, with the optional constraints $\boldsymbol{\Delta}_{k}=\boldsymbol{\Delta}$ or $\boldsymbol{\Delta}_{k}=\boldsymbol{\mathrm{I}}_{p}$ for $k=1,\ldots,K$. Depending on whether a particular constraint is present or not, a set of 12 possible models arises which is referred to as the expanded parsimonious Gaussian mixture models (EPGMM) family. A detailed description is provided in Table 2 of [10]. These models are estimated by the alternating expectation-conditional maximization (AECM) algorithm [34].

EMMIXmfa [35] is a freely available software in the form of an R [36] package, implementing the approach used by [5] and allows both constrained and uncostrained error variance per component. The pgmm package [37] is available from the Comprehensive R Archive Network and estimates MFA models using the EPGMM family [8, 10]. Since pgmm offers substantially greater flexibility than EMMIXmfa, we only report results based on the former package.

2.4 Bayesian formulation

This section introduces the Bayesian framework for the MFA model, given $K$, $q$ and a random sample size of $n$ observations $\boldsymbol{x}=(\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})$ from Equation (10). Our aim is to estimate the model parameters: $(\boldsymbol{w},\boldsymbol{\theta})$, where $\boldsymbol{\theta}=\{(\boldsymbol{\mu}_{k},\boldsymbol{\Lambda}_{k},\boldsymbol{\Sigma}_{k});k=1,\ldots,K\}$. The cluster assignments as well as the latent factors $(\boldsymbol{z},\boldsymbol{y})$ are treated as missing data. Let $\mathcal{D}(\cdots)$ denote the Dirichlet distribution and $\mathcal{G}(\alpha,\beta)$ denote the Gamma distribution with mean $\alpha/\beta$. Let also $\boldsymbol{\Lambda}_{kr\cdot}$ denote the $r$-th row of the matrix of factor loadings $\boldsymbol{\Lambda}_{k}$; $k=1,\ldots,K$; $r=1,\ldots,p$. The following prior assumptions are imposed on the model parameters:

[TABLE]

where all variables are assumed mutually independent and $\nu_{r}=\min\{r,q\}$; $r=1,\ldots,p$; $\ell=1,\ldots,q$; $j=1,\ldots,K$. In Equation (15) $\boldsymbol{\Omega}=\mbox{diag}(\omega_{1}^{2},\ldots,\omega_{q}^{2})$ denotes a $q\times q$ diagonal matrix, where the diagonal entries are distributed independently according to Equation (17). The prior assumptions are the same as the ones introduced in [11] for fixed $K$ and $q$.

Given the fixed set of hyperparameters $(K,q,\alpha,\beta,\boldsymbol{\xi},\boldsymbol{\Psi},\gamma,g,h)$, the joint probability density function of the model is written as:

[TABLE]

and its graphical representation is shown in Figure 1. In order to estimate the model parameters we use the Gibss sampler to approximately sample from the posterior distribution $f(\boldsymbol{y},\boldsymbol{z},\boldsymbol{w},\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{\Lambda},\boldsymbol{\Omega}|\boldsymbol{x})$. The reader is referred to A.1 for details.

2.5 Overfitted mixture model

Assume that the observed data has been generated from a mixture model with $K_{0}$ components

[TABLE]

where $f_{k}\in\mathcal{F}_{\Theta}=\{f(\cdot|\boldsymbol{\theta}):\boldsymbol{\theta}\in\Theta\}$; $k=1,\ldots,K_{0}$ denotes a member of a parametric family of distributions. Consider that an overfitted mixture model $f_{K}(\boldsymbol{x})$ with $K>K_{0}$ components is fitted to the data. It has been shown [18] that the asymptotic behaviour of the posterior distribution of the $K-K_{0}$ redudant components depends on the prior distribution of mixing proportions $(\boldsymbol{w})$. Let $d$ denote the dimension of free parameters of the emission distribution $f_{k}$. For the case of a Dirichlet prior distribution in Equation (13), if $\max\{\gamma_{k};k=1,\ldots,K\}<d/2$ then the posterior weight of the extra components converges to zero (Theorem 1 of [18]).

Let $f_{K}(\boldsymbol{\theta},\boldsymbol{z}|\boldsymbol{x})$ denote the joint posterior distribution of model parameters and latent allocation variables for a model with $K$ components. When using an overfitted mixture model, the inference on the number of clusters reduces to (a): choosing a sufficiently large value of mixture components ($K$), (b): running a typical MCMC sampler for drawing samples from the posterior distribution $f_{K}(\boldsymbol{\theta},\boldsymbol{z}|\boldsymbol{x})$ and (c) inferring the number of “alive” mixture components. Note that at MCMC iteration $t=1,2,\ldots$ (c) reduces to keeping track of the number of elements in the set $\boldsymbol{K_{0}}^{(t)}=\{k=1,\ldots,K:\sum_{i=1}^{n}I(z_{i}^{(t)}=k)>0\}$, where $z_{i}^{(t)}$ denotes the simulated allocation of observation $i$ at iteration $t$. At this point, we underline the fact that the MCMC sampler always operates on a space with $K$ components. The empty components at a given iteration may contain allocated observations in subsequent iterations.

In our case the dimension of free parameters in the emission distribution is equal to $d=2p+pq-\frac{q(q-1)}{2}$. We set $\gamma_{1}=\ldots=\gamma_{K}=\frac{\gamma}{K}$, thus the distribution of mixing proportions in Equation (13) becomes

[TABLE]

where $0<\gamma<d/2$ denotes a pre-specified positive number. Such a value is chosen for two reasons. At first, it is smaller than $d/2$ so the asymptotic results of [18] ensure that extra components will be emptied as $n\rightarrow\infty$. Second, this choice can be related to standard practice when using Bayesian non-parametric clustering methods where the parameters of a mixture are drawn from a Dirichlet process [38], that is, a Dirichlet process mixture model [39]. This approach is used in a variety of applications that involve clustering of subjects into an unknown number of groups, see for example [40] where they cluster gene expression profiles using infinite mixture models. Under this point of view, the prior distribution in Equation (19) represents a finite-sum random probability measure that approximates the Dirichlet process with concentration parameter $\gamma$ [41].

2.6 Selecting the Number of Factors

Recall that the MFA model has been defined assuming that the number of factors $q$ is fixed. However, the dimensionality of the latent factor space is not known and should be estimated. In [11] the number of factors is treated as a random variable and a Birth-Death MCMC sampler [12] is used to draw samples from the joint posterior distribution of model parameters and $q$. Treating $q$ as a random variable is preferable from a Bayesian point of view, however it is noted that the sampler may exhibit poor mixing, especially in high dimensions. In our setup, we consider a simpler approach where we produce MCMC samples conditionally on each value of $q$ in a pre-specified range and choose the best one according to model selection criteria.

The following penalized likelihood criteria were used: Akaike’s information Criterion (AIC) [42], Bayesian Information Criterion (BIC) [43], and Deviance Information Criterion (DIC) [44]. Let $D(\theta)=-2\log f(x;\theta)$ denote the deviance. Then,

[TABLE]

where $d_{K_{0},q}$ is the number of free parameters of a model with $K_{0}$ clusters and $q$ factors, whereas $p_{D}=\overline{D(\boldsymbol{\theta})}-D(\widehat{\boldsymbol{\theta}})$ denotes the effective dimension of the model. It is well-known that DIC tends to overfit [45], so we also consider an alternative version which increases the penalty term, namely $\mbox{DIC}_{2}=D(\widehat{\boldsymbol{\theta}})+3p_{D}$ [46, 47]. For AIC and BIC, $\widehat{\boldsymbol{\theta}}$ corresponds to the Maximum Likelihood estimate of $\boldsymbol{\theta}$. Note that under certain conditions [48], the quantity $-\mbox{BIC}/2$ achieves asymptotic consistency by roughly approximating the logarithm of the Bayes factor [49] between two models. The definition of $\widehat{\boldsymbol{\theta}}$ is not unique for DIC (see the discussion in [50]). We have chosen $\widehat{\boldsymbol{\theta}}$ as the simulated parameter values that maximizes the observed log-likelihood across the MCMC run.

It should be noted that in order to compute the observed likelihood as well as the number of free parameters $d_{K_{0},q}$ in AIC and BIC we used only the parameter values of “alive” components, corresponding to the most probable number of clusters (after rescaling the weights in order to sum to 1). It does not make sense to include the redundant parameters in the computation of the number of parameters simply because the results will be inconsistent when considering various values for the total number of clusters in the overfitted mixture ($K$). Also recall here that asymptotically the redundant components are assigned zero posterior weight. This means that the contribution of the extra components to the observed log-likelihood is zero.

2.7 Prior Parallel Tempering

It is well known that the posterior surface of mixture models can exhibit many local modes [51, 52]. In such cases simple MCMC algorithms may become trapped in minor modes and demand a very large number of iterations to sufficiently explore the posterior distribution. In order to produce a well-mixing MCMC sample and improve the convergence of our algorithm we utilize ideas from parallel tempering schemes [53, 54, 55], where different chains are running in parallel and they are allowed to switch states. Each chain corresponds to a different posterior distribution, and usually each one represents a “heated” version of the target posterior distribution. This is achieved by raising the original target to a power $T$ with $0\leqslant T\leqslant 1$, which flattens the posterior surface, thus, easier to explore when using an MCMC sampler.

In the context of overfitting mixture models, [19] introduced a prior parallel tempering scheme. Under this approach, each heated chain corresponds to a model with identical likelihood as the original, but with a different prior distribution. Although the prior tempering can be imposed on any subset of parameters, it is only applied to the Dirichlet prior distribution of mixing proportions. Let us denote by $f_{i}(\boldsymbol{\varphi}|\boldsymbol{x})$ and $f_{i}(\boldsymbol{\varphi})$; $i=1,\ldots,J$, the posterior and prior distribution of the $i$-th chain, respectively. Obviously, $f_{i}(\boldsymbol{\varphi}|\boldsymbol{x})\propto f(\boldsymbol{x}|\boldsymbol{\varphi})f_{i}(\boldsymbol{\varphi})$. Let $\boldsymbol{\varphi}^{(t)}_{i}$ denote the state of chain $i$ at iteration $t$ and assume that a swap between chains $i$ and $j$ is proposed. The proposed move is accepted with probability $\min\{1,A\}$ where

[TABLE]

and $\widetilde{f}_{i}(\cdot)$ corresponds to the probability density function of the Dirichlet prior distribution related to chain $i=1,\ldots,J$. According to Equation (19), this is

[TABLE]

for a pre-specified set of parameters $\gamma_{(j)}>0$ for $j=1,\ldots,J$.

2.8 Label Switching Problem

The likelihood of a mixture model with $K$ components is invariant with respect to permutations of the parameters $(w_{1},\theta_{1}),\ldots,(w_{K},\theta_{K})$. Typically, the same invariance property holds for the prior distribution of $(\boldsymbol{w},\boldsymbol{\theta})$. Therefore, the posterior distribution $f(\boldsymbol{w},\boldsymbol{\theta}|\boldsymbol{x})$ will be invariant to permutations of the labels and it will exhibit (a multiple of) $K!$ symmetric modes. The label switching problem [56, 57, 58] refers to the fact that an MCMC sample that has sufficiently explored the posterior surface will be switching among those symmetric areas.

Note that the symmetry of the likelihood is not of great practical importance under a frequentist point of view, because the EM algorithm will converge to a mode of the likelihood surface. However, under a Bayesian approach it burdens the estimation procedure since the marginal posterior distribution of $(w_{k},\theta_{k})$ will be the same for all $k=1,\ldots,K$. In order to derive meaningful estimates of the marginal posterior distributions as well as estimates of the posterior means the inference should take into account the label switching problem [59, 52, 60, 61, 62]. A variety of relabelling algorithms is available in the label.switching package [25], more specifically we have used the ECR reordering algorithm [60].

Finally we underline that reordering the MCMC sample will only deal with the label switching problem but not the sign switching problem related to the rows of $\boldsymbol{\Lambda}_{k}$; $k=1,\ldots,K$, as shown in Figure 2. Figure 2.(a) illustrates both label and sign switching. After successfully dealing with label switching, the MCMC trace of factor loadings shown in Figure 2.(b) is identifiable up to a sign indeterminacy. But as mentioned earlier, for this particular set of parameters we restrict our attention to sign-invariant parametric functions.

2.9 Details of the implementation

Number of overfitted mixture components In our applications we have considered overfitted mixtures with $K=20$ components. In all cases the sampled values of the number of “alive” components was strictly smaller than this upper bound, at least for models where the number of factors was equal to or larger than the “true” one. If this is not the case for any of the fitted models, the user should consider increasing this upper bound.

Data normalization and prior parameters Before running the sampler, the raw data is standardized by applying the $z$-transformation

[TABLE]

where $\bar{x}_{r}=\frac{\sum_{i=1}^{n}x_{ir}}{n}$ and $s^{2}_{r}=\frac{1}{n-1}\sum_{i=1}^{n}\left(x_{ir}-\bar{x}_{r}\right)^{2}$. The main reason for using standardized data is that the sampler mixes better. Furthermore, it is easier to choose prior parameters that are not depending on the observed data, that is, using the data twice. In any other case, one could use empirical prior distributions as reported in [11], see also [15]. For example, a sensible choice for selecting the prior parameter $\boldsymbol{\xi}$ in case of non-standardized data is to set it equal to the vector containing the mean of each variable. The parameters $\alpha$ and $\beta$ could be selected in a way so that the mode of the prior distribution of $\boldsymbol{\sigma}^{2}_{kr}$ is equal to the variance of each variable, with a large standard deviation. However, in cases where the data contains variables which are measured in very different scales, the algorithm may be quite sensitive to the selection of these prior parameters. Therefore, we advise to use standardized data at the cost of losing interpretability of the parameter estimates.

For the case of standardized data, the prior parameters are specified in Table 1.

Prior parallel tempering We used a total of $J=8$ parallel chains where the prior distribution of mixing proportions for chain $j$ in Equation (21) is selected as

[TABLE]

where $\delta>0$. Since $\gamma=1$ (as shown in Table 1) and $K=20$, it follows from Equation (19) that the parameter vector of the Dirichlet prior of mixture weights which corresponds to the target posterior distribution ($j=1$) is equal to $(0.05,\ldots,0.05)$. Also in our examples we have used $0.25\leqslant\delta\leqslant 5$, but in general we strongly suggest to tune this parameter until a reasonable acceptance rate is achieved. Each chain runs in parallel and every 10 iterations we randomly select two adjacent chains $(j,j+1)$, $j\in\{1,\ldots,J-1\}$ and propose to swap their current states. A proposed swap is accepted with probability $A$ in Equation (20).

Initialization strategy The sampler may require a large number of iterations to discover the high posterior probability areas when initialized from totally random starting values. Of course this behaviour is alleviated as the number of parallel chains ($J$) increases, for example [19] obtained good mixing when $J=30$ for mixtures of univariate normal distributions. In our case, the number of parameters is dramatically larger and an even larger number of parallel chains would be required to obtain similar levels of mixing. However we would rather to test our method in cases where the number of available cores takes smaller values, e.g. $J=8$ or even $J=4$, which are typical numbers of threads in modern-day laptops. Under this point of view, the following two-stage initialization scheme is adopted.

We used an initial period of 100 MCMC iterations where each chain is initialized from totally random starting values, but under a Dirichlet prior distribution with large prior parameter values. These values were chosen in a way that the asymptotic results of [18] guarantee that the redundant mixture components will have non-negligible posterior weights. More specifically for chain $j$ we assume $w\sim\mathcal{D}(\gamma^{\prime}_{j},\ldots,\gamma^{\prime}_{j})$ with $\gamma^{\prime}_{(j)}=\frac{d}{2}+(j-1)\frac{d}{2(J-1)}$, for $j=1,\ldots,J$. We observed that such an initialization quickly reaches to a state where the true clusters are split into many subgroups. Then, we initialize the actual model by this state. In this case, the sampler will spend some time combining the split clusters to homogeneous groups. As soon as all similar subgroups are merged into homogeneous clusters, the sampler will start exploring the stationary distribution of the chain. A comparison with a simpler random starting scheme is presented in A.4 (see also A.3).

Number of parallel chains, MCMC iterations and post-processing Our MCMC sampler ran for $20000$ iterations using $8$ heated chains that run on parallel using the same number of threads. The first $5000$ iterations were discarded as burn-in period and then each chain was thinned by keeping every $10^{\mbox{th}}$ iteration. The MCMC sample corresponding to the retained iterations was reordered according to the ECR algorithm [60] available in the label.switching package [25].

3 Results

At first, we use an extended simulation study in order to evaluate the ability of our method to estimate the correct clustering structure. Our results are benchmarked in terms of clustering estimation accuracy against state-of-the-art software for fitting MFA models with the EM algorithm as implemented in the R package pgmm (version 1.2) [37]. Finally, we illustrate the applicability of our method to three publicly available datasets with known ground-truth classifications. In what follows, we use the acronym fabMix to refer to the proposed method.

3.1 Simulation study

Figure 3.(a) visualizes a simulated dataset with 10 clusters using the procedure in A.2 for $p=40$ variables. The sample size is set to $n=500$ (only a randomly sampled subset of 200 observations is actually displayed in Figure 3). On panels (b) and (c), we juxtapose two similar datasets with identical cluster labels and marginal means, but different covariance structure. It is evident that the clusters exhibit varying levels of overlapping between different scenarios. The underlying latent factor space consists of $q=4$ dimensions for the first two scenarios and $q=1$ for scenario 3.

Both methods estimated MFA models with a number of factors between $q=1,\ldots,8$. The EM-based method fitted a series of mixture models with $K\in\{1,\ldots,20\}$ and selected the optimal number of clusters and factors using the BIC. We also used BIC in order to choose the number of factors in our Bayesian MFA model. The reader is referred to A.5 for a detailed comparison between BIC, AIC and DIC for selecting the number of factors in our model.

It is well known that likelihood methods may be highly sensitive to starting values and that it is quite challenging to adopt optimal schemes for the initialization of the EM algorithm in mixture models (see for example [63]). As with all EM-type algorithms, pgmm may converge to local maxima. We used six starting values: five runs were initialized from random starting values and one run was initialized using the $k$-means clustering algorithm. Regarding the parameterization we have considered all 12 models available in pgmm and the best one is selected according to BIC. The wallclock run-time for fitting all 12 pgmm models with six different starts for $1\leqslant K\leqslant 20$ and $1\leqslant q\leqslant 8$ for Scenario 1 is 62 hours, using 12 parallel threads (one thread per model). The corresponding computing time for fitting the two parameterizations of the overfitted MFA with $K=20$ components for the same range of factors and $20000$ MCMC iterations was equal to 71 hours, based on 8 cores (one core per heated chain).

The results are summarized in Table 2, in terms of the estimation of the number of clusters, factors and clustering accuracy based on the Adjusted Rand Index (ARI) [64]. Note that fabMix exhibits excellent performance in all cases. The best models selected from pgmm correspond to the parameterizations: "UCU", (scenarios 1 and 2) and "CCU" (scenario 3). In all cases, the model with the same variance of errors is selected from fabMix.

Next, we repeat scenario 1 for a total of 300 synthetic datasets, considering that the true number of clusters ranges in the set $K_{0}\in\{1,\ldots,10\}$. For each value of $K_{0}$ we generated 10 datasets consisting of $n=500$ observations, 10 datasets with $n=1000$ and 10 datasets with $n=2000$. As shown in Figure 4 the proposed method is able to infer the true number of clusters as well as the correct clustering structure in terms of the ARI. Under six different starts, the EM algorithm implementation in pgmm gives quite accurate results but note that sometimes the number of clusters may be overestimated when $K\geqslant 9$ and $n=500,1000$ (as shown in the first row of Table 2). Both methods always select the parameterization that corresponds to the one used to generate the data. More specifically, fabMix selects the model with the same variance of errors per cluster and pgmm selects the model corresponding to the "UCU" parameterization (excluding a single case where the "UUU" model was selected).

3.2 Wave dataset

We selected a randomly sampled subset of 1500 observations from the wave dataset [65], available from the UCI machine learning repository [66]. According to the available “ground-truth” classification of the dataset, there are 3 equally weighted underlying classes of $21$-dimensional continuous data, where each one is generated from a combination of 2 of 3 “base” waves, shown in Figure 5. The lower panel displays the sample correlation matrix per cluster, using the corrplot package [67].

We applied the proposed method considering that the number of factors ranges in the set $q\in\{1,\ldots,4\}$. According to BIC, the selected number of factors is equal to $q=1$. Conditionally on $q=1$, our method selects the true value of clusters as the most probable number of “alive components”. In particular the corresponding estimated posterior probability is $\widehat{\rm{P}}(K_{0}=3|\boldsymbol{x})=0.64$. The same model is also chosen by pgmm. The estimated clusters with respect to the true allocation of each observation is shown in Table 3. It is evident that the resulting estimates are very similar to each other (notice that pgmm and fabMix achieve the same ARI) and exhibit strong agreement with the true clustering of the data. fabMix selects the parameterization with the same variance of errors per component and pgmm selects the "UCU" parameterization.

As far as factor interpretation is concerned, we estimated the “regularized score” of each variable to the (single) factor per cluster. For this purpose, we only considered the simulated values that correspond to three alive components and reorder the MCMC values to undo the label switching. Then, we take ergodic averages in order to estimate $\zeta_{krj}=\zeta_{kr}$ (single factor $j=1$) in (12) for $k=1,2,3$ and $r=1,\ldots,21$, as shown in Figure 6. For the first cluster ($k=1$) we note that $\widehat{\zeta}_{1r}$ takes positive values for $2\leqslant r\leqslant 10$ and negative values for $12\leqslant r\leqslant 20$, while $\widehat{\zeta}_{1r}\approx 0$ for $r=1,11,21$. Note the strong agreement when comparing these values with the patterns in the correlation matrix in Figure 5 for true cluster 1. It is evident that the block of variables $\{2,\ldots,10\}$ has different correlation behaviour than $\{12,\ldots,20\}$. Moreover, observe that the block of variables $\{1,11,21\}$ is not correlated to any of the other variables. Similar conclusions can be drawn when comparing the regularized scores for clusters 2 and 3 with their correlation structure.

3.3 Italian wines

The wine dataset [68], available at the pgmm package [37], contains $p=27$ variables measuring chemical and physical properties of $n=178$ wines. There are three types of wine (Barolo, Grignolino, Barbera) from the Piedmont region of Italy which were collected over the ten year period 1970–1979. We applied our method considering $K=20$ and $1\leqslant q\leqslant 5$. For the same number of factors we also applied pgmm, considering that the number of components ranges in the set $1\leqslant K\leqslant 5$. The number of initializations and heated chains are the same as in Section 3.2. As shown in Table 4, fabMix selects a model with 4 clusters and 2 factors and pgmm selects a model with 3 clusters and 4 factors. The best clustering performance corresponds to pgmm, which achieves an almost perfect classification. Furthermore, we report that fabMix selects the parameterization with different variance of errors per component, while pgmm selects the "CUU" parameterization.

3.4 Coffee dataset

The coffee dataset [69] consists of $n=43$ coffee samples collected from beans corresponding to the Arabica and Robusta species. For each sample thirteen variables are observed: water, pH value, fat, chlorogenic acid, bean weight, free acid, caffeine, neochlorogenic acid, extract yield, mineral content, trigonelline, isochlorogenic acid and total chlorogenic acid. The dataset has been previously analyzed using mixtures of factor analyzers by [8] and is available at the pgmm package [37]. Following [8], the total chlorogenic acid is excluded from the analysis since it is the sum of the chlorogenic, neochlorogenic and isochlorogenic acid values.

We applied our method considering $K=20$ and $1\leqslant q\leqslant 5$. For the same number of factors we also applied pgmm, considering that the number of components ranges in the set $1\leqslant K\leqslant 5$. As shown in Table 5, fabMix selects a model with 2 clusters and 1 factor and pgmm selects a model with 5 clusters and 3 factors. The best clustering performance in terms of the adjusted and raw Rand indexes corresponds to fabMix which achieves a perfect classification. Furthermore, we report that fabMix selects the parameterization with different variance of errors per component and pgmm selects the "CCUU" parameterization.

We observed that in this case the estimates obtained by the EM algorithm are quite sensitive to the selection of starting values, a behaviour which typically indicates a highly multimodal likelihood surface. Indeed, based on 100 independent calls to pgmm consisting of six initializations under different values, the algorithm selected a model with $K=2,3,4,5$ clusters with frequencies equal to $41,1,30,28$, respectively. For this reason, the results reported in Table 5 are based on a very large number (1000) of random starting values for pgmm (as well as one run based on $K$-means starting values). We note that a perfect classification was achieved when pgmm selected a model with $K=2$ clusters, as reported in Table 15 of [8]. For each of the 41 runs which selected a model with $K=2$, the "CCUU" parameterization was returned with a BIC value equal to $1292.821$. But the pgmm results in Table 5 correspond to a smaller BIC value ($840.27$) (thus, a better model). Finally, we mention that the same model is still ranked first when using the ICL criterion [70].

4 Discussion and further remarks

This study presented a solution to the complex problem of clustering multivariate data using Bayesian mixtures of factor analyzers. The proposed model builds upon the prior assumptions of [11] assuming a fixed number of clusters and factors. Extending the Bayesian framework of [19] we demonstrated that estimating an overfitting mixture model is a straightforward and efficient approach to the problem at hand. We also used prior parallel tempering schemes in order to improve the mixing of the algorithm. The posterior inference conditionally on a specific model (number of alive clusters) is possible after applying suitable algorithms to deal with label switching.

According to our simulation study, the proposed method can accurately infer the number and composition of the underlying clusters. Our results were benchmarked against state-of-the-art software for estimating MFA models via the EM algorithm, that is, the R package pgmm. Based on six starts of pgmm, we concluded that our method is quite competitive with the EM algorithm in estimation accuracy and can lead to improved inference as the number of clusters grows large. Starting values and the usage of multiple runs are very important for maximum likelihood methods. At this point recall that our benchmarking simulation study generated 100 datasets with dimension $500\times 40$, 100 datasets with dimension $1000\times 40$, and 100 datasets with dimension $2000\times 40$. Using a single EM run for each dataset, the number of compared models for a given parameterization is equal to $K_{\mbox{max}}\times Q_{\mbox{max}}$, with $K_{\mbox{max}}=20$ (maximum number of components) and $Q_{\mbox{max}}$ denoting the maximum number of factors, hence this number should be multiplied by the total number of runs in the case of initializing the EM from multiple starting values (as done in pgmm). On the contrary, the total number of compared models for each parameterization under our Bayesian approach is equal to $Q_{\mbox{max}}$.

We are currently extending the method to the case where the multivariate data contains missing values. This is straightforward to implement under our Bayesian set-up by adding one extra Gibbs sampling step that generates samples from the full conditional distribution in Equation (9). Thus, our future work will focus on presenting, benchmarking and applying this extra functionality to real datasets, as well as improving the practical implementation with the addition of user-friendly output summaries and plot methods. Another interesting generalization of our model would be to propose Bayesian estimation of overfitted mixture models using the whole set of parameterizations introduced in the family of EPGMMs. In this case the dimensionality of the parameter space is reduced further, thus, it is expected that our method will provide more flexible results in certain applications (such as the Italian Wine dataset where pgmm performed better than fabMix).

The source code of the proposed algorithm is hosted online at https://github.com/mqbssppe/overfittingFABMix in the form of an R package, together with scripts that reproduce the results of Section 3. Future versions of the package will be also available on CRAN.

5 Supplementary Material

A.1 contains some technical details regarding the implementation of the Gibbs sampler. A detailed description of the generation of synthetic data can be found in A.2. A.3 compares the proposed overfitting Bayesian MFA model against standard Bayesian MFA models. A.4 discusses convergence diagnostics and comparison with a simpler initialization scheme. Detailed results on the estimation of the number of factors is given in A.5.

6 Acknowledgments

Research was funded by MRC award MR/M02010X/1. The author would like to acknowledge the assistance given by IT services and use of the Computational Shared Facility of the University of Manchester. The suggestions of two anonymous reviewers helped to improve the findings of this study.

Appendix A Supplementary Material

A.1 Gibbs sampling updates

In the following, $(x|y)$ denotes the distribution of $x$ conditionally on $y$. Given $K$ and $q$ the Gibbs sampler updates each parameter according to the following scheme.

Give some initial values $(\boldsymbol{\Omega}^{(0)},\boldsymbol{\Lambda}^{(0)},\boldsymbol{\mu}^{(0)},\boldsymbol{z}^{(0)},\boldsymbol{\Sigma}^{(0)},\boldsymbol{w}^{(0)},\boldsymbol{y}^{(0)})$. 2. 2.

At iteration $t=1,2,\ldots$

1. (a)

update $\boldsymbol{\Omega}^{(t)}\sim\left(\boldsymbol{\Omega}|\boldsymbol{\Lambda}^{(t-1)}\right)$ 2. (b)

update $\boldsymbol{\Lambda}^{(t)}\sim\left(\boldsymbol{\Lambda}|\boldsymbol{\Omega}^{(t)},\boldsymbol{\Sigma}^{(t-1)},\boldsymbol{x},\boldsymbol{y}^{(t-1)},\boldsymbol{z}^{(t-1)}\right)$ 3. (c)

update $\boldsymbol{\mu}^{(t)}\sim\left(\boldsymbol{\mu}|\boldsymbol{\Lambda}^{(t)},\boldsymbol{\Sigma}^{(t-1)},\boldsymbol{x},\boldsymbol{y}^{(t-1)},\boldsymbol{z}^{(t-1)}\right)$ 4. (d)

update $\boldsymbol{z}^{(t)}\sim\left(\boldsymbol{z}|\boldsymbol{w}^{(t-1)},\boldsymbol{\mu}^{(t)},\boldsymbol{\Lambda}^{(t)},\boldsymbol{\Sigma}^{(t-1)},\boldsymbol{x}\right)$ 5. (e)

update $\boldsymbol{\Sigma}^{(t)}\sim\left(\boldsymbol{\Sigma}|\boldsymbol{x},\boldsymbol{z}^{(t)},\boldsymbol{\mu}^{(t)}\right)$ 6. (f)

update $\boldsymbol{w}^{(t)}\sim\left(\boldsymbol{w}|\boldsymbol{z}^{(t)}\right)$ 7. (g)

update $\boldsymbol{y}^{(t)}\sim\left(\boldsymbol{y}|\boldsymbol{x},\boldsymbol{z}^{(t)},\boldsymbol{\mu}^{(t)},\boldsymbol{\Sigma}^{(t)},\boldsymbol{\Lambda}^{(t)}\right)$.

All full conditional distributions are detailed in [11]. However, in some cases, it is easy to analytically integrate some parameters and accelerate convergence. Note that in steps (2.b) and (2.d) we have integrated out $\boldsymbol{\mu}$ and $\boldsymbol{y}$, respectively. Thus, in steps (2.b) and (2.d), the parameters are updated from the corresponding collapsed conditional distributions, rather than the full conditional distributions.

In particular, after integrating out $\boldsymbol{\mu}=(\boldsymbol{\mu}_{1},\ldots,\boldsymbol{\mu}_{K})$, it follows that

[TABLE]

independent for $k=1,\ldots,K$; $r=1,\ldots,p$. The parameters on the right-hand side of Equation (22) are defined as

[TABLE]

where

[TABLE]

for $k=1,\ldots,K$. Finally, from Equations (7) and (10) it immediately follows that

[TABLE]

independent for $i=1,\ldots,n$, where $f(\cdot;\boldsymbol{\mu},\boldsymbol{\Sigma})$ denotes the probability density function of the multivariate normal distribution with mean $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$. Note that in order to compute the right hand side of the last equation, inversion of the $p\times p$ matrix $\boldsymbol{\Lambda}_{k}\boldsymbol{\Lambda}_{k}^{T}+\boldsymbol{\Sigma}_{k}$ is required. Since $\boldsymbol{\Sigma}_{k}$ is diagonal, this matrix can be efficiently inverted using the Sherman–Morrison–Woodbury formula [71]:

[TABLE]

for $k=1,\ldots,K$.

A.2 Simulation study details

The synthetic data in Section 3.1 was simulated according to mixtures of multivariate normal distributions, as shown in Equation (10). Given the number of mixture components ($K$) and number of factors ($q$) for each dataset, the true values of parameters were generated according to the following scheme.

[TABLE]

independent for $r=1,\ldots,p$; $k=1,\ldots,K$. Observe that the diagonal matrix containing the variance of errors ($\sigma^{2}_{krr}$) is the same for each mixture component. Under a slight abuse of notation in the $q\times p$ matrix $\boldsymbol{\Lambda}_{k}^{T}$, a star ($\bigstar$) denotes independent realizations from a random variable following a

[TABLE]

distribution. Moreover, a square ($\square$) denotes independent realizations from a $\mathcal{N}(0,1)$ random variable. Each block matrix containing a star consists of $q$ rows and $\lfloor p/q\rfloor$ columns, where $\lfloor\cdot\rfloor$ denotes the integer part of a positive number. In case that $p/q$ is not integer, we fill the remaining $p-q\lfloor p/q\rfloor$ columns with squares (corresponding to the last block of $\boldsymbol{\Lambda}_{k}^{T}$). The parameters of the distribution (24) are generated as follows:

[TABLE]

and

[TABLE]

for $k=1,\ldots,K$. In all cases $\phi_{k}$ and $\rho_{k}$ are assumed independent for $k=1,\ldots,K$.

A.3 Benchmarking against standard initialization and Bayesian MFA models

In this section we compare against standard Bayesian MFA models. We considered $30$ simulated datasets with $n=500$ and $p=20$ with a fixed number of factors $q=2$, which is assumed known. Then, we estimated the number of clusters using the BIC by fitting a series of standard Bayesian MFA models with a number of components equal to $K=1,\ldots,20$. The results are compared against the proposed method using $J=8$ parallel chains. Note that no parallel chains are considered in the standard Bayesian MFA MCMC sampler. In order to make the comparison as fair as possible, we considered that the number of MCMC iterations is equal to $Jm$ and $m$ for standard and overfitting MFA, respectively, with $m=20000$ as in the previous sections.

As shown in Figure 7(a), the proposed method gives more accurate results than standard Bayesian MFA models, especially as the number of components grows large. Figure 7(b) displays the behaviour of the BIC values of the standard Bayesian MFA model for two specific datasets where the true number of components is equal to $K=9$. Observe that the standard MCMC sampler selects $K=10$ clusters for dataset 1 and $K=9$ clusters for dataset 2. Furthermore, note that the BIC values, as a function of the number of components, are not smooth for both datasets. Such a behaviour is typical in cases where the algorithm has converged to minor modes (see [63]).

A.4 Convergence diagnostics

In order to evaluate the convergence of our algorithm, we have ran the sampler for the same number of iterations ($m=20000$) based on different starting values for a total of 10 runs using a simulated dataset with $K=10$ clusters and $q=4$ factors. The number of observations is equal to $n=500$ for a total of $p=40$ variables generated according to scenario 1. Using the coda package [72], the Gelman-Rubin diagnostic criterion [73] varies between $(1.00,1.012)$ (mixing proportions), $(1.00,1.04)$ (means) and $(1.00,1.13)$ (covariance matrices). Note that the clusters have been relabelled across all runs so that they agree to each other. This range of values is quite smaller than $1.2$, so according to [74] we can be fairly confident that convergence has been reached. Figure 8 illustrates the evolution of the Gelman-Rubin diagnostic criterion for the mixing proportions as the number of iterations increases.

Furthermore, we compare the performance of our MCMC sampler when a different initialization scheme is applied. Recall that the default initialization scheme is based on a first stage of 100 iterations where the mixing proportions are a-priori distributed according to a Dirichlet distribution with parameters that lead to overfitting, as discussed in Section 2.9. Now we consider a simpler scheme that does not include this stage and all parameters are randomly generated from the default prior distributions. Under the same number of iterations and independent runs, we calculated that the Gelman-Rubin diagnostic criterion [73] varies between $(1.00,1.45)$ (mixing proportions), $(1.00,1.86)$ (means) and $(1.15,2.31)$ (covariance matrix). This means that the MCMC sampler has not converged when using this simpler initialization scheme, under the same number of iterations.

A.5 Estimating the number of factors

In this section we assess the ability of our algorithm to infer the number of factors. We generated synthetic data from mixtures of factor analyzers where the true number of factors varies in the set $q\in\{1,\ldots,5\}$ and for each value 15 datasets were simulated. The true number of clusters is uniformly distributed on the set $\{1,\ldots,10\}$. For each case we considered that $n=1000$ and $p=40$.

A typical output of the proposed method for a single dataset is shown in Figure 9. Note that we have also considered the case $q=0$, which is explicitly taken into account into our method. Observe that when the number of factors is less than the true one ($q=4$) the posterior distribution of the number of alive components supports larger values than the true number of clusters ($K=4$). As soon as the true number of factors is reached, then the posterior distribution of the number of alive clusters remains concentrated at the true value.

For each dataset, the number of factors was estimated by running our algorithm considering that $q\in\{0,1,\ldots,10\}$ and selecting the best model according to AIC, BIC, DIC and $\mbox{DIC}_{2}$. As shown in Fig 9 (down), both AIC and BIC are able to correctly estimate the number of factors. On the other hand, DIC and $\mbox{DIC}_{2}$ tend to overfit, as often concluded in the literature (see e.g. [45]). However, even in the case that the Deviance Information Criterion is used for choosing the number of factors, the corresponding estimate of the underlying clusters (as well as the number of “alive” clusters) is the same as the one returned by BIC or AIC. This behaviour is also illustrated in Figure 9 (up) for a single dataset. Consequently, given that the MCMC sampler has converged, using the DIC or $\mbox{DIC}_{2}$ does not necessarily means that the inferred clustering is “wrong”.