\section{Introduction}
\label{sec:intro}

Finite mixture models have seen increasing use over the last decades \citep[for a recent survey, see][]{mcnicholas16b}.
Because of their flexibility, they are a suitable statistical tool for modeling a wide range of phenomena characterized by unobserved heterogeneity, and constitute a powerful device for clustering and classification.
Specifically, in their ``direct application'', each mixture component represents a group (or cluster) within data, and the scope is to identify these groups and estimate the parameters of the conditional-group distributions \citep{mclachlan2000finite, mcnicholas2016mixture}.
If no exogenous variables explain the means and the variances of each component, they are also called unconditional mixture models.
However, when there is a linear relationship between some variables, important insight can be gained by accounting for functional dependencies between them.
For this reason, finite mixtures of regression models with fixed covariates (FMR) have been proposed in the literature \citep[see][for examples]{desarbo1988, fruhwirth2006finite}.
Finite mixtures of regression models with concomitant variables (FMRC; \citealp{dayton1988concomitant}) are an extension of FMR where the mixing weights depend on some concomitant variables (which are often the same covariates) and are usually modeled by a multinomial logistic model (see \citealp{ingrassia2016decision, ingrassia2019cluster} and \citealp{Mazz:Batt:Ingr:Punz:Mode:2019} for details).
Unfortunately, these methodologies do not explicitly use the distribution of the covariates for clustering, i.e.~the assignment of data points to clusters does not directly utilize any information from the distribution of the covariates.
Differently from these approaches, finite mixtures of regressions with random covariates (\citealp{gershenfeld1997, gershenfeld1999}), also known as cluster-weighted models (CWMs), allow for such functional dependency.
This occurs because, for each mixture component, CWMs decompose the joint distribution of responses and covariates into the product between the marginal distribution of the covariates and the conditional distribution of the responses given the covariates.

Several CWMs have been introduced in the univariate and multivariate literature.
Most of them consider a univariate response variable, along with a set of covariates, modeled by a univariate and a multivariate distribution, respectively \citep[see][for examples]{ingrassia2012local, ingrassia2014model, punzo2014flexible}.
Instead, fewer CWMs exist in the case of a multivariate response \citep{punzo2017robust, dang2017multivariate}.
However, over the years, there has been an increasing interest in applications involving matrix-variate (three-way) data. 
This data structure can occur in several and different application domains, such as longitudinal data on multiple response variables, spatial multivariate data, multivariate repeated measures or spatio-temporal data.
Nevertheless, there exists a limited number of contributions involving matrix-variate regression models.
First introduced by \citet{viroli2012matrix}, these models are a natural generalization of the multivariate-multiple regression (see also \citealp{anderlucci2014matrix}).
In the mixture framework, finite mixtures of matrix-variate regressions (MN-FMR) have been recently considered by \citet{melnykov2019studying}.
There are no matrix-variate CWMs in the literature and this paper aims to fill this gap by introducing and discussing a matrix-variate CWM in which the cluster-specific marginal distribution of the covariates, and the cluster-specific conditional distribution of the responses given the covariates, are assumed to be matrix normal.

The remainder of the paper is organized as follows.
In Section~\ref{sec:back}, some preliminary aspects are described.
The matrix normal cluster-weighted model (MN-CWM) and the expectation-conditional maximization (ECM) algorithm for parameter estimation are discussed in Section~\ref{sec:model}. 
Computational and operational aspects are laid out in Section~\ref{sec:comp}.
In the simulation study outlined in Section~\ref{sec:artif}, the parameter recovery and the classification performance of the MN-CWM are investigated as well as the capability of the Bayesian information criterion (BIC; \citealp{schwarz1978estimating}) to detect the underlying group structure.
The MN-CWM is also therein compared to the MN-FMR and to the multivariate-multiple normal CWM (MMN-CWM).
The application of the MN-CWM to two real datasets concerning educational indicators and the Italian non-life insurance market is therefore analyzed in Section~\ref{sec:data}, whereas some conclusions and ideas for future developments are drawn in Section~\ref{sec:conclude}.

\section{Background}
\label{sec:back}

\subsection{Matrix normal distribution}
\label{subsec:MXVN}

