# Finite Mixtures of Multivariate Poisson-Log Normal Factor Analyzers for
Clustering Count Data

Andrea Payne111Co-first authors 222School of Mathematics and Statistics,
Carleton University, Ottawa, ON, Canada. Anjali Silva11footnotemark: 1
333Princess Margaret Cancer Centre - University Health Network, Toronto, ON,
Canada; University of Toronto Libraries, Toronto, ON, Canada.  Steven J.
Rothstein444Department of Molecular and Cellular Biology, University of
Guelph, Guelph, ON, Canada.
Paul D. McNicholas555Department of Mathematics and Statistics, McMaster
University, Hamilton, ON, Canada. Sanjeena Subedi22footnotemark: 2

###### Abstract

A mixture of multivariate Poisson-log normal factor analyzers is introduced by
imposing constraints on the covariance matrix, which resulted in flexible
models for clustering purposes. In particular, a class of eight parsimonious
mixture models based on the mixtures of factor analyzers model are introduced.
Variational Gaussian approximation is used for parameter estimation, and
information criteria are used for model selection. The proposed models are
explored in the context of clustering discrete data arising from RNA
sequencing studies. Using real and simulated data, the models are shown to
give favourable clustering performance. The GitHub R package for this work is
available at https://github.com/anjalisilva/mixMPLNFA and is released under
the open-source MIT license.

##  1 Introduction

Model-based clustering is a technique that utilizes finite mixture models to
cluster data (Wolfe, 1965; Hand, 1989; McLachlan and Peel, 2000; McNicholas,
2016). The general distribution function for mixture models can be given as

$\displaystyle
f(\mathbf{y}|\pi_{1},\ldots,\pi_{G},\boldsymbol{\vartheta}_{1},\ldots,\boldsymbol{\vartheta}_{G})=\sum_{g=1}^{G}\pi_{g}f_{g}(\mathbf{y}|\boldsymbol{\vartheta}_{g}),$

where $G$ is the total number of clusters, $f_{g}(\cdot)$ is the distribution
function with parameters $\boldsymbol{\vartheta}_{g}$, and $\pi_{g}>0$ is the
mixing weight of the $g^{\text{th}}$ component such that
$\sum_{g=1}^{G}\pi_{g}=1$. Mixture model-based clustering methods can be over-
parameterized in high-dimensional spaces, especially as the number of clusters
increases. Subspace clustering allows for clustering data in low-dimensional
subspaces while utilizing information from all the dimensions and by
introducing restrictions to mixture parameters (Bouveyron and Brunet, 2012).
Restrictions are introduced to the model parameters to obtain parsimonious
models, which provide various possible covariance structures while reducing
the number of parameters to fit.

Clustering multivariate count data is of interest in a variety of fields,
including, but not limited to, bioinformatics, text mining, and sports
analytics (Subedi and Browne, 2020). Specifically within bioinformatics, this
form of clustering has become increasingly popular as researchers aim to sort
data into homogeneous subgroups so differing characteristics can be identified
and analyzed (Subedi and Browne, 2020). Many novel clustering methods have
been proposed in recent years; de Souto et al. (2008) compare these methods to
more traditional clustering methods, finding that a mixture model-based
clustering approach outperformed the rest in the context of clustering cancer
tissues for subtype identification. In the years since this publication,
further research has proposed advancements to this approach and new methods
for mixture model clustering for variety of data types (Bouveyron and Brunet,
2012; Gollini and Murphy, 2014; Subedi and McNicholas, 2014; Vrbik and
McNicholas, 2014; Subedi et al., 2015; Browne and McNicholas, 2015; Dang et
al., 2015; Kosmidis and Karlis, 2016; Tortora et al., 2019; Subedi and
McNicholas, 2020; Subedi and Browne, 2020). These proposed methods focus on a
variety of distributions and methods of parameter estimation.

RNA sequencing (RNA-seq) is a specific area of interest within the field of
bioinformatics for applications of model-based clustering. RNA-seq is used to
determine the transcriptional dynamics of a biological system by measuring the
expression levels of many genes (Wang et al., 2009; Roberts et al., 2011).
This technique provides counts of reads that can be mapped back to a
biological entity. Di et al. (2011) note that negative binomial-based models
can account for variability in RNA-seq data while also requiring less
statistical power. As such, negative binomial models have been utilized to
cluster RNA-seq data as well (Si et al., 2014). However, as stated by Silva et
al. (2019), clustering methods based on a univariate distribution such as
negative binomial distribution can struggle with clustering multivariate RNA-
seq data as such models assume independence between variables. Hence, in this
paper, we will develop a method based on multivariate Poisson-log normal
distributions (Aitchison and Ho, 1989; Silva et al., 2019) which can account
for overdispersion and also accommodate positive and negative correlations.

The factor analysis model by Spearman (1904) assumes that a $d$-dimensional
vector of observed variables, $\mathbf{\mathbf{X}}$, can be modelled by a
$K$-dimensional vector of latent factors, where $K<d$. As a result, the factor
analysis model is useful in modelling the covariance structure of high-
dimensional data using a small number of latent variables. The mixture of
factor analyzers model was later introduced by Ghahramani et al. (1996) and
this model can perform clustering and local dimensionality reduction within
each cluster concurrently. Consider $n$ independent $d$-dimensional continuous
variables $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}$, which come from a
heterogeneous population with $G$ subgroups. In the mixture of factor
analyzers framework, $\mathbf{X}_{i}$ is modelled as

$\mathbf{X}_{i}=\mbox{\boldmath$\mu$}_{g}+\mbox{\boldmath$\Lambda$}_{g}\mathbf{U}_{ig}+\mbox{\boldmath$\epsilon$}_{ig}\text{
with probability }\pi_{g},$

for $i=1,\ldots,n$ and $g=1,\ldots,G$. Here, $\mbox{\boldmath$\mu$}_{g}$ is a
$d\times 1$ vector of $g$th component mean, $\mbox{\boldmath$\Lambda$}_{g}$ is
a $d\times K$ matrix of $g$th component factor loadings,
$\mathbf{U}_{ig}\sim\mathcal{N}_{K}(\mathbf{0},\mathbf{I}_{K})$ is a $K\times
1$ vector of $g$th component latent factors, and
$\mbox{\boldmath$\epsilon$}_{ig}\sim\mathcal{N}_{d}(\mathbf{0},\mbox{\boldmath$\Psi$}_{g})$
is a $d\times 1$ vector of $g$th component errors with
$\mbox{\boldmath$\Psi$}_{g}=\text{diag}(\psi_{g1},\ldots,\psi_{gd})$. Note
that the $\mathbf{U}_{ig}$ are independently distributed and are independent
of the $\mbox{\boldmath$\epsilon$}_{ig}$, which are also independently
distributed. Under this model, the density of $\mathbf{X}_{i}$ from a mixture
of factor analyzer models is

$f(\mathbf{x}_{i}\mid\boldsymbol{\Theta})=\sum_{g=1}^{G}\pi_{g}f_{g}(\mathbf{x}_{i}\mid\mbox{\boldmath$\mu$}_{g},\mbox{\boldmath$\Lambda$}_{g}\mbox{\boldmath$\Lambda$}_{g}^{\prime}+\mbox{\boldmath$\Psi$}_{g}),$

where,
$\boldsymbol{\Theta}=(\pi_{1},\ldots,\pi_{G},\mbox{\boldmath$\mu$}_{1},\ldots,\mbox{\boldmath$\mu$}_{G},\mathbf{\Lambda}_{1},\ldots,\mathbf{\Lambda}_{G},\mathbf{\Psi}_{1},\ldots,\mathbf{\Psi}_{G})$.
Further,
$\mathbf{X}_{i}|\mathbf{u}_{ig}\sim\mathcal{N}_{d}(\mbox{\boldmath$\mu$}_{g}+\mbox{\boldmath$\Lambda$}_{g}\mathbf{u}_{ig},\mbox{\boldmath$\Psi$}_{g})$.
It should be noted that the $\mathbf{\Lambda}_{g}$ is not uniquely defined for
$K>1$. Therefore the $K$-dimensional space in which the factors lie can be
determined, but the directions of these factors cannot be determined. However,
this does not affect the clustering algorithm, because
$\mathbf{\Lambda}_{g}\mathbf{\Lambda}_{g}^{\prime}$ is unique. The number of
free $g^{th}$ component covariance parameters that are reduced using the
factor analysis model is

$\frac{1}{2}d(d+1)-\Big{[}dK+d-\frac{1}{2}K(K-1)\Big{]}=\frac{1}{2}\Big{[}(d-K)^{2}-(d+K)\Big{]},$

given that $(d-K)^{2}>(d+K)$ (Lawley and Maxwell, 1962; McNicholas, 2016).

In 2008, this work was extended and a family of eight parsimonious Gaussian
mixture models (PGMMs; McNicholas and Murphy, 2008) were introduced with
parsimonious covariance structures. The PGMM family arises by considering the
general mixture of factor analyzers model
($\mathbf{\Sigma}_{g}=\mathbf{\Lambda}_{g}\mathbf{\Lambda}_{g}^{\prime}+\mathbf{\Psi}_{g}$)
and by allowing the constraints $\mathbf{\Lambda}_{g}=\mathbf{\Lambda}$,
$\mathbf{\Psi}_{g}=\mathbf{\Psi}$, and the isotropic constraint
$\mathbf{\Psi}_{g}=\psi_{g}\mathbf{I}_{d}$. These covariance structures can
have as few as $dK-K(K-1)/2+1$ free parameters or as many as
$G[dK-K(K-1)/2+d]$ free parameters. With the introduction of the factor
analysis structure, the number of covariance parameters is linear in data
dimensionality, thus making this family well suited for analysis of high-
dimensional data (McNicholas and Murphy, 2010). The constraints allow for
assuming a common structure in the component covariance matrix
$\mathbf{\Sigma}_{g}$ between clusters, if appropriate, and this enables a
parsimonious model.

Previously, a model-based clustering methodology using a mixture of
multivariate Poisson-log normal distributions (MPLN; Aitchison and Ho, 1989)
was developed to analyze multivariate count measurements from RNA-seq studies
(Silva et al., 2019). A $d$-dimensional random variable following a
$G$-component mixtures of MPLN distribution is said to have a total of $G-1$
\+ $Gd$ \+ $Gd(d+1)/2$ free parameters. Here, $G-1$ parameters are contributed
by the mixing proportions, $Gd$ from the means and $Gd(d+1)/2$ from the
covariance matrices. Since the largest contribution is through the covariance
matrices, it is a natural focus for the introduction of parsimony.

In this work, a family of mixtures of MPLN factor analyzers that is analogous
to the PGMM family is developed, by considering the constraints
$\mathbf{\Lambda}_{g}=\mathbf{\Lambda}$ and $\mathbf{\Psi}_{g}=\mathbf{\Psi}$.
This family is referred to as the parsimonious mixtures of the MPLN factor
analyzers family (MPLNFA). The proposed model simultaneously performs factor
analysis and cluster analysis, by assuming that the discrete observed data
have been generated by a factor analyzer model in the continuous latent
variables. Details of parameter estimation are provided, and both real and
simulated data illustrations are used to demonstrate the clustering ability.

##  2 Methodology

###  2.1 Mixtures of Factor Analyzer

Suppose $\mathbf{Y}_{ij}$ denotes the observed counts of the $i^{th}$ sample
for the $j^{th}$ variable and $\mathbf{X}_{ij}$ denotes the associated
underlying latent variable. To develop a mixture of factor analyzers, we
assume the following hierarchical structure with the $K$-dimensional latent
factor $\mathbf{U}_{i}$:

$\displaystyle Y_{ij}$ $\displaystyle\mid
X_{ij}=x_{ij}\sim\text{Poisson}\left(\exp\\{x_{ij}+\log{C_{i}}\\}\right),$
$\displaystyle\mathbf{X}_{i}$
$\displaystyle\mid\mathbf{U}_{i}=\mathbf{u}_{i}\sim
N(\boldsymbol{\mu}_{g}+\boldsymbol{\Lambda}_{g}\mathbf{u}_{i},\boldsymbol{\Psi}_{g}),$
$\displaystyle\mathbf{U}_{i}$ $\displaystyle\sim
N(\mathbf{0}_{K},\mathbf{I}_{K}).$

Here, $C_{i}$ is a fixed known constant representing the normalized library
sizes added to account for the differences in the library sizes of the
$i^{th}$ sample.

In model-based clustering, we also have unobserved component membership
indicator variable $\mathbf{Z}$ such that $Z_{ig}=1$ if the observation
$i^{th}$ belongs to group $g$ and $Z_{ig}=0$ otherwise. Hence, the complete
data now comprises observed expression levels $\mathbf{y}$, underlying latent
variable $\mathbf{X}$, and unknown group membership $\mathbf{Z}$. The
complete-data likelihood can therefore be written as:

$\displaystyle L(\boldsymbol{\vartheta})$
$\displaystyle=\prod_{g=1}^{G}\prod_{i=1}^{n}\left[\pi_{g}\left\\{\prod_{j=1}^{d}f(y_{ij}\mid
x_{ij},C_{i})\right\\}\leavevmode\nobreak\
f(\mathbf{x}_{i}\mid\mathbf{u}_{i},\boldsymbol{\mu}_{g},\boldsymbol{\Lambda}_{g},\boldsymbol{\Psi}_{g})\leavevmode\nobreak\
f(\mathbf{u}_{i})\right]^{z_{ig}},$

where $\boldsymbol{\vartheta}$ denotes all the model parameters. The complete
data log-likelihood can therefore be written as:

$\displaystyle l_{c}$ $\displaystyle(\boldsymbol{\vartheta})=$
$\displaystyle\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\left[\log\pi_{g}+\left\\{\sum_{j=1}^{d}\log
f(y_{ij}\mid x_{ij},k_{j})\right\\}+\log
f(\mathbf{x}_{i}\mid\mathbf{u}_{i},\boldsymbol{\mu}_{g},\boldsymbol{\Lambda}_{g},\boldsymbol{\Psi}_{g})+\log
f(\mathbf{u}_{i})\right],$

To utilize the expectation-maximization (EM) framework for parameter
estimation, we require
$\operatorname{\mathbb{E}}(Z_{ig}\mathbf{X}_{i}\mid\mathbf{y}_{i},\mathbf{u}_{i},\boldsymbol{\mu}_{g},\boldsymbol{\Lambda}_{g},\boldsymbol{\Psi}_{g})$,
$\operatorname{\mathbb{E}}(Z_{ig}\mathbf{X}_{i}\mathbf{X}_{i}^{\prime}\mid\mathbf{y}_{i},\mathbf{u}_{i},\boldsymbol{\mu}_{g},\boldsymbol{\Lambda}_{g},\boldsymbol{\Psi}_{g})$,
$\operatorname{\mathbb{E}}(Z_{ig}\mathbf{U}_{i}\mid\mathbf{y}_{i},\mathbf{x}_{i},\boldsymbol{\mu}_{g},\boldsymbol{\Lambda}_{g},\boldsymbol{\Psi}_{g})$
and
$\operatorname{\mathbb{E}}(Z_{ig}\mathbf{U}_{i}\mathbf{U}_{i}^{\prime}\mid\mathbf{y}_{i},\mathbf{x}_{i},\boldsymbol{\mu}_{g},\boldsymbol{\Lambda}_{g},\boldsymbol{\Psi}_{g})$.
To compute these quantities, we need the posterior distribution of
$\mathbf{X}_{i}\mid\mathbf{y}_{i},\mathbf{u}_{i}$ which does not have a closed
form. Thus, parameter estimation relies on the Markov chain Monte Carlo
expectation-maximization (MCMC-EM) approach which can be computationally
intensive.

###  2.2 Variational approximation of mixtures of MPLN factor analyzers

Subedi and Browne (2020) proposed a framework for parameter estimation
utilizing variational Gaussian approximation (VGA) for MPLN-based mixture
models. Variational approximation (Wainwright et al., 2008) is an approximate
inference technique that uses a computationally convenient approximating
density in place of a more complex but ‘true’ posterior density obtained by
minimizing the Kullback-Leibler (KL) divergence between the true and the
approximating densities. Markov chain Monte Carlo expectation-maximization
(MCMC-EM) has also been used for parameter estimation of MPLN-based mixture
models, but VGA was shown to be computationally efficient (Silva et al.,
2023). Here, we propose a mixture of Poisson log-normal distributions with
factor analyzers that utilize a variational EM framework for parameter
estimation.

Parameter estimation for a mixture of factor analyzers is typically done using
an alternating expectation conditional maximization (AECM) framework that
assumes different specifications of missing data at different stages. Here, we
will adopt a similar notion resulting in a two-stage iterative EM-type
approach.

###  2.3 Stage 1

In Stage 1, we treat $\mathbf{X}$ and $\mathbf{Z}$ as missing such that

$\displaystyle Y_{ij}$ $\displaystyle\mid
X_{ij}=x_{ij}\sim\text{Poisson}\left(\exp\\{x_{ij}+\log{C_{i}}\\}\right)$
$\displaystyle\mathbf{X}_{i}$ $\displaystyle\sim
N(\boldsymbol{\mu}_{g},\boldsymbol{\Sigma}_{g}),$

where
$\boldsymbol{\Sigma}_{g}=\boldsymbol{\Lambda}_{g}\boldsymbol{\Lambda}_{g}^{T}+\boldsymbol{\Psi}_{g}$.Therefore,
the component-specific marginal density of the observed data $\mathbf{y}_{i}$
can be written as

$\displaystyle f(\mathbf{y}_{i}\mid
Z_{ig}=1,\boldsymbol{\mu}_{g},\boldsymbol{\Sigma}_{g})=\int_{{}^{d}}\left[\prod_{j=1}^{d}f_{p}(y_{ij}\mid
x_{ij},C_{i},Z_{ig}=1)\right]\leavevmode\nobreak\ f_{N}(\mathbf{x}_{i}\mid
Z_{ig}=1,\boldsymbol{\mu}_{g},\boldsymbol{\Sigma}_{g})\leavevmode\nobreak\
d\mathbf{x}_{i},$

where $x_{ij}$ and $y_{ij}$ are the $j^{th}$ element of $\mathbf{x}_{i}$ and
$\mathbf{y}_{i}$ respectively, $f_{p}(\cdot)$ is the probability mass function
of the Poisson distribution with mean $e^{x_{ij}+\log C_{j}}$ and
$f_{N}(\cdot)$ is the probability density function of $d$-dimensional Gaussian
distribution with mean $\boldsymbol{\mu}_{g}$ and covariance
$\boldsymbol{\Sigma}_{g}$. Note that the marginal distribution of $\mathbf{Y}$
involves multiple integrals and cannot be further simplified. The complete-
data log-likelihood using the marginal density of $\mathbf{y}_{i}\mid
Z_{ig}=1$ can be written as:

$l(\boldsymbol{\vartheta})=\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\log\pi_{g}+\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\log
f(\mathbf{y}_{i}\mid Z_{ig}=1,\boldsymbol{\mu}_{g},\boldsymbol{\Sigma}_{g}).$

As the marginal of $\mathbf{Y}_{i}\mid Z_{ig}=1$ cannot be simplified, we use
variational Gaussian approximation to approximate the marginal of
$\mathbf{Y}_{i}$ from component $g$ in the mixtures of MPLN distributions.
Suppose, we have an approximating density $q(\mathbf{x}_{ig})$, we can write

$\displaystyle\log f(\mathbf{y}_{i}\mid Z_{ig}=1)$
$\displaystyle=F(q_{ig},\mathbf{y}_{i})+D_{KL}(q_{ig}\|f_{ig}),$

where
$D_{KL}(q_{ig}\|f_{ig})=\int_{{}^{d}}q(\mathbf{x}_{ig})\log\frac{q(\mathbf{x}_{ig})}{f(\mathbf{x}_{i}\mid\mathbf{y}_{i},Z_{ig}=1)}d\mathbf{x}_{ig}$
is the Kullback-Leibler (KL) divergence between
$f(\mathbf{x}_{i}\mid\mathbf{y}_{i},Z_{ig}=1)$ and approximating distribution
$q(\mathbf{x}_{ig})$, and

$F(q_{ig},\mathbf{y}_{i})=\int_{{}^{d}}\left[\log
f(\mathbf{y}_{i},\mathbf{x}_{i}\mid Z_{ig}=1)-\log
q(\mathbf{x}_{ig})\right]q(\mathbf{x}_{ig})d\mathbf{x}_{ig},$

is our evidence lower bound (ELBO) for each observation $\mathbf{y}_{i}$.
Details can be found in Appendix A.1. The complete data log-likelihood of the
mixtures of MPLN distributions can be written as:

$\displaystyle l_{c}(\boldsymbol{\vartheta}\mid\mathbf{y})$
$\displaystyle=\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\log\pi_{g}+\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\left[F(q_{ig},\mathbf{y}_{i})+D_{KL}(q_{ig}\|f_{ig})\right].$

In VGA, $q(\mathbf{x}_{ig})$ is assumed to be a Gaussian distribution.
Assuming $q(\mathbf{x}_{ig})=\mathscr{N}(\mathbf{m}_{ig},\mathbf{S}_{ig})$,
the ELBO for each observation $\mathbf{y}_{i}$ becomes

$\displaystyle F(q_{ig},\mathbf{y}_{i})$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{g})^{\prime}\boldsymbol{\Sigma}_{g}^{-1}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{g})-\frac{1}{2}\,\mbox{tr}(\boldsymbol{\Sigma}_{g}^{-1}\mathbf{S}_{ig})+\frac{1}{2}\log|\mathbf{S}_{ig}|-\frac{1}{2}\log|\boldsymbol{\Sigma}_{g}|+\frac{d}{2}$
$\displaystyle+\mathbf{m}_{ig}^{\prime}\mathbf{y}_{i}+\sum_{j=1}^{d}\left(\log
C_{i}\right)y_{ij}-\sum_{j=1}^{d}\left\\{e^{\log
C_{i}+m_{igj}+\frac{1}{2}S_{ig,jj}}+\log(y_{ij}!)\right\\}.$

Thus, to minimize the KL divergence between the true and approximating
density, we find the values of the variational parameters $\mathbf{m}_{ig}$
and $\mathbf{S}_{ig}$ that maximize our ELBO. Similar to Arridge et al. (2018)
and Subedi and Browne (2020), estimation of the mean ($\mathbf{m}_{ig}$) and
variance ($\mathbf{S}_{ig}$) can be obtained via Newton’s method and fixed-
point method.

Therefore, in Stage 1, we update
$\operatorname{\mathbb{E}}(Z_{ig}\mid\mathbf{y}_{i})$, variational parameters
$\mathbf{m}_{ig}$ and $\mathbf{S}_{ig}$, and model parameters $\pi_{g}$ and
$\boldsymbol{\mu}_{g}$.

  1. 1.

Conditional on the variational parameters $\mathbf{m}_{ig},\mathbf{S}_{ig}$
and on $\boldsymbol{\mu}_{g}$ and $\boldsymbol{\Sigma}_{g}$, the
$\operatorname{\mathbb{E}}(Z_{ig})$ is computed. Given $\boldsymbol{\mu}_{g}$
and $\boldsymbol{\Sigma}_{g}$,

$\displaystyle\operatorname{\mathbb{E}}(Z_{ig}\mid\mathbf{y}_{i})=\frac{\pi_{g}f(\mathbf{y}\mid\boldsymbol{\mu}_{g},\boldsymbol{\Sigma}_{g})}{\sum_{h=1}^{G}\pi_{h}f(\mathbf{y}\mid\boldsymbol{\mu}_{h},\boldsymbol{\Sigma}_{h})}.$

Note that this involves the marginal distribution of $\mathbf{Y}$ which is
difficult to compute. Hence, similar to Subedi and Browne (2020), we use an
approximation of $\operatorname{\mathbb{E}}(Z_{ig})$ where we replace the
marginal density of the exponent of ELBO such that

$\displaystyle\widehat{Z}_{ig}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{\pi_{g}\exp\left[F\left(q_{ig},\mathbf{y}_{i}\right)\right]}{\sum_{h=1}^{G}\pi_{h}\exp\left[F\left(q_{ih},\mathbf{y}_{i}\right)\right]}.$

  2. 2.

Update the variational parameters $\mathbf{m}_{ig}$ and $\mathbf{S}_{ig}$ as
following:

    1. (a)

Fixed-point method for updating $\mathbf{S}_{ig}$ is

$\displaystyle\mathbf{S}_{ig}^{(t+1)}=\left\\{\boldsymbol{\Sigma}_{g}^{-1}+\mathbf{I}\odot\mathbf{exp}\left[\log
C_{i}+\mathbf{m}_{ig}^{(t)}+\frac{1}{2}\,\mbox{diag}\left(\mathbf{S}_{ig}^{(t)}\right)\right]\mathbf{1}_{d}^{\prime}\right\\}^{-1}$

where the vector function
$\mathbf{exp}\left[\mathbf{a}\right]=(e^{a_{1}},\ldots,e^{a_{d}})^{\prime}$ is
a vector of exponential each element of the $d$-dimensional vector
$\mathbf{a}$,
$\mbox{diag}(\mathbf{S})=(\mathbf{S}_{11}\ldots,\mathbf{S}_{dd})$ puts the
diagonal elements of the $d\times d$ matrix $\mathbf{S}$ into a d-dimensional
vector, $\odot$ the Hadmard product and $\mathbf{1}_{d}$ is a d-dimensional
vector of ones ;

    2. (b)

Newton’s method to update $\mathbf{m}_{ig}$ is

