Infinite Mixtures of Infinite Factor Analysers

TL;DR

The paper introduces IMIFA, a flexible Bayesian model for clustering high-dimensional data that automatically infers the number of clusters and factors, eliminating the need for pre-specification and extensive model selection.

Contribution

IMIFA is the first model to combine infinite mixtures with infinite factor analysis, allowing automatic inference of cluster and factor numbers with improved performance.

Findings

IMIFA automatically infers the number of clusters and factors.

IMIFA outperforms traditional models in clustering accuracy.

IMIFA reduces computational burden and quantifies uncertainty.

Abstract

Factor-analytic Gaussian mixture models are often employed as a model-based approach to clustering high-dimensional data. Typically, the numbers of clusters and latent factors must be specified in advance of model fitting, and remain fixed. The pair which optimises some model selection criterion is then chosen. For computational reasons, models in which the number of latent factors differ across clusters are rarely considered. Here the infinite mixture of infinite factor analysers (IMIFA) model is introduced. IMIFA employs a Pitman-Yor process prior to facilitate automatic inference of the number of clusters using the stick-breaking construction and a slice sampler. Furthermore, IMIFA employs multiplicative gamma process shrinkage priors to allow cluster-specific numbers of factors, automatically inferred via an adaptive Gibbs sampler. IMIFA is presented as the flagship of a family of factor-analytic mixture models, providing flexible approaches to clustering high-dimensional data. Applications to a benchmark data set, metabolomic spectral data, and a manifold learning handwritten digit example illustrate the IMIFA model and its advantageous features. These include obviating the need for model selection criteria, reducing the computational burden associated with the search of the model space, improving clustering performance by allowing cluster-specific numbers of factors, and quantifying uncertainty in the numbers of clusters and cluster-specific factors.

Tables

Table 1. Table 1: Hyperparameter specifications for the IMIFA model. Note that the specification of the beta distribution in the prior for d 𝑑 d amounts to a standard uniform.

Parameter(s)	Hyperparameter(s)	Value(s)
${𝝁}_g$	$\phi$	$0.01$
${𝚿}_g$	$({\alpha }_0,{\beta }_0)$	$(2.5,3)$
${\varphi }_{jkg}$	$({\nu }_1,{\nu }_2)$	$(3,2)$
${\delta }_{1g}$	$({\alpha }_1,{\beta }_1)$	$(2.1,1)$
${\delta }_{kg}$	$({\alpha }_2,{\beta }_2)$	$(3.1,1)$
${\sigma }_g$	$({\varrho }_1,{\varrho }_2)$	$(3,2)$
$\alpha$	$(a,b)$	$(2,4)$
$d$	$(a^{\prime },b^{\prime },\kappa )$	$(1,1,0.5)$

Table 2. Table 2: Results of fitting a range of models, including the full IMIFA family, to the Italian olive oil data, detailing the number of candidate models explored, the run-time relative to the IMIFA run (approx. 782 782 782 seconds), the posterior mean or fixed value of α 𝛼 \alpha , the posterior mean of d 𝑑 d , modal estimates of G 𝐺 G and 𝐐 𝐐 \mathbf{Q} , and the ARI and misclassification rate as evaluated against the known area labels, under the optimal or modal model as appropriate.

Model	# Models	Relative Time	$\alpha$	$d$	$G$	$𝐐$	ARI	Error (%)
IMIFA	$1$	$1.00$	$0.48$	$0.01$	$4$	$6$, $3$, $6$, $2$	$0.94$	$8.39$
IMFA	$7$	$4.14$	$0.62$	$0.01$	$5$	$6$, $6$, $6$, $6$, $6$	$0.91$	$14.86$
OMIFA	$1$	$1.19$	$0.02$	—	$4$	$6$, $3$, $6$, $4$	$0.93$	$9.97$
OMFA	$7$	$5.11$	$0.03$	—	$5$	$6$, $6$, $6$, $6$, $6$	$0.85$	$15.56$
MIFA	$9$	$3.41$	$1$	—	$5$	$6$, $3$, $6$, $6$, $4$	$0.92$	$10.31$
MFA	$63$	$13.86$	$1$	—	$2$	$5$, $5$	$0.82$	$17.13$
IFA	$1$	$0.11$	—	—	$1$	$6$	—	—
FA	$7$	$0.37$	—	—	$1$	$6$	—	—
mclust\tnotextnA	$115$	$0.01$	—	—	$6$	—	$0.56$	$38.64$
MFMA\tnotextnA	$1,350$	$4.68$	—	—	$4$	$5$, $5$, $5$, $5$	$0.68$	$20.28$
pgmm\tnotextnA  ,\tnotextnB	$588$	$4.46$	—	—	$5$	$6$, $6$, $6$, $6$, $6$	$0.53$	$35.84$

Table 3. Table 3: Confusion matrices of the MAP IMIFA clustering of the Italian olive oils against ( 3a ) the known 3 3 3 area labels and ( 3b ) the new labelling in which northern Italy is split into its constituent sub-regions.

	1	2	3	4
Southern Italy	$323$	$0$	$0$	$0$
Sardinia	$0$	$98$	$0$	$0$
Northern Italy	$0$	$0$	$103$	$48$

Table 4. a 3 area cross tabulation

	1	2	3	4
Southern Italy	$323$	$0$	$0$	$0$
Sardinia	$0$	$98$	$0$	$0$
Northern Italy	$0$	$0$	$103$	$48$

Table 5. b 4 area cross tabulation

	1	2	3	4
Southern Italy	$323$	$0$	$0$	$0$
Sardinia	$0$	$98$	$0$	$0$
East Liguria & Umbria	$0$	$0$	$100$	$0$
West Liguria	$0$	$0$	$3$	$48$

Table 6. Table 4: Cross tabulation of the IMIFA model’s MAP clustering (rows) against true digit labels (columns) for the USPS data. Cells that are 0 0 are blank for clarity. Posterior means π ^ g subscript ^ 𝜋 𝑔 \smash{\widehat{\pi}_{g}} and modal estimates q ^ g subscript ^ 𝑞 𝑔 \smash{\widehat{q}_{g}} , with associated 95 % percent 95 95\% credible intervals, are also given.

	0	1	2	3	4	5	6	7	8	9	${\hat{\pi }}_g$	${\hat{q}}_g$
1	$359$										$0.05$	$4[2,8]$
2	$58$		$12$			$3$	$2$				$0.01$	$3[2,7]$
3	$108$										$0.01$	$2[1,4]$
4	$9$										$0.00$	$16[3,16]$
5	$95$										$0.01$	$4[1,8]$
6	$308$					$3$					$0.04$	$7[4,10]$
7		$844$			$2$						$0.12$	$2[0,4]$
8		$133$							$1$		$0.02$	$1[0,4]$
9		$2$	$392$	$10$		$1$					$0.05$	$7[5,12]$
10	$59$		$121$	$93$	$19$	$91$	$13$	$2$	$25$	$4$	$0.06$	$12[9,16]$
11				$136$		$64$					$0.03$	$5[2,9]$
12					$38$	$1$		$1$			$0.01$	$2[0,8]$
13	$25$		$3$	$7$	$98$	$51$	$2$	$36$	$59$	$28$	$0.04$	$8[5,12]$
14	$48$		$73$	$61$	$62$	$135$	$32$	$1$	$16$	$6$	$0.06$	$8[6,12]$
15	$1$						$83$				$0.01$	$3[1,7]$
16	$1$						$74$				$0.01$	$2[1,5]$
17		$2$			$4$	$19$	$381$		$2$		$0.06$	$2[1,6]$
18								$207$			$0.03$	$4[1,8]$
19	$123$	$8$	$129$	$348$	$247$	$184$	$77$	26	$420$	$84$	$0.23$	$6[3,9]$
20		$16$	$1$	$3$	$120$	$1$		$338$	$19$	$451$	$0.13$	$2[1,6]$
21					$62$	$3$		$34$		$71$	$0.02$	$3[1,6]$

Table 7. Table B.1 : Aggregated simulation study results for the IMIFA model under different dimensionality scenarios and settings of the concentration and discount parameters α 𝛼 \alpha and d 𝑑 d (posterior mean estimates thereof in parentheses where appropriate). The modal estimates of G 𝐺 G and associated estimates of q g ​ ∀ g subscript 𝑞 𝑔 for-all 𝑔 q_{g}\>\>\forall\>\>g are reported (with 95 % percent 95 95\% credible intervals in brackets). Clustering performance is assessed through the average percentage error rate against the known cluster labels.

Dimension	$\alpha$	$d$	$G$	$q_1$	$q_2$	$q_3$	Error (%)
$N=25$	$0.5$	$0$	$3[3,3]$	$5[3,9]$	$5[3,9]$	$5[3,9]$	$0$
	$1$	$0$	$3[3,3]$	$5[3,9]$	$5[3,9]$	$5[3,9]$	$0$
	$5$	$0$	$3[3,4]$	$5[3,9]$	$5[3,9]$	$5[3,9]$	$6.4$
	$\left(0.57\right)$	$0$	$3[3,3]$	$5[3,9]$	$5[3,9]$	$5[3,9]$	$0$
	$\left(0.51\right)$	$\left(0.05\right)$	$3[3,3]$	$5[3,9]$	$5[3,9]$	$5[3,9]$	$0$
$N=50$ $\left(N=p\right)$	$0.5$	$0$	$3[3,3]$	$5[4,7]$	$5[4,7]$	$5[4,7]$	$0$
	$1$	$0$	$3[3,3]$	$5[4,7]$	$5[4,7]$	$5[4,7]$	$0$
	$5$	$0$	$3[3,3]$	$5[4,7]$	$5[4,7]$	$5[4,7]$	$0$
	$\left(0.52\right)$	$0$	$3[3,3]$	$5[4,7]$	$5[4,7]$	$5[4,7]$	$0$
	$\left(0.48\right)$	$\left(0.03\right)$	$3[3,3]$	$5[4,7]$	$5[4,7]$	$5[4,7]$	$0$
$N=300$ $\left(N>p\right)$	$0.5$	$0$	$3[3,3]$	$5[4,6]$	$5[4,6]$	$5[4,6]$	$0$
	$1$	$0$	$3[3,3]$	$5[4,6]$	$5[4,6]$	$5[4,6]$	$0$
	$5$	$0$	$3[3,3]$	$5[4,6]$	$5[4,6]$	$5[4,6]$	$0$
	$\left(0.42\right)$	$0$	$3[3,3]$	$5[4,6]$	$5[4,6]$	$5[4,6]$	$0$
	$\left(0.39\right)$	$\left(0.02\right)$	$3[3,3]$	$5[4,6]$	$5[4,6]$	$5[4,6]$	$0$

Table 8. Table C.1 : Clustering performance of the IMIFA model on expanded noisy versions of the Italian olive oil data. The run-time relative to running IMIFA on the original data, posterior mean of the PYP parameters α 𝛼 \alpha and d 𝑑 d , modal estimates of G 𝐺 G and 𝐐 𝐐 \mathbf{Q} , ARI, and percentage error rate are all given.