A $p \times r$ continuous random matrix $\bY \in \real^{p \times r}$ has a matrix normal (MN) distribution, denoted by $\pazocal{N}_{p \times r} \left(\bM,\bPhi,\bPsi\right)$, if its density is 
\begin{equation}
\resizebox{.94\textwidth}{!}{
\bphi_{p \times r}\left(\bY;\bM,\bPhi,\bPsi\right)=\left(2\pi\right)^{-\frac{pr}{2}}\left|\bPhi\right|^{-\frac{r}{2}}\left|\bPsi\right|^{-\frac{p}{2}}\exp\left\{-\frac{1}{2}\tr\left[\bPhi^{-1}\left(\bY-\bM\right)\bPsi^{-1}\left(\bY-\bM\right)^\top\right]\right\},
$}
\label{eq:MatNorm}
\end{equation}
where $\bM$ is the $p\times r$ mean matrix, and $\bPhi$ and $\bPsi$ are the $p \times p$ and $r \times r$ covariance matrices associated with the $p$ rows and $r$ columns, respectively.
An equivalent definition specifies the $\left(p\times r\right)$-matrix normal distribution as a special case of the $pr$-multivariate normal distribution.
Specifically,
\begin{equation}
\bY \sim \pazocal{N}_{p \times r} \left(\bM,\bPhi,\bPsi\right)\Leftrightarrow \vect\left(\bY\right) \sim \pazocal{N}_{pr}\left(\vect\left(\bM\right),\bPsi \otimes \bPhi \right),
\label{eq:vectorize}
\end{equation}
where $\pazocal{N}_{pr}\left(\cdot\right)$ denotes the $pr$-variate normal distribution, $\vect\left(\cdot\right)$ is the vectorization operator and $\otimes$ denotes the Kronecker product.
However, a MN distribution has the desirable feature of simultaneously modeling and identify the between and within-variables variabilities as well as reducing the number of free parameters from $pr\left(pr+1\right)/2$ to $\left[r\left(r+1\right)/2\right]+\left[p\left(p+1\right)/2\right]-1$.

\subsection{The matrix-variate regression model}
\label{subsec:MatReg}

Let $\bY \in \real^{p \times r}$ be a continuous random matrix of dimension $p \times r$, containing $p$ responses measured in $r$ occasions.
Suppose we observe a set of $q$ covariates for each occasion, inserted in a matrix $\bX$ of dimension $q \times r$.
A generic matrix-variate regression model for $\bY$ has the form
\begin{equation}
\bY = \bbeta\bw^{\top} + \bB \bX + \bU ,
\label{eq:regModel1}
\end{equation}
where $\bbeta$ is the $p \times 1 $ vector consisting in the parameters related with the intercept, $\bw$ is a $r \times 1$ vector of ones, $\bB$ is the $p \times q$ matrix containing the parameters related to the $q$ covariates and $\bU$
is the $p \times r$ error term matrix.
Model~\eqref{eq:regModel1} can be expressed in compact notation as 
\begin{equation}
\bY = \bB^{*} \bX^{*} + \bU ,
\label{eq:regModel2}
\end{equation}
where $\bB^{*}$ is the $p \times \left(1+q\right)$ matrix involving all the parameters to be estimated and $\bX^{*}$ is the $\left(1+q\right) \times r$ matrix containing the information about the intercept and $q$ covariates \citep{viroli2012matrix}.
If we assume $\bU \sim \pazocal{N}_{p \times r}\left(\b0,\bPhi,\bPsi\right)$, then $\bY|\bX^{*} \sim \pazocal{N}_{p \times r}\left(\bB^{*} \bX^{*},\bPhi,\bPsi\right)$.
Therefore, a matrix-variate regression can be viewed as an encompassing framework containing as special cases the multivariate-multiple regression, when $r = 1$, and the univariate-multiple regression when $r = 1$ and $p = 1$.

\section{Methodology}
\label{sec:model}

\subsection{The matrix normal CWM}
\label{subsec:MatCWM}

Let $\left(\bX,\bY\right)$ be a pair of random matrices, defined as in Section~\ref{subsec:MatReg}, with joint distribution $p\left(\bX,\bY\right)$.
Then, a general matrix CWM has the following joint distribution
\begin{equation}
p\left(\bX,\bY\right) = \sum\limits_{g=1}^{G} p_{g}\left(\bY|\bX^{*}\right) p_{g}\left(\bX\right) \pi_{g},
\label{eq:CWM}
\end{equation}
where $p_{g}\left(\bY|\bX^{*}\right)$ is the cluster-specific conditional distribution of the responses, $p_{g}\left(\bX\right)$ is the cluster-specific marginal distribution of the covariates, and $\pi_{g}$ is the mixing weight (with $\pi_{g} > 0 $ and $\sum\limits_{g=1}^{G} \pi_{g} =1$).
Furthermore, we assume that in each group the conditional expectation $\E\left(\bY|\bX^{*}\right)$ is a linear function of $\bX^{*}$ depending on some parameters. 