$\displaystyle\mathbf{m}_{ig}^{(t+1)}$
$\displaystyle=\mathbf{m}_{ig}^{(t)}\\!-\\!\mathbf{S}_{ig}^{(t+1)}\\!\left\\{\\!\mathbf{y}_{i}-\mathbf{exp}\left[\log
C_{i}+\mathbf{m}_{ig}^{(t)}\\!+\\!\frac{1}{2}\,\mbox{diag}\left(\mathbf{S}_{ig}^{(t+1)}\right)\right]\\!\right.$
$\displaystyle\left.-\\!\boldsymbol{\Sigma}_{g}^{-1}\left(\mathbf{m}_{ig}^{(t)}-\boldsymbol{\mu}_{g}\right)\\!\right\\}\\!.$

  3. 3.

Given $\widehat{Z}_{ig}$ and the variational parameters $\mathbf{m}_{ig}$ and
$\mathbf{S}_{ig}$, the parameters $\boldsymbol{\pi}$ and
$\boldsymbol{\mu}_{g}$ are updated as:

$\displaystyle\widehat{\pi}_{g}$
$\displaystyle=\frac{\sum_{i=1}^{n}\widehat{Z}_{ig}}{n},$
$\displaystyle\widehat{\boldsymbol{\mu}}_{g}$
$\displaystyle=\frac{\sum_{i=1}^{n}\widehat{Z}_{ig}\mathbf{m}_{ig}^{(t+1)}}{\sum_{i=1}^{n}\widehat{Z}_{ig}}.$

###  2.4 Stage 2

In Stage 2, we treat $\mathbf{X}$, $\mathbf{U}$ and $\mathbf{Z}$ as missing
such that

$\displaystyle Y_{ij}$ $\displaystyle\mid
X_{ij}=x_{ij}\sim\text{Poisson}\left(\exp\\{x_{ij}+\log{C_{i}}\\}\right)$
$\displaystyle\mathbf{X}_{i}$
$\displaystyle\mid\mathbf{U}_{i}=\mathbf{u}_{i}\sim
N(\boldsymbol{\mu}_{g}+\boldsymbol{\Lambda}_{g}\mathbf{u}_{i},\boldsymbol{\Psi}_{g}),$
$\displaystyle\mathbf{U}_{i}$ $\displaystyle\sim
N(\mathbf{0}_{K},\mathbf{I}_{K}).$

Here, we use $d$ for the dimensionality of the observed data and $K$ as the
dimensionality of the latent variable. Suppose, we have an approximating
density $q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})$, the
marginal of the log of the probability mass function of $\mathbf{Y}_{i}\mid
Z_{ig}=1$ can be written as:

$\displaystyle\log f(\mathbf{y}_{i}\mid Z_{ig}=1)$
$\displaystyle=F(q_{ig},\mathbf{y}_{i})+D_{KL}(q_{ig}\|f_{ig}),$

where
$D_{KL}(q_{ig}\|f_{ig})=\int_{{}^{p}}\int_{{}^{d}}q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})\log\frac{q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})}{f(\mathbf{x}_{i},\mathbf{u}_{i}\mid\mathbf{y}_{i},Z_{ig}=1)}\leavevmode\nobreak\
d\mathbf{x}_{ig}\leavevmode\nobreak\ d\mathbf{u}_{ig}$ is the Kullback-Leibler
(KL) divergence between
$f(\mathbf{x}_{i},\mathbf{u}_{i}\mid\mathbf{y}_{i},Z_{ig}=1)$ and
approximating distribution $q(\mathbf{x}_{ig},\mathbf{u}_{ig})$, and

$F(q_{ig},\mathbf{y}_{i})=\int_{{}^{p}}\int_{{}^{d}}\Big{[}\log
f(\mathbf{y}_{i},\mathbf{x}_{i},\mathbf{u}_{i}\mid Z_{ig}=1)-\log
q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})\Big{]}q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})\leavevmode\nobreak\
d\mathbf{x}_{ig}\leavevmode\nobreak\ d\mathbf{u}_{ig},$

is our evidence lower bound (ELBO) for each observation $\mathbf{y}_{i}$.

If we assume $q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})$ to
be factorable such that
$q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})=q_{\mathbf{x}}(\mathbf{x}_{ig})q_{\mathbf{u}}(\mathbf{u}_{ig})$,

$\displaystyle F$ $\displaystyle(q_{ig},\mathbf{y}_{i})=$
$\displaystyle\int_{{}^{p}}\Bigg{[}\left\\{\int_{{}^{d}}\Big{[}\log
f(\mathbf{y}_{i}\mid\mathbf{x}_{i},\mathbf{u}_{i},Z_{ig}=1)+\log
f(\mathbf{x}_{i}\mid\mathbf{u}_{i},Z_{ig}=1)-\log
q_{\mathbf{x}}(\mathbf{x}_{i})\Big{]}\leavevmode\nobreak\
q_{\mathbf{x}}(\mathbf{x}_{ig})\leavevmode\nobreak\ d\mathbf{x}_{ig}\right\\}$
$\displaystyle+\log f(\mathbf{u}_{i}\mid Z_{ig}=1)-\log
q_{\mathbf{u}}(\mathbf{u}_{ig})\Bigg{]}q_{\mathbf{u}}(\mathbf{u}_{ig})\leavevmode\nobreak\
d\mathbf{u}_{ig}.$

Using the approximating density $q_{\mathbf{x}}(\mathbf{x}_{ig})$ from Stage 1
(i.e. $q_{\mathbf{x}}(\mathbf{x}_{ig})=N(\mathbf{m}_{ig},\mathbf{S}_{ig})$)
and assuming
$q_{\mathbf{u}}(\mathbf{u}_{ig})=N(\mathbf{P}_{ig},\mathbf{Q}_{ig})$, we get

$\displaystyle F(q_{ig},\mathbf{y}_{i})$
$\displaystyle=\mathbf{m}_{ig}^{\prime}\mathbf{y}_{i}+\sum_{j=1}^{d}\left(\log
C_{i}\right)y_{ij}-\sum_{j=1}^{d}\left\\{e^{\log
C_{i}+m_{igj}+\frac{1}{2}S_{ig,jj}}+\log(y_{ij}!)\right\\}+\frac{1}{2}\log|\mathbf{S}_{ig}|+\frac{d}{2}$
$\displaystyle-\frac{1}{2}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})^{T}\boldsymbol{\Psi}_{g}^{-1}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})+(\mathbf{m}_{ig}-\boldsymbol{\mu}_{g})^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}-\frac{1}{2}\mathbf{P}_{ig}^{T}\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}$
$\displaystyle-\frac{1}{2}\,\mbox{tr}(\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{Q}_{ig})-\frac{1}{2}\,\mbox{tr}\left(\boldsymbol{\Psi}_{g}^{-1}\mathbf{S}_{ig}\right)-\frac{1}{2}\log|\boldsymbol{\Psi}_{g}|-\frac{1}{2}\mathbf{P}_{ig}^{T}\mathbf{P}_{ig}+\frac{1}{2}\log|\mathbf{Q}_{ig}^{-1}|$
$\displaystyle-\frac{1}{2}\,\mbox{tr}(\mathbf{Q}_{ig})+\frac{K}{2}.$

The variational parameters $\mathbf{P}_{ig}$ and $\mathbf{Q}_{ig}$ that
maximize this ELBO will also minimize the KL divergence. Therefore, in Stage
2, we update variational parameters $\mathbf{P}_{ig}$ and $\mathbf{Q}_{ig}$,
and model parameters $\boldsymbol{\Lambda}_{g}$ and $\boldsymbol{\Psi}_{g}$.

  1. 1.

The update for the variational parameters $\mathbf{P}_{ig}$ and
$\mathbf{Q}_{ig}$ are:

$\displaystyle\widehat{\mathbf{Q}}_{ig}^{(t+1)}$
$\displaystyle=(\mathbf{I}_{p}+\hat{\boldsymbol{\Lambda}}_{g}^{T}\hat{\boldsymbol{\Psi}}_{g}^{-1}\hat{\boldsymbol{\Lambda}}_{g})^{-1}$
$\displaystyle=\mathbf{I}_{p}-\hat{\boldsymbol{\beta}}_{g}\hat{\boldsymbol{\Lambda}}_{g}(\text{after
a bit of simplification where
}\hat{\boldsymbol{\beta}}_{g}=\hat{\boldsymbol{\Lambda}}_{g}^{T}(\hat{\boldsymbol{\Psi}}_{g}+\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\Lambda}}_{g}^{T})^{-1}).$
$\displaystyle\widehat{\mathbf{P}}_{ig}^{(t+1)}$
$\displaystyle=\widehat{\mathbf{Q}}_{ig}^{(t+1)}\hat{\boldsymbol{\Lambda}}_{g}^{T}\hat{\boldsymbol{\Psi}}^{-1}_{g}(\mathbf{m}_{ig}^{(t+1)}-\hat{\boldsymbol{\mu}}_{g})$
$\displaystyle=\hat{\boldsymbol{\beta}}_{g}(\mathbf{m}_{ig}^{(t+1)}-\hat{\boldsymbol{\mu}}_{g}).$

  2. 2.

The update for $\boldsymbol{\Psi}_{g}$ and $\boldsymbol{\Lambda}_{g}$ are
obtained by repeating the following steps until convergence:

$\displaystyle\widehat{\boldsymbol{\Theta}}_{g}$
$\displaystyle=\mathbf{I}_{p}-\hat{\boldsymbol{\beta}}_{g}\hat{\boldsymbol{\Lambda}}_{g}+\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}_{g}^{T},$
$\displaystyle\hat{\boldsymbol{\beta}}_{g}$
$\displaystyle=\hat{\boldsymbol{\Lambda}}_{g}^{T}(\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\Lambda}}_{g}^{T}+\hat{\boldsymbol{\Psi}}_{g})^{-1},$
$\displaystyle\hat{\boldsymbol{\Lambda}}_{g}$
$\displaystyle=\mathbf{W}_{g}\hat{\boldsymbol{\beta}}^{T}\widehat{\boldsymbol{\Theta}}_{g}^{-1},$
(1) $\displaystyle\hat{\boldsymbol{\Psi}}_{g}$
$\displaystyle=\,\mbox{diag}\left(\mathbf{W}_{g}-\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{\sum_{i=1}^{n}\widehat{Z}_{ig}\mathbf{S}_{ig}}{\sum_{i=1}^{n}\widehat{Z}_{ig}}\right),$
(2)

where
$\mathbf{W}_{g}=\frac{\sum_{i=1}^{n}\widehat{Z}_{ig}(\mathbf{m}_{ig}^{(t+1)}-\hat{\boldsymbol{\mu}}_{g})(\mathbf{m}_{ig}^{(t+1)}-\hat{\boldsymbol{\mu}}_{g})^{T}}{\sum_{i=1}^{n}\widehat{Z}_{ig}}$Note:
convergence is assumed when the difference in Frobenius norm of
$\hat{\boldsymbol{\Lambda}}_{g}$ and $\hat{\boldsymbol{\Psi}}_{g}$ from
successive iterations are both less than $10^{-6}$.

##  3 Family of models

To obtain a family of parsimonious models, constraints can be imposed on
$\boldsymbol{\Lambda}_{g}$ and $\boldsymbol{\Psi}_{g}$ as described in Table 1
that results in a family of 8 models.

Table 1: Family of models based on the constraints on
$\boldsymbol{\Lambda}_{g}$ and $\boldsymbol{\Psi}_{g}$.

Model ID |  Loading Matrix $\boldsymbol{\Lambda}_{g}$ |  Error Variance $\boldsymbol{\Psi}_{g}$ |  Isotropic $\psi_{g}\mathbf{I}$ | Covariance parameters
---|---|---|---|---
UUU | unconstrained | unconstrained | unconstrained | $G\left[dK-K\left(K-1\right)/2\right]+Gd$
UUC | unconstrained | unconstrained | constrained | $G\left[dK-K\left(K-1\right)/2\right]+G$
UCU | unconstrained | constrained | unconstrained | $G\left[dK-K\left(K-1\right)/2\right]+d$
UCC | unconstrained | constrained | constrained | $G\left[dK-K\left(K-1\right)/2\right]+1$
CUU | constrained | unconstrained | unconstrained | $\left[dK-K\left(K-1\right)/2\right]+Gd$
CUC | constrained | unconstrained | constrained | $\left[dK-K\left(K-1\right)/2\right]+G$
CCU | constrained | constrained | unconstrained | $\left[dK-K\left(K-1\right)/2\right]+d$
CCC | constrained | constrained | constrained | $\left[dK-K\left(K-1\right)/2\right]+1$

The estimation of $\boldsymbol{\Psi}_{g}$ and $\boldsymbol{\Lambda}_{g}$ under
various constraints are provided below:

  * •

Under the isotropic constraint (i.e.
$\boldsymbol{\Psi}_{g}=\psi_{g}\mathbf{I}_{d}$),

$\hat{\psi}_{g}=\frac{1}{d}\,\mbox{tr}\left(\mathbf{W}_{g}-\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}}{n_{g}}\right).$
(3)

  * •