Scenario	Relative Time	$\alpha$	$d$	$G$	$𝐐$	ARI	Error (%)
$N=572,p=18$	$1.86$	$0.48$	$0.01$	$4$	$3$, $4$, $4$, $3$	$0.85$	$12.59$
$N=572,p=58$	$3.14$	$0.47$	$0.01$	$4$	$1$, $2$, $2$, $2$	$0.74$	$14.69$
$N=572,p=108$	$5.64$	$0.46$	$0.02$	$4$	$0$, $1$, $0$, $2$	$0.73$	$17.66$
$N=582,p=8$	$1.10$	$0.57$	$0.01$	$5$	$6$, $2$, $2$, $2$, $0$	$0.94$	$6.87$
$N=622,p=8$	$1.09$	$0.56$	$0.01$	$5$	$4$, $1$, $1$, $2$, $0$	$0.95$	$6.59$
$N=672,p=8$	$1.07$	$0.53$	$0.01$	$5$	$4$, $1$, $2$, $2$, $0$	$1.00$	$0.45$

Equations

$$f\negthinspace\left(\mathbf{x}_{i}\,|\,\bm{\theta}\right)=\sum_{g=1}^{G}\pi_{g}\textrm{N}_{p}\big{(}\mathbf{x}_{i};\bm{\mu}_{g},\bm{\Sigma}_{g}=\bm{\Lambda}{\mathstrut}_{g}\bm{\Lambda}{\mathstrut}_{g}^{\top}+\bm{\Psi}_{g}\big{)},\vspace{-0.25em}$$

$${\mathbf{\Lambda }}_{jg}=\left({\lambda }_{j1g},\dots ,{\lambda }_{jqg}\right)\sim {\text{N}}_q\left(0,{\mathcal{I}}_q\right).\text{vspace}-0.25em$$

$${\mathbf{\mu }}_g\sim {\text{N}}_p\left(\tilde{\mathbf{\mu }},{\phi }^{-1}{\mathcal{I}}_p\right),\text{vspace}-0.25em$$

$$\displaystyle\begin{split}\lambda_{jkg}\,|\,\phi_{jkg},\tau_{kg},\sigma_{g}&\sim\textrm{N}_{1}\negthinspace\left(0,\phi_{jkg}^{-1}\tau_{kg}^{-1}\sigma_{g}^{-1}\right),&\quad\phi_{jkg}&\sim\textrm{Ga}\negthinspace\left(\nu_{1},\nu_{2}\right),\\[-0.96873pt] \tau_{kg}&=\prod_{h=1}^{k}\delta_{hg}\,,&\quad\sigma_{g}&\sim\textrm{Ga}\negthinspace\left(\varrho_{1},\varrho_{2}\right),\\[-0.96873pt] \delta_{1g}&\sim\textrm{Ga}\negthinspace\left(\alpha_{1},\beta_{1}\right),&\quad\delta_{hg}&\sim\textrm{Ga}\negthinspace\left(\alpha_{2},\beta_{2}\right)\>\>\forall\>\>h\geq 2,\end{split}$$

$$f\negthinspace\left(\mathbf{x}_{i}\,|\,\bm{\theta}\right)=\sum_{g=1}^{\infty}\pi_{g}\textrm{N}_{p}\big{(}\mathbf{x}_{i};\bm{\mu}_{g},\bm{\Lambda}{\mathstrut}_{g}\bm{\Lambda}{\mathstrut}_{g}^{\top}+\bm{\Psi}_{g}\big{)}.\vspace{-0.33em}$$

$$\displaystyle\begin{split}\upsilon_{g}&\sim\textrm{Beta}\negthinspace\left(1-d,\alpha+gd\right),&\quad\bm{\theta}_{g}&\sim H_{0}\,,\\ \pi_{g}&=\upsilon_{g}\prod_{l=1}^{g-1}\left(1-\upsilon_{l}\right),&\quad H&=\sum_{g=1}^{\infty}\pi_{g}\delta_{\bm{\theta}_{g}}\sim\textrm{PYP}\negthinspace\left(\alpha,d,H_{0}\right),\end{split}$$

$$f\left(\mathbf{X},\mathbf{\eta },\mathbf{Z},\mathbf{u},\mathbf{{\rm Y}},\mathbf{\theta }\right)$$

$$\displaystyle=\Bigg{\{}\negthinspace\prod_{i=1}^{N}\prod_{g\,\in\,\bm{\mathcal{A}}_{\bm{\xi}}\negthinspace\left(u_{i}\right)}\negthinspace\textrm{N}_{p}\big{(}\mathbf{x}_{i};\bm{\mu}_{g}+\bm{\Lambda}_{g}\bm{\eta}{\mathstrut}_{i},\bm{\Psi}_{g}\big{)}^{z_{ig}}\negthinspace\Bigg{\}}\Bigg{\{}\negthinspace\prod_{i=1}^{N}\textrm{N}_{q}\big{(}\bm{\eta}{\mathstrut}_{i};\mathbf{0},\bm{\mathcal{I}}_{q}\big{)}\negthinspace\Bigg{\}}$$

$$\displaystyle\phantom{=\,\,}\left\{\prod_{i=1}^{N}\prod_{g=1}^{\infty}\negthinspace\left(\dfrac{\pi_{g}}{\xi_{g}}\mathds{1}{{\big{(}u_{i}<\xi_{g}\big{)}}}\negthinspace\right)^{z_{ig}}\negthinspace\right\}\negmedspace\left\{\prod_{g=1}^{\infty}\dfrac{\left(1-\upsilon_{g}\right)^{\alpha+gd-1}}{\upsilon_{g}^{d}\leavevmode\nobreak\ \mathrm{B}\negthinspace\left(1-d,\alpha+gd\right)}\right\}\negthinspace f\negthinspace\left(\bm{\theta}\right)\negthinspace,$$

$$\Big{\lvert}{\big{\lVert}\bm{\mathcal{H}}\big{\rVert}}_{\mathcal{F}}-{\big{\lVert}\bm{\mathcal{H}}^{\left(r\right)}\big{\rVert}}_{\mathcal{F}}\Big{\rvert}\leq{\big{\lVert}\bm{\mathcal{H}}-\bm{\mathcal{H}}^{\left(r\right)}\big{\rVert}}_{\mathcal{F}}\leq{\big{\lVert}\bm{\mathcal{H}}\big{\rVert}}_{\mathcal{F}}+{\big{\lVert}\bm{\mathcal{H}}^{\left(r\right)}\big{\rVert}}_{\mathcal{F}}\,.$$

$${\mathbf{\mu }}_g|\dots$$

$$\mathbf{\eta }{}_i|z_{ig}=1,\dots$$

$${\psi }_{jg}|\dots$$

$${\mathbf{\Lambda }}_{jg}|\dots$$

$${\varphi }_{jkg}|\dots$$

$${\delta }_{1g}|\dots$$

$${\delta }_{kg}|\dots$$

$${\sigma }_g|\dots$$

$${\upsilon }_g|\dots$$

$$u_i|z_{ig}=1,\dots$$

$${\mathbf{\Omega }}_{\mathbf{\mu }}$$

$${\mathbf{\Omega }}_{\mathbf{\eta }}$$

$${\mathbf{\Omega }}_{\mathbf{\lambda }}$$

$${\tau }_{hg}^{\left(k\right)}$$

$${\pi }_g$$

$${\mathcal{S}}_{jg}$$

$${\psi }_g|\dots$$

$${\psi }_j|\dots$$

$$\psi |\dots$$

$$p_{ig}=\Pr\negthinspace\left(z_{ig}=1\,|\,\mathbf{x}_{i},\ldots\right)=\frac{\pi_{g}\textrm{N}_{p}\big{(}\mathbf{x}_{i};\bm{\mu}_{g},\bm{\Lambda}{\mathstrut}_{g}\bm{\Lambda}{\mathstrut}_{g}^{\top}+\bm{\Psi}_{g}\big{)}}{\sum_{g=1}^{\widetilde{G}}\pi_{g}\textrm{N}_{p}\big{(}\mathbf{x}_{i};\bm{\mu}_{g},\bm{\Lambda}{\mathstrut}_{g}\bm{\Lambda}{\mathstrut}_{g}^{\top}+\bm{\Psi}_{g}\big{)}},$$

$$p_{ig}\propto\textrm{N}_{p}\big{(}\mathbf{x}_{i};\bm{\mu}_{g},\bm{\Lambda}{\mathstrut}_{g}\bm{\Lambda}{\mathstrut}_{g}^{\top}+\bm{\Psi}_{g}\big{)}\frac{\pi_{g}}{\xi_{g}}\mathds{1}{{\big{(}u_{i}<\xi_{g}\big{)}}}.$$

$$\underset{g\,\in\,\left\{1,\ldots,\widetilde{G}\right\}}{\operatorname{arg}\,\operatorname{max}}\;\big{(}\negthinspace\ln\negthinspace\left(p_{ig}\right)+s_{ig}\big{)}.$$

$$\alpha |d,\dots$$

$$d|\alpha ,\dots$$

$$\alpha\,|\,G_{0},\chi,\ldots\sim\omega_{\chi}\textrm{Ga}\big{(}a+G_{0},b-\ln\negthinspace\left(\chi\right)\big{)}+\left(1-\omega_{\chi}\right)\textrm{Ga}\big{(}a+G_{0}-1,b-\ln\negthinspace\left(\chi\right)\big{)},\vspace{-0.25ex}$$

$$\left(\chi |\alpha ,G_0\right)$$

$$\frac{{\omega }_{\chi }}{1-{\omega }_{\chi }}$$

$$\min\big{(}1,\left(\pi_{h}/\pi_{g}\right)^{n_{g}-n_{h}}\negmedspace\big{)}.$$

$$\min\big{(}1,\left(1-\upsilon_{l+1}\right)^{n_{l}}\negmedspace/\left(1-\upsilon_{l}\right)^{n_{l+1}}\negmedspace\big{)};$$

$$\alpha\,|\,\mathbf{Z},\widetilde{G},\ldots\propto\frac{\Gamma\big{(}\alpha\widetilde{G}\big{)}}{\Gamma\big{(}N+\alpha\widetilde{G}\big{)}}\Bigg{\{}\prod_{g\colon n_{g}>0}\frac{\Gamma\left(n_{g}+\alpha\right)}{\Gamma\left(\alpha\right)}\Bigg{\}}p\left(\alpha\right).$$

$$f\big{(}\mathbf{x}_{i}\,|\,\bm{\eta}{\mathstrut}_{i},\bm{\theta}\big{)}=\sum_{g=1}^{G}\pi_{g}\textrm{N}_{p}\big{(}\mathbf{x}_{i};\bm{\mu}_{g}+\bm{\Lambda}_{g}\bm{\eta{\mathstrut}}_{i},\bm{\Psi}_{g}\big{)}.$$

Full text

Here the infinite mixture of infinite factor analysers (IMIFA) model is introduced. IMIFA employs a Pitman-Yor process prior to facilitate automatic inference of the number of clusters using the stick-breaking construction and a slice sampler. Automatic inference of the cluster-specific numbers of factors is achieved using multiplicative gamma process shrinkage priors and an adaptive Gibbs sampler. IMIFA is presented as the flagship of a family of factor-analytic mixtures.

