Clustering Multivariate Data using Factor Analytic Bayesian Mixtures with an Unknown Number of Components

TL;DR

This paper introduces an advanced Bayesian mixture model for clustering multivariate data that automatically determines the number of clusters and covariance structures using overfitting models, MCMC sampling, and model selection criteria.

Contribution

It extends overfitting Bayesian mixture models by incorporating multiple covariance parameterizations and a robust sampling scheme, improving clustering flexibility and model selection.

Findings

Effective in estimating the number of clusters automatically.

Performs well on simulated and real datasets.

Provides an accessible R package for implementation.

Abstract

Recent work on overfitting Bayesian mixtures of distributions offers a powerful framework for clustering multivariate data using a latent Gaussian model which resembles the factor analysis model. The flexibility provided by overfitting mixture models yields a simple and efficient way in order to estimate the unknown number of clusters and model parameters by Markov chain Monte Carlo (MCMC) sampling. The present study extends this approach by considering a set of eight parameterizations, giving rise to parsimonious representations of the covariance matrix per cluster. A Gibbs sampler combined with a prior parallel tempering scheme is implemented in order to approximately sample from the posterior distribution of the overfitting mixture. The parameterization and number of factors is selected according to the Bayesian Information Criterion. Identifiability issues related to label switching are dealt by post-processing the simulated output with the Equivalence Classes Representatives algorithm. The contributed method and software are demonstrated and compared to similar models estimated using the Expectation-Maximization algorithm on simulated and real datasets. The software is available online at https://CRAN.R-project.org/package=fabMix.

Tables

Table 1. Table 1: Arguments of the fabMix() function.

Argument	Description
model	Any non-empty subset of c("UUU", "CUU", "UCU", "CCU", "UCC", "UUC", "CUC", "CCC"), indicating the fitted models. By default, all models are fitted.
Kmax	Number of components in the overfitted mixture (integer, at least equal to two). Default: 20.
nChains	Number of parallel (heated) chains. When dirPriorAlphas is supplied, this argument can be ignored.
dirPriorAlphas	vector of length nChains in the form of an increasing sequence of positive scalars. Each entry contains the (common) prior Dirichlet parameter for the corresponding chain. Default: dirPriorAlphas = c(1, 1 + dN*(2:nChains - 1))/Kmax, where dN = 1, for nChains > 1. Otherwise, dirPriorAlphas = 1/Kmax.
rawData	The observed data in the form of an $n\times p$ matrix. Clustering is performed on the rows of the matrix.
outDir	Name of the output folder. An error is thrown if the directory already exists inside the current working directory. Note: it should not correspond to an absolute path, e.g.: outDir = ‘example‘ is acceptable, but outDir = ‘C:\User\Documents\example‘ is not.
mCycles	Number of MCMC cycles. Each cycle consists of nIterPerCycle MCMC iterations. At the end of each cycle a swap of the state of two randomly chosen adjacent chains is attempted.
burnCycles	Number of cycles that will be discarded as burn-in period.
g	Prior parameter $g$. Default value: g = 0.5.
h	Prior parameter $h$. Default value: g = 0.5.
alpha_sigma	Prior parameter $\alpha$. Default value: alpha_sigma = 0.5.
beta_sigma	Prior parameter $\beta$. Default value: beta_sigma = 0.5.
q	A vector of strictly positive integers, containing the number of factors to be fitted.
normalize	Logical value indicating whether the observed data will be normalized. Default value: TRUE (recommended)
nIterPerCycle	Number of iteration per MCMC cycle. Default value: 10.
warm_up_overfitting	Number of iterations for the overfitting initialization scheme. Default value: 500.
warm_up	Number of iterations that will be used to initialize the models before starting proposing switchings. Default value: 5000.
overfittingInitialization	Logical value indicating whether the chains are initialized via the overfitting initialization scheme. Default: TRUE (recommended).
rmDir	Logical value indicating whether to delete the outDir directory. Default: TRUE.
parallelModels	Model-level parallelization: An optional integer specifying the number of cores that will be used in order to fit in parallel each member of model. Default: NULL (no model-level parallelization).

Table 2. Table 2: Output returned to the user of the fabMix() function.

Object	Description
bic	Bayesian Information Criterion per model and number of factors.
class	The estimated single best clustering of the observations according to the selected model.
n_Clusters_per_model	The most probable number of clusters (number of non-empty components of the overfitting mixture) per model and number of factors
posterior_probability	The posterior probability of the estimated allocations according to the selected model
covariance_matrix	The estimated posterior mean of the covariance matrix per cluster according to the selected model
mu	The estimated posterior mean of the mean per cluster according to the selected model
weights	The estimated posterior mean of the mixing proportions according to the selected model
mcmc	A list containing the MCMC draws for the parameters of the selected model.
Kmap_prob	The posterior probability of the Maximum A Posteriori number of alive clusters for each parameterization and factor level.

Table 3. Table 3: Selected number of clusters, factors, parameterization and adjusted rand index for simulated data 1 and 2.

Data $(K,q)$	fabMix				pgmm
	$\hat{K}$	$\hat{q}$	model	ARI	$\hat{K}$	$\hat{q}$	model	ARI
1 $(6,2)$	6	2	UUC	1	7	2	UUC	.95
2 $(2,3)$	2	2	UUC	.98	2	2	CUC	.88

Table 4. Table 4: Selected number of clusters, factors, parameterization and adjusted rand index for the publicly available data.

Data (K)	fabMix				pgmm
	$\hat{K}$	$\hat{q}$	model	ARI	$\hat{K}$	$\hat{q}$	model	ARI
Coffee (2)	2	1	CUU	1	4	4	CUU	.29
Wave (3)	3	1	UCU	.61	3	1	UCU	.61
Wine (3)	5	4	CUU	.83	3	4	CUU	.97
Yeast (5)	5	6	CUU	.50	20	10	CUC	.20

Table 5. Table B.5: 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 6. Table B.6: Setup for our simulation scenarios with unequal cluster sizes. In all cases the dimensionality of the multivariate data is equal to p = 30 𝑝 30 p=30 .

	Scenario 1	Scenario 2
True number of clusters $\left(K\right)$	$5$	$2$
True number of factors $\left(q\right)$	$2$	$3$
True mixing proportions $\left(𝒘\right)$	$(\frac{1}{15},\frac{2}{15},\frac{3}{15},\frac{4}{15},\frac{5}{15})$	$(\frac{1}{20},\frac{19}{20})$
Sample size ($n$)	$50,100,200,300,400,500$	$50,100,200,300,\dots ,1500$

Equations

$${\mathbf{x}}_i$$

$$({\mathbf{y}}_i,{\mathbf{\epsilon }}_i)$$

$$\mathbf{\Sigma }$$

$${\mathbf{x}}_i|{\mathbf{y}}_i$$

$${\mathbf{x}}_i\sim {\mathcal{N}}_p(\mathbf{\mu },\mathbf{\Lambda }{\mathbf{\Lambda }}^T+\mathbf{\Sigma }),\text{mbox}iidfori=1,\dots ,n.$$

$${\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),\text{mbox}iidi=1,\dots ,n$$

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

$$({\mathbf{y}}_i,{\mathbf{\epsilon }}_i|z_i=k)\sim {\mathcal{N}}_q(0,{\mathrm{I}}_q){\mathcal{N}}_p(0,{\mathbf{\Sigma }}_k),$$

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

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

$${\mathbf{\Lambda }}_1=\dots ={\mathbf{\Lambda }}_K=\mathbf{\Lambda }.$$

$${\mathbf{\Sigma }}_1=\dots ={\mathbf{\Sigma }}_K=\mathbf{\Sigma }.$$

$${\mathbf{\Sigma }}_k$$