Under the equal variance constraint (i.e.,
$\boldsymbol{\Psi}_{g}=\boldsymbol{\Psi}$),

$\displaystyle\hat{\boldsymbol{\Psi}}$
$\displaystyle=\,\mbox{diag}\left(\sum_{g=1}^{G}\frac{n_{g}}{n}\mathbf{W}_{g}-\sum_{g=1}^{G}\frac{n_{g}}{n}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right).$

  * •

Under the equal Variance and isotropic constraint (i.e.,
$\boldsymbol{\Psi}_{g}=\psi\mathbf{I}_{d}$),

$\displaystyle\hat{\psi}$
$\displaystyle=\frac{1}{d}\,\mbox{tr}\left(\sum_{g=1}^{G}\frac{n_{g}}{n}\mathbf{W}_{g}-\sum_{g=1}^{G}\frac{n_{g}}{n}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right).$
(4)

  * •

Under the equal loading matrices constraint (i.e.
$\boldsymbol{\Lambda}_{g}=\boldsymbol{\Lambda}$), the loading matrix cannot be
solved directly and must be solved in a row-by-row manner as suggested by
McNicholas and Murphy (2008). Hence,

$\displaystyle\hat{\lambda}_{i}$
$\displaystyle=\mathbf{r}_{i}\left(\sum_{g=1}^{G}\frac{n_{g}}{\hat{\Psi}_{g,(i)}}\hat{\boldsymbol{\Theta}}_{g}\right)^{-1},$
(5)

where $\lambda_{i}$ is the $i^{th}$ row of the matrix $\boldsymbol{\Lambda}$,
$\Psi_{g,(i)}$ is the $i^{th}$ diagonal element of $\boldsymbol{\Psi}_{g}$ and
$\mathbf{r}_{i}$ is the $i^{th}$ row of the matrix
$\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Psi}}_{g}^{-1}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}_{g}^{T}$.

Full details of the parameter updates for all models can be found in Appendix
A.3.

##  4 Results

###  4.1 Simulation Study

Simulation studies were conducted to illustrate the ability to recover the
true underlying parameters by the MPLNFA. Three simulation settings were
utilized, each from a different model in the family of models. Within the
three settings, three different numbers of clusters and latent factors were
utilized and 100 datasets were generated with 1,000 observations per dataset
under each setting. The exact specifications of the simulation settings are
summarized in Table 2.

Table 2: Summary of the simulation settings.

Setting | Number of | Number of $n$ | Dimensions, $d$ | True Model
---|---|---|---|---
| Simulated Datasets | Observations, $n$ |  |
1 | 100 | 1000 | 8 | UCC
2 | 100 | 1000 | 10 | CCC
3 | 100 | 1000 | 10 | UUU

Setting | Number of | Mixing Proportion, $\pi$ | Number of | True Model
---|---|---|---|---
| Clusters, $G$ |  | Latent factors, $K$ |
1 | 4 | (0.11, 0.43, 0.24, 0.22) | 2 | UCC
2 | 2 | (0.32, 0.68) | 3 | CCC
3 | 3 | (0.23, 0.44, 0.33) | 4 | UUU

Table 3: Comparison of true and recovered parameters for all three simulation
settings.

Simulation Setting 1 - UCC
---
Group 1 | $\boldsymbol{\mu}$ | (6, 3, 3, 6, 3, 6, 3, 3)
$\boldsymbol{\hat{\mu}}$ | (5.99, 3.00, 3.00, 6.00, 3.00, 6.00, 3.00, 3.00)
Standard Deviation | (0.07, 0.10, 0.08, 0.13, 0.08, 0.07, 0.10, 0.08)
MSE($\boldsymbol{\Sigma}$) | 0.01
Group 2 | $\boldsymbol{\mu}$ | (1, 3, 5, 1, 3, 5, 3, 5)
$\boldsymbol{\hat{\mu}}$ | (1.44, 2.85, 5.14, 1.42, 2.86, 5.15, 3.14, 4.86)
Standard Deviation | (1.05, 0.36, 0.36, 1.04, 0.35, 0.35, 0.36, 0.35)
MSE($\boldsymbol{\Sigma}$) | 0.07
Group 3 | $\boldsymbol{\mu}$ | (4, 2, 6, 4, 2, 6, 4, 4)
$\boldsymbol{\hat{\mu}}$ | (3.58, 2.16, 5.86, 3.59, 2.14, 5.86, 3.87, 4.13)
Standard Deviation | (1.03, 0.34, 0.36, 1.05, 0.36, 0.36, 0.36, 0.35)
MSE($\boldsymbol{\Sigma}$) | 0.07
Group 4 | $\boldsymbol{\mu}$ | (5, 3, 5, 3, 5, 3, 3, 5)
$\boldsymbol{\hat{\mu}}$ | (5.00, 2.99, 4.99, 3.00, 5.00, 3.01, 3.00, 5.00)
Standard Deviation | (0.04, 0.05, 0.05, 0.05, 0.03, 0.04, 0.03, 0.06)
MSE($\boldsymbol{\Sigma}$) | 0.00
Simulation Setting 2 - CCC
Group 1 | $\boldsymbol{\mu}$ | (6, 3, 3, 6, 3, 6, 3, 3, 6, 3)
$\boldsymbol{\hat{\mu}}$ | (6.00, 3.00, 3.00, 6.00, 3.00, 6.01, 3.00, 3.00, 6.00, 3.00)
Standard Deviation | (0.05, 0.07, 0.05, 0.07 0.05, 0.05, 0.07, 0.06, 0.05, 0.06)
MSE($\boldsymbol{\Sigma}$) | 0.00
Group 2 | $\boldsymbol{\mu}$ | (5, 3, 5, 3, 5, 5, 3, 5, 3, 5)
$\boldsymbol{\hat{\mu}}$ | (4.99, 3.01, 5.00, 3.00, 5.00, 5.00, 3.01, 5.00, 3.00, 5.00)
Standard Deviation | (0.03, 0.05, 0.04, 0.04, 0.04, 0.04, 0.04, 0.05, 0.03, 0.05)
MSE($\boldsymbol{\Sigma}$) | 0.00
Simulation Setting 3 - UUU
Group 1 | $\boldsymbol{\mu}$ | (4, 6, 4, 2, 2, 4, 6, 4, 6, 2)
$\boldsymbol{\hat{\mu}}$ | (4.00, 6.00, 4.00, 1.99, 2.00, 4.00, 5.98, 4.00, 6.00, 2.00)
Standard Deviation | (0.05, 0.10, 0.06, 0.12, 0.07, 0.07, 0.09, 0.08, 0.09, 0.08)
MSE($\boldsymbol{\Sigma}$) | 0.01
Group 2 | $\boldsymbol{\mu}$ | (5, 5, 3, 3, 7, 5, 3, 3, 7, 7)
$\boldsymbol{\hat{\mu}}$ | (5.00, 5.00, 2.99, 3.00, 7.00, 5.00, 3.00, 3.00, 7.00, 7.00)
Standard Deviation | (0.04, 0.07, 0.07, 0.07, 0.04, 0.03, 0.08, 0.06, 0.06, 0.07)
MSE($\boldsymbol{\Sigma}$) | 0.02
Group 3 | $\boldsymbol{\mu}$ | (2, 4, 4, 7, 2, 4, 7, 2, 7, 4)
$\boldsymbol{\hat{\mu}}$ | (1.99, 3.99, 3.99, 6.99, 1.98, 4.00, 7.01, 1.97, 6.99, 4.00)
Standard Deviation | (0.09, 0.06, 0.08, 0.08, 0.09, 0.07, 0.09, 0.07, 0.10, 0.07)
MSE($\boldsymbol{\Sigma}$) | 0.01

For the simulation data analyses, the normalization factors representing
library size estimate for samples were set to 1. All simulation analyses were
performed with R version 4.2.2 Patched “Innocent and Trusting”, released on
November 11, 2022, (R Core Team, 2022).

We fitted all the proposed eight models from the MPLNFA family to all datasets
for $G=1\text{ to }G_{o}+1$, where $G_{o}$ is the actual number of clusters
and the number of latent factors ranged $K=1\text{ to }K_{o}+1$, where $K_{o}$
is the actual number of latent factors. All models were initialized with
$k$-means initialization with three random starts.

The simulations demonstrated that our proposed method can recover the correct
parameters for MPLNFA (see Table 3). The clustering results obtained using the
BIC for model selection and the corresponding average ARI values are
summarized in Table 4. As seen in Table 4, BIC selected the correct number of
clusters, latent factors, and component scale matrix for the first simulation
setting, with $M=$ UCC, with an average ARI value of 0.99. BIC provides
similar results for settings 2 and 3 ($M=$ CCC and UUU, respectively), with an
average ARI value of 1.00.

Table 4: Model selection results of the clusters (average ARI, standard
deviation), latent factors and component scale matrices for each simulation
setting compared to true values.

|  | Cluster, $G$ |  |  | Latent factor, $K$ |  |  | Component scale matrix, $M$ |  |
---|---|---|---|---|---|---|---|---|---|---
Setting |  | Original | MPLNFA | mPoisson | mNB |  |  | Original | MPLNFA | mPoisson | mNB |  |  | Original | MPLNFA | mPoisson | mNB |  |
1 |  | 4 | 4 | 5 | 2 |  |  | 2 | 2 | NA | NA |  |  | UCC | UCC | NA | NA |  |
|  |  | (0.99, 0.01) | (0.51, 0.08) | (0.36, 0.19) |  |  |  |  |  |  |  |  |  |  |  |  |  |
2 |  | 2 | 2 | 3 | 2 |  |  | 3 | 3 | NA | NA |  |  | CCC | CCC | NA | NA |  |
|  |  | (1.00, 0.00) | (1.00, 0.00) | (1.00, 0.00) |  |  |  |  |  |  |  |  |  |  |  |  |  |
3 |  | 3 | 3 | 4 | 3 |  |  | 4 | 4 | NA | NA |  |  | UUU | UUU | NA | NA |  |
|  |  | (1.00, 0.00) | (0.69, 0.02) | (0.52, 0.11) |  |  |  |  |  |  |  |  |  |  |  |  |  |

Additionally, in Table 4, we clustered the same data with two other common
approaches, also using BIC for model selection. The first, labelled ‘mPoisson’
in Table 4, is a mixture of Poisson distributions via the ${\tt HTSCluster}$ R
package (version 2.0.10) (Rau et al., 2015). This approach performed well in
the fully constrained setting (i.e., simulation setting 2), with an average
ARI of 1.00. However, it faltered with settings 1 and 3, where unconstrained
elements were introduced. For settings 1 and 3, the average ARI of this
approach was 0.51 and 0.69, respectively. Further, it consistently selected
models with more groups than the actual number in the simulated datasets.
These results were anticipated as a known limitation of this approach is that
it cannot account for overdispersion, so to capture the larger variability, it
will use additional components. The second approach is a mixture of Negative
Binomial distributions, implemented via the ${\tt MBCluster.Seq}$ R package
(version 1.0) and labelled ‘mNB’ in Table 4 (Si et al., 2014). This approach
works very well on simulation setting 2, with an average ARI of 1.00 and
identifies the correct number of groups, but it does not perform well on the
other two settings. In Simulation setting 1, it has an average ARI of 0.36 and
two fewer groups than the true number of groups, and in setting 3, it has an
average ARI of 0.52 despite selecting the correct number of groups. A
limitation of this approach is that it assumes independence between variables,
and therefore cannot account for correlation, which greatly hinders its
performance. As seen in these results, our proposed MPLNFA model performed the
best for all simulation settings due to its ability to account for
overdispersion and both positive and negative correlation, both of which are
known limitations of the other two compared approaches.

###  4.2 Transcriptome Data Analysis

To illustrate the applicability of the MPLNFA model, it was applied to a
publicly available, curated RNA sequencing dataset of breast invasive
carcinoma (BRCA_RNASeqGene-20160128) from The Cancer Genome Atlas (TCGA)
Program (Ramos et al., 2020). Breast cancer is a commonly diagnosed cancer and
a leading cause of cancer-related deaths in women worldwide (Wilkinson and
Gathani, 2022). Clustering has been applied to understand breast cancer, by
clustering of genes to identify co-expression networks, or clustering of
tissues to identify subtypes (Sørlie et al., 2001; Kreike et al., 2007). The
original dataset explored here had 878 patient tissues and 20,502 genes. The
status of the estrogen receptor (ER), progesterone receptor (PR), and human
epidermal growth factor receptor 2 (HER2) play a crucial role in breast cancer
and its clinical management (Chen et al., 2022). Here, we focused only on a
subset of tissues for which hormone receptor status was available and were
females, resulting in a sample size of 712 tissues. Typically, only a subset
of genes from the experiment are used for analysis purposes in order to reduce
noise. The original number of 20,502 genes in the dataset was filtered to
remove genes with missing expression values, resulting in 17,586 genes. The
top 50 genes with the most variable gene expression across the 712 tissues
were selected. The final dataset had a size of 712 tissues with 50 genes.