Applications to benchmark data, metabolomic spectral data, and a handwritten digit example illustrate the IMIFA model’s advantageous features. These include obviating the need for model selection criteria, reducing the computational burden associated with the search of the model space, improving clustering performance by allowing cluster-specific numbers of factors, and uncertainty quantification.

Keywords: adaptive Markov chain Monte Carlo, factor analysis, model-basedclustering, Pitman-Yor process, multiplicative gamma process.

1 Introduction

In cases where the number of variables $p$ is comparable to or greater than the number of observations $N$, many clustering techniques tend to perform poorly or be intractable. Factor analysis (FA; Knott & Bartholomew, 1999) is a well-known approach to parsimoniously modelling data. Bai & Li (2012) outline some computational difficulties which arise when $N\ll p$. Model-based clustering methods which rely on latent factor models have long been successfully utilised to cluster high-dimensional data. Ghahramani & Hinton (1996) propose a mixture of factor analysers model (MFA) with cluster-specific parsimonious covariance matrices and estimate it via an expectation-maximisation algorithm; McLachlan & Peel (2000) provide a succinct overview. Estimation of MFA models has also been considered in a Bayesian framework (Diebolt & Robert, 1994; Richardson & Green, 1997). McNicholas & Murphy (2008, 2010) develop a suite of similar parsimonious Gaussian mixture models. Other related developments in this area include Baek et al. (2010) and Viroli (2010), among others.

Clustering using a MFA model typically requires specifying the number of clusters and factors in advance of model fitting. Generally, a range of MFA models with different numbers of clusters and factors are fitted and then compared through the use of information criteria, such as the Bayesian Information Criterion (BIC; Kass & Raftery, 1995) or the Deviance Information Criterion (Spiegelhalter et al., 2002, 2014). Within a Bayesian framework Fokoué & Titterington (2003) use a stochastic model selection approach but do not simultaneously choose the optimal number of clusters and factors. Conducting an exhaustive search of the model space is computationally expensive; the cost is typically reduced by only considering models in which the number of factors is common across clusters. Regardless, even searching the reduced model space can be computationally onerous. The problem of identifying the optimal model is exacerbated by the fraught task of choosing among the range of model selection tools available, which often suggest different optimal models. Moreover, enforcing a common number of factors across clusters may lead to poor clustering performance due to a lack of flexibility.

The infinite mixture of infinite factor analysers (IMIFA) model is introduced here. It theoretically allows infinitely many components and infinitely many factors within each component. The need to select a model selection criterion is obviated and quantification of the uncertainty in the optimal numbers of non-empty clusters and cluster-specific factors is facilitated. IMIFA relies on an infinite mixture model through the use of a nonparametric Pitman-Yor process (PYP) prior (Perman et al., 1992; Pitman & Yor, 1997), of which the well-known Dirichlet process (DP; Ferguson, 1973) is a special case. The infinite mixture model framework allows the number of clusters present to be automatically inferred; here the stick-breaking construction (Pitman, 1996) and an independent slice-efficient sampler (Kalli et al., 2011) are employed to facilitate this.

By allowing infinitely many factors within each cluster, IMIFA addresses the difficulty in choosing the optimal number of factors. This facilitates fitting factor-analytic models which are more flexible, in the sense that the number of factors may be cluster-specific, thereby potentially improving clustering performance. This is achieved by assuming multiplicative gamma process (MGP) shrinkage priors (Bhattacharya & Dunson, 2011; Durante, 2017) on the cluster-specific factor loading matrices, thus generalising the MGP prior to the mixture setting. Such a prior allows the degree of shrinkage of the factor loadings towards zero to increase as the factor number tends towards infinity. The number of factors with non-negligible loadings can be considered as the ‘active’ number of factors within each cluster. Following Bhattacharya & Dunson (2011), a computationally efficient adaptive Gibbs sampling algorithm is employed for estimation. Thus, the choice of the numbers of active factors in different clusters is automated.

The IMIFA model with its PYP-MGP priors thus offers a single-pass and therefore computationally efficient approach to clustering high-dimensional data. It can be viewed as the most flexible model at the head of a family of Bayesian factor-analytic mixture models. Section 2 develops the hierarchy of the IMIFA model family, beginning with the MFA model and concluding with the flagship IMIFA model. Between these extremes the novel finite mixture of infinite factor analysers model (MIFA) is introduced. Overfitted factor-analytic mixtures (Papastamoulis, 2018) also belong to the IMIFA family; the overfitted mixture of infinite factor analysers (OMIFA) model is also introduced here.

Section 3 considers implementation of the IMIFA family of models. A benchmarking experiment is conducted on the well-known Italian olive oil data set. A real data application follows through the cluster analysis of spectral metabolomic data from an epilepsy study. Finally an illustrative application is provided through clustering United States Postal Service handwritten digit data, a setting for which fitting sub-models of the IMIFA family is practically infeasible. Comparisons against other clustering methods are provided throughout. Simulation studies demonstrating the performance of IMIFA under different scenarios are deferred to Appendix B. Section 4 concludes the article with a discussion of IMIFA and thoughts on future research directions.

A software implementation of IMIFA and its family of sub-models is provided by the associated R package IMIFA (Murphy et al., 2021), with which all results were generated, which is freely available from https://www.r-project.org (R Core Team, 2021).

2 The IMIFA Model Family

The hierarchy of the IMIFA family of models is delineated herein, including a review of extant methodologies, the introduction of novel sub-models, and concluding with the flagship IMIFA model. Prior specifications, Markov chain Monte Carlo (MCMC) inferential procedures, approaches to posterior predictive model checking, and model-specific implementation issues that arise in practice are addressed.

2.1 Mixtures of Factor Analysers

Mixtures of factor analysers are Gaussian latent variable models often used for clustering high-dimensional data. For each of $G$ clusters in these finite mixtures, the cluster-specific FA model in cluster $g$ is given by $\smash{\mathbf{x}_{i}-\bm{\mu}_{g}=\bm{\Lambda}_{g}\bm{\eta}{\mathstrut}_{i}+\bm{\varepsilon}_{ig}}$. The observed feature vector $\smash{\mathbf{x}_{i}=\left(x_{i1},\ldots,x_{ip}\right)^{\top}}$ with mean $\smash{\bm{\mu}_{g}}$ and covariance matrix $\smash{\bm{\Sigma}_{g}}$ is assumed to linearly depend on a $q$-vector $\left(q\ll p\right)$ of latent common factor scores $\smash{\bm{\eta}{\mathstrut}_{i}}$ and additional sources of variation called specific factors $\smash{\bm{\varepsilon}_{ig}}$. It is assumed that $\smash{\bm{\eta}{\mathstrut}_{i}}$ has a $q$-variate Gaussian distribution $\smash{\textrm{N}_{q}\negthinspace\left(\mathbf{0},\bm{\mathcal{I}}_{q}\right)},$ where $\smash{\bm{\mathcal{I}}_{q}}$ denotes the $q\times q$ identity matrix, and that $\smash{\bm{\varepsilon}_{ig}\sim\textrm{N}_{p}\negthinspace\left(\mathbf{0},\bm{\Psi}_{g}\right)}$, where $\smash{\bm{\Psi}_{g}}$ is a diagonal matrix with non-zero elements $\smash{\psi_{1g},\ldots,\psi_{pg}}$ known as uniquenesses. Here, $\smash{\bm{\Lambda}_{g}}$ denotes the $p\times q$ factor loadings matrix of cluster $g$ and notably $q=0$ is permitted.

To facilitate estimation, a latent cluster indicator vector $\smash{\mathbf{z}_{i}=(z_{i1},\ldots,z_{iG})^{\top}}$ is introduced such that $\smash{z_{ig}=1}$ if observation $i$ belongs to cluster $g$ and $\smash{z_{ig}=0}$ otherwise. Hence, $\smash{\mathbf{z}_{i}}$ has a $\smash{\textrm{Mult}\negthinspace\left(1,\bm{\pi}\right)}$ distribution where $\smash{\bm{\pi}=\left(\pi_{1},\ldots,\pi_{G}\right)^{\top}}$ are the cluster mixing proportions which sum to $1$. A symmetric uniform Dirichlet prior $\bm{\pi}\sim\textrm{Dir}_{G}\negthinspace\left(\bm{\alpha}=\left(\alpha,\ldots,\alpha\right)=\mathbf{1}\right)$ is assumed. Upon marginalising out $\smash{\mathbf{z}_{i}}$ and $\smash{\bm{\eta}{\mathstrut}_{i}}$, MFA yields a parsimonious finite sum covariance structure for the observed data

[TABLE]

where $\textrm{N}_{p}\negthinspace\left(\mathbf{x}_{i};\cdot,\cdot\right)$ denotes the density of a $p$-variate Gaussian distribution evaluated at $\mathbf{x}_{i}$ and $\smash{\bm{\theta}_{g}=\left\{\bm{\mu}_{g},\bm{\Lambda}_{g},\bm{\Psi}_{g}\right\}}$ are the cluster-specific FA parameters for which inference is straightforward under a Gibbs sampling scheme. Imposing constraints on $\smash{\mathbf{\Psi}_{g}}$ (McNicholas & Murphy, 2008, 2010) and/or fixing $\smash{\pi_{g}=\nicefrac{{1}}{{G}}\>\>\forall\>\>g}$ may be useful in some settings.

2.1.1 Prior Specification and Practical Issues

The conditionally conjugate nature of the various prior distributions detailed below facilitates MCMC sampling via straightforward Gibbs updates. A multivariate Gaussian prior is assumed for the factor loadings of the variable $j$ across the $q$ factors of cluster $g$:

[TABLE]

Similarly, a diffuse multivariate Gaussian prior is assumed for the component means,

[TABLE]

where $\widetilde{\bm{\mu}}$ is the overall sample mean and the scalar $\varphi$ controls the level of diffusion.

An inverse gamma prior $\psi_{jg}\sim\textrm{IG}\negthinspace\left(\alpha_{0},\beta_{j}\right)$ is assumed for the uniquenesses of variable $j$ in cluster $g$. Guided by Frühwirth-Schnatter & Lopes (2010), hyperparameters are chosen to ensure $\psi_{jg}$ is bounded away from [math], thereby avoiding Heywood problems. With a sufficiently large shape $\alpha_{0}$, variable-specific scales are derived from the sample precision matrix $\mathbf{S}^{\star}=\mathbf{S}^{-1}$ via $\beta_{j}=\left(\alpha_{0}-1\right)\negthinspace/S^{\star}_{jj}$. However, when $\nicefrac{{N}}{{p}}$ is close to or less than $1$, or when $\mathbf{S}^{-1}$ is otherwise unavailable, $\mathbf{S}^{\star}$ is replaced by a ridge-type estimator \widehat{\mathbf{S}^{-1}}=\big{(}\beta_{0}+\nicefrac{{N}}{{2}}\big{)}\big{(}\beta_{0}\bm{\mathcal{I}}_{p}+0.5\sum_{i=1}^{N}\mathbf{x}{\mathstrut}_{i}\mathbf{x}{\mathstrut}_{i}^{\top}\big{)}^{-1}, which combines the the inverse Wishart prior $\mathbf{S}^{-1}\sim\textrm{W}_{p}\negthinspace\left(\beta_{0},\beta_{0}\bm{\mathcal{I}}_{p}\right)$ with the sample information, where $\beta_{0}$ is a hyperparameter (Frühwirth-Schnatter & Lopes, 2018). For unstandardised data, this estimator is constructed for the inverse correlation matrix and then appropriately scaled using the diagonal entries of $\mathbf{S}$ (Wang et al., 2015). When the variances are roughly balanced, constraining $\bm{\Psi}_{g}$ to $\psi_{g}\bm{\mathcal{I}}_{p}\,$, and/or using $\beta_{j}=\beta=\left(\alpha_{0}-1\right)\negthinspace/{\max\negthinspace\left(\operatorname{diag}\negthinspace\left(\mathbf{S}^{\star}\right)\right)}$, provides additional parsimony. Notably, the isotropic constraint provides the link between factor analysis and probabilistic principal component analysis (Tipping & Bishop, 1999).