In this paper, we focus on model~\eqref{eq:CWM} by assuming that both $p_{g}\left(\bY|\bX^{*}\right)$ and $p_{g}\left(\bX\right)$ are matrix normal densities, and $\E\left(\bY|\bX^{*}\right) = \bB^{*} \bX^{*}$, as described in Section~\ref{subsec:MatReg}.
Thus, model~\eqref{eq:CWM} can be rewritten as
\begin{equation}
\resizebox{.94\textwidth}{!}{
p\left(\bX,\bY;\bTheta\right) = \sum\limits_{g=1}^{G} \bphi_{p \times r}\left(\bY|\bX^{*};\bB_{g}^{*}\;\bX^{*},\bPhi_{\bY_{g}},\bPsi_{\bY_{g}}\right) \bphi_{q \times r}\left(\bX;\bM_{g},\bPhi_{\bX_{g}},\bPsi_{\bX_{g}}\right) \pi_{g},
$}
\label{eq:NOR_CWM}
\end{equation}
where $\bTheta$ denotes the set of all parameters.
For ease of understanding, a subscript with the variable name is added to the respective covariance matrices.
Notice that there is an identifiability issue concerning the covariance matrices.
Indeed, for a MN distribution, $\bPsi \otimes \bPhi = \bPsi^{*} \otimes \bPhi^{*}$ if $\bPhi^{*}= a\bPhi$ and $\bPsi^{*}= a^{-1}\bPsi$.
As a result, matrices $\bPhi$ and $\bPsi$ are identifiable up to a multiplicative constant $a$ \citep{melnykov2019studying}.
According to \citet{gallaugher17,gallaugher2018finite}, a way to obtain a unique solution is to fix the first diagonal element of the row covariance matrix to 1.
Therefore, in model~\eqref{eq:NOR_CWM} we adopt this approach by setting the first diagonal element of $\bPhi_{\bY_{g}}$ and $\bPhi_{\bX_{g}}$ to $1$.

If the MN-CWM was not available, a possible approach would be to vectorize the matrices and consider the MMN-CWM, of which it is a special case; see Equation~\eqref{eq:vectorize}.
However, such a procedure leads to two main concerns.
The first one is the overparameterization of the vectorized model.
Secondly, this increased number of parameters may have direct effects on model selection, as will be better analyzed in Section~\ref{subsec:simu4}.

\subsection{Parameter estimation}
\label{sec:EM}

Parameter estimation is carried out via the expectation-conditional maximization (ECM) algorithm \citep{meng1997algorithm}.
The only difference with respect to the expectation-maximization (EM) algorithm \citep{dempster1977maximum} is that the M-step is replaced by a sequence of simpler and computationally convenient CM-steps.
The EM algorithm cannot be directly implemented since there is no closed form solution for the covariance matrices of the MN distribution, i.e.~one of the two depends on the value of the other at the previous iteration \citep{dutilleul1999mle}.