In terms of patient characterization, cohort ages ranged from 26 to 90, with
no age data available for 9%. In terms of the pathologic stage, tissues from
11 stages (stages i, ia, ib, ii, iia, iib, iiia, iiib, iiic, iv, x) were
present, with no data for 0.98%. Staging quantifies the amount of cancer in
the body, with a higher stage number indicating more cancer spread (Cuthrell
and Tzenios, 2023). Stage x indicates that the status cannot be determined. A
majority of the cohort was identified as white (74%), followed by black (8%),
Asian (7%), American Indian or Alaska Native (0.15%), with no data available
for 10.85%. In the cohort, 76.5% were ER-positive and 66.7% were PR-positive.
In terms of HER2, 80.7% were negative, 13.8% were positive, 1.1% were
equivocal and 1.8% did not have data. All patients were female.

In order to identify subtypes of breast invasive carcinoma, the tissues were
clustered for a clustering range from $G=1\text{ to }8$, and a latent factor
range of $k=1\text{ to }10$, for all eight models listed in Table 1. All data
analyses were performed using $k$-means initialization with $3$ runs. Both BIC
and ICL selected a model with $G=5,k=6$ and $M=\text{CUU}$. The 5 clusters
analyzed here are interchangeably referred to as the 5 subtypes of breast
invasive carcinoma. It was observed that the tissues were assigned to clusters
with high posterior probabilities via the MPLNFA model (see Additional Figure
1). In this model, Clusters 1–5 were composed of
$201(28\%),118(17\%),83(12\%),61(8\%),\text{and }249(35\%)$ tissues,
respectively. See Additional File 1 for the tissue composition of each
cluster. In simulation studies, where the true labels are available, the ARI
can be used to assess clustering performance. Because no true labels are
available in real data, the predicted clusters were assessed in terms of
overall survival, clinical characteristics and gene expression patterns. The
distinct expression patterns between the 5 clusters are provided in Figure 1
(Gu, 2022). The genetic subtypes were associated with a divergent overall
survival (OS) in the overall cohort (P <0.0001) as shown in Figure 2. The most
favourable OS was observed for Cluster 2 and the most unfavourable OS was
observed in Clusters 4, followed by 3. In terms of race, age and pathologic
stage, no significant associations were identified with the clusters. However,
it was interesting to note that all races were represented in Cluster 4 and
Cluster 4 also had the least proportion of tissues from stage i.

Figure 1: Heatmap showing, for Cluster 1 (C1) through Cluster 5 (C5), log-
transformed gene expression patterns. This is for the $G=5,k=6$ and
$M=\text{CUU}$ model selected by both BIC and ICL for the breast invasive
carcinoma RNA-seq dataset (n = 712). The red and blue colours represent the
expression levels, where red represents high expression and blue represents
low expression. The rows represent the tissues and the columns represent the
top $50$ most variable genes included in the study. ER is estrogen receptor,
PR is progesterone receptor, and HER2 is human epidermal growth factor
receptor 2.

Figure 2: Kaplan-Meier plot of overall survival for Cluster 1 (C1) through
Cluster 5 (C5) for the $G=5,k=6$ and $M=\text{CUU}$ model selected by both BIC
and ICL for the breast invasive carcinoma RNA-seq dataset.

The ER and PR positive status has shown to be associated with favourable
patient outcomes (Davey et al., 2021). Interestingly, Cluster 4 and Cluster 2,
both contained the highest number of ER and PR negative cases. However, among
all clusters, Cluster 4 contained the highest number of HER2 positive cases,
while Cluster 2 had the most HER2 negative cases. HER2 overexpression or HER2
positive status in breast cancer has shown to be associated with poor
prognosis and decreased survival (Chen et al., 2022). The PR is an estrogen-
regulated gene, so ER-positive tumours are generally PR positive, and ER-
negative tumours are PR negative (Bae et al., 2015). This is observed across
clusters in the current analysis. Further, the ER positive, PR positive status
is known to have a better prognosis compared to ER negative, PR negative
status (Bae et al., 2015). This is evident for Cluster 1 and 5, where the
majority of cases have both receptors positive.

To explore expression patterns, 50 genes in the dataset were divided into non-
exclusive, three groups: developmental-related, metabolism and enzymatic-
related or neurological-related genes. Cluster 1-3, and 5 showed variable
expression levels among the three gene groups. However, Cluster 4 stood out as
having high expression in developmental-related genes, followed by a mixed
expression in metabolism-related genes, and low expression in neurological-
related genes. The expression patterns observed here were supported in the
literature. We selected some genes to illustrate this. For example, within
development-related genes, NKAIN1 has high expression across all clusters.
This gene encodes the sodium/potassium transporting ATPase interacting with 1
protein in humans. An association between high expression of NKAIN1 and breast
cancer has been shown and NKAIN1 mRNA levels have been identified as a
potential biomarker in the diagnosis of breast cancer (Ali-Hassanzadeh et al.,
2022). Further, HS6ST3 is the Heparan Sulfate 6-O-Sulfotransferase 3 protein-
coding gene. It has been shown to be highly expressed in breast cancer cell
lines and its silencing has diminished the tumour growth and progression
(Iravani et al., 2017). Interestingly, Cluster 4 showed a high expression of
HS6ST3, Cluster 3 showed an intermediate level and Cluster 2 showed the lowest
expression, in accordance with OS differences observed.

In terms of metabolism-related genes, GSTA1 is the Glutathione S-Transferase
Alpha 1 protein-coding gene. Reduced GSTA1 expression levels have been shown
to be associated with increased breast cancer (Ahn et al., 2006), while higher
expression has shown longer OS and disease-free survival (Liu et al., 2020).
Here, Cluster 2 showed low expression of GSTA1 and Cluster 4 showed higher
expression, agreeing with OS observed in the current study. The C1orf64 gene
encodes the steroid receptor-associated and regulated protein. It has been
shown to have a high expression in breast tumours (Naderi, 2017). Here,
Cluster 4 showed high expression of C1orf64 and Cluster 2 showed lower
expression. In terms of neurological-related genes, ZIC1 is the zinc finger of
the cerebellum 1 gene. It is known to have a tumour-suppressive role in breast
cancer (Han et al., 2018). An increase in ZIC1 expression has been shown to
suppress the growth of breast cancer cells and xenograft tumours (Han et al.,
2018). Interestingly, Cluster 4 showed low expression of ZIC1 and Cluster 2
showed a high expression, supporting the OS differences observed.

##  5 Discussion

A mixture of factor analyzers model for MPLN distribution as well as a family
of mixture models based theorem is introduced. This is the first use of a
mixture of MPLN factor analyzer distributions within the literature. To our
knowledge, this is also the first use of a mixture of discrete factor
analyzers within the literature. The proposed models are well-suited to high-
dimensional applications as the number of scale parameters is linear in data
dimensionality for all eight models, as opposed to in traditional MPLN, where
the parameters grow quadratically. Further, a single MPLN factor analysis
model can be obtained as a special case of the mixture of MPLNFA, i.e., with
$G=1$. Extensions applicable to factor analyzers model, such as the mixture of
common factor analyzers (Baek et al., 2010) can be applied to the MPLNFA
family in future. A mixture of common factor analyzers is a restrictive form
of the mixture of factor analyzers model and can be useful when the number of
clusters and dimensionality are very large.

Overall, applying MPLNFA identified 5 subtypes of breast invasive carcinoma
that differ from each other in terms of patient survival outcomes, clinical
covariates and gene expression. Connecting gene expression patterns with
clinical covariates and patient outcomes is important in understanding disease
diversity. Findings from our study provide a deeper understanding of the
molecular and clinical heterogeneity of breast cancer and these may inform the
development of subtype-specific treatment strategies to improve patient
outcomes in future.

### Acknowledgments

This research was supported by the Postdoctoral Fellowship award from the
Canadian Institutes of Health Research (Silva), the National Sciences and
Engineering Research Council of Canada grant (2021-03812 - Subedi, 2023-06030
- McNicholas), and the Canada Research Chair Program (22020-00303 - Subedi,
2022-00494 - McNicholas).

###  5.1 Software availability

The GitHub R package for this work is available at:
https://github.com/anjalisilva/mixMPLNFA and is released under the open source
MIT license.

##  Appendix A Appendix

###  A.1 Variational Approximations

The marginal log-density of $\mathbf{Y}_{i}\mid Z_{ig}=1$ in Stage 1 can be
written as:

$\displaystyle\log f(\mathbf{y}_{i}\mid Z_{ig}=1)$
$\displaystyle=\int_{{}^{d}}\log f(\mathbf{y}_{i}\mid
Z_{ig}=1)\leavevmode\nobreak\ q(\mathbf{x}_{ig})\leavevmode\nobreak\
d\mathbf{x}_{ig}$
$\displaystyle=\int_{{}^{d}}\log\frac{f(\mathbf{y}_{i},\mathbf{x}_{i}\mid
Z_{ig}=1)/q(\mathbf{x}_{ig})}{f(\mathbf{x}_{i}\mid\mathbf{y}_{i},Z_{ig}=1)/q(\mathbf{x}_{ig})}\leavevmode\nobreak\
q(\mathbf{x}_{ig})\leavevmode\nobreak\ d\mathbf{x}_{ig}$
$\displaystyle=\int_{{}^{d}}\left[\log f(\mathbf{y}_{i},\mathbf{x}_{i}\mid
Z_{ig}=1)-\log q(\mathbf{x}_{ig})\right]q(\mathbf{x}_{ig})\leavevmode\nobreak\
d\mathbf{x}_{ig}+D_{KL}(q_{ig}\|f_{ig})$
$\displaystyle=F(q_{ig},\mathbf{y}_{i})+D_{KL}(q_{ig}\|f_{ig}),$

The marginal log-density of $\mathbf{Y}_{i}\mid Z_{ig}=1$ in Stage 2 can be
written as:

$\displaystyle\log f(\mathbf{y}_{i}\mid Z_{ig}=1)$
$\displaystyle=\int_{{}^{K}}\int_{{}^{d}}\log f(\mathbf{y}_{i}\mid
Z_{ig}=1)\leavevmode\nobreak\
q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})\leavevmode\nobreak\
d\mathbf{x}_{ig}\leavevmode\nobreak\ d\mathbf{u}_{ig}$
$\displaystyle=\int_{{}^{K}}\int_{{}^{d}}\log\frac{f(\mathbf{y}_{i},\mathbf{x}_{i},\mathbf{u}_{i}\mid
Z_{ig}=1)/q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})}{f(\mathbf{x}_{i},\mathbf{u}_{i}\mid\mathbf{y}_{i},Z_{ig}=1)/q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})}\leavevmode\nobreak\
q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})\leavevmode\nobreak\
d\mathbf{x}_{ig}$ $\displaystyle=\int_{{}^{K}}\int_{{}^{d}}\Big{[}\log
f(\mathbf{y}_{i},\mathbf{x}_{i},\mathbf{u}_{i}\mid Z_{ig}=1)-\log
q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})\Big{]}$
$\displaystyle\times
q_{\mathbf{x},\mathbf{u}}(\mathbf{x}_{ig},\mathbf{u}_{ig})\leavevmode\nobreak\
d\mathbf{x}_{ig}\leavevmode\nobreak\ d\mathbf{u}_{ig}+D_{KL}(q_{ig}\|f_{ig})$
$\displaystyle=F(q_{ig},\mathbf{y}_{i})+D_{KL}(q_{ig}\|f_{ig}),$

ELBO in Stage 2 under the assumption that $q_{x,u}(\mathbf{x},\mathbf{u})$ is
factorable:

$\displaystyle F(q_{ig},\mathbf{y}_{i})$
$\displaystyle=\int_{{}^{K}}\int_{{}^{d}}\Big{[}\log
f(\mathbf{y}_{i}\mid\mathbf{x}_{i},\mathbf{u}_{i},Z_{ig}=1)+\log
f(\mathbf{x}_{i}\mid\mathbf{u}_{i},Z_{ig}=1)+\log f(\mathbf{u}_{i}\mid
Z_{ig}=1)$ $\displaystyle-\log q_{\mathbf{x}}(\mathbf{x}_{ig})-\log
q_{\mathbf{u}}(\mathbf{u}_{ig})\Big{]}\leavevmode\nobreak\
q_{\mathbf{x}}(\mathbf{x}_{ig})\leavevmode\nobreak\
q_{\mathbf{u}}(\mathbf{u}_{ig})\leavevmode\nobreak\
d\mathbf{x}_{ig}\leavevmode\nobreak\ d\mathbf{u}_{ig}$
$\displaystyle=\int_{{}^{K}}\Bigg{[}\left\\{\int_{{}^{d}}\Big{[}\log
f(\mathbf{y}_{i}\mid\mathbf{x}_{i},\mathbf{u}_{i},Z_{ig}=1)+\log
f(\mathbf{x}_{i}\mid\mathbf{u}_{i},Z_{ig}=1)\right.$ $\displaystyle\left.-\log
q_{\mathbf{x}}(\mathbf{x}_{i})\Big{]}\leavevmode\nobreak\
q_{\mathbf{x}}(\mathbf{x}_{ig})\leavevmode\nobreak\
d\mathbf{x}_{ig}\right\\}+\log f(\mathbf{u}_{i}\mid Z_{ig}=1)-\log
q_{\mathbf{u}}(\mathbf{u}_{ig})\Bigg{]}q_{\mathbf{u}}(\mathbf{u}_{ig})\leavevmode\nobreak\
d\mathbf{u}_{ig}.$