The rotational invariance property which makes FA models non-identifiable is well known: most covariance matrices $\bm{\Sigma}$ cannot be uniquely factored as $\smash{\bm{\Lambda}\bm{\Lambda}^{\top}+\bm{\Psi}}$ when $q>1$. Though identifiability of $\bm{\Lambda}$ is not strictly necessary for the purposes of clustering or inferring $\bm{\Sigma}$, addressing the identifiability problem offline using the parameter expanded approach of Ghosh & Dunson (2008) in tandem with Procrustean methods, as in McParland et al. (2014) and Aßmann et al. (2016), yields interpretable posterior summaries. Another practical issue is the label switching phenomenon (Frühwirth-Schnatter, 2010) which is addressed offline using the cost-minimising permutation given by the square assignment algorithm (Carpaneto & Toth, 1980). Finally, optimal FA and MFA models are chosen using the BIC-MCMC criterion (Frühwirth-Schnatter, 2011) where necessary in what follows.

2.2 Mixtures of Infinite Factor Analysers

To overcome the requirement to specify $q$, infinite factor analysis (IFA) models are employed (Bhattacharya & Dunson, 2011). The IFA model is a factor analysis model which assumes a multiplicative gamma process (MGP) shrinkage prior on the loadings matrix. This prior allows the degree of shrinkage towards zero to increase as the column index $k\rightarrow\infty$, mitigating against the factor splitting phenomenon. Here the IFA model is generalised to the mixture setting, leading to the novel mixture of infinite factor analysers (MIFA) model. Under MIFA, the MGP prior is placed on each parameter expanded $\smash{\bm{\Lambda}_{g}}$ matrix with no restrictions on its entries, thereby making the induced prior on $\smash{\bm{\Sigma}_{g}}$ invariant to the ordering of the variables. The MGP prior is conditionally conjugate, facilitating block Gibbs updates of the loadings and hence rapid mixing. Thus, the MGP prior in mixture settings is given by

[TABLE]

where $\tau_{kg}$ is a column shrinkage parameter for the $k$-th column in the $g$-th cluster’s loadings matrix $\bm{\Lambda}_{g}\>\>\forall\>\>k=1,\ldots,\infty$, and $\textrm{Ga}\negthinspace\left(\alpha,\beta\right)$ denotes the gamma distribution with mean $\nicefrac{{\alpha}}{{\beta}}$. The role of the local shrinkage parameters $\phi_{1kg},\ldots,\phi_{pkg}$ for the $p$ elements in column $k$ of $\bm{\Lambda}_{g}$ is to favour sparsity while also preserving the signal of non-zero loadings. Lastly, the cluster shrinkage parameter $\sigma_{g}$ reflects the belief that the degree of shrinkage is cluster-specific. A schematic illustration of the MGP prior is given in Figure 1; note that loadings can shrink arbitrarily close, but not exactly, to zero.

Bhattacharya & Dunson (2011) fix $\smash{\beta_{1}=\beta_{2}=1}$ and recommend that $\smash{\alpha_{2}>\beta_{2}}$. However, Durante (2017) elaborates on the cumulative shrinkage properties and roles played by hyperparameters, showing in particular that $\smash{\alpha_{2}>\beta_{2}+1}$ is necessary in order to have column-specific variances $\smash{\tau_{kg}^{-1}}$ that decrease in expectation with growing $k$. It is also recommended that $\smash{\alpha_{2}}$ be moderately large relative to $\smash{\alpha_{1}}$ (to ensure that the cumulative shrinkage property for which the prior was developed holds) and to avoid excessively high values for $\smash{\alpha_{1}}$ (to avoid over-shrinking to increasingly low-dimensional factorisations). While Bhattacharya & Dunson (2011) assume $\textrm{Ga}\negthinspace\left(\nu,\nu\right)$ priors for the local shrinkage parameters, here more general settings are used to allow control over prior non-informativity. In the spirit of Durante (2017), the expectation $\smash{\nu_{2}/\negthinspace\left(\nu_{1}-1\right)}$ of the induced inverse gamma prior on $\smash{\phi_{jkg}^{-1}}$ is suggested to be $\leq 1$ to induce sparsity on average. Furthermore, following the guidelines of Durante (2017), it is generally advisable that all MGP hyperparameters are chosen such that the first two moments of the associated hyperprior are defined, as this leads to superior performance in terms of the expected deviation between the true and estimated covariance matrices. In the mixture setting, $\smash{\alpha_{1}}$ and $\smash{\alpha_{2}}$ may need to be higher than the values suggested by Durante (2017) to enforce a greater degree of shrinkage in clusters with few units; this aspect is highlighted in simulation studies in Appendix B.

2.2.1 The Adaptive Gibbs Sampler

An adaptive Gibbs sampler (AGS) is employed when performing inference for MIFA. This dynamically shrinks the loadings matrices (and the infinite scores matrix $\bm{\eta}$) to have finite numbers of columns, by selecting the number of ‘active’ factors. This practically facilitates posterior computation while closely approximating the IFA model, without requiring specification of $\smash{\mathbf{Q}=\left(q_{1},\ldots,q_{G}\right)^{\top}}$. However, a strategy is required for choosing appropriate truncation levels, $\smash{\widehat{q}_{g}}$, that strike a balance between missing important factors and wasting computational effort. For computational reasons, a conservatively high upper bound is used, such that q_{g}^{\star}=\min\negthinspace\big{(}\left\lfloor 3\negthinspace\left(p\right)\right\rceil,N-1,p-1\big{)}\>\>\forall\>\>g. The number of factors in each $\bm{\Lambda}_{g}$ is then adaptively tuned as the MCMC chain progresses. Adaptation can be made to occur only after the burn-in period, in order to ensure the true posterior distribution is being sampled from before truncating the loadings matrices.

At the $t$-th iteration, adaptation occurs with probability $\mathrm{p}\negthinspace\left(t\right)=\exp\negthinspace\left(-b_{0}-b_{1}t\right)$, with $b_{0}\geq 0\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\ b_{1}>0$ chosen so that adaptation occurs often at the beginning of the chain but then decreases exponentially fast in frequency. Here $b_{0}=0.1$ and $\smash{b_{1}=5\times 10^{-5}}$ are used. With probability $\mathrm{p}\negthinspace\left(t\right)$, loadings columns having some pre-specified proportion of elements $\varsigma$ in a small neighbourhood $\epsilon$ of zero are monitored. If there are no such columns, an additional column is added by simulation from the MGP prior. Otherwise redundant columns are discarded and the AGS proceeds with all parameters corresponding to non-redundant columns retained. Choice of $\varsigma\leavevmode\nobreak\ \mbox{and}\leavevmode\nobreak\ \epsilon$ can be delicate, as there is an implicit trade-off between these two fixed tuning parameters; smaller $\varsigma$ and larger $\epsilon$ speed up the algorithm by favouring the discarding of factors during the adaptation step, and vice versa. Typically, $\varsigma$ should be kept close to $1$ and $\epsilon$ should be kept small, relative to the scale of the data. Here, $\varsigma=\lfloor 0.7\times p\rfloor/p$ and $\epsilon=0.1$ are found to strike an appropriate balance. The dimension of the matrix $\bm{\eta}$ of factor scores at a given iteration are set to $p\times\overline{q}=p\times\max\left(\mathbf{Q}\left(t\right)\right)$; rows corresponding to observations currently assigned to a cluster with fewer latent factors than $\overline{q}$ are padded with zeros. Notably, here $\smash{\widehat{q}_{g}}$ may shrink to [math] thus allowing diagonal covariance structure within a component. If this occurs, the decision to simulate a new column is based on a binary trial with probability $1-\varsigma$ as there are no loadings columns to monitor.

The numbers of active factors in each cluster for each retained posterior sample can be used to construct a barchart approximation to the posterior distribution of $q_{g}$. The posterior mode is used to estimate each $q_{g}$, with credible intervals quantifying uncertainty. Another strategy, which circumvents the need to pre-specify $\varsigma$ and $\epsilon$ is to forego adaptation (provided the computational burden of doing so is tolerable) and estimate $\widehat{q}_{g}$ from the number of non-redundant columns in the posterior mean loadings matrices. However, this approach is not considered further here.

In any case, the main advantages of MIFA are that different clusters can be modelled by different numbers of factors and that the model search is reduced to one for $G$ only, as $q_{g}$ is estimated automatically during model fitting. Here, for MIFA models, the optimal $G$ is chosen via the BICM (BIC-Monte (Carlo)) proposed by Raftery et al. (2007), with \smash{\mbox{BICM}=2\ln\big{(}\overline{\mathcal{L}}\big{)}-2s_{l}^{2}\left(\ln\negthinspace\left(N\right)-1\right)}, where $\overline{\mathcal{L}}$ and $s_{l}^{2}$ are the sample mean and sample variance, respectively, of the log-likelihood values calculated for each retained posterior sample. This criterion is particularly useful in the context of nonparametric models where the number of free parameters is difficult to quantify, though we caution that it may be biased in favour of $G=1$ models, under which the log-likelihoods tend to exhibit less variability, and that a large number of posterior samples are required to ensure stable estimation of $s_{l}^{2}$.

2.2.2 Other Infinite Factor Models

This work offers an extension of the MGP prior and its related AGS routine to the mixture modelling context. Wang et al. (2016) develop a related model employing a multiplicative exponential process prior. Other nonparametric approaches to inferring the number of factors include Knowles & Ghahramani (2007), in which a two-parameter Indian Buffet Process (IBP) prior is assumed on an infinite binary matrix underlying the factor scores, thus selecting features of interest, with associated standard Gaussian weights. A closely related approach using the beta process (BP) is provided by Paisley & Carin (2009). In Knowles & Ghahramani (2011) and Ročková & George (2016), an IBP prior is instead assumed for sparsifying the loadings. These models assume a single sparse infinite factor model for the whole data set. However, embedding them in a mixture modelling setting, similar to the IMIFA framework, is intuitively feasible.