$${\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{x}}_i\sim {\sum }_{k=1}^K{w}_k{\mathcal{N}}_p({\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k{\mathbf{\Lambda }}_k^T+\mathbf{\Sigma })$$

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

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

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

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

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

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

$$\mathbf{\Lambda }=\left(\begin{array}{cccc}{\lambda }_{11}& 0& \cdots & 0\\ {\lambda }_{21}& {\lambda }_{22}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\lambda }_{q1}& {\lambda }_{q2}& \cdots & {\lambda }_{qq}\\ ⋮& ⋮& \ddots & ⋮\\ {\lambda }_{p1}& {\lambda }_{p2}& \cdots & {\lambda }_{pq}\end{array}\right).$$

$$\mathbf{w}$$

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

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

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

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

$$n_k$$

$${\mathbf{A}}_k$$

$${\mathbf{B}}_k$$

$${\mathbf{\tau }}_{kr}$$

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

$$s_{kr}$$

$$\mathbf{T}$$

$${\mathbf{M}}_k$$

$$\mathrm{P}(z_i=k|\cdots )\propto w_{k}f\left({\mathbf{x}}_i;{\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k{\mathbf{\Lambda }}_k^T+{\mathbf{\Sigma }}_k\right),$$

$$\mathbf{w}|\cdots \sim \mathcal{D}(\gamma +n_1,\dots ,\gamma +n_K).$$

$${\mathbf{\mu }}_k|\cdots \sim {\mathcal{N}}_p\left({\mathbf{A}}_k^{-1}{\mathbf{B}}_k,{\mathbf{A}}_k^{-1}\right),$$

$${\mathbf{\Lambda }}_{kr\cdot }|\cdots \sim {\mathcal{N}}_{{\nu }_r}\left({\left[{\mathbf{\Omega }}^{-1}+{\mathbf{\Delta }}_{kr}\right]}^{-1}{\mathbf{\tau }}_{kr},{\left[{\mathbf{\Omega }}^{-1}+{\mathbf{\Delta }}_{kr}\right]}^{-1}\right),$$

$${\mathbf{\Lambda }}_{r\cdot }|\cdots \sim {\mathcal{N}}_{{\nu }_r}\left({\left[{\mathbf{\Omega }}^{-1}+{\mathbf{\Delta }}_{\bullet r}\right]}^{-1}{\mathbf{\tau }}_{\bullet r},{\left[{\mathbf{\Omega }}^{-1}+{\mathbf{\Delta }}_{\bullet r}\right]}^{-1}\right),$$

Full text

1 Introduction

Factor Analysis (FA) explains relationships among a set of observed variables using a set of latent variables. This is typically achieved by expressing the observed multivariate data as a linear combination of a smaller set of unobserved and uncorrelated variables known as factors. 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$. Let $\mathcal{N}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ denotes the $p$-dimensional normal distribution with mean $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$ and also denote by $\boldsymbol{\mathrm{I}}_{p}$ the $p\times p$ identity matrix. The following equations summarize the typical FA model.

[TABLE]

Before proceeding note that we are not differentiating the notation between random variables and their corresponding realizations. Bold upper-case letters are used for matrices, bold lower-case letters are used for vectors and normal text for scalars.

In Equation (1) we assume that $\boldsymbol{x}_{i}$ is expressed as a linear combination of a latent vector of factors $\boldsymbol{y}_{i}\in\mathbb{R}^{q}$. The $p\times q$ dimensional matrix $\boldsymbol{\Lambda}=(\lambda_{rj})$ contains the factor loadings, while $\boldsymbol{\mu}=(\mu_{1},\ldots,\mu_{p})$ contains the marginal mean of $\boldsymbol{x}_{i}$. The unobserved 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}$ as shown in Equation (2), where $\boldsymbol{0}$ denotes a vector of zeros. Note that the error terms $\boldsymbol{\varepsilon}_{i}$ are independent from $\boldsymbol{y}_{i}$. Furthermore, the errors are consisting of independent random variables $\varepsilon_{i1},\ldots,\varepsilon_{ip}$, as implied by the diagonal covariance matrix $\boldsymbol{\Sigma}$ in Equation (3). 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. The previous assumptions lead to

[TABLE]

According to Equation (5), the covariance matrix of the marginal distribution of $\boldsymbol{x}_{i}$ is equal to $\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{T}+\boldsymbol{\Sigma}$. This is the crucial characteristic of factor analytic models, where they aim to explain high-dimensional dependencies using a set of lower-dimensional uncorrelated factors (Kim and Mueller, 1978; Bartholomew et al., 2011).

Mixtures of Factor Analyzers (MFA) are generalizations of the typical FA model, by assuming that Equation (5) becomes

[TABLE]

where $K$ denotes the number of mixture components. The vector of mixing proportions $\boldsymbol{w}:=(w_{1},\ldots,w_{K})$ contains the weight of each component, with $0\leqslant w_{k}\leqslant 1$; $k=1,\ldots,K$ and $\sum_{k=1}^{K}w_{k}=1$. Note that the mixture components are characterized by different parameters $\boldsymbol{\mu}_{k},\boldsymbol{\Lambda}_{k},\boldsymbol{\Sigma}_{k}$, $k=1,\ldots,K$. Thus, MFAs are particularly useful when the observed data exhibits unusual characteristics such as heterogeneity. That being said, this approach aims to capture the behaviour of each cluster within a component of the mixture model. A comprehensive perspective on the history and development of MFA models is given in Chapter 3 of the monograph by McNicholas (2016).

Early works applying the Expectation-Maximization (EM) algorithm (Dempster et al., 1977) for estimating MFA are the ones from Ghahramani et al. (1996); Tipping and Bishop (1999); McLachlan and Peel (2000). McNicholas and Murphy (2008, 2010) introduced the family of parsimonious Gaussian mixture models (PGMM) by considering the case where the factor loadings and/or error variance may be shared or not between the mixture components. These models are estimated by the alternating expectation-conditional maximization algorithm (Meng and Van Dyk, 1997) and have superior performance compared to other approaches (McNicholas and Murphy, 2008). Under a Bayesian setup, Fokoué and Titterington (2003) estimate the number of mixture components and factors by simulating a continuous-time stochastic birth-death point process using a Birth-Death MCMC algorithm (Stephens, 2000). More recently, Papastamoulis (2018b) estimated Bayesian MFA models with an unknown number of components using overfitting mixtures.

In recent years there is a growing progress on the usage of overfitting mixture models in Bayesian analysis (Rousseau and Mengersen, 2011; van Havre et al., 2015). An overfitting mixture model consists of a number of components which is much larger than its true (and unknown) value. Under suitable prior assumptions (see Appendix A) introduced by Rousseau and Mengersen (2011), 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.

Other Bayesian approaches to estimate the number of components in a mixture model include the Reversible jump MCMC (RJMCMC) (Green, 1995; Richardson and Green, 1997; Dellaportas and Papageorgiou, 2006; Papastamoulis and Iliopoulos, 2009), Birth-death MCMC (BDMCMC) (Stephens, 2000) and allocation sampling (Nobile and Fearnside, 2007; Papastamoulis and Rattray, 2017) algorithms. However, overfitting mixture models are straightforward to implement, while the rest approaches require either careful design of various move types that bridge models with different number of clusters, or analytical integration of parameters.

The overall message is that there is a need for developing an efficient Bayesian method that will combine the previously mentioned frequentist advances on parsimonious representations of MFAs and the flexibility provided by the Bayesian viewpoint. This study aims at filling this gap by extending the Bayesian method of Papastamoulis (2018b) to the family of parsimonious Gaussian mixtures of McNicholas and Murphy (2008). Furthermore, we illustrate the proposed method using the R (Ihaka and Gentleman, 1996; R Core Team, 2016) package fabMix (Papastamoulis, 2018a) available as a contributed package from the Comprehensive R Archive Network at https://CRAN.R-project.org/package=fabMix. The proposed method efficiently deals with many inferential problems (see e.g. Celeux et al. (2000a)) related to mixture posterior distributions, such as (i) inferring the number of non-empty clusters using overfitting models, (ii) efficient exploration of the posterior surface by running parallel heated chains and (iii) incorporating advanced techniques that succesfully deal with the label switching issue (Papastamoulis, 2016).