Assuming that the approximating distributions (i.e.
$q_{\mathbf{x}}(\mathbf{x}_{ig})=N(\mathbf{m}_{ig},\mathbf{S}_{ig})$ and
$q_{\mathbf{u}}(\mathbf{u}_{ig})=N(\mathbf{P}_{ig},\mathbf{Q}_{ig})$), the
ELBO in Stage 2 becomes:

$\displaystyle F(q_{ig},\mathbf{y}_{i})$
$\displaystyle=\int_{{}^{K}}\Bigg{[}\mathbf{m}_{ig}^{\prime}\mathbf{y}_{i}+\sum_{j=1}^{d}\left(\log
k_{j}\right)y_{ij}-\sum_{j=1}^{d}\left\\{e^{\log
k_{j}+m_{igj}+\frac{1}{2}S_{ig,jj}}+\log(y_{ij}!)\right\\}$
$\displaystyle+\frac{1}{2}\log|\mathbf{S}_{ig}|+\frac{d}{2}-\frac{1}{2}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{g}-\boldsymbol{\Lambda}_{g}\mathbf{u}_{i})^{\prime}\boldsymbol{\Psi}_{g}^{-1}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{g}-\boldsymbol{\Lambda}_{g}\mathbf{u}_{i})-\frac{1}{2}\,\mbox{tr}\left(\boldsymbol{\Psi}_{g}^{-1}\mathbf{S}_{ig}\right)$
$\displaystyle-\frac{1}{2}\log|\boldsymbol{\Psi}_{g}|+\log
f(\mathbf{u}_{i}\mid Z_{ig}=1)-\log
q_{\mathbf{u}}(\mathbf{u}_{ig})\Bigg{]}q_{\mathbf{u}}(\mathbf{u}_{ig})\leavevmode\nobreak\
d\mathbf{u}_{ig}$
$\displaystyle=\mathbf{m}_{ig}^{\prime}\mathbf{y}_{i}+\sum_{j=1}^{d}\left(\log
k_{j}\right)y_{ij}-\sum_{j=1}^{d}\left\\{e^{\log
k_{j}+m_{igj}+\frac{1}{2}S_{ig,jj}}+\log(y_{ij}!)\right\\}+\frac{1}{2}\log|\mathbf{S}_{ig}|+\frac{d}{2}$
$\displaystyle-\frac{1}{2}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})^{T}\boldsymbol{\Psi}_{g}^{-1}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})+(\mathbf{m}_{ig}-\boldsymbol{\mu}_{g})^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}-\frac{1}{2}\mathbf{P}_{ig}^{T}\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}$
$\displaystyle-\frac{1}{2}\,\mbox{tr}(\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{Q}_{ig})-\frac{1}{2}\,\mbox{tr}\left(\boldsymbol{\Psi}_{g}^{-1}\mathbf{S}_{ig}\right)-\frac{1}{2}\log|\boldsymbol{\Psi}_{g}|-\frac{1}{2}\mathbf{P}_{ig}^{T}\mathbf{P}_{ig}+\frac{1}{2}\log|\mathbf{Q}_{ig}^{-1}|$
$\displaystyle-\frac{1}{2}\,\mbox{tr}(\mathbf{Q}_{ig})+\frac{K}{2}.$

###  A.2 Family of Models

####  A.2.1 Isotropic assumption:
$\boldsymbol{\Psi}_{g}=\psi_{g}\mathbf{I}_{d}$

$\displaystyle l(\boldsymbol{\vartheta})$
$\displaystyle=\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\log\pi_{g}+\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\log
f(\mathbf{y}_{i}\mid
Z_{ig}=1,\boldsymbol{\mu}_{g},\boldsymbol{\Psi}_{g},\boldsymbol{\Lambda}_{g}).$
$\displaystyle=\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\log\pi_{g}+\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\left\\{F(q_{ig},\mathbf{y}_{i})+D_{KL}(q_{ig}\mid\mid
f_{ig})\right\\}$

Under the isotropic constraint (i.e.
$\boldsymbol{\Psi}_{g}=\psi_{g}\mathbf{I}_{d}$),

$\displaystyle\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}F(q_{ig},\mathbf{y}_{i})$
$\displaystyle=\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\left[-\frac{\psi_{g}^{-1}}{2}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})^{T}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})+\psi_{g}^{-1}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{g})^{T}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}\right.$
$\displaystyle\left.+\frac{d}{2}\log\psi_{g}^{-1}-\frac{\psi_{g}^{-1}}{2}\mathbf{P}_{ig}^{T}\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}-\frac{\psi_{g}^{-1}}{2}\,\mbox{tr}(\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Lambda}_{g}\mathbf{Q}_{g})-\frac{\psi_{g}^{-1}}{2}\,\mbox{tr}\left(\mathbf{S}_{ig}\right)\right]+C,$
$\displaystyle=\sum_{g=1}^{G}\left[-\frac{\psi_{g}^{-1}}{2}n_{g}\,\mbox{tr}(\mathbf{W}_{g})+\psi_{g}^{-1}n_{g}\,\mbox{tr}(\mathbf{W}_{g}\boldsymbol{\Lambda}_{g}\boldsymbol{\beta}_{g})+\frac{dn_{g}}{2}\log\psi_{g}^{-1}\right.$
$\displaystyle\left.-\frac{\psi_{g}{{}^{-1}}n_{g}}{2}\,\mbox{tr}(\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Lambda}_{g}\mathbf{Q}_{g})-\frac{\psi_{g}{{}^{-1}}}{2}n_{g}\,\mbox{tr}(\boldsymbol{\beta}^{T}\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Lambda}_{g}\boldsymbol{\beta}_{g}\mathbf{W}_{g})\right.$
$\displaystyle\left.-\frac{\psi_{g}^{-1}}{2}\,\mbox{tr}\left(\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)\right]+C,$

where $C$ is a constant with respect to $\psi_{g}^{-1}$ and
$n_{g}=\sum_{i=1}^{n}z_{ig}$. Taking derivative with respect to
$\psi^{-1}_{g}$ and setting it to 0, we get

$\displaystyle d\hat{\psi}_{g}$
$\displaystyle-\,\mbox{tr}(\mathbf{W}_{g})+2\,\mbox{tr}(\mathbf{W}_{g}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g})-\,\mbox{tr}(\hat{\boldsymbol{\beta}}^{T}\hat{\boldsymbol{\Lambda}}_{g}^{T}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g})$
$\displaystyle-\,\mbox{tr}(\hat{\boldsymbol{\Lambda}}_{g}^{T}\hat{\boldsymbol{\Lambda}}_{g}\mathbf{Q}_{g})-\,\mbox{tr}\left(\frac{\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}}{n_{g}}\right)=0.$

Simplifying this, we get

$\hat{\psi}_{g}=\frac{1}{d}\,\mbox{tr}\left(\mathbf{W}_{g}-\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}}{n_{g}}\right).$
(6)

####  A.2.2 Equal Variance: $\boldsymbol{\Psi}_{g}=\boldsymbol{\Psi}$

Under the equal variance constraint (i.e.
$\boldsymbol{\Psi}_{g}=\boldsymbol{\Psi}$),

$\displaystyle\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}F(q_{ig},\mathbf{y}_{i})$
$\displaystyle=\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\left[-\frac{1}{2}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})^{T}\boldsymbol{\Psi}^{-1}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})\right.$
$\displaystyle\left.+(\mathbf{m}_{ig}-\boldsymbol{\mu}_{g})^{T}\boldsymbol{\Psi}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}+\frac{1}{2}\log|\boldsymbol{\Psi}^{-1}|-\frac{1}{2}\mathbf{P}_{ig}^{T}\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Psi}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}\right.$
$\displaystyle\left.-\frac{1}{2}\,\mbox{tr}(\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Psi}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{Q}_{g})-\frac{1}{2}\,\mbox{tr}\left(\boldsymbol{\Psi}^{-1}\mathbf{S}_{ig}\right)\right]+C,$
$\displaystyle=\sum_{g=1}^{G}\left[-\frac{n_{g}}{2}\,\mbox{tr}(\boldsymbol{\Psi}^{-1}\mathbf{W}_{g})+n_{g}\,\mbox{tr}(\mathbf{W}_{g}\boldsymbol{\Psi}^{-1}\boldsymbol{\Lambda}_{g}\boldsymbol{\beta}_{g})-\frac{n_{g}}{2}\,\mbox{tr}(\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Psi}^{-1}\boldsymbol{\Lambda}_{g}\mathbf{Q}_{g})\right.$
$\displaystyle\left.+\frac{n_{g}}{2}\log|\boldsymbol{\Psi}^{-1}|-\frac{\psi_{g}{{}^{-1}}}{2}n_{g}\,\mbox{tr}(\boldsymbol{\beta}^{T}\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Psi}^{-1}\boldsymbol{\Lambda}_{g}\boldsymbol{\beta}_{g}\mathbf{W}_{g})\right.$
$\displaystyle\left.-\frac{1}{2}\,\mbox{tr}\left(\boldsymbol{\Psi}^{-1}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)\right]+C,$

where $C$ is a constant with respect to $\boldsymbol{\Psi}$ and
$n_{g}=\sum_{i=1}^{n}z_{ig}$. Taking derivative with respect to
$\boldsymbol{\Psi}^{-1}$ and setting it to 0, we get

$\displaystyle n\hat{\boldsymbol{\Psi}}$
$\displaystyle-\sum_{g=1}^{G}n_{g}\mathbf{W}_{g}+2\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}-\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\beta}}^{T}\hat{\boldsymbol{\Lambda}}_{g}^{T}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}$
$\displaystyle-\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Lambda}}_{g}^{T}\hat{\boldsymbol{\Lambda}}_{g}\mathbf{Q}_{g}-\,\mbox{tr}\left(\frac{\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}}{n_{g}}\right)=0.$

Simplifying this, we get

$\displaystyle\hat{\boldsymbol{\Psi}}$
$\displaystyle=\,\mbox{diag}\left(\sum_{g=1}^{G}\frac{n_{g}}{n}\mathbf{W}_{g}-\sum_{g=1}^{G}\frac{n_{g}}{n}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right).$
$\displaystyle=\,\mbox{diag}\left(\sum_{g=1}^{G}\frac{n_{g}}{n}\mathbf{W}_{g}-\sum_{g=1}^{G}\frac{n_{g}}{n}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)$

####  A.2.3 Equal Variance and isotropic constraint
$\boldsymbol{\Psi}_{g}=\psi\mathbf{I}_{d}$

$\displaystyle\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}F(q_{ig},\mathbf{y}_{i})$
$\displaystyle=\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\left[-\frac{\psi^{-1}}{2}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})^{T}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{ig})+\psi^{-1}(\mathbf{m}_{ig}-\boldsymbol{\mu}_{g})^{T}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}\right.$
$\displaystyle\left.+\frac{d}{2}\log\psi^{-1}-\frac{\psi^{-1}}{2}\mathbf{P}_{ig}^{T}\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Lambda}_{g}\mathbf{P}_{ig}-\frac{\psi^{-1}}{2}\,\mbox{tr}(\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Lambda}_{g}\mathbf{Q}_{g})-\frac{\psi^{-1}}{2}\,\mbox{tr}\left(\mathbf{S}_{ig}\right)\right]+C,$
$\displaystyle=\sum_{g=1}^{G}\left[-\frac{\psi^{-1}}{2}n_{g}\,\mbox{tr}(\mathbf{W}_{g})+\psi^{-1}n_{g}\,\mbox{tr}(\mathbf{W}_{g}\boldsymbol{\Lambda}_{g}\boldsymbol{\beta}_{g})+\frac{dn_{g}}{2}\log\psi^{-1}\right.$
$\displaystyle\left.-\frac{\psi{{}^{-1}}n_{g}}{2}\,\mbox{tr}(\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Lambda}_{g}\mathbf{Q}_{g})-\frac{\psi{{}^{-1}}}{2}n_{g}\,\mbox{tr}(\boldsymbol{\beta}^{T}\boldsymbol{\Lambda}_{g}^{T}\boldsymbol{\Lambda}_{g}\boldsymbol{\beta}_{g}\mathbf{W}_{g})\right.$
$\displaystyle\left.-\frac{\psi^{-1}}{2}\,\mbox{tr}\left(\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)\right]+C,$