Indeed, Chen et al. (2010) employ the BP prior, coupled with a Dirichlet process prior, to perform clustering in a manifold learning setting. While the BP and IBP priors achieve exact sparsity, which may be advantageous in certain applications, the MGP prior has a weaker notion of sparsity by virtue of cumulatively shrinking an infinite series arbitrarily close to zero, thereby preserving small signals. The block updates of each row of $\bm{\Lambda}_{g}$ facilitated by the MGP prior and parameter expansion mean the AGS approach is a simpler, more computationally efficient alternative to the BP and IBP priors.

2.3 Overfitted Mixtures of (Infinite) Factor Analysers

While MIFA obviates the need to pre-specify $\mathbf{Q}$, the issue of model choice is not yet fully resolved. Overfitted mixtures (Rousseau & Mengersen, 2011; van Havre et al., 2015) are one means of extending MIFA; indeed Papastamoulis (2018) proposes an overfitted mixture of factor analysers (OMFA), albeit with finite factors. Here, the overfitted mixture of infinite factor analysers (OMIFA) model is introduced.

In overfitted mixtures the symmetric Dirichlet prior on $\bm{\pi}$ plays an important role. Estimation is approached by initially overfitting the number of clusters expected to be present. Small values of the hyperparameter $\bm{\alpha}$ encourage emptying out excess components in the posterior distribution; the uniform prior with $\bm{\alpha}=\mathbf{1}$ is rather indifferent in this respect. The sampler is initialised with a conservatively high number of components: \smash{G^{\star}=\max\negthinspace\big{(}\left\lceil 3\ln\negthinspace\left(N\right)\right\rceil,25,N-1\big{)}}, though this may be too high if it is close to $N$. While $\widetilde{G}=G^{\star}$ remains fixed throughout the MCMC chain, the number of non-empty clusters is recorded at each iteration of the sampler as G_{0}=\widetilde{G}-\sum_{g=1}^{\widetilde{G}}\mathds{1}{{\big{(}\sum_{i=1}^{N}z_{ig}=0\big{)}}} where $\mathds{1}{{\left(\cdot\right)}}$ is the indicator function. The true $G$ is estimated by $\widehat{G}$, the $G_{0}$ value visited most often. Cluster-specific inference is conducted only on samples corresponding to those visits. For the OMIFA model, the AGS is modified to handle empty components: the MGP-related parameters are simulated from the relevant priors and each corresponding $\bm{\Lambda}_{g}$ matrix is restricted to having $\overline{q}$ factors, i.e. the same number of columns currently in the matrix of factor scores $\bm{\eta}$, either by truncation or by padding with zeros, as required.

2.4 Infinite Mixtures of (Infinite) Factor Analysers

Embedding MFA and MIFA in an infinite mixture setting leads, respectively, to the infinite mixture of finite factor analysers model (IMFA) and the flagship infinite mixture of infinite factor analysers model (IMIFA). These models employ a nonparametric Pitman-Yor process (PYP) prior which is easily incorporated into the MCMC sampling scheme.

The PYP is a stochastic process whose draws are discrete probability measures, whereby $H\sim\textrm{PYP}\negthinspace\left(\alpha,d,H_{0}\right)$ denotes a PYP probability distribution $H$, with base distribution $H_{0}$ interpreted as the mean of the PYP, discount parameter $d\in\left[0,1\right)$, and concentration parameter $\alpha>-d$. For the PYP mixture model IMFA and the PYP-MGP mixture model IMIFA $H_{0}$ comes from the factor-analytic mixture (1), hence

[TABLE]

The stick-breaking representation of the PYP (Pitman, 1996) is used as a prior process for generating the mixing proportions in (2). This construction views $\left\{\pi_{1},\pi_{2},\ldots\right\}$ as pieces of a unit-length stick that is sequentially broken in an infinite process, with stick-breaking proportions $\bm{\Upsilon}=\left\{\upsilon_{1},\upsilon_{2},\ldots\right\}$, summarised as

[TABLE]

where $\smash{\delta_{\bm{\theta}_{g}}}$ is the Dirac delta centred at $\smash{\bm{\theta}_{g}}$, such that draws are composed of a sum of infinitely many point masses. The PYP reduces to the DP when $d=0$, in which case mass shifts to the right with increasing dispersion as $\alpha$ increases, implying an a priori larger number of components. However, some important distributional features fundamentally differ when $d\neq 0$ (De Blasi et al., 2015). The PYP exhibits heavier tail behaviour and allows the stick-breaking distribution to vary according to the component index $g$, without sacrificing much in the way of tractability. In particular, increasing $d$ values have the effect of flattening the prior, controlling its degree of non-informativity.

Slice sampling (Walker, 2007; Kalli et al., 2011) is used here to yield samples from the PYP by adaptively truncating the number of components needed to be sampled at each iteration. By introducing an auxiliary variable $u_{i}>0$ which preserves the marginal distribution of the data, and denoting by $\bm{\xi}=\left\{\xi_{1},\xi_{2},\ldots\right\}$ a positive sequence of infinite quantities which sum to $1$, the joint density of $\left(\mathbf{x},\mathbf{u}\right)$ is given by $f\negthinspace\left(\mathbf{x},\mathbf{u}\,|\,\bm{\theta},\bm{\xi}\right)=\sum_{g=1}^{\infty}\pi_{g}\textrm{Unif}\negthinspace\left(\mathbf{u};0,\xi_{g}\right)f\negthinspace\left(\mathbf{x}\,|\,\bm{\theta}_{g}\right)$. Since only a finite number of $\xi_{g}$ are greater than $\mathbf{u}$, the conditional density of $\mathbf{x}\,|\,\mathbf{u}$ can be written as a finite mixture with $\widetilde{G}=\max_{i\,\in\,\left\{1,\ldots,N\right\}}\left(\lvert\bm{\mathcal{A}}_{\bm{\xi}}\negthinspace\left(u_{i}\right)\rvert\right)$ ‘active’ components at each iteration, where $|\cdot|$ denotes cardinality and $\bm{\mathcal{A}}_{\bm{\xi}}\negthinspace\left(\mathbf{u}\right)=\{g\colon\mathbf{u}<\xi_{g}\}$. Though $G$ is infinite in theory, $\widetilde{G}$ can be at most equal to $N$. Thus, the infinite mixture of (infinite) factor analysers models can be sampled from.

Typical implementations of the slice sampler arise when $\xi_{g}=\pi_{g}$ (Walker, 2007) but independent slice-efficient sampling (Kalli et al., 2011) allows for a deterministic decreasing sequence, e.g. geometric decay, given by $\xi_{g}=\left(1-\rho\right)\rho^{g-1}$ where $\rho\in\left[0,1\right)$ is a fixed value to be chosen with care. Higher values generally lead to better mixing but longer run-times, as the average cardinality of $\bm{\mathcal{A}}_{\bm{\xi}}\negthinspace\left(\mathbf{u}\right)$ increases, and vice versa. Setting $\rho=0.75$, in line with the recommendations of Kalli et al. (2011), appears to strike an appropriate balance in the applications considered here.

2.4.1 Inference for Infinite Mixtures of Factor Analysers Models

For clarity, what follows focuses on the IMIFA model where inference proceeds via the independent slice-efficient sampler with geometric decay. Inference for other models in the IMIFA family is closely related. The joint density of the IMIFA model is

[TABLE]

where $\mathrm{B}\negthinspace\left(\cdot\right)$ is the Beta function and $f(\bm{\theta})$ is the product of the previously defined collection of conditionally conjugate priors with additional layers for hyperparameters. Only the parameters of the $\widetilde{G}$ active components are sampled at each iteration. The algorithm is initialised with the same $G^{\star}$ value detailed in Section 2.3, typically above the anticipated number to which the algorithm will converge, in the spirit of Hastie et al. (2014). Here, however, $\widetilde{G}$ can theoretically exceed this value. For computational reasons, a finite upper limit is placed on $\widetilde{G}$ with $\max\negthinspace\left(G^{\star},\min\negthinspace\left(N-1,50\right)\right)$ found to be sufficiently large. However, $\widetilde{G}$ is only regarded as a set of proposals as to where to allocate observations; as in Section 2.3, it is the subset of non-empty clusters $G_{0}$ that is of inferential interest.

Bayesian approaches to clustering are known to be sensitive to initial cluster allocations. While starting values for $\mathbf{z}_{i}$ can be obtained by any means, model-based agglomerative hierarchical clustering (Scrucca et al., 2016) is used here. Though this is fast and intuitive given that IMIFA models are initialised at a conservatively high number of components, which are then merged as the sampler proceeds, heavily imbalanced initial cluster sizes are cautioned against. By extension, initial cluster means and mixing proportions are computed empirically. Other parameter starting values are simulated from their relevant prior distributions. The adaptive inferential algorithm for IMIFA then proceeds mostly via Gibbs updates. For those which are multivariate Gaussian, using the Cholesky factor of the covariance matrices and employing block updates speeds up the algorithm (Rue & Held, 2005). The allocations $\mathbf{z}_{i}$ are sampled in a fast, numerically stable fashion, using the Gumbel-Max trick (Yellott, 1977). Finally, state spaces for applications of IMIFA to real data can be highly multimodal with well-separated regions of high posterior probability coexisting, corresponding to clusterings with different numbers of components. Thus, label switching moves (Papaspiliopoulos & Roberts, 2008) are incorporated in order to improve mixing. Details of the Gibbs updates, Gumbel-Max trick, and label switching moves are provided in Appendix A.

2.4.2 Assessing Model Fit and Mixing

As is good statistical practice, posterior predictive model checking (Gelman et al., 2004) is employed. Sampled model parameters from the MCMC chain are used to generate replicate data from the posterior predictive distribution. Valid posterior samples, after conditioning on $\smash{\widehat{G}}$, are those for which $\max\negthinspace\left(\mathbf{Q}\left(t\right)\right)\negthinspace\geq\max\negthinspace\left(\widehat{q}_{1},\ldots,\widehat{q}_{\widehat{G}}\right)$ such that the dimension of the estimated scores matrix $\smash{\widehat{\bm{\eta}}}$ is preserved. To assess model fit, histograms of the modelled data $\mathbf{X}$ are compared to histograms of the replicate data in a global sense using the Posterior Predictive Reconstruction Error (PPRE), calculated as follows:

Gather the histogram bin counts of each variable in $\mathbf{X}$ into the $h\times p$ matrix $\bm{\mathcal{H}}$, where $h$ is the maximum number of bins across all variables and $\bm{\mathcal{H}}$ is padded with zeros as required. 2. 2.

Generate $r\in\left\{1,\ldots,R\right\}$ data sets $\bm{\mathcal{X}}^{\left(r\right)}$ from the posterior predictive distribution. 3. 3.

Create a similar matrix of histogram bin counts $\bm{\mathcal{H}}^{\left(r\right)}$ for each $\bm{\mathcal{X}}^{\left(r\right)}$ using the same break-points with which $\bm{\mathcal{H}}$ was constructed (with endpoint bins extended to $\pm\,\infty$). 4. 4.