The rest of the paper is organized as follows. Section 2 reviews the basic concepts of parsimonious MFAs. Identifiability problems and corresponding treatments are detailed in Section 2.1. The Bayesian model is introduced in Section 2.2. Section 3 presents the full conditional posterior distributions of the model. The MCMC algorithm is described in Section 3.2. A detailed presentation of the main function of the contributed R package is given in Section 4. Our method is illustrated and compared to similar models estimated by the EM algorithm in Sections 5.1 and 5.2 using an extended simulation study and 4 publicly available datasets, respectively. We conclude in Section 6 with a summary of our findings and directions for further research. An Appendix contains further discussion on overfitting mixture models (Appendix A), details of the MCMC sampler (Appendix B) and additional simulation results (Appendix C).

2 Parsimonious Mixtures of Factor Analyzers

Consider the latent allocation variables $z_{i}$ which assign observation $\boldsymbol{x}_{i}$ to a component $k=1,\ldots,K$ for $i=1,\ldots,n$. A-priori each observation is generated from component $k$ with probability equal to $w_{k}$, that is,

[TABLE]

independent for $i=1,\ldots,n$. Note that the allocation vector $\boldsymbol{z}:=(z_{1},\ldots,z_{n})$ is not observed, so it should be treated as missing data. We assume that $\boldsymbol{z}_{i}$ and $\boldsymbol{y}_{i}$ are independent, thus Equation (2) is now written as:

[TABLE]

and conditional on the cluster membership and latent factors we obtain that

[TABLE]

Consequently,

[TABLE]

independent for $i=1,\ldots,n$. From Equations (7) and (10) we derive that the marginal distribution of $\boldsymbol{x}_{i}$ is the finite mixture model in Equation (6).

Following McNicholas and Murphy (2008), the factor loadings and error variance per component may be common or not among the $K$ components in Equation (6). If the factor are constrained, then:

[TABLE]

If the error variance is constrained, then:

[TABLE]

Furthermore, the error variance may be isotropic (i.e. proportional to the identity matrix) or not and depending on whether contstraint (12) is disabled or enabled:

[TABLE]

We note that under constraint (13), the model is referred to as a mixture of probabilistic principal component analyzers (Tipping and Bishop, 1999).

Depending on whether a particular constraint is present or not, the following set of 8 parameterizations arises.

[TABLE]

independent for $i=1,\ldots,n$. Following the pgmm nomenclature (McNicholas and Murphy, 2008): the first, second and third letter denotes whether $\boldsymbol{\Lambda}_{k}$, $\boldsymbol{\Sigma}_{k}=\mathrm{diag}(\sigma^{2}_{k1},\\ \ldots,\sigma^{2}_{kp})$ and $\sigma^{2}_{kj}$, $k=1,\ldots,K$; $j=1,\ldots,p$, are constrained (C) or unconstrained (U), respectively. A novelty of the present study is to offer a Bayesian framework for estimating the whole family of the previous parameterizations (note that Papastamoulis (2018b) estimated the UUU and UCU parameterizations).

2.1 Label switching and other identifiability problems

Let $L(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{\phi}|\boldsymbol{x})=\prod_{i=1}^{n}\sum_{k=1}^{K}w_{k}f(x_{i}|\theta_{k},\boldsymbol{\phi})$, $(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{\phi})\in\mathcal{P}_{K-1}\times\Theta^{K}\times\Phi$ denote the likelihood function of a mixture of $K$ densities, where $\mathcal{P}_{K-1}$ denotes the parameter space of the mixing proportions $\boldsymbol{w}$, $\boldsymbol{\theta}=(\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{K})$ are the component-specific parameters and $\boldsymbol{\phi}$ denotes a (possibly empty) collection of parameters that are common between all components. For instance, consider the UCU parameterization where $\boldsymbol{\theta}_{k}=(\boldsymbol{\mu}_{k},\boldsymbol{\Lambda}_{k})$ for $k=1,\ldots,K$ and $\boldsymbol{\phi}=\boldsymbol{\Sigma}$. For any permutation $\tau=(\tau_{1},\ldots,\tau_{K})$ of the set $\{1,\ldots,K\}$, the likelihood of mixture models is invariant to permutations of the component labels: $L(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{\phi}|\boldsymbol{x})=L(\tau\boldsymbol{w},\tau\boldsymbol{\theta},\boldsymbol{\phi}|\boldsymbol{x})$. Thus, the likelihood surface of a mixture model with $K$ components will exhibit $K!$ symmetric areas. If $(\boldsymbol{w}^{*},\boldsymbol{\theta}^{*},\boldsymbol{\phi}^{*})$ corresponds to a mode of the likelihood, the same will hold for any permutation $(\tau\boldsymbol{w}^{*},\tau\boldsymbol{\theta}^{*},\boldsymbol{\phi}^{*})$.

Label switching (Redner and Walker, 1984) is the commonly used term to describe this phenomenon. Under a Bayesian point of view, in the case that the prior distribution is also invariant to permutations (which is typically the case, see e.g. Marin et al. (2005); Papastamoulis and Iliopoulos (2013)), the same invariance property will also hold for the posterior distribution $f(\boldsymbol{w},\boldsymbol{\theta},\boldsymbol{\phi}|\boldsymbol{x})$. Consequently, the marginal posterior distributions of mixing proportions and component-specific parameters will be coinciding, i.e.: $f(w_{1}|\boldsymbol{x})=\ldots=f(w_{K}|\boldsymbol{x})$ and $f(\theta_{1}|\boldsymbol{x})=\ldots=f(\theta_{K}|\boldsymbol{x})$. Thus, when approximating the posterior distribution via MCMC sampling, the standard practice of ergodic averages for estimating quantities of interest (such as the mean of the marginal posterior distribution for each parameter) becomes meaningless. In order to deal with this identifiability problem we post-process the simulated MCMC output using a deterministic relabelling algorithm, that is, the Equivalence Classes Representatives (ECR) algorithm (Papastamoulis and Iliopoulos, 2010; Papastamoulis, 2014), as implemented in the R package label.switching (Papastamoulis, 2016).

A second source of identifiability problems is related to orthogonal transformations of the matrix of factor loadings. A popular practice (Geweke and Zhou, 1996; Fokoué and Titterington, 2003; Mavridis and Ntzoufras, 2014; Papastamoulis, 2018b) to overcome this issue, is to preassign values to some entries of $\boldsymbol{\Lambda}$, in particular we set the entries of the upper diagonal of the first $q\times q$ block matrix of $\boldsymbol{\Lambda}$ equal to zero:

[TABLE]

Another problem is related to the so-called “sign switching” phenomenon, see e.g. Conti et al. (2014). Simultaneously switching the signs of a given row $r$ of $\boldsymbol{\Lambda}$; $r=1,\ldots,p$ and $\boldsymbol{y}_{i}$ does not alter the likelihood. Thus, $\boldsymbol{\Lambda}$ and $\boldsymbol{y}_{i}$; $i=1,\ldots,n$ are not marginally identifiable due to sign-switching across the MCMC trace. However, this is not a problem in our implementation, since all parameters of the marginal density of $\boldsymbol{x}_{i}$ in (6) are identified (see also the discussion for sign-invariant parametric functions in Papastamoulis (2018b)).

Parameter expanded approaches are preferred in the recent literature (Bhattacharya and Dunson, 2011; McParland et al., 2017), because the mixing of the MCMC sampler is improved. In our implementation, we are able to obtain excellent mixing using the popular approach of restricting elements of $\boldsymbol{\Lambda}$: the reader is referred to Figure 2 of Papastamoulis (2018b), where it is obvious that our MCMC sampler has the ability to rapidly move between the multiple modes of the target posterior distribution of $\boldsymbol{\Lambda}$ (more details on convergence diagnostics are also presented in Appendix A.4 of Papastamoulis (2018b)).

2.2 Prior assumptions