where $C$ is a constant with respect to $\psi^{-1}$ and
$n_{g}=\sum_{i=1}^{n}z_{ig}$. Taking derivative with respect to
$\psi^{-1}_{g}$ and setting it to 0, we get

$\displaystyle dn\hat{\psi}$
$\displaystyle-\,\mbox{tr}(\sum_{g=1}^{G}n_{g}\mathbf{W}_{g})+2\,\mbox{tr}(\sum_{g=1}^{G}n_{g}\mathbf{W}_{g}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g})-\,\mbox{tr}(\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\beta}}^{T}\hat{\boldsymbol{\Lambda}}_{g}^{T}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g})$
$\displaystyle-\,\mbox{tr}(\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Lambda}}_{g}^{T}\hat{\boldsymbol{\Lambda}}_{g}\mathbf{Q}_{g})-\,\mbox{tr}\left(\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)=0.$

And, simplifying this, we get

$\displaystyle\hat{\psi}$
$\displaystyle=\frac{1}{d}\,\mbox{tr}\left(\sum_{g=1}^{G}\frac{n_{g}}{n}\mathbf{W}_{g}-\sum_{g=1}^{G}\frac{n_{g}}{n}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)$
(7)

####  A.2.4 Equal Loading Matrices:
$\boldsymbol{\Lambda}_{g}=\boldsymbol{\Lambda}$

Under the equal loading matrices constraint (i.e.
$\boldsymbol{\Lambda}_{g}=\boldsymbol{\Lambda}$),

$\displaystyle\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}F(q_{ig},\mathbf{y}_{i})$
$\displaystyle=\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\left((\mathbf{m}_{ig}-\boldsymbol{\mu}_{g})^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}\mathbf{P}_{ig}-\frac{1}{2}\mathbf{P}_{ig}^{T}\boldsymbol{\Lambda}^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}\mathbf{P}_{ig}\right.$
$\displaystyle\left.-\frac{1}{2}\,\mbox{tr}(\boldsymbol{\Lambda}^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}\mathbf{Q}_{g})\right)$
$\displaystyle=\sum_{g=1}^{G}\left(n_{g}\,\mbox{tr}\left(\boldsymbol{\beta}_{g}\mathbf{W}_{g}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}\right)-\frac{n_{g}}{2}\,\mbox{tr}\left(\boldsymbol{\beta}_{g}^{T}\boldsymbol{\Lambda}^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}\boldsymbol{\beta}_{g}\mathbf{W}_{g}\right)\right.$
$\displaystyle\left.-\frac{n_{g}}{2}tr(\boldsymbol{\Lambda}^{T}\boldsymbol{\Psi}_{g}^{-1}\boldsymbol{\Lambda}\mathbf{Q}_{g})\right).$

Taking derivative with respect to $\boldsymbol{\Lambda}$ and setting it equal
to 0, we get

$\displaystyle 0$
$\displaystyle=\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Psi}}_{g}^{-1}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}_{g}^{T}-\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Psi}}_{g}^{-1}\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}_{g}^{T}-\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Psi}}_{g}^{-1}\hat{\boldsymbol{\Lambda}}\mathbf{Q}_{g}$
$\displaystyle=\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Psi}}_{g}^{-1}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}_{g}^{T}-\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Psi}}_{g}^{-1}\hat{\boldsymbol{\Lambda}}\left(\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}_{g}^{T}-\mathbf{I}_{p}-\hat{\boldsymbol{\beta}}_{g}\hat{\boldsymbol{\Lambda}}\right)$
$\displaystyle=\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Psi}}_{g}^{-1}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}_{g}^{T}-\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Psi}}_{g}^{-1}\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\Theta}}_{g}.$

In this case, the loading matrix cannot be solved directly and must be solved
in a row-by-row manner as suggested by McNicholas and Murphy (2008). Hence,

$\displaystyle\hat{\lambda}_{i}$
$\displaystyle=\mathbf{r}_{i}\left(\sum_{g=1}^{G}\frac{n_{g}}{\hat{\Psi}_{g,(i)}}\hat{\boldsymbol{\Theta}}_{g}\right)^{-1},$
(8)

where $\lambda_{i}$ is the $i^{th}$ row of the matrix $\boldsymbol{\Lambda}$,
$\Psi_{g,(i)}$ is the $i^{th}$ diagonal element of $\boldsymbol{\Psi}_{g}$ and
$\mathbf{r}_{i}$ is the $i^{th}$ row of the matrix
$\sum_{g=1}^{G}n_{g}\hat{\boldsymbol{\Psi}}_{g}^{-1}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}_{g}^{T}$.

###  A.3 Parameter Updates

Table 5: Parameter Updates for the family of models based on constraints on $\boldsymbol{\Lambda}_{g}$ and $\boldsymbol{\Psi}_{g}$. Model |  Parameters |  Estimates
---|---|---
UUU |  $\boldsymbol{\beta}_{g}$ |  $\hat{\boldsymbol{\beta}}_{g}=\hat{\boldsymbol{\Lambda}}_{g}^{T}(\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\Lambda}}_{g}^{T}+\boldsymbol{\Psi}_{g})^{-1}$
|  $\boldsymbol{\Theta}_{g}$ |  $\hat{\boldsymbol{\Theta}}_{g}=\mathbf{I}_{K}-\hat{\boldsymbol{\beta}_{g}}\hat{\boldsymbol{\Lambda}}_{g}+\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}\hat{\boldsymbol{\beta}_{g}}^{T}$
|  $\boldsymbol{\Lambda}_{g}$ |  $\hat{\boldsymbol{\Lambda}}_{g}=\mathbf{W}_{g}\hat{\boldsymbol{\beta}}^{T}\widehat{\boldsymbol{\Theta}}_{g}^{-1}$
|  $\boldsymbol{\Psi}_{g}$ |  $\hat{\boldsymbol{\Psi}}_{g}=\,\mbox{diag}\left(\mathbf{W}_{g}-\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{\sum_{i=1}^{n}\widehat{Z}_{ig}\mathbf{S}_{ig}}{\sum_{i=1}^{n}\widehat{Z}_{ig}}\right)$
UUC |  $\boldsymbol{\beta}_{g}$ |  $\hat{\boldsymbol{\beta}}_{g}=\hat{\boldsymbol{\Lambda}}_{g}^{T}(\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\Lambda}}_{g}^{T}+\psi_{g}\mathbf{I}_{d})^{-1}$
|  $\boldsymbol{\Theta}_{g}$ |  $\hat{\boldsymbol{\Theta}}_{g}=\mathbf{I}_{K}-\hat{\boldsymbol{\beta}_{g}}\hat{\boldsymbol{\Lambda}}_{g}+\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}\hat{\boldsymbol{\beta}_{g}}^{T}$
|  $\boldsymbol{\Lambda}_{g}$ |  $\hat{\boldsymbol{\Lambda}}_{g}=\mathbf{W}_{g}\hat{\boldsymbol{\beta}}^{T}\widehat{\boldsymbol{\Theta}}_{g}^{-1}$
|  $\boldsymbol{\Psi}_{g}=\psi_{g}\mathbf{I}_{d}$ |  $\hat{\psi}_{g}=\frac{1}{d}\,\mbox{tr}\left(\mathbf{W}_{g}-\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}}{n_{g}}\right)$
UCU |  $\boldsymbol{\beta}_{g}$ |  $\hat{\boldsymbol{\beta}}_{g}=\hat{\boldsymbol{\Lambda}}_{g}^{T}(\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\Lambda}}_{g}^{T}+\boldsymbol{\Psi})^{-1}$
|  $\boldsymbol{\Theta}_{g}$ |  $\hat{\boldsymbol{\Theta}}_{g}=\mathbf{I}_{K}-\hat{\boldsymbol{\beta}_{g}}\hat{\boldsymbol{\Lambda}}_{g}+\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}\hat{\boldsymbol{\beta}_{g}}^{T}$
|  $\boldsymbol{\Lambda}_{g}$ |  $\hat{\boldsymbol{\Lambda}}_{g}=\mathbf{W}_{g}\hat{\boldsymbol{\beta}}^{T}\widehat{\boldsymbol{\Theta}}_{g}^{-1}$
|  $\boldsymbol{\Psi}_{g}=\boldsymbol{\Psi}$ |  $\hat{\boldsymbol{\Psi}}=\,\mbox{diag}\left(\sum_{g=1}^{G}\frac{n_{g}}{n}\mathbf{W}_{g}-\sum_{g=1}^{G}\frac{n_{g}}{n}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)$
UCC |  $\boldsymbol{\beta}_{g}$ |  $\hat{\boldsymbol{\beta}}_{g}=\hat{\boldsymbol{\Lambda}}_{g}^{T}(\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\Lambda}}_{g}^{T}+\psi\mathbf{I})^{-1}$
|  $\boldsymbol{\Theta}_{g}$ |  $\hat{\boldsymbol{\Theta}}_{g}=\mathbf{I}_{K}-\hat{\boldsymbol{\beta}_{g}}\hat{\boldsymbol{\Lambda}}_{g}+\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}\hat{\boldsymbol{\beta}_{g}}^{T}$
|  $\boldsymbol{\Lambda}_{g}$ |  $\hat{\boldsymbol{\Lambda}}_{g}=\mathbf{W}_{g}\hat{\boldsymbol{\beta}}^{T}\widehat{\boldsymbol{\Theta}}_{g}^{-1}$
|  $\boldsymbol{\Psi}_{g}=\psi\mathbf{I}$ |  $\hat{\psi}=\frac{1}{d}\,\mbox{tr}\left(\sum_{g=1}^{G}\frac{n_{g}}{n}\mathbf{W}_{g}-\sum_{g=1}^{G}\frac{n_{g}}{n}\hat{\boldsymbol{\Lambda}}_{g}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)$
CUU |  $\boldsymbol{\beta}_{g}$ |  $\hat{\boldsymbol{\beta}}_{g}=\hat{\boldsymbol{\Lambda}}^{T}(\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\Lambda}}^{T}+\boldsymbol{\Psi}_{g})^{-1}$
|  $\boldsymbol{\Theta}_{g}$ |  $\hat{\boldsymbol{\Theta}}_{g}=\mathbf{I}_{K}-\hat{\boldsymbol{\beta}_{g}}\hat{\boldsymbol{\Lambda}}+\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}\hat{\boldsymbol{\beta}_{g}}^{T}$
|  $\boldsymbol{\Lambda}_{g}=\boldsymbol{\Lambda}$ |  Updated using Equation 5.
|  $\boldsymbol{\Psi}_{g}$ |  $\hat{\boldsymbol{\Psi}}_{g}=\,\mbox{diag}\left(\mathbf{W}_{g}-\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{\sum_{i=1}^{n}\widehat{Z}_{ig}\mathbf{S}_{ig}}{\sum_{i=1}^{n}\widehat{Z}_{ig}}\right)$
CUC |  $\boldsymbol{\beta}_{g}$ |  $\hat{\boldsymbol{\beta}}_{g}=\hat{\boldsymbol{\Lambda}}^{T}(\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\Lambda}}^{T}+\psi_{g}\mathbf{I}_{d})^{-1}$
|  $\boldsymbol{\Theta}_{g}$ |  $\hat{\boldsymbol{\Theta}}_{g}=\mathbf{I}_{K}-\hat{\boldsymbol{\beta}_{g}}\hat{\boldsymbol{\Lambda}}+\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}\hat{\boldsymbol{\beta}_{g}}^{T}$
|  $\boldsymbol{\Lambda}_{g}=\boldsymbol{\Lambda}$ |  Updated using Equation 5.
|  $\boldsymbol{\Psi}_{g}=\psi_{g}\mathbf{I}_{d}$ |  $\hat{\psi}_{g}=\frac{1}{d}\,\mbox{tr}\left(\mathbf{W}_{g}-\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\beta}}_{g}\mathbf{W}_{g}+\frac{\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}}{n_{g}}\right)$
CCU |  $\boldsymbol{\beta}_{g}=\boldsymbol{\beta}$ |  $\hat{\boldsymbol{\beta}}=\hat{\boldsymbol{\Lambda}}^{T}(\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\Lambda}}^{T}+\boldsymbol{\Psi})^{-1}$
|  $\boldsymbol{\Theta}_{g}$ |  $\hat{\boldsymbol{\Theta}}_{g}=\mathbf{I}_{K}-\hat{\boldsymbol{\beta}}\hat{\boldsymbol{\Lambda}}+\hat{\boldsymbol{\beta}}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}^{T}$
|  $\boldsymbol{\Lambda}_{g}=\boldsymbol{\Lambda}$ |  Updated using Equation 5.
|  $\boldsymbol{\Psi}_{g}=\boldsymbol{\Psi}$ |  $\hat{\boldsymbol{\Psi}}=\,\mbox{diag}\left(\sum_{g=1}^{G}\frac{n_{g}}{n}\mathbf{W}_{g}-\sum_{g=1}^{G}\frac{n_{g}}{n}\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\beta}}\mathbf{W}_{g}+\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)$
CCC |  $\boldsymbol{\beta}_{g}=\boldsymbol{\beta}$ |  $\hat{\boldsymbol{\beta}}=\hat{\boldsymbol{\Lambda}}^{T}(\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\Lambda}}^{T}+\psi\mathbf{I}_{d})^{-1}$
|  $\boldsymbol{\Theta}_{g}$ |  $\hat{\boldsymbol{\Theta}}_{g}=\mathbf{I}_{K}-\hat{\boldsymbol{\beta}}\hat{\boldsymbol{\Lambda}}+\hat{\boldsymbol{\beta}}\mathbf{W}_{g}\hat{\boldsymbol{\beta}}^{T}$
|  $\boldsymbol{\Lambda}_{g}=\boldsymbol{\Lambda}$ |  Updated using Equation 5.
|  $\boldsymbol{\Psi}_{g}=\psi\mathbf{I}$ |  $\hat{\psi}=\frac{1}{d}\,\mbox{tr}\left(\sum_{g=1}^{G}\frac{n_{g}}{n}\mathbf{W}_{g}-\sum_{g=1}^{G}\frac{n_{g}}{n}\hat{\boldsymbol{\Lambda}}\hat{\boldsymbol{\beta}}\mathbf{W}_{g}+\frac{1}{n}\sum_{g=1}^{G}\sum_{i=1}^{n}z_{ig}\mathbf{S}_{ig}\right)$