Compute the Frobenius norm ${\left\lVert\cdot\right\rVert}_{\mathcal{F}}$ between $\bm{\mathcal{H}}$ and $\bm{\mathcal{H}}^{\left(r\right)}$, standardising to the 0-1 scale using the triangle inequality:

[TABLE]

The distribution of PPRE values can be visualised using boxplots and summarised by the median, with credible intervals quantifying uncertainty. This discrepancy measure is well-suited to assessing model adequacy for mixtures of multivariate data: it accounts for inherent multimodality and gives a global quantitative measure of agreement between the distributions of the observed variables and their posterior predictive counterparts.

Convergence of the MCMC chains is assessed using the potential scale reduction factor (PSRF; Brooks & Gelman, 1998; Plummer et al., 2006). Random allocations of the initial cluster labels, resulting in different draws from the relevant priors for parameter initialisation, are used to construct the multiple overdispersed chains required. The MAP labels of each chain are matched to the main chain prior to computing the diagnostics; $\bm{\Lambda}_{g}$ matrices are also rotated to a common template for each cluster. Good convergence is indicated by upper PSRF $95\%$ confidence interval limits close to $1$; this is a stricter requirement than the PSRF values themselves being near $1$.

2.4.3 Comparing the IMIFA Family Models

Though IMIFA and OMIFA come with the computational complexities inherent in nonparametric methods, diminishing adaptation, and extra tuning parameters, their advantages over other models in the IMIFA family are numerous: i) flexibility, in the sense that models where $q_{g}\neq q_{g}^{\prime}$ can be fitted, ii) computational efficiency, in the sense that the burden is reduced relative to searching over a range of fitted MFA or MIFA models, iii) removing the need for model selection criteria, and iv) the ability to quantify the uncertainty in $\smash{\widehat{G}}$ and $\smash{\widehat{q}_{g}}$. Both methods offer simpler alternatives to reversible jump MCMC (Richardson & Green, 1997) and birth-death MCMC (Stephens, 2000). Hence, among the IMIFA family, the infinite factor models are recommended over the finite factor models and the infinite and overfitted mixtures are recommended over the finite mixtures. However, the MIFA model is appropriate if one wishes to fix $G$ but infer $q_{g}$.

While infinite mixtures are often used for density estimation, they are also employed to infer the number of components in cluster analyses (e.g. Kim et al. 2006; Xing et al. 2006; Yerebakan et al. 2014). However, Miller & Harrison (2013, 2014) raise concerns about the guarantee of posterior consistency for the number of non-empty clusters, showing the number uncovered is typically greater than or equal to the truth, often with several vanishingly small clusters inferred. These concerns highlight the need for practitioners to pay due consideration to the uncertainty in the number of clusters offered by IMIFA models. Relatedly, Frühwirth-Schnatter & Malsiner-Walli (2019) compare infinite mixtures to overfitted (‘sparse finite’) mixtures. They highlight that overfitted mixtures are useful for applications in which the data arise from a moderate number of clusters, even as the sample size increases, whereas infinite mixtures are suited to cases where the number of clusters also increases. However, they show that clustering results are driven less by the assumption of whether the data arose from a finite or infinite mixture, but by the hyperprior on the DP parameters or the sparseness of the Dirichlet prior in the overfitted setting. Indeed, they show that overfitted and infinite mixtures yield comparable clustering performance on the observed data when these hyperpriors are matched. This matching leads to ‘sparse’ infinite mixtures that avoid overfitting the number of clusters. Similar behaviour is observed for the PYP prior in the applications in Section 3, where the IMIFA and OMIFA models, with matched hyperpriors, give comparable results.

The issue of choosing $\alpha$ can make implementing overfitted models challenging. With fixed $\alpha=\nicefrac{{\gamma}}{{G^{\star}}}$, the prior approximates a DP with concentration parameter $\gamma$ as $G^{\star}$ tends to infinity (Green & Richardson, 2001). Here, following Frühwirth-Schnatter & Malsiner-Walli (2019), a $\textrm{Ga}\negthinspace\left(a,bG^{\star}\right)$ hyperprior is assumed for $\alpha$. This favours small values and allows $\alpha$ to be updated via Metropolis-Hastings. In the infinite mixture setting, learning the PYP parameters (which also requires Metropolis-Hastings steps) and adopting the label-switching moves enables accurate inference on $G_{0}$. A joint hyperprior $\mathrm{p}\negthinspace\left(\alpha,d\right)=\mathrm{p}\negthinspace\left(\alpha\,|\,d\right)\mathrm{p}\negthinspace\left(d\right)$ is assumed (Carmona et al., 2019) where $\mathrm{p}\negthinspace\left(\alpha\,|\,d\right)=\textrm{Ga}\negthinspace\left(\alpha+d;a,b\right)$; choosing a large $b$ encourages clustering (Müller & Mitra, 2013). A spike-and-slab hyperprior $d\sim\kappa\delta_{0}+\left(1-\kappa\right)\textrm{Beta}\negthinspace\left(a^{\prime},b^{\prime}\right)$ is assumed. The estimated proportion $\smash{\widehat{\kappa}}$ can then be used to assess whether the data arose from a DP or a PYP at little extra computational cost. See Appendix A for further details.

3 Illustrative Applications

The flexibility and performance of the IMIFA model and its related model family are demonstrated below through application to benchmark and real data sets. All results are obtained through the IMIFA R package; code to reproduce many of the results is available in the associated vignette111 https://cran.r-project.org/web/packages/IMIFA/vignettes/IMIFA.html.. Appendix B reports on simulation studies demonstrating the performance of IMIFA under different scenarios, including effects of the $\nicefrac{{N}}{{p}}$ ratio, the PYP parameters, imbalanced cluster sizes, uncommon $q_{g}$, and the degree of loadings sparsity. Appendix C explores the robustness of IMIFA.

MCMC chains were run for $50,000$ iterations, except for Section 3.3 in which $20,000$ were run. Every 2nd sample was thinned and the first $20\%$ of iterations were discarded as burn-in. All computations were performed on a Dell Latitude 5491 laptop, equipped with a 6-core 2.60 GHz Intel Core i7-8850H processor and 16 GB of RAM. Where necessary, the optimal finite and infinite factor models are chosen by the BIC-MCMC and BICM criteria, respectively. Throughout, $\smash{\widehat{\cdot}}$ denotes the posterior mode, posterior mean, or relevant optimal value. Unless otherwise stated, data were mean-centred and unit-scaled and no constraints were imposed on the uniquenesses. Hyperprior specifications are detailed in Table 1. While there are many hyperparameters to select, the choices are all reasonably standard. However, poor settings may introduce additional factors or clusters to maintain flexibility and so care in specifying hyperparameters is advised.

3.1 Benchmark Data: Italian Olive Oils

The Italian olive oil data (Forina & Tiscornia, 1982; Forina et al., 1983) is often clustered using factor-analytic models, e.g. McNicholas (2010). The data detail the percentage composition of $8$ fatty acids in $572$ Italian olive oils, known to originate from three areas: southern and northern Italy and Sardinia. Each area is composed of different regions: southern Italy comprises north Apulia, Calabria, south Apulia, and Sicily; Sardinia is divided into inland and coastal Sardinia; and northern Italy comprises Umbria and east and west Liguria. Hence, the true number of clusters is hypothesised to correspond to either $3$ areas or $9$ regions.

The full family of IMIFA models is fitted to the olive oil data with results detailed in Table 2. Models relying on pre-specification of finite ranges of $G$ and/or $q$ are based on $G=1,\ldots,9$ and $q=0,\ldots,6$. Clustering performance is evaluated using the adjusted Rand index (ARI; Hubert & Arabie, 1985) and the misclassification rate, compared to the $3$ area labels. The $\alpha$ parameter is reported as its fixed value or posterior mean, as appropriate. Table 2 shows the flexibility and accuracy of the developed model family, and of the IMIFA model in particular which has the best clustering performance. Additionally, IMIFA is the most computationally efficient model considered, among those in the IMIFA family achieving clustering, as it requires only one run. This speed improvement would be exacerbated with larger data sets. However, methods requiring fitting of multiple models were run here in series; parallel implementations would reduce runtimes. Finally, models with different numbers of cluster-specific factors show improved clustering performance compared to the corresponding finite factor model in every case.

Table 2 also shows that the performance of the IMIFA model compares favourably to the best parsimonious Gaussian mixture model, fit via the pgmm R package (McNicholas et al., 2019) and the best mixture of factor mixture analysers (MFMA) model (Viroli, 2010), evaluated with $1,\ldots,5$ components in both layers. Models with zero factors were not considered in either case. IMIFA also outperforms the best constrained Gaussian mixture model fitted using mclust (Scrucca et al., 2016). These finite mixtures are fit via maximum likelihood and use the BIC for model selection after fitting a large number of candidate models.

It is also notable that within the set of IMIFA models relying on information criteria, those deemed optimal were not necessarily optimal in a clustering sense. For instance, the $4$-cluster MIFA model yields an ARI of $0.94$ and a misclassification rate of $6.99\%$, with respect to the $3$ area labels, despite its sub-optimal BICM. Similarly, the BICM and BIC-MCMC criteria suggest different optimal MFA models. For the IMIFA model $\smash{\widehat{\kappa}\approx 0.89}$, suggesting similar inference would have resulted under a DP prior. Indeed, the results obtained by the OMIFA and OMFA models are similar to those of their infinite mixture counterparts, though the latter provide a better fit to the data (see Figure 5).

Figure 3 shows a barchart approximation to the posterior distribution of $G$ under the IMIFA model. The modal value of $4$, visited in $\approx 90\%$ of posterior samples, is used as the estimate of the true number of clusters (with $95\%$ credible interval $[4,5]$). Table 3a tabulates the MAP clustering against the $3$ area labels and suggests this solution makes geographic sense, in that northern oils are cleanly split into two sub-clusters. Cluster 1 contains all of the $323$ southern Italy oils: this large cluster requires the largest number of factors ($\widehat{q}_{1}=6\leavevmode\nobreak\ [5,6]$, with $95\%$ credible intervals in brackets). Some of the other clusters require notably fewer ($\widehat{q}_{2}=3\leavevmode\nobreak\ [1,6]$, $\widehat{q}_{3}=6\leavevmode\nobreak\ [3,6]$, and $\widehat{q}_{4}=2\leavevmode\nobreak\ [1,4]$). Table 3b gives the confusion matrix with oils from the north labelled by their associated region(s), yielding an ARI of $0.994$ and a misclassification rate of $0.52\%$. Figure 3 shows the uncertainty in the allocations to these clusters. Only three oils have large probability of belonging to a cluster other than the one to which they were assigned by the IMIFA model.

To assess sensitivity to starting values, the IMIFA model was re-fitted using multiple random initial allocations, implying also different random draws from the priors for parameter starting values. These runs led to identical inference about $\widehat{G}$ and $\widehat{\mathbf{Q}}$ and equivalent clustering performance. These overdispersed chains were used to compute the upper $95\%$ PSRF confidence limits depicted in Figure 5, which indicate good convergence. The PPRE boxplots in Figure 5 demonstrate the superior fit of the IMIFA model (with a median PPRE of $0.10$) to the olive oil data, compared to the other IMIFA family models. Histograms comparing the bin counts between the modelled and replicate data sets for each variable, under the IMIFA model, are given in Appendix D.