We assume that the number of mixture components ($K$) has a sufficiently large value so that it overestimates the “true” number of clusters. Unless otherwise stated, the default choice is $K=20$. All prior assumptions of the overfitting mixture models are discussed in detail in Papastamoulis (2018b). For ease of presentation we repeat them in this section. 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 (17) $\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 (19). A graphical representation of the hierarchical model is given in Figure 1 of Papastamoulis (2018b). The default values of the remaining fixed hyper-parameters are given in Appendix B.

The previous assumptions refer to the case of the unconstrained parameter space, that is, the UUU paramaterization. Clearly, they should be modified accordingly when a constrained model is used. Under constraint (11), the prior distribution in Equation (17) becomes $\boldsymbol{\Lambda}_{r\cdot}\sim\mathcal{N}_{\nu_{r}}(\boldsymbol{0},\boldsymbol{\Omega})$, independent for $r=1,\ldots,p$. Under constraints (12) and (13), the prior distribution in Equation (18) becomes $\sigma_{r}^{-2}\sim\mathcal{G}(\alpha,\beta)$, independent for $r=1,\ldots,p$. Finally, under constraints (12) and (14), the prior distribution in Equation (18) becomes $\sigma^{-2}\sim\mathcal{G}(\alpha,\beta)$.

3 Inference

This section describes the full conditional posterior distributions of model parameters and the corresponding MCMC sampler. Due to conjugacy, all full conditional posterior distributions are available in closed forms.

3.1 Full conditional posterior distributions

Let us define the following quantities:

[TABLE]

for $k=1,\ldots,K$; $r=1,\ldots,p$. For a generic sequence of the form $\{G_{rc};r\in\mathcal{R},c\in\mathcal{C}\}$ we also define $G_{\bullet c}=\sum_{r}G_{rc}$ and $G_{r\bullet}=\sum_{c}G_{rc}$. Finally, $(x|\cdots)$ denotes the conditional distribution of $x$ given the value of all remaining variables.

From Equations (6) and (7) it immediately follows that for $k=1,\ldots,K$