## References

  * Ahn et al. (2006) Ahn, J., M. D. Gammon, R. M. Santella, M. M. Gaudet, J. A. Britton, S. L. Teitelbaum, M. B. Terry, A. I. Neugut, S. M. Eng, Y. Zhang, C. Garza, and C. B. Ambrosone (2006, 04).  Effects of glutathione S -transferase A1 ( GSTA1 ) genotype and potential modifiers on breast cancer risk .  Carcinogenesis 27(9), 1876–1882.
  * Aitchison and Ho (1989) Aitchison, J. and C. Ho (1989).  The multivariate Poisson-log normal distribution.  Biometrika 76(4), 643–653.
  * Ali-Hassanzadeh et al. (2022) Ali-Hassanzadeh, M., A. Golestan, A. Tahmasebi, and A. Ramezani (2022, May).  In silico identification of nkain1 gene as a potential breast cancer associated biomarker and validation using real time pcr technique.  Journal of Jiroft University of Medical Sciences.
  * Arridge et al. (2018) Arridge, S. R., K. Ito, B. Jin, and C. Zhang (2018).  Variational gaussian approximation for poisson data.  Inverse Problems 34(2), 025005.
  * Bae et al. (2015) Bae, S. Y., S. Kim, J. H. Lee, H.-c. Lee, S. K. Lee, W. H. Kil, S. W. Kim, J. E. Lee, and S. J. Nam (2015).  Poor prognosis of single hormone receptor-positive breast cancer: similar outcome as triple-negative breast cancer.  BMC Cancer 15(1), 1–9.
  * Baek et al. (2010) Baek, J., G. J. McLachlan, and L. K. Flack (2010).  Mixtures of factor analyzers with common factor loadings: Applications to the clustering and visualization of high-dimensional data.  IEEE Transactions on Pattern Analysis and Machine Intelligence 32(7), 1298–1309.
  * Bouveyron and Brunet (2012) Bouveyron, C. and C. Brunet (2012).  Simultaneous model-based clustering and visualization in the fisher discriminative subspace.  Statistics and Computing 22(1), 301–324.
  * Browne and McNicholas (2015) Browne, R. P. and P. D. McNicholas (2015).  A mixture of generalized hyperbolic distributions.  Canadian Journal of Statistics 43(2), 176–198.
  * Chen et al. (2022) Chen, M., J. Wu, D. Liu, W. Chen, C. Lin, L. Andriani, S. Ding, O. Huang, J. He, X. Chen, W. Chen, Y. Li, K. Shen, and L. Zhu (2022).  Combined estrogen receptor and progesterone receptor level can predict survival outcome in human epidermal growth factor receptor 2-positive early breast cancer.  Clinical Breast Cancer 22(2), e147–e156.
  * Cuthrell and Tzenios (2023) Cuthrell, K. M. and N. Tzenios (2023).  Breast cancer: Updated and deep insights.  International Research Journal of Oncology 6(1), 104–118.
  * Dang et al. (2015) Dang, U. J., R. P. Browne, and P. D. McNicholas (2015).  Mixtures of multivariate power exponential distributions.  Biometrics 71(4), 1081–1089.
  * Davey et al. (2021) Davey, M. G., É. J. Ryan, P. J. Folan, N. O’Halloran, M. R. Boland, M. K. Barry, K. J. Sweeney, C. M. Malone, R. J. McLaughlin, M. J. Kerin, et al. (2021).  The impact of progesterone receptor negativity on oncological outcomes in oestrogen-receptor-positive breast cancer.  BJS Open 5(3), zrab040.
  * de Souto et al. (2008) de Souto, M. C., I. G. Costa, D. S. de Araujo, T. B. Ludermir, and A. Schliep (2008).  Clustering cancer gene expression data: a comparative study.  BMC Bioinformatics 9(1), 1.
  * Di et al. (2011) Di, Y., D. W. Schafer, J. S. Cumbie, and J. H. Chang (2011).  The nbp negative binomial model for assessing differential gene expression from rna-seq.  Statistical Applications in Genetics and Molecular Biology 10(1).
  * Ghahramani et al. (1996) Ghahramani, Z., G. E. Hinton, et al. (1996).  The EM algorithm for mixtures of factor analyzers.  Technical report, Technical Report CRG-TR-96-1, University of Toronto.
  * Gollini and Murphy (2014) Gollini, I. and T. B. Murphy (2014).  Mixture of latent trait analyzers for model-based clustering of categorical data.  Statistics and Computing 24(4), 569–588.
  * Gu (2022) Gu, Z. (2022).  Complex heatmap visualization.  iMeta.
  * Han et al. (2018) Han, W., F. Cao, X.-J. Gao, H.-B. Wang, F. Chen, S.-J. Cai, C. Zhang, Y.-W. Hu, J. Ma, X. Gu, and et al. (2018, Sep).  Zic1 acts a tumor suppressor in breast cancer by targeting survivin.  International Journal of Oncology.
  * Hand (1989) Hand, D. J. (1989).  Journal of the Royal Statistical Society. Series C (Applied Statistics) 38(2), 384–385.
  * Iravani et al. (2017) Iravani, O., B.-H. Bay, and G. W.-C. Yip (2017).  Silencing hs6st3 inhibits growth and progression of breast cancer cells through suppressing igf1r and inducing xaf1.  Experimental Cell Research 350(2), 380–389.
  * Kosmidis and Karlis (2016) Kosmidis, I. and D. Karlis (2016).  Model-based clustering using copulas with applications.  Statistics and Computing 26(5), 1079–1099.
  * Kreike et al. (2007) Kreike, B., M. van Kouwenhove, H. Horlings, B. Weigelt, H. Peterse, H. Bartelink, and M. J. van de Vijver (2007).  Gene expression profiling and histopathological characterization of triple-negative/basal-like breast carcinomas.  Breast Cancer Research 9, 1–14.
  * Lawley and Maxwell (1962) Lawley, D. N. and A. E. Maxwell (1962).  Factor analysis as a statistical method.  Journal of the Royal Statistical Society. Series D (The Statistician) 12(3), 209–229.
  * Liu et al. (2020) Liu, X., X. Sui, C. Zhang, K. Wei, Y. Bao, J. Xiong, Z. Zhou, Z. Chen, C. Wang, H. Zhu, and et al. (2020, Jan).  Glutathione s-transferase a1 suppresses tumor progression and indicates better prognosis of human primary hepatocellular carcinoma.  Journal of Cancer 11.
  * McLachlan and Peel (2000) McLachlan, D. and D. Peel (2000).  Finite Mixture Models.  Wiley.
  * McNicholas (2016) McNicholas, P. D. (2016).  Mixture Model-Based Classification (1st ed.).  Chapman and Hall/CRC Press.
  * McNicholas and Murphy (2008) McNicholas, P. D. and T. B. Murphy (2008).  Parsimonious gaussian mixture models.  Statistics and Computing 18(3), 285–296.
  * McNicholas and Murphy (2010) McNicholas, P. D. and T. B. Murphy (2010, 08).  Model-based clustering of microarray expression data via latent Gaussian mixture models.  Bioinformatics 26(21), 2705–2712.
  * Naderi (2017) Naderi, A. (2017, May).  C1orf64 is a novel androgen receptor target gene and coregulator that interacts with 14-3-3 protein in breast cancer.  Oncotarget.
  * R Core Team (2022) R Core Team (2022).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Ramos et al. (2020) Ramos, M., L. Geistlinger, S. Oh, L. Schiffer, R. Azhar, H. Kodali, I. de Bruijn, J. Gao, V. J. Carey, M. Morgan, and L. Waldron (2020).  Multiomic integration of public oncology databases in bioconductor.  JCO Clinical Cancer Informatics (4), 958–971.  PMID: 33119407.
  * Rau et al. (2015) Rau, A., C. Maugis-Rabusseau, M.-L. Martin-Magniette, and G. Celeux (2015).  Co-expression analysis of high-throughput transcriptome sequencing data with Poisson mixture models.  Bioinformatics 31(9), 1420–1427.
  * Roberts et al. (2011) Roberts, A., H. Pimentel, C. Trapnell, and L. Pachter (2011).  Identification of novel transcripts in annotated genomes using RNA-seq.  Bioinformatics 27, 2325–2329.
  * Si et al. (2014) Si, Y., P. Liu, P. Li, and T. P. Brutnell (2014).  Model-based clustering for RNA-seq data.  Bioinformatics 30(2), 197–205.
  * Silva et al. (2023) Silva, A., X. Qin, S. J. Rothstein, P. D. McNicholas, and S. Subedi (2023).  Finite mixtures of matrix variate Poisson-log normal distributions for three-way count data.  Bioinformatics 39(5), btad167.
  * Silva et al. (2019) Silva, A., S. J. Rothstein, P. D. McNicholas, and S. Subedi (2019).  A multivariate poisson-log normal mixture model for clustering transcriptome sequencing data.  BMC Bioinformatics 20(1), 394.
  * Sørlie et al. (2001) Sørlie, T., C. M. Perou, R. Tibshirani, T. Aas, S. Geisler, H. Johnsen, T. Hastie, M. B. Eisen, M. Van De Rijn, S. S. Jeffrey, et al. (2001).  Gene expression patterns of breast carcinomas distinguish tumor subclasses with clinical implications.  Proceedings of the National Academy of Sciences 98(19), 10869–10874.
  * Spearman (1904) Spearman, C. (1904).  The proof and measurement of association between two things.  The American Journal of Psychology 15(1), 72–101.
  * Subedi and Browne (2020) Subedi, S. and R. Browne (2020).  A parsimonious family of multivariate poisson-lognormal distributions for clustering multivariate count data.  arXiv preprint arXiv:2004.06857.
  * Subedi and McNicholas (2014) Subedi, S. and P. D. McNicholas (2014).  Variational bayes approximations for clustering via mixtures of normal inverse gaussian distributions.  Advances in Data Analysis and Classification 8(2), 167–193.
  * Subedi and McNicholas (2020) Subedi, S. and P. D. McNicholas (2020).  A variational approximations-dic rubric for parameter estimation and mixture model selection within a family setting.  Journal of Classification, 1–20.
  * Subedi et al. (2015) Subedi, S., A. Punzo, S. Ingrassia, and P. D. McNicholas (2015).  Cluster-weighed $t$-factor analyzers for robust model-based clustering and dimension reduction.  Statistical Methods & Applications 24(4), 623–649.
  * Tortora et al. (2019) Tortora, C., B. C. Franczak, R. P. Browne, and P. D. McNicholas (2019).  A mixture of coalesced generalized hyperbolic distributions.  Journal of Classification 36(1), 26–57.
  * Vrbik and McNicholas (2014) Vrbik, I. and P. D. McNicholas (2014).  Parsimonious skew mixture models for model-based clustering and classification.  Computational Statistics and Data Analysis 71, 196–210.
  * Wainwright et al. (2008) Wainwright, M. J., M. I. Jordan, et al. (2008).  Graphical models, exponential families, and variational inference.  Foundations and Trends® in Machine Learning 1(1–2), 1–305.
  * Wang et al. (2009) Wang, Z., M. Gerstein, and M. Snyder (2009).  RNA-Seq: a revolutionary tool for transcriptomics.  Nature Review Genetics 10, 57–63.
  * Wilkinson and Gathani (2022) Wilkinson, L. and T. Gathani (2022, Feb).  Understanding breast cancer as a global health concern.  The British Journal of Radiology.
  * Wolfe (1965) Wolfe, J. (1965, 05).  A computer program for the maximum-likelihood analysis of types.  pp. 57.