3.2 Spectral Metabolomic Data

IMIFA is employed to cluster spectral metabolomic data for which $N\ll p$ (Figure 6). The data are nuclear magnetic resonance spectra consisting of $p=189$ spectral peaks from urine samples of $N=18$ participants, half of which are known to have epilepsy (Carmody & Brennan, 2010; Nyamundanda et al., 2010). Interest lies in uncovering any underlying clustering structure given the $N\ll p$ setting.

Data were mean-centred and Pareto scaled (van den Berg et al., 2006). Although $N\ll p$, no restrictions are imposed on the uniquenesses as the sample variances are quite imbalanced. Fitting MIFA models for $G=1,\ldots,5$ is feasible as $N$ is small. The BICM criterion chooses $\smash{\widehat{G}=2}$ as optimal and one participant is misclassified. IMIFA, however, unanimously visits a 2-cluster model and perfectly uncovers the group structure.

The modal estimates of the number of factors in each IMIFA cluster are $\widehat{q}_{1}=3\leavevmode\nobreak\ [2,9]$ and $\widehat{q}_{2}=5\leavevmode\nobreak\ [4,13]$ (see Figure 7). Cluster 1 corresponds to the control group and Cluster 2 to the epileptic participants. Figure 8 illustrates the $\smash{p\times\widehat{q}_{g}}$ posterior mean loadings matrices, based on retained samples with $\smash{\widehat{q}_{g}}$ or more factors, after Procrustes rotation to a common template for both clusters. The sparsity and shrinkage induced by the MGP prior is apparent, as is the greater complexity in Cluster 2, given the greater variation in colour and larger number of factors. For instance, many elevated loadings are visible for chemical shift values between $8$ and $10$ for the first two factors in Cluster 2; this activity is not present for other factors in either cluster. In general, the distributions of the loadings within a factor exhibit narrow spread around zero, particularly for the cluster of control participants, with the exception of the regions of the spectrum corresponding to the large peaks between chemical shifts of $3$ and $5$ in Figure 6.

IMIFA outperforms the optimal $\smash{\widehat{G}=3}$ mclust model and the optimal $\smash{\widehat{G}=2}$, $\smash{\widehat{q}=5}$ pgmm model, with respective ARI values of $0.73$ and $0.27$. The clustering performance of the optimal MFMA model is identical to the optimal MIFA model described above. Given the $N\ll p$ nature of the data, spectral clustering with the Gaussian kernel (Ng et al., 2001) is also considered. The eigengap heuristic suggests $\smash{\widehat{G}=2}$ and a perfect clustering is achieved almost instantaneously. However, the approach does not characterise the uncovered clusters in an interpretable manner, nor provide estimates of cluster membership uncertainty as given by model-based clustering approaches such as IMIFA.

The median PPRE for the IMIFA model of $0.21\leavevmode\nobreak\ [0.18,0.24]$ shows good model fit, given the size and dimensionality of the data. The median PSRF upper $95\%$ confidence limits, using three randomly initialised auxiliary chains, for the cluster means, uniquenesses, loadings, and mixing proportions of $1.01\,(0.01)$, $1.00\,(<0.01)$, $1.01\,(0.08)$, and $1.00\,(<0.01)$ respectively, show good mixing also (standard deviations in parentheses). Notably, all chains yield the same inference about $\smash{\widehat{G}}$ and $\smash{\widehat{\mathbf{Q}}}$. So too, again, does the OMIFA model, although its model fit is inferior (median PPRE=$0.26$).

3.3 Handwritten Digit Data

A final illustration of IMIFA is given through its application to handwritten digit data from the United States Postal Service (USPS; Le Cun et al., 1990; Hastie et al., 2001). Here $N=7,291$ images of the digits $0,\ldots,9$ are considered, taken from handwritten zip codes. The data are not balanced in terms of digit labels. Each image is a $16\times 16$ grayscale grid concatenated into a $p=256$-dimensional vector; data were mean-centred but not scaled. Such data are often considered in the context of manifold learning, positing that the data dimensionality is artificially high. Given $N$ and $p$, fitting a range of MFA or MIFA models is practically infeasible. Results of a single IMIFA run are presented here. For these data, it is reasonable to expect the number of components to grow as the sample size grows. It is anticipated that the flexibility afforded by having cluster-specific numbers of factors will help characterise digits with different geometric features.

The IMIFA model visited a $\widehat{G}=21$ cluster solution in all posterior samples; Table 4 cross-tabulates the MAP clustering against the known digit labels and achieves an ARI of $0.33$. The median PPRE of $0.05$ $[0.04,0.06]$ indicates good model fit. The overdispersed chains used to compute the PSRF diagnostics lead to identical inference about the number of clusters but slightly different inference about the modal numbers of cluster-specific factors. The ARI values between each resulting pair of MAP partitions all exceed $0.93$. As before, good mixing is indicated by median PSRF upper $95\%$ confidence limits for the cluster means, uniquenesses, and mixing proportions of $1.01\,(0.01)$, $1.01\,(0.01)$, and $1.01\,(<0.01)$, respectively. In computing the diagnostic for the loadings ($1.14\,(0.35)$), only the first factor (common to all loadings matrices across all clusters in all chains) was considered, for reasons of fairness and computational resource constraints.

Generally, IMIFA assigns images of the same digit, albeit written differently, to different clusters. Posterior mean images for each cluster are shown in Figure 9, ordered, as is Table 4, from 0 to 9 according to the digit most frequently assigned to the related cluster. Cluster 7 and the smaller cluster 8 capture the digit 1 written in a straight and slanted fashion, respectively. Clusters 15–17 represent the digit 6 written with extended, medium, and compact loop curvature, respectively. Notably, cluster 15 requires more factors than clusters 16 and 17. A similar interpretation follows for clusters 20 and 21 ($\smash{\widehat{q}_{20}=2,\widehat{q}_{21}=3}$), mostly capturing the digit 9 with small and large loops, respectively. Cluster 19 appears to mostly represent the digit 8 and has a large number of factors ($\smash{\widehat{q}_{19}=6}$) in comparison, say, to clusters 7 and 8 ($\smash{\widehat{q}_{7}=2,\widehat{q}_{8}=1}$) which capture the digit 1. This is intuitive, as 8 is a more geometrically complex digit than 1. However, some clusters appear to be diluted by the confusion of the so-called ‘closed’ 4, in contrast to the ‘open’ 4 in cluster 12, with the digits 3, 5, and 8 (cluster 19) and the digits 7 (written with a horizontal bar) and 9 (clusters 20 and 21). Many clusters capture the most common digit 0, with differing degrees of elongation and border thickness. Of concern here is cluster 4, containing just 9 units; the fact that $\smash{\widehat{q}_{4}=16}$, the upper AGS limit, suggests that the model struggles to shrink the number of factors in poorly populated clusters. This difficulty is highlighted further in the simulation studies in Appendix B. Finally, Table 4 indicates that clusters 10, 13, and 14 also capture several other digits, all of which are reflected in the blurriness of the resulting posterior mean images and in $\smash{\widehat{q}_{10}}$, $\smash{\widehat{q}_{13}}$, and $\smash{\widehat{q}_{14}}$ being quite large. The cluster-membership uncertainties are visualised in Appendix D.

It is computationally infeasible to run mclust, pgmm, or MFMA on these large data, as an exhaustive model search would be too vast. For comparative purposes, a DP-BP model (Chen et al., 2010) is fitted; this approach also simultaneously assumes infinitely many components and factors. It finds $43$ clusters, each with around $14$ factors, and achieves an ARI of $0.32$. Cross tabulating this clustering against the $21$ clusters of the IMIFA model shows that some of the DP-BP clusters are encapsulated by the larger IMIFA clusters. IMIFA is thus the more parsimonious approach and affords greater cluster-specific factor flexibility. Additionally, a finite mixture of matrix-normal distributions (Viroli, 2011) is also fitted. This approach accounts for the grid nature of the data, but is computationally infeasible for $G>15$ and requires a model selection strategy. The optimal model according to BIC yields $\smash{\widehat{G}=12}$ and $\mbox{ARI}=0.38$. While neither IMIFA nor the DP-BP model account for the spatial structure in the data, they demonstrate comparative performance without the need for a computationally expensive model search.

4 Discussion

The proposed IMIFA model is a Bayesian nonparametric approach clustering high-dimensional data using factor-analytic mixture models. By extending the MGP prior (Bhattacharya & Dunson, 2011) to the PYP-MGP setting, the model sidesteps the fraught and computationally intensive task of determining the optimal number of clusters and factors using model selection criteria. Thus, the IMIFA model is recommended when fitting factor-analytic mixtures in settings where an exhaustive model search is computationally infeasible. Though IMIFA is not entirely choice-free, it achieves improved clustering results by allowing factor-analytic models of different dimensions in different clusters. If small clusters are inferred, one may wish to prune or merge small clusters with the larger clusters (West et al., 1994) or assess whether the small clusters are in fact of domain-specific interest. While comparative performance can be achieved by the IMIFA and OMIFA models, one may wish to fit a MIFA or OMIFA model when the expectation is that the number of clusters is fixed or unlikely to grow with $N$, respectively.

Future research directions are varied and plentiful. Incorporating covariates, in the spirit of Bayesian factor regression models (West, 2003; Carvalho et al., 2008), would allow for direct inclusion of the weight and urine pH covariates available with the metabolomic data, for example. Furthermore, the models could be extended to the (semi-)supervised model-based classification setting where all (or some) of the data are labelled. While constraints on the uniquenesses across variables and/or clusters are allowed, there is scope for also constraining the loadings across clusters. Though the number of factors would no longer be cluster-specific, the common number of loadings columns would be estimated in a similarly automatic fashion. However, incorporating covariance matrix constraints in the IMIFA model family problematically reintroduces the need for model selection strategies, in order to choose between them, though the BICM criterion could feasibly be used for this purpose also.

As proposed by Bhattacharya & Dunson (2011), the MGP hyperparameters could be learned via Metropolis-Hastings, and thus also be made cluster-specific. This could help combat some difficulties identified in the simulation studies in Appendix B. For example, learning those related to local shrinkage may help when loadings are notably dense. Learning those related to column shrinkage may help in settings with many small clusters, where IMIFA struggles to adaptively truncate loadings columns. In principle, a further global shrinkage parameter $\varpi$ could be added to the MGP prior to borrow information across clusters, i.e. \smash{\lambda_{jkg}\,|\,\ldots\sim\textrm{N}_{1}\big{(}0,\phi_{jkg}^{-1}\tau_{kg}^{-1}\sigma_{g}^{-1}\varpi^{-1}\big{)}}. Alternatively, the infinite factor prior of Legramanti et al. (2020) could be employed, which decouples control over the cumulative shrinkage rate and the active loadings terms. The IBP prior (Knowles & Ghahramani, 2007, 2011; Ročková & George, 2016) could be also extended to the infinite mixture setting, as per the DP-BP model of Chen et al. (2010). Thus, the IMIFA family can in fact be considered to be wider than the range of models presented here, with potential for further expansion.