[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. Using the Sherman–Morrison–Woodbury formula (see e.g. Hager (1989)), the inverse matrix is equal to $\boldsymbol{\Sigma}_{k}^{-1}-\boldsymbol{\Sigma}_{k}^{-1}\boldsymbol{\Lambda}_{k}\boldsymbol{M}_{k}^{-1}\boldsymbol{\Lambda}_{k}^{T}\boldsymbol{\Sigma}^{-1}_{k}$, for $k=1,\ldots,K$. The full conditional posterior distribution of mixing proportions is a Dirichlet distribution with parameters

[TABLE]

The full conditional posterior distribution of the marginal mean per component is

[TABLE]

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

The full conditional posterior distribution of the factor loadings without any restriction is

[TABLE]

independent for $k=1,\ldots,K$; $r=1,\ldots,p$. Under constraint (11) we obtain that

[TABLE]

independent for $r=1,\ldots,p$.

The full conditional distribution of error variance without any restriction is

[TABLE]

independent for $k=1,\ldots,K$; $r=1,\ldots,p$. Under constraint (12) we obtain that

[TABLE]

independent for $r=1,\ldots,p$. Under constraints (12) and (13) we obtain that

[TABLE]

independent for $k=1,\ldots,K$. Under constraints (12) and (14) we obtain that

[TABLE]

The full conditional distribution of latent factors is given by

[TABLE]

independent for $i=1,\ldots,n$. Finally, the full conditional distribution for $\omega_{\ell}$ is

[TABLE]

while under constraint (11) we obtain that

[TABLE]

independent for $\ell=1,\ldots,q$.

3.2 MCMC sampler

Given the number of factors ($q$) and a model parameterization, a Gibbs sampler (Geman and Geman, 1984; Gelfand and Smith, 1990) coupled with a prior parallel tempering scheme (Geyer, 1991; Geyer and Thompson, 1995; Altekar et al., 2004) is used in order to produce a MCMC sample from the joint posterior distribution. Each heated chain ($j=1,\ldots,\mbox{\tt nChains}$) 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 (van Havre et al., 2015). The inference is based on the output of the first chain ($j=1$) of the prior parallel tempering scheme (van Havre et al., 2015). The number of factors and model parameterization is selected according to the Bayesian Information Criterion (BIC) (Schwarz, 1978), conditional on the most probable number of alive clusters per model (see Papastamoulis (2018b) for a detailed comparison of BIC with other alternatives).

Let $\mathcal{M}$ and $\mathcal{Q}$ denote the set of model parameterizations and number of factors. In the following pseudocode, $x\leftarrow[y|z]$ denotes that $x$ is updated from a draw from the distribution $f(y|z)$ and $\theta_{j}^{(t)}$ denotes the value of $\theta$ at the $t$-th iteration of the $j$-th chain.

For $(m,q)\in\mathcal{M}\times\mathcal{Q}$

1. (a)

Obtain initial values $(\boldsymbol{\Omega}_{j}^{(0)}$, $\boldsymbol{\Lambda}_{m;j}^{(0)}$, $\boldsymbol{\mu}^{(0)}_{j}$, $\boldsymbol{z}^{(0)}_{j}$, $\boldsymbol{\Sigma}_{m;j}^{(0)}$, $\boldsymbol{w}^{(0)}_{j}$, $\boldsymbol{y}^{(0)}_{j})$ by running the overfitting initialization scheme, for $j=1,\ldots,\mbox{\tt nChains}$. 2. (b)

For MCMC iteration $t=1,2,\ldots$ update

1. i.

For chain $j=1,\ldots,\mbox{\tt nChains}$

1. A.

$\boldsymbol{\Omega}^{(t)}_{j}\leftarrow\left[\boldsymbol{\Omega}|\boldsymbol{\Lambda}_{mj}^{(t-1)}\right]$.

If $m\in\{\mbox{UUU,UCU,UUC, UCC}\}$ use (30)

else use (31). 2. B.

$\boldsymbol{\Lambda}_{m;j}^{(t)}\leftarrow\left[\boldsymbol{\Lambda}|\boldsymbol{\Omega}_{j}^{(t)},\boldsymbol{\mu}_{j}^{(t-1)},\boldsymbol{\Sigma}_{m;j}^{(t-1)},\boldsymbol{x},\boldsymbol{y}_{j}^{(t-1)},\boldsymbol{z}_{j}^{(t-1)}\right]$

If $m\in\{\mbox{UUU,UCU,UUC, UCC}\}$ use (23)

else use (24). 3. C.

$\boldsymbol{\mu}_{j}^{(t)}\leftarrow\left[\boldsymbol{\mu}|\boldsymbol{\Lambda}_{m}^{(t)},\boldsymbol{\Sigma}_{m}^{(t-1)},\boldsymbol{x},\boldsymbol{y}^{(t-1)},\boldsymbol{z}_{j}^{(t-1)}\right]$ according to (22). 4. D.

$\boldsymbol{z}_{j}^{(t)}\leftarrow\left[\boldsymbol{z}|\boldsymbol{w}_{j}^{(t-1)},\boldsymbol{\mu}_{j}^{(t)},\boldsymbol{\Lambda}_{m;j}^{(t)},\boldsymbol{\Sigma}_{m;j}^{(t-1)},\boldsymbol{x}\right]$ according to (20). 5. E.

$\boldsymbol{w}_{j}^{(t)}\leftarrow\left[\boldsymbol{w}|\boldsymbol{z}_{j}^{(t)}\right]$ according to (21) with prior parameter $\gamma=\gamma_{(j)}$. 6. F.

$\boldsymbol{\Sigma}_{m;j}^{(t)}\leftarrow\left[\boldsymbol{\Sigma}|\boldsymbol{x},\boldsymbol{z}_{j}^{(t)},\boldsymbol{\mu}_{j}^{(t)},\boldsymbol{\Lambda}_{m;j}^{(t)},\boldsymbol{y}_{j}^{(t-1)}\right]$

If $m\in\{\mbox{UUU,CUU}\}$ use (25)

else if $m\in\{\mbox{UCU,CCU}\}$ use (26)

else if $m\in\{\mbox{UUC,CUC}\}$ use (27)

else use (28). 7. G.

$\boldsymbol{y}_{j}^{(t)}\leftarrow\left[\boldsymbol{y}|\boldsymbol{x},\boldsymbol{z}_{j}^{(t)},\boldsymbol{\mu}_{j}^{(t)},\boldsymbol{\Sigma}_{m;j}^{(t)},\boldsymbol{\Lambda}_{m;j}^{(t)}\right]$ according to (29). 2. ii.

Select randomly $1\leqslant j^{*}\leqslant\mbox{\tt nChains}-1$ and propose to swap the states of chains $j^{*}$ and $j^{*}+1$. 3. (c)

For chain $j=1$ compute BIC conditionally on the most probable number of alive clusters. 2. 2.

Select the best $(m,q)$ model corresponding to chain $j=1$ according to BIC and reorder the simulated output of the selected model according to ECR algorithm, conditional on the most probable number of alive clusters.

The MCMC algorithm is initialized using random starting values arising from the “overfitting initialization” procedure introduced by Papastamoulis (2018b). For further details on steps 1.(a) (MCMC initialization) and 1.(b).ii (prior parallel tempering scheme) the reader is referred to Appendix B (see also Sections 2.6, 2.7 and 2.9 of Papastamoulis (2018b)).

4 Using the fabMix package

The main function of the fabMix package is fabMix(), with its arguments shown in Table 1. This function takes as input a matrix rawData of observed data where rows and columns correspond to observations and variables of the dataset, respectively. The parameters of the Dirichlet prior distribution ($\gamma_{(j)};j=1,\ldots,\mbox{nChains}$) of the mixing proportions are controlled by dirPriorAlphas. The range for the number of factors is specified in the q argument. Valid input for q is any positive integer vector between 1 and the Leddermann bound (Ledermann, 1937) implied by the number of variables in the dataset. By default, all 8 parameterizations are fitted, however the user can specify in model any non-empty subset of them.

The fabMix() function simulates a total number of $\mbox{{\tt nChains}}\times\mbox{{\tt length(models)}}\times\mbox{{\tt length(q)}}$ MCMC chains. For each parameterization and number of factors, the (nChains) heated chains are processed in parallel while swaps between pairs of chains are proposed. Parallelization is possible in the parameterization level as well, using the argument parallelModels. This means that parallelModels are running in parallel where each one of them runs nChains chains in parallel, provided that the number of available threads is at least equal to $\mbox{{\tt nChains}}\times\mbox{{\tt parallelModels}}$. In order to parallelize our code, the doParallel (Revolution Analytics and Steve Weston, 2015), foreach (Revolution Analytics and Steve Weston, 2014) and doRNG (Gaujoux, 2018) packages are imported.

The prior parameters $g,h,\alpha,\beta$ in Equations (18) and (19) correspond to g, h,alpha_sigma and beta_sigma, respectively, with a (common) default value equal to $0.5$. It is suggested to run the algorithm using normalize = TRUE, in order to standardize the data before running the MCMC sampler. The default behaviour of our method is to normalize the data, thus, all reported estimates refer to the standardized dataset. In the case that the most probable number of mixture components is larger than 1, the ECR algorithm is applied in order to undo the label switching problem. Otherwise, the output is post-processed so that the generated parameters of the (single) alive component are switched to the first component of the overfitting mixture.

The sampler will first run for warm_up iterations before starting to propose swaps between pairs of chain. By default, this stage consists of 5000 iterations. After that, each chain will run for a series of mCycles MCMC cycles, each one consisting of nIterPerCycle MCMC iterations (steps A, B, $\ldots$, G of the pseudocode). The updates of factors loadings according to (23) and (24) at step B of the pseudocode are implemented using object-oriented programming using the Rcpp and RcppArmadillo libraries (Eddelbuettel and François, 2011; Eddelbuettel and Sanderson, 2014). At the end of each cycle, a swap between a pair of chains is proposed.

Obviously, the total number of MCMC iterations is equal to $\mbox{{\tt warm\_up}}+\mbox{{\tt mCycles}}\times\mbox{{\tt nIterPerCycle}}$ and the first $\mbox{{\tt warm\_up}}+\mbox{{\tt burnCycles}}\times\mbox{{\tt nIterPerCycle}}$ iterations are discarded as burn-in. Given the default values of nIterPerCycle, warm_up and overfittingInitialization, choices between $50\leqslant\mbox{{\tt burnCycles}}\leqslant 500<\mbox{{\tt mCycles}}\leqslant 1500$ are typical in our implementation (see also the convergence analysis in Papastamoulis (2018b)).

While the function runs, some basic information is printed either on the screen (if parallelModels is not enabled) or in separate text files inside the output folder (in the opposite case), such as the progress of the sampler as well as the acceptance rate of proposed swaps between chains. The output which is returned to the user is detailed in Table 2. The full MCMC output of the selected model is returned as a list (named as mcmc) consisting of mcmc objects, a class imported from the coda package (Plummer et al., 2006). In particular, mcmc consists of the following:

• $\bullet$

y: object of class mcmc containing the simulated factors.

• $\bullet$

w: object of class mcmc containing the simulated mixing proportions of the alive components, reordered according to ECR algorithm.

• $\bullet$

Lambda: list containing objects of class mcmc with the simulated factor loadings of the alive components, reordered according to ECR algorithm. Note that this particular parameter is not identifiable due to sign-switching across the MCMC trace.

• $\bullet$

mu: list containing objects of class mcmc with the simulated marginal means of the alive components, reordered according to ECR algorithm.

• $\bullet$

z: matrix of the simulated latent allocation variables of the mixture model, reordered according to ECR algorithm.

• $\bullet$

Sigma: list containing objects of class mcmc with the simulated variance of errors of the alive components, reordered according to ECR algorithm.

• $\bullet$

K_all_chains: matrix of the simulated values of the number of alive components per chain.

The user can call the print, summary and plot methods of the package in order to easily retrieve and visualize various summaries of the output, as exemplified in the next section.

5 Examples

This section illustrates our method. At first we demonstrate a typical implementation on two single simulated datasets and explain in detail the workflow. Then we perform an extensive simulation study for assesing the ability of the proposed method to recover the correct clustering and compare our findings to the pgmm package (McNicholas and Murphy, 2008, 2010; McNicholas et al., 2010, 2015). Application to four publicly available datasets is provided next.

5.1 Simulation study

We simulated synthetic data of $p=30$ variables consisting of $n=300$ observations and $K=6$ clusters (dataset 1) and $n=200$, $K=2$ (dataset 2), as shown in Figure 1. Both of them were generated using MFA models with $q=2$ (dataset 1) and $q=3$ (dataset 2) factors. The two datasets exhibit different characteristics: the variance of errors per cluster ($\boldsymbol{\Sigma}_{k}$) is significantly larger in dataset 2 compared to dataset 1. In addition, the selection of factor loadings in dataset 2 result to more complex covariance structure. The generating mechanism, described in detail in Papastamoulis (2018b), is available in the fabMix package via the simData() and simData2() functions, as shown below.

library(’fabMix’)

dataset 1

set.seed(1) n = sample(100*(1:10), 1) # sample size q = sample(1:3, 1) # number of factors K = sample(1:10, 1) # number of clusters

results to n = 300, q = 2, K = 6

p = 30 # number of variables

inverse variance of errors

sINV <- array(data = NA, dim = c(K,p)) for(k in 1:K){sINV[k,] <- 1/(1+20*log(k+1))} dataset1 <- simData(sameSigma=FALSE, K.true=K,

• n = n, q = q, p = p, sINV_values=sINV)

synthetic dataset 2

set.seed(30) n = 200; q = 3; K = 2; p = 30

inverse variance of errors

sINV <- array(data = NA, dim = c(K,p)) for(k in 1:K){sINV[k,] <- 1/(1+1000*log(k+1))} dataset2 <- simData2(sameSigma=FALSE, K.true=K,

• n = n, q = q, p = p, sINV_values = sINV)

Next we estimate the 8 overfitting Bayesian MFA models with $K_{\mbox{max}}=20$ mixture components assuming that the number of factors ranges in the set $1\leqslant q\leqslant 5$. The MCMC sampler runs nChains = 4 heated chains, each one consisting of mCycles = 700 cycles, while the first burnCycles = 100 are discarded. Recall that each MCMC cycle consists of nIterPerCycle = 10 usual MCMC iterations and that there is an additional warm up period of the MCMC sampler (before starting to propose chain swaps) corresponding to 5000 usual MCMC iterations.

Kmax <- 20 # number of components nChains <- 4 # number of chains qRange <- 1:5 # number of factors

Run fabMix() for dataset 1

set.seed(1)

fm1 <- fabMix(nChains = nChains,

• rawData = dataset1$data, outDir = "tmp1",
• Kmax = Kmax, mCycles = 700, burnCycles = 100,
• q = qRange, parallelModels = 4)

Run fabMix() for dataset 2

set.seed(1)

fm2 <- fabMix(nChains = nChains,

• rawData = dataset2$data, outDir = "tmp2",
• Kmax = Kmax, mCycles = 700, burnCycles = 100,
• q = qRange, parallelModels = 4)

The argument parallelModels = 4 implies that 4 parameterizations will be processed in parallel. In addition, each model will use nChains = 4 threads to run in parallel the specified number of chains. Our jobscript used 16 threads so in this case the $\mbox{{\tt parallelModels}}\times{{\tt nChains}}=16$ jobs are efficiently allocated.

5.1.1 Methods for printing, summarizing and plotting the output

The print method for a fabMix.object displays some basic information for a given run of the fabMix function. The following output corresponds to the first dataset.

print(fm1)

• Run information: Number of fitted models: (5 factor levels) x (8 parameterizations) = 40 models. Selected model: UUC model with K = 6 clusters and q = 2 factors.

• Maximum A Posteriori (MAP) number of ‘alive’ clusters and selected number of factors (BIC) per model: model K_MAP K_MAP_prob q BIC_q chain_swap 1 UUU 4 1.00 3 4295.0 8.43% 2 CUU 5 0.70 3 4248.9 13.29% 3 UCU 7 0.71 2 3345.8 25.86% 4 CCU 13 0.70 2 3410.5 93.71% 5 UCC 7 0.55 2 3270.7 22.71% 6 UUC 6 1.00 2 2274.5 19.43% 7 CUC 10 0.41 3 2868.3 78.43% 8 CCC 13 0.52 2 3378.3 92.43%

• Estimated number of observations per cluster (selected model): label 4 7 13 14 15 17 60 55 41 72 50 22

The following output corresponds to the print method for the fabMix function for the second dataset.

print(fm2)

• Run information: Number of fitted models: (5 factor levels) x (8 parameterizations) = 40 models. Selected model: UUC model with K = 2 clusters and q = 2 factors.

• Maximum A Posteriori (MAP) number of ‘alive’ clusters and selected number of factors (BIC) per model: model K_MAP K_MAP_prob q BIC_q chain_swaps 1 UUU 2 1.00 2 14776.6 4.86% 2 CUU 2 1.00 2 14676.0 3% 3 UCU 2 0.64 3 15043.1 4.57% 4 CCU 3 0.48 3 14836.2 12.71% 5 UCC 3 0.64 2 15141.9 4.29% 6 UUC 2 0.95 2 14558.5 3.14% 7 CUC 3 0.85 3 14598.7 4.14% 8 CCC 4 0.44 3 14851.6 11.14%

• Estimated number of observations per cluster (selected model): label 6 20 113 87

We conclude that the selected models correspond to the UUC parameterization with $K=6$ clusters and $q=2$ factors for dataset 1 and $K=2$, $q=2$ for dataset 2. The selected number of clusters and factors for the whole range of 8 models is displayed next, along with the estimated posterior probability of the number of alive clusters per model (K_MAP_prob), the value of the BIC for the selected number of factors (BIC_q) as well as the proportion of the accepted swaps between the heated MCMC chains in the last column. The frequency table of the estimated single best clustering of the datasets is displayed in the last field. We note that the labels of the frequency table correspond to the labels of the alive components of the overfitting mixture model, that is, components 4, 7, 13, 14, 15, and 17 for dataset 1 and components 6 and 20 for dataset 2. Clearly, these labels can be renamed to $1,2,3,4,5,6$ and $1,2$ respectively, but we prefer to retain the raw output of the sampler as a reminder of the fact that it corresponds to the alive components of the overfitted mixture model.

The summary method of the fabMix package summarizes the MCMC output for the selected model by calculating posterior means and quantiles for the mixing proportions, marginal means and the covariance matrix per (alive) cluster. A snippet of the output for dataset 2 is shown below.

s <- summary(fm2)

• ‘Alive’ cluster labels: [1] "6" "20"

• Posterior mean of the mixing proportions: 6 20 0.58 0.42

• Posterior mean of the marginal means: Cluster label Variable 6 20 V1 -0.06 0.08 V2 -0.02 0.02 ............... V30 -0.01 0.00

• Posterior mean of the covariance matrix:

  Covariance matrix for cluster ‘6’: V1 V2 ... V30 V1 1.12 0.53 ... -0.18 V2 0.53 1.11 ... -0.14 .......................... V30 -0.18 -0.14 ... 1.27

  Covariance matrix for cluster ‘20’: V1 V2 ... V30 V1 0.66 0.39 ... -0.05 V2 0.39 0.88 ... -0.03 .......................... V30 -0.05 -0.03 ... 0.57

Quantiles for each parameter: quantile parameter 2.5% 25% 50% 75% 97.5% weight_6 0.51 0.55 0.58 0.51 0.65 weight_20 0.35 0.40 0.42 0.45 0.50 mean_6_V1 -0.27 -0.13 -0.06 0.01 0.13 mean_20_V1 -0.09 0.03 0.08 0.14 0.25 mean_6_V2 -0.22 -0.09 -0.02 0.04 0.15 mean_20_V2 -0.17 -0.05 0.02 0.08 0.19 ........................................... mean_6_V30 -0.21 -0.08 -0.01 0.07 0.20 mean_20_V30 -0.16 -0.06 0.00 0.06 0.17 cov_6_V1_V1 0.89 1.03 1.10 1.19 1.46 cov_20_V1_V1 0.50 0.59 0.64 0.72 0.85 cov_6_V1_V2 0.37 0.46 0.52 0.59 0.78 cov_20_V1_V2 0.24 0.32 0.38 0.44 0.56 ........................................... cov_6_V1_V30 -0.37 -0.24 -0.18 -0.12 0.01 cov_20_V1_V30 -0.16 -0.08 -0.05 -0.01 0.06 ........................................... cov_6_V2_V30 -0.34 -0.21 -0.13 -0.08 0.04 cov_20_V2_V30 -0.17 -0.07 -0.03 0.03 0.12 ........................................... cov_6_V30_V30 1.03 1.18 1.26 1.37 1.62 cov_20_V30_V30 0.45 0.51 0.56 0.61 0.73

The printed output is also returned to the user via s$posterior_means and s$quantiles.

The plot() method of the package generates the following types of graphics output:

1. (1)

Plot of the BIC values per factor level and parameterization. 2. (2)

Plot of the posterior means of marginal means ($\boldsymbol{\mu}_{k}$) per (alive) cluster and Highest Density Intervals of the corresponding normal distribution along with its assigned data. 3. (3)

The coordinate projection plot of the mclust package (Fraley and Raftery, 2002; Fraley et al., 2012), that is, a scatterplot of the assigned data per cluster for each pair of variables. 4. (4)

Visualization of the posterior mean of the correlation matrix per cluster using the corrplot package. 5. (5)

The MAP estimate of the factor loadings ($\boldsymbol{\Lambda}_{k}$) per (alive) cluster.

The following commands produce plot (1) for datasets 1 and 2.

plot(fm1, what = ’BIC’) plot(fm2, what = ’BIC’)

The produced plots are shown in Figure 2. Note that each point in the plot is labeled by an integer, which corresponds to the MAP number of alive components for the specific combination of factors and parameterization.

The following commands produce plot (2) for datasets 1 and 2.

plot(fm1, what = ’classification_matplot’,

• class_mfrow = c(3,2), confidence = 0.95)

plot(fm2, what = ’classification_matplot’,

• class_mfrow = c(2,1), confidence = 0.95)

The created plots are shown in Figure 3. The class_mfrow arguments control the rows and columns of the layout and it should consists of 2 integers with their product equal to the selected number of (alive) clusters. In addition, a legend is placed on the bottom of the layout. The value(s) in the confidence argument draws the Highest Density Interval(s) of the estimated normal distribution. Note that these plots display the original and not the scaled dataset which is used in the MCMC sampler. Therefore, the central curve and confidence limits displayed in the specific plot correspond to the mean and variance (multiplied by the appropriate quantile of the standard Normal distribution) of the random variables arising by applying the inverse of the z-transformation on the MCMC estimates reported by the fabMix function.

Figure 4 visualizes the correlation matrix for the first cluster of each dataset, using the corrplot package. The argument $\mbox{{\tt sig\_correlation}}=\alpha$ is used for marking cases where the equally tailed $(1-\alpha)$ Bayesian credible interval contains zero. The following commands generate the plots in Figure 4.

plot(fm1, what = ’correlation’, +Ψsig_correlation = 0.05) plot(fm2, what = ’correlation’,

• sig_correlation = 0.05)

5.1.2 Assessing clustering accuracy and comparison with pgmm

In this section we compare our findings against the ground-truth in simulated datasets and also compare against the pgmm package, considering the same range of clusters and factors per dataset. For each combination of number of factors, components and parameterization, the pgmmEM() algorithm was initialized using 3 random starting values as well as the K-means clustering algorithm, that is, 4 different starts in total. Note that the number of different starts of the EM algorithm is set equal to number of parallel chains in the MCMC algorithm. The input data is standardized in both algorithms.

As shown in Table 3, the adjusted Rand index (ARI) (Rand, 1971) between fabMix and the ground-truth classification is equal to 1 and 0.98 for simulated dataset 1 and 2, respectively. The corresponding ARI for pgmm equals to 0.98 and 0.88, respectively. In both cases our method finds the correct number of clusters, however pgmm overestimates $K$ in dataset 1. Both methods select the UUC parameterization in dataset 1, but in dataset 2 different models are selected (UUC by fabMix and CUC by pgmm).

The selected number of factors equals 2, however in dataset 2 the “true” number of factors equals 3. The underestimation of the number of factors in dataset 2 remains true for a wide range of similar data: in particular we generated synthetic datasets with identical parameter values as the ones in dataset 2 but each time the sample size was increasing by 200 observations. We observed that the correct number of factors is returned when $n\geqslant 1600$ for fabMix and $n\geqslant 1800$ for pgmm.

Next we replicate the two distinct simulation procedures (according to the simData() and simData2() functions of the package) used to generate the previously described datasets, but considering that $1\leqslant K\leqslant 10$ (true number of clusters) and $1\leqslant q\leqslant 3$ (true number of factors). The number of variables remains the same as before, that is, $p=30$ and the sample size is drawn uniformly at random in the set $\{100,200,\ldots,1000\}$. We will use the terms ’Scenario 1’ and ’Scenario 2’ to label the two different simulation procedures. In Scenario 1 the diagonal of the variance of errors is generated as $\sigma_{kr}^{2}=1+20\log(k+1)$, $r=1,\ldots,p$, whereas in Scenario 2: $\sigma_{kr}^{2}=1+u_{r}\log(k+1)$, where $u_{r}\sim\mbox{Uniform}(500,1000)$, $r=1,\ldots,p$; $k=1,\ldots,K$. In general, Scenario 1 generates datasets with well separated clusters. On the other hand, the amount of error variance in Scenario 2 makes the clusters less separated. For a given simulated dataset with $K_{\mbox{true}}$ clusters and $q_{\mbox{true}}$ factors, we are considering that the total number of components in the overfitting mixture model (fabMix) as well as the maximum number of components fitted from pgmm is set equal to $K_{\mbox{max}}=K_{\mbox{true}}+6$ and that the number of factors ranges between $1\leqslant q\leqslant q_{\mbox{true}}+2$. These bounds are selected in order to speed up computation time without introducing any bias in the resulting inference (as confirmed by a smaller pilot study). For each scenario 500 datasets were simulated.

The main findings of the simulation study are illustrated in Figure 5. Note that in Scenario 1 fabMix almost always finds the correct clustering structure: the boxplots of the adjusted Rand Index are centered at 1 and, on the second row, the boxplots of the estimated number of clusters are centered at the corresponding true value. On the other hand, observe that for $K\geqslant 6$ pgmm has the tendency to overestimate the number of clusters. In the more challenging Scenario 2 the estimates of the number of cluster exhibit larger variability. However note that for $K=8,9,10$ the number of clusters selected by fabMix is closer to the true value than pgmm, a fact which is also reflected in the ARI where fabMix tends to have larger values than pgmm. For both scenarios, the estimation of the number of factors is in strong agreement between the two methods, as shown in the third row of Figure 5. In the last row, the selected parameterization is shown. Observe that the results are fairly consistent between the two methods.

Finally, we note that in the presented simulation study, the generated clusters have equal sizes (on average). The reader is referred to Appendix C for exploring the performance of the compared methods in the presence of small and large clusters with respect to the size of the available data ($n$).

5.2 Publicly available datasets

In this section we analyze 4 publicly available datasets: a subset of the wave dataset (Breiman et al., 1984; Lichman, 2013) available at the fabMix package, the wine dataset (Forina et al., 1986) available at the pgmm package, the coffee dataset (Streuli, 1973) available at the pgmm package, and the standardized yeast cell cycle data (Cho et al., 1998) available at http://faculty.washington.edu/kayee/model/. Note that Papastamoulis (2018b) analyzed the first three datasets but only considering the UUU and UCU parameterizations for fabMix.

The coffee dataset consists of $n=43$ coffee samples of $p=12$ variables, collected from beans corresponding to the Arabica and Robusta species (thus, $K=2$). The wave dataset consists of a randomly sampled subset of 1500 observations from the wave dataset (Breiman et al., 1984), available from the UCI machine learning repository (Lichman, 2013). According to the available ground-truth classification of the dataset, there are 3 equally weighted underlying classes of $21$-dimensional continuous data. The wine dataset (Forina et al., 1986), available at the pgmm package (McNicholas et al., 2015), contains $p=27$ variables measuring chemical and physical properties of $n=178$ wines, grouped in three types (thus, $K=3$). The reader is referred to McNicholas and Murphy (2008); Papastamoulis (2018b) for more detailed descriptions of the the data.

The yeast cell cycle data (Cho et al., 1998) quantifies gene expression levels over two cell cycles (17 time points). The dataset has previously been used for evaluating the effectiveness of model-based clustering techniques (Yeung et al., 2001). We used the standardized subset of the 5-phase criterion, containing $n=384$ genes measured at $p=17$ time points. The expression levels of the $n=384$ genes peak at different time points corresponding to the five phases of cell cycle, so this five class partition of the data is used as the ground-truth classification.

We applied our method using the 8 parameterizations of overfitting mixtures with Kmax = 20 components for $1\leqslant q\leqslant q_{\mbox{max}}$ factors using nChains = 4 heated chains. We set $q_{\mbox{max}}=5$ for the coffee, wave and wine datasets, while $q_{\mbox{max}}=10$ for the yeast cell cycle dataset.The number of MCMC cycles was set to mCycles = 1100, while the first burnCycles = 100 were discarded as burn-in. The 8 parameterizations are processed in parallel on parallelModels = 4 cores, while each heated chain of a given parameterization is also running in parallel. All other prior parameters were fixed at their default values.

We have also applied pgmm considering the same range of clusters and factors per dataset. For each combination of number of factors, components and parameterization, the EM algorithm was initialized using 5 random starting values as well as the K-means clustering algorithm, that is, 6 different starts in total. For the coffee dataset a larger number of different starts is required as discussed in Papastamoulis (2018b).

Table 4 summarizes the results for each of the publicly available data. We conclude that fabMix performs better than pgmm at the coffee and yeast datasets. In the wine dataset, on the other hand, pgmm performs better than fabMix, but we underline the improved performance of our method compared to the one reported by Papastamoulis (2018b) where only the UUU and UCU parameterizations were fitted. The two methods are in agreement on the wave dataset. The plot command of the fabMix package displays the estimated clusters according to the CUU model with 6 factors for the yeast dataset, as shown in Figure 6.

6 Discussion and further remarks

This study offered an efficient Bayesian methodology for model-based clustering of multivariate data using mixtures of factor analyzers. The proposed model extended the ideas of Papastamoulis (2018b) building upon the previously introduced set of parsimonious Gaussian mixture models (McNicholas and Murphy, 2008; McNicholas et al., 2010). The additional parameterizations improved the performance of the proposed method compared to Papastamoulis (2018b) where only two out of eight parameterizations were available. Furthermore, our contributed R package makes the proposed method available to a wider audience of researchers.

The computational cost of our MCMC method is larger than the EM algorithm, as shown in Figure 7. But of course, when a point estimate is required, the EM algorithm is the quickest solution. When a point estimate is not sufficient, our method offers an attractive Bayesian treatment of the problem. Clearly, the Bayesian approach does show further advantages (as in the simulated datasets according to Scenario 1, as well as in the coffee and yeast datasets), where the multimodality of the likelihood potentially causes the EM to converge to local maxima.

A direction for future research is to generalize the method in order to automatically detect the number of factors in a fully Bayesian manner. This is possible by e.g. treating the number of factors as a random variable and implementing a reversible jump mechanism in order to update it inside the MCMC sampler. Another possibility would be to incorporate strategies for searching the space of sparse factor loading matrices allowing posterior inference for factor selection (Bhattacharya and Dunson, 2011; Mavridis and Ntzoufras, 2014; Conti et al., 2014). Recent advances on infinite mixtures of infinite factor models (Murphy et al., 2019) also allow for direct inference of the number of clusters and factors and could boost the flexibility of our modelling approach.

Acknowledgements.

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 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. Rousseau and Mengersen (2011) showed that the asymptotic behaviour of the posterior distribution of the $K-K_{0}$ redundant components depends on the prior distribution of mixing proportions $(\boldsymbol{w})$. Let $d$ denote the dimension of free parameters of the distribution $f_{k}$. For the case of a Dirichlet prior distribution,

[TABLE]

if

[TABLE]

then the posterior weight of the extra components converges to zero (Theorem 1 of Rousseau and Mengersen (2011)).

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$.

In our case the dimension of free parameters in the $k$-th mixture component is equal to $d=2p+pq-\frac{q(q-1)}{2}$. Following Papastamoulis (2018b), we set $\gamma_{1}=\ldots=\gamma_{K}=\frac{\gamma}{K}$, thus the distribution of mixing proportions in Equation (A.32) 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 Rousseau and Mengersen (2011) 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 (Ferguson, 1973), that is, a Dirichlet process mixture model (Neal, 2000).

Appendix B Details of the MCMC sampler

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 Fokoué and Titterington (2003), see also Dellaportas and Papageorgiou (2006). For the case of standardized data, the prior parameters are specified in Table B.5. Standardized data is also used as input to pgmm.

Prior parallel tempering It is well known that the posterior surface of mixture models can exhibit many local modes (Celeux et al., 2000b; Marin et al., 2005). 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 Geyer (1991); Geyer and Thompson (1995); Altekar et al. (2004), 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, van Havre et al. (2015) introduced a prior parallel tempering scheme, which is also applied by Papastamoulis (2018b). 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 (van Havre et al., 2015). 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 (A.33), this is

[TABLE]

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

In our examples we used a total of $J=4$ parallel chains where the prior distribution of mixing proportions for chain $j$ in Equation (B.35) is selected as

[TABLE]

where $\delta>0$. For example, when the overfitting mixture model uses $K=20$ components and $\gamma=1$ (the default value shown in Table B.5), it follows from Equation (A.33) 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 $\delta=1$, 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 (B.34).

“Overfitting initialization” strategy We briefly describe the “overfitting initialization” procedure introduced by Papastamoulis (2018b). We used an initial period of 500 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 Rousseau and Mengersen (2011) 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$. Then, we initialize the actual model by this state. According to Papastamoulis (2018b), this specific scheme was found to outperform other initialization procedures.

Appendix C Additional simulations

In the simulation section of the manuscript, the weights of the simulated datasets have been randomly generated from a Dirichlet distribution with mean equal to $1/K$, conditional on the number of clusters ($K$). Thus, on average, the true cluster sizes are equal. In this section we examine the performance of the proposed method in the presence of unequal cluster sizes with respect to the size ($n$) of the observed data.

We replicate the simulation mechanism for scenarios 1 and 2 presented in the main text, but now we consider unequal (true) cluster sizes, as detailed in Table B.6. For each case, the sample size is increasing (as shown in the last column of Table B.6) while keeping all others parameters (that is, the true values of marginal means and factor loadings) constant. As shown in Table B.6, in scenario 1 there are $5$ clusters and $2$ factors, whereas in scenario 2 there are $2$ clusters and $3$ factors. In total 3 different examples per scenario are considered: for a given scenario, the component-specific parameters are different in each example but the weights are the same. An instance of our three examples (per scenario) using $n=200$ simulated observations is shown at Figure B.8. Observe that in all cases the “true clusters” are not easily distinguishable, especially in scenario 2 where there is a high degree of cluster overlapping.

We applied fabMix and pgmm using the same number (4) of parallel chains (for fabMix) and different starts (for pgmm) as in the simulations presented in the main paper. The results are summarized in Figures C.9 and C.10 for scenarios 1 and 2, respectively. The adjusted Rand Index is displayed in the first line of each Figure, where the horizontal axis denotes the sample size ($n$) of the synthetic data. The dotted black line corresponds to the adjusted Rand index between the ground-truth and the cluster assignments arising when applying the Maximum A Posteriori rule using the true parameter values that generated the data, that is,

[TABLE]

where $K^{*}$, $(w_{1}^{*},\ldots,w^{*}_{K^{*}})$ and $(\theta_{1}^{*},\ldots,\theta^{*}_{K^{*}})$ denote the values of number of components, mixing proportions and parameters of the multivariate normal densities of the mixture model used to generate the data. Observe that in all three examples of Scenario 1 the dotted black line is always equal to 1, but this is not the case in the more challenging Scenario 2 due to enhanced levels of cluster overlapping.

The adjusted Rand index between the ground truth clustering and the estimated cluster assignments arising from fabMix and pgmm are shown in the the first row of Figures C.9 and C.10. Clearly, the compared methods have similar performance as the sample size increases, but for smaller values of $n$ the proposed method outperforms the pgmm package.

The estimated number of clusters, shown at the second row of Figures C.9 and C.10 agree (in most cases) with the true number of clusters, but note that our method is capable of detecting the right value earlier than pgmm. Two exceptions occur at $n=200$ for example 2 of scenario 1 where fabMix (red line at 2nd row of Figure C.9) inferred 6 alive clusters instead of 5, as well as at $n=800$ for example 1 of scenario 2 where fabMix (red line at 2nd row of Figure C.10) inferred 3 alive clusters instead of 2.

Finally, the last row in Figures C.9 and C.10 displays the inferred number of factors for Scenarios 1 and 2, respectively. In every single case, the estimate arising from fabMix is at least as close as the estimate arising from pgmm to the corresponding true value. Note however that in example 1 of scenario 2 both methods detect a smaller number of factors (2 instead of 3 factors). In all other cases we observe that as the sample size increases both methods infer the “true” number of factors.