Let $\Sa = \left\{\left(\bX_{i},\bY_{i}\right)\right\}_{i=1}^{N}$ be a sample of $N$ independent observations from model~\eqref{eq:NOR_CWM}.
Then, the incomplete-data likelihood function is
\begin{equation}
\resizebox{.94\textwidth}{!}{
  \begin{aligned}
	L\left(\bTheta|\Sa\right) = & \prod_{i=1}^{N} p\left(\bX_{i},\bY_{i};\bTheta\right)\\
                            = & \prod_{i=1}^{N}\left[\sum\limits_{g=1}^{G} \bphi_{p \times r}\left(\bY_{i}|\bX_{i}^{*};\bB_{g}^{*}\;\bX_{i}^{*},\bPhi_{\bY_{g}},\bPsi_{\bY_{g}}\right) \bphi_{q \times r}\left(\bX_{i};\bM_{g},\bPhi_{\bX_{g}},\bPsi_{\bX_{g}}\right) \pi_{g}\right].
\end{aligned}
}										
\label{eq:lik}
\end{equation}
Within the formulation of mixture models, $\Sa$ is viewed as being incomplete because, for each observation, we do not know its component membership.
Let $\bz_{i} = \left( z_{i1},\ldots,z_{iG}\right)^{\top}$ be the component membership vector such that $z_{ig}= 1$ if $\left(\bX_{i},\bY_{i}\right)$ comes from group $g$ and $z_{ig}= 0$ otherwise.
Now, the complete-data are $\Sa_{c} = \left\{\left(\bX_{i},\bY_{i},\bz_{i}\right)\right\}_{i=1}^{N}$, and the complete-data likelihood is
\begin{equation}
\resizebox{.94\textwidth}{!}{
\begin{aligned}
L_{c}\left(\bTheta|\Sa_{c}\right) = \prod_{i=1}^{N} \prod_{g=1}^{G} \left[\bphi_{p \times r}\left(\bY_{i}|\bX_{i}^{*}; \bB_{g}^{*}\;\bX_{i}^{*},\bPhi_{\bY_{g}},\bPsi_{\bY_{g}}\right) \bphi_{q \times r}\left(\bX_{i};\bM_{g},\bPhi_{\bX_{g}},\bPsi_{\bX_{g}}\right) \pi_{g}\right]^{z_{ig}}.
\end{aligned}
}
\label{eq:likcomp}
\end{equation}
Therefore, the corresponding complete-data log-likelihood can be written as
\begin{align}
l_{c}\left(\bTheta|\Sa_{c}\right) = & \sum\limits_{i=1}^{N} \sum\limits_{g=1}^{G} z_{ig} \ln\left(\pi_{g}\right) + \sum\limits_{i=1}^{N} \sum\limits_{g=1}^{G} z_{ig} \ln\left[\bphi_{p \times r}\left(\bY_{i}|\bX_{i}^{*};\bB_{g}^{*}\;\bX_{i}^{*},\bPhi_{\bY_{g}},\bPsi_{\bY_{g}}\right)\right]\nonumber\\ 
& + \sum\limits_{i=1}^{N} \sum\limits_{g=1}^{G} z_{ig} \ln\left[\bphi_{q \times r}\left(\bX_{i};\bM_{g},\bPhi_{\bX_{g}},\bPsi_{\bX_{g}}\right)\right].\label{eq:llk}
\end{align
In the following, by adopting the notation used in \citet{tomarchio2020two}, the parameters marked with one dot correspond to the updates at the previous iteration and those marked with
two dots represent the updates at the current iteration.

\paragraph{E-Step} 
\label{subsec:E-step}

The E-step requires calculation of the conditional expectation of \eqref{eq:llk}, given the observed data and the current estimate of the parameters $\dot{\bTheta}$.
To do this, we need to calculate 
\begin{align}
\ddot{z}_{ig} & = \E_{\dot{\bTheta}} \left[Z_{ig}|\bX_{i},\bY_{i}\right]\nonumber\\
                          &= \frac{\dot{\pi}_{g} \bphi_{p \times r}(\bY_{i}|\bX_{i}^{*};\dot{\bB}_{g}^{*}\bX_{i}^{*},\dot{\bPhi}_{\bY_{g}},\dot{\bPsi}_{\bY_{g}}) \bphi_{q \times r}(\bX_{i};\dot{\bM}_{g},\dot{\bPhi}_{\bX_{g}},\dot{\bPsi}_{\bX_{g}})}  {\sum\limits_{j=1}^{G} \dot{\pi}_{j} \bphi_{p \times r}(\bY_{i}|\bX_{i}^{*};\dot{\bB}_{j}^{*}\bX_{i}^{*},\dot{\bPhi}_{\bY_{j}},\dot{\bPsi}_{\bY_{j}}) \bphi_{q \times r}(\bX_{i};\dot{\bM}_{j},\dot{\bPhi}_{\bX_{j}},\dot{\bPsi}_{\bX_{j}})}, \label{eq:E-step}
\end{align}
which corresponds to the posterior probability that the unlabeled observation $\left(\bX_{i},\bY_{i}\right)$ belongs to the $g$th component of the mixture.

\paragraph{CM-Step 1}
\label{subsec:cm1}

At the first CM-step, we maximize the expectation of the complete-data log-likelihood with respect to $\bTheta_{1}=\left\{\pi_g,\bM_{g},\bPhi_{\bX_{g}},\bB_{g}^{*},\bPhi_{\bY_{g}}\right\}_{g=1}^{G}$, fixing $\bTheta_{2}=\left\{\bPsi_{\bX_{g}},\bPsi_{\bY_{g}}\right\}_{g=1}^{G}$ at $\dot\bTheta_{2}$.
In particular, we obtain
\begin{equation}
\ddot{\pi}_g   = \frac{\sum_{i=1}^N \ddot{z}_{ig}}{N}, 
\label{eq:cwm_cm1.1} 
\end{equation}
\begin{equation}
\ddot{\bM}_g    = \frac{\sum_{i=1}^N \ddot{z}_{ig}\bX_i}{\sum_{i=1}^N \ddot{z}_{ig}}, 
\label{eq:cwm_cm1.2} 
\end{equation}
\vspace{-8mm}
\begin{align}
\ddot{\bB}_{g}^{*} & = \left[\sum_{i=1}^N \ddot{z}_{ig} \bY_{i} (\dot{\bPsi}_{\bY_{g}})^{-1} \bX_{i}^{*\top}\right] \left[\sum_{i=1}^N \ddot{z}_{ig} \bX_{i}^{*} (\dot{\bPsi}_{\bY_{g}})^{-1} \bX_{i}^{*\top}\right]^{-1}, \label{eq:cwm_cm1.3} \\
\ddot{\bPhi}_{\bX_{g}} & = \frac{\sum_{i=1}^N \ddot{z}_{ig}\left(\bX_i-\ddot{\bM}_g\right)\left(\dot\bPsi_{\bX_{g}}\right)^{-1}\left(\bX_i-\ddot{\bM}_g\right)^{\top}}{r \sum_{i=1}^N \ddot{z}_{ig}}, \label{eq:cwm_cm1.4} \\
\ddot{\bPhi}_{\bY_{g}} & = \frac{\sum_{i=1}^N \ddot{z}_{ig}\left(\bY_i-\ddot{\bB}^{*}_{g}\bX_{i}^{*}\right)\left(\dot\bPsi_{\bY_{g}}\right)^{-1}\left(\bY_i-\ddot{\bB}^{*}_{g}\bX_{i}^{*}\right)^{\top}}{r \sum_{i=1}^N \ddot{z}_{ig}}. \label{eq:cwm_cm1.5} 
\end{align}

\paragraph{CM-Step 2}

At the second CM-step, we maximize the expectation of the complete-data log-likelihood with respect to $\bTheta_{2}$, keeping fixed $\bTheta_{1}$ at $\ddot\bTheta_{1}$.
Therefore, we have
\begin{align}
\ddot{\bPsi}_{\bX_{g}} & = \frac{\sum_{i=1}^N \ddot{z}_{ig}\left(\bX_i-\ddot{\bM}_g\right)^{\top}\left(\ddot\bPhi_{\bX_{g}}\right)^{-1}\left(\bX_i-\ddot{\bM}_g\right)}{q \sum_{i=1}^N \ddot{z}_{ig}},\label{eq:cwm_cm2.1} \\
\ddot{\bPsi}_{\bY_{g}} & = \frac{\sum_{i=1}^N \ddot{z}_{ig}\left(\bY_i-\ddot{\bB}^{*}_{g}\bX_{i}^{*}\right)^{\top}\left(\ddot\bPhi_{\bY_{g}}\right)^{-1}\left(\bY_i-\ddot{\bB}^{*}_{g}\bX_{i}^{*}\right)}{p \sum_{i=1}^N \ddot{z}_{ig}}.\label{eq:cwm_cm2.2} 
\end{align}

\section{Computational and operational aspects}
\label{sec:comp}

The code for the ECM algorithm previously described, and for the analysis in the following sections, is written within the R computing environment \citep{R_software}.

\subsection{ECM initialization}
\label{subsec:init}

The choice of the starting values is an important aspect for EM-based algorithms (see, e.g., \citealp{maitra2010simulating, michael2016effective}).
A common way to start an EM-based algorithm consists in providing initial values for the $z_{ig}$ in~\eqref{eq:E-step}, during the first E-step of the algorithm \citep{mclachlan2000finite}.
In our case, we need to provide an initial value also for $\bPsi_{\bX_{g}}$ and $\bPsi_{\bY_{g}}$ in~\eqref{eq:cwm_cm1.4} and~\eqref{eq:cwm_cm1.5}, respectively, during the first CM-step 1 of the algorithm.
Therefore, the following initialization strategy is implemented:
\begin{enumerate}
\item Generate $G$ random positive definite matrices for both $\bPsi_{\bX_{g}}$ and $\bPsi_{\bY_{g}}$.
This is done via the \texttt{genPositiveDefMat()} function of the \textbf{clusterGeneration} package, by using the ``eigen'' method. For further details see \citet{clusterG}. 
\item generate $N$ random vectors $z_{i}=\left(z_{i1},\ldots,z_{iG}\right)^{\top}$, $i=1,\ldots,N$.
This is done by using the following three approaches:
\begin{enumerate}[label*=\arabic*.]
    \item in a ``soft'' way, by generating $G$ positive random values from a uniform distribution on [0,1] for each observation, that are subsequently normalized to have a unitary sum.
		Being purely random, this procedure is repeated 15 times, and the solution maximizing the observed-data log-likelihood among these runs is considered;
		\item in a ``hard'' way, by using the classification produced by the $k$-means algorithm on the vectorized and merged data. 
		Specifically, after computing $\left\{\vect\left(\bX_{i}\right)\right\}_{i=1}^{N}$ and $\left\{\vect\left(\bY_{i}\right)\right\}_{i=1}^{N}$, the data are merged so that for each 
		observation we have a vector of dimension $\left(pr+qr\right) \times 1$;
		\item in a ``hard'' way, by using the classification produced by mixtures of matrix-normal distributions, computed on the merged data.
		In detail, for each observation we have a $\left(p+q\right) \times r$ matrix.
  \end{enumerate}
\end{enumerate}
The approach producing the highest observed-data log-likelihood is finally selected.

\subsection{Spurious clusters}
\label{subsec:spur}

A well-known issue in mixture models is related to the possibility for EM-based algorithms to converge to spurious local maximizers.
This is not a failing of these algorithms, and such solutions reflect a random pattern in the data rather than an underlying group structure.
They can be usually detected by the presence of data groups with very few observations or small eigenvalues compared to the other ones.
Consequently, these solutions often have a high likelihood, but are of little practical use or real-world interpretation \citep{mclachlan2000finite}.
Various approaches for mitigate this problem have been proposed in the literature.
For example, \citet{leisch2004flexmix} and \citet{dang2015families} impose a minimum size for the clusters of $\pi_{g}=0.05$.
\citet{melnykov2019studying} fought with this issue by excluding all the solutions that involved clusters consisting of less than 5 points. 
We faced this problem by removing the solutions that included clusters with estimated ${\pi}_{g}\approx 0.05$ or less, and with eigenvalues close to zero, in one or more of its estimated covariance matrices.

\subsection{Model selection and performance evaluation}
\label{subsec:classperf}

It is often the case that the number of groups $G$ is not known in advance, and model selection criteria are commonly used for estimating it.
Among them, the BIC is one of the most popular, and will be used in the following.
Furthermore, when the true classification of the data is known, the adjusted Rand index (ARI; \citealp{hubert1985comparing}) can be adopted to evaluate the clustering performance of a model.
Specifically, the ARI calculates the agreement between the true classification and the one predicted by the model.
An ARI of 1 indicates perfect agreement between the two partitions, whereas the expected value of the ARI under a random classification is 0.
It will be used in the manuscript along with the misclassification percentage rate $\eta$, which is the percentage of units classified in the wrong classes.

\section{Simulation studies}
\label{sec:artif}


\subsection{Simulation 1: A focus on the matrix-normal CWM}
\label{subsec:simu1}

In this study, several aspects related to our model are analyzed.
First of all, since the ECM algorithm is used to fit the model, it is desirable to evaluate its parameter recovery, i.e.~whether it can recover the generating parameters accurately. 
For this reason, data are generated from a four-component MN-CWM with $p=q=r=3$.
Two scenarios are then evaluated, according to the different level of overlap of the mixture components.
In the first scenario (labeled as ``Scenario A$_{1}$''), the mixture components are well-separated both in $\bX$, by assuming relatively distant mean matrices, and in $\bY|\bX^{*}$, by using different intercepts and slopes.
On the contrary, in the second scenario (labeled as ``Scenario B$_{1}$''), there is a certain amount of overlap, since the intercepts are all equal among the mixture components, while the slopes and the mean matrices assume approximately the same values among the mixture components.
The parameters used for Scenario A$_{1}$ are displayed in~