For applied problems, a mismatch between the assumed model and the data distribution will impact inference. Miller & Harrison (2013, 2014) highlight that posterior consistency for the number of non-empty clusters in infinite mixtures is contingent on correct specification of the component distributions. While they do not discourage the use of infinite mixtures for clustering, they show that a few tiny extra clusters are typically fitted and suggest robustifying inference. If the data distribution is close to but not exactly a finite mixture of Gaussians, an infinite Gaussian mixture will introduce more components as the amount of data increases. Potential avenues of exploration thus include considering the IMIFA model with the heavy tailed multivariate $t$-distribution (Peel & McLachlan, 2000). Similarly, modelling of complex component distributions can be achieved by considering the MFMA approach in the context of infinite factor models. Defining robust inference functions as in Lee & MacEachern (2014) or using nonparametric unimodal component distributions as in Rodriguez & Walker (2014) may also prove fruitful. Another means of robustifying inference is to explicitly include a noise component with zero factors to capture outliers which depart from the component multivariate normality assumption. Finally, a ‘coarsened’ posterior (Miller & Dunson, 2018) could be used for addressing misspecification, by conditioning on the event that the model generates data close to the observed data in a distributional sense.

Acknowledgements

This research was supported by Science Foundation Ireland (SFI/12/RC/2289_P2). The authors thank the members of the UCD Working Group in Statistical Learning and Prof. Adrian Raftery’s Working Group in Model-based Clustering and Prof. David Dunson for helpful discussion. The authors also thank Prof. Lorraine Brennan (UCD), for the metabolomic data, and the anonymous reviewers for constructive feedback from which this work greatly benefitted.

Appendices

A Posterior Conditional Distributions:

Technical details for sampling from the IMIFA model

The structure of the Metropolis-within-Gibbs sampler to conduct inference for the IMIFA model and the exact forms of the required conditional distributions are detailed below. Note that $\textrm{Ga}\negthinspace\left(\alpha,\beta\right)$ refers throughout to the gamma distribution with mean $\nicefrac{{\alpha}}{{\beta}}$. The number of observations in a component is denoted by $n_{g}$, where $\smash{\mathbf{n}=\left(n_{1},\ldots,n_{\widetilde{G}}\right)^{\top}}$ sums to $N$, where $\widetilde{G}$ is the current number of active components (of which some may be empty), and $\widetilde{q}_{g}$ is the current sample of the cluster-specific number of active factors. Algorithms for sampling other models in the IMIFA family can all be considered as special cases of what follows. The algorithm is implemented in the associated R package IMIFA (Murphy et al., 2021).

For $g=1,\ldots,\widetilde{G}$:

[TABLE]

Here $\mathbf{x}^{(j)}$ denotes the $j$-th column of the data matrix, $\lambda_{jkg}^{2}$ denotes a single squared loading, and $\smash{\tau_{kg}=\prod_{h=1}^{k}\delta_{hg}}$ is updated after every update of $\delta_{hg}$.

Parsimonious parameterisations of the component covariance matrices are easily incorporated. Uniquenesses can be constrained to be isotropic, with $\bm{\Psi}_{g}=\operatorname{diag}\negthinspace\left(\psi_{g},\ldots,\psi_{g}\right)$, leading to a model that corresponds to an infinite mixture and infinite-dimensional extension of probabilistic principal components analysers (Tipping & Bishop, 1999). Uniquenesses can also be constrained across clusters, with or without the isotropic constraint across variables. These restrictions define the models in the pgmm family (McNicholas & Murphy, 2008) named UUC, UCU, and UCC, respectively, to which the following Gibbs updates are Bayesian analogues

[TABLE]

In the contexts of finite and overfitted mixtures (i.e. MFA, MIFA, OMFA, and OMIFA) $\mathbf{z}_{i}\,|\,\mathbf{x}_{i},\ldots\sim\textrm{Mult}\negthinspace\left(1,p_{i1},\ldots,p_{i\widetilde{G}}\right)$, with

[TABLE]

whereas under the IMIFA and IMFA models

[TABLE]

The allocations $\mathbf{z}_{i}$ are sampled in a fast, numerically stable fashion, using the unnormalised log-probabilities and independent draws from the standard Gumbel distribution (Yellott, 1977) via $s_{ig}=-\ln\negthinspace\left(m_{ig}\right)$, with $m_{ig}\sim\textrm{Exp}\negthinspace\left(\lambda=1\right)$. Observation $i$ is assigned the label $g$ satisfying

[TABLE]

For the IMIFA and IMFA models, the sampler need only find the maximum over, and only draw Gumbel noise for, log-probabilities for which the indicator function in (3) evaluates to $1$.

Sampling the parameters of the PYP for non-zero $d$ values requires the introduction of Metropolis-Hastings steps within the Gibbs sampler. A joint hyperprior of the form $\mathrm{p}\negthinspace\left(\alpha,d\right)=\mathrm{p}\negthinspace\left(d\right)\mathrm{p}\negthinspace\left(\alpha\,|\,d\right)$ is assumed, as per Jara et al. (2010). Firstly, the hyperprior for the discount parameter $d$ is similar to the one assumed by Carmona et al. (2019); a mixture of a point-mass at zero and a continuous beta distribution, in order to consider the DP special case with $d=0$ with positive probability, i.e. $d\sim\kappa\delta_{0}+\left(1-\kappa\right)\textrm{Beta}\negthinspace\left(a^{\prime},b^{\prime}\right)$. This facilitates explicit comparison between DP models and encompassing PYP alternatives. Secondly, the hyperprior for $\alpha$ is given conditionally on $d$, s.t. $\left(\alpha\,|\,d\right)\sim\textrm{Ga}\negthinspace\left(\alpha+d;a,b\right)$, and includes the constraint $\alpha>-d$ by shifting the support of the gamma density to the interval $\left(-d,\infty\right)$; choosing a large $b$ value is particularly relevant as it encourages clustering (Müller & Mitra, 2013).

The likelihood for $\alpha$ and $d$ is given by the exchangeable partition probability function induced by the PYP (Pitman, 1995). Thus, the required conditional posterior distributions are

[TABLE]

Sampling from the distributions in (4) and (5), while always considering the support $\alpha>-d$, proceeds as per Carmona et al. (2019); a Metropolis-Hastings step is implemented for the discount parameter with independent proposal distribution $0.5\delta_{0}+0.5\textrm{Beta}\negthinspace\left(d;1,1\right)$, and a random walk Metropolis-Hastings step with proposal distribution given by $\alpha^{\star}\,|\,\alpha\sim\textrm{Unif}\negthinspace\left(\alpha-\zeta,\alpha+\zeta\right)$ is implemented for the concentration parameter, where $\zeta$ (= $2$ in our implementation) is used to control the acceptance rate. For $d$, the mutation rate is considered rather than the acceptance rate, whereby a move is only considered accepted if the proposal differs from the current value.

However, when the DP prior is assumed, or when the sampled value of $d$ is exactly zero under the PYP prior, $\alpha$ is updated according to the auxiliary variable routine of West (1992), with Gibbs updates by simulation from a weighted mixture of two gamma distributions, via

[TABLE]

where $G_{0}$ denotes the current number of non-empty clusters,

[TABLE]

The complementary label switching moves of Papaspiliopoulos & Roberts (2008), which are effective at swapping similar and unequal clusters, respectively, are also incorporated. Firstly, the labels of two randomly chosen non-empty clusters $g$ and $h$ are swapped with probability

[TABLE]

Secondly, the labels of neighbouring active components $l$ and $l+1$ are swapped with probability

[TABLE]

if accepted, $\upsilon_{l}$ and $\upsilon_{l+1}$ are also swapped. Cluster-specific parameters are reordered accordingly after each accepted move. Finally, for updating $\alpha$ under the sparse finite OMIFA or OMFA models, a random walk Metropolis-Hastings step is implemented, with a Gaussian proposal distribution, where

[TABLE]

Further details of this update can be found in Malsiner-Walli et al. (2016) and Frühwirth-Schnatter & Malsiner-Walli (2019).

B Simulation Studies

The performance of the novel IMIFA model with its PYP-MGP priors, in terms of inferring both the number of clusters and the cluster-specific numbers of factors, is assessed here through simulation studies. Section B.1 explores sensitivity to the PYP parameters in a range of dimensionality scenarios, with balanced cluster sizes and a common number of factors. The simulation study in Section B.2 is more challenging; a larger number of clusters (many of which are small) are simulated for $N<p$ data, with different numbers of cluster-specific factors (some of which are large). The final simulation study in Section B.3 mirrors the design in Section B.2, only here the true $\bm{\Lambda}_{g}$ matrices used to generate the data are sparse.

B.1 Simulation Study 1

Firstly, data with $G=3$ clusters and $p=50$ variables are simulated with $q_{g}=4\>\>\forall\>\>g$, and with $\bm{\pi}=\left(\nicefrac{{1}}{{3}},\nicefrac{{1}}{{3}},\nicefrac{{1}}{{3}}\right)$ so that clusters are roughly equally sized. Other model parameters are simulated, with $\bm{\eta}{\mathstrut}_{i}\sim\textrm{N}_{q}\negthinspace\left(\mathbf{0},\bm{\mathcal{I}}_{q}\right)$, $\psi_{jg}\sim\textrm{IG}\negthinspace\left(2,1\right)$, and $\bm{\Lambda}_{jg}\sim\textrm{N}_{q}\negthinspace\left(\mathbf{0},\bm{\mathcal{I}}_{q}\right)$. Notably, loadings are not drawn from the MGP prior (Bhattacharya & Dunson, 2011) underpinning the IMIFA model. To ensure clusters are reasonably closely located, $\bm{\mu}_{g}\sim\textrm{N}_{p}\negthinspace\left(\left(2g-G-1\right)\mathbf{1},\bm{\mathcal{I}}_{p}\right)$. The data are then simulated according to the conditional mixture model

[TABLE]

To evaluate performance in different settings, sample sizes less than, equal to, and greater than $p$ are considered, i.e. $N=25,50$, and $300$. Sensitivity to the PYP and DP parameters is explored by firstly assuming a DP prior with various values of $\alpha$ less than, equal to, and greater than $1$, and by allowing $\alpha$ to be learned as per West (1992), and secondly by incorporating Metropolis-Hastings steps to learn both $\alpha$ and $d$, assuming a PYP prior.

Results, provided in Table B.1, are based on $10$ replicate data sets, standardised prior to model fitting, for each scenario. MCMC chains were run for $25,000$ iterations, with every 2nd sample thinned and the first $20\%$ of iterations discarded as burn-in. Cluster labels were initialised using mclust (Scrucca et al., 2016), as hierarchical clustering gave poor, heavily imbalanced starting values. As the cluster-specific $\bm{\Lambda}_{g}$ and $\bm{\Psi}_{g}$ parameters could still induce separation among clusters, pairwise scatterplots from one randomly chosen raw replicate data set under the $N>p$ scenario are shown in Figure B.1 to demonstrate the extent of the overlap; for visual clarity, only $5$ randomly chosen variables are depicted.