# Directed Clustering of Multivariate Data Based on Linear or Quadratic Latent Variable Models
## 3 Methodology  
### Modeling  
Consider multivariate data \(x_{i}\in\mathbb{R}^{m}\), \(m>1\), \(i=1,2,\ldots,n\). The model proposed by [16] considers the data generated by a latent variable plus Gaussian noise as follows:  
\[x_{i}=\alpha+\beta z_{i}+\varepsilon_{i}, \tag{1}\]  
where the univariate random effects \(z_{i}\) are realizations of a random variable \(Z\) with density function \(\phi\) (where no distributional assumption is made), \(\alpha\in\mathbb{R}^{m}\) and \(\beta\in\mathbb{R}^{m}\) are \(m\)-variate parameter vectors, and \(\varepsilon_{i}\in\mathbb{R}^{m}\) are independent Gaussian errors with a diagonal variance matrix \(\Sigma=\text{diag}(\sigma_{j}^{2})_{\{1\leq j\leq m\}}\in\mathbb{R}^{m\times m}\).  
This model is limited by its linear nature, particularly in handling multivariate data structures with the curvature. Given the existence of a non-linear structure in the motivating example data, it becomes apparent that a non-linear model is needed to capture the curvature effectively. In this paper, we propose a quadratic extension of the linear random effect model to allow the fitting of non-linear multivariate data. This extended model can be written as follows:  
\[x_{i}=\alpha+\beta z_{i}+\eta z_{i}{}^{2}+\varepsilon_{i}, \tag{2}\]  
where \(z_{i}\in\mathbb{R}\) is again a univariate random effect, which can be considered the curve parameterization of a quadratic curve through \(m\)-variate space, which is defined by \(m\)-variate parameters \(\alpha\), \(\beta\) and \(\eta\). As before, \(\varepsilon_{i}\sim N_{m}(0,\Sigma)\) are independent Gaussian errors. For a simplified presentation, we combine both models into a single formulation, as follows:  
\[x_{i}=g(z_{i},\theta)+\varepsilon_{i}, \tag{3}\]  
where \(g(z,\theta)=\alpha+\beta z\) with \(\theta=(\alpha^{T},\beta^{T})^{T}\), or \(g(z,\theta)=\alpha+\beta z+\eta z^{2}\) with \(\theta=(\alpha^{T},\beta^{T},\eta^{T})^{T}\) according to Equations (1) or (2), respectively, and proceed (as long as possible) with a unified presentation of the methodology.  
Under formulation (3), the conditional probability density function of \(x_{i}\) is a multivariate normal distribution, which can be expressed as follows:  
\[f(x_{i}|z_{i},\theta)=\] \[(2\pi)^{-m/2}|\Sigma|^{-1/2}\exp\biggl{\{}-\frac{1}{2}(x_{i}-g( z_{i},\theta))^{T}\Sigma^{-1}(x_{i}-g(z_{i},\theta))\biggr{\}}.\]  
To find the likelihood function, we need to first obtain the marginal probability density function \(f(x_{i}|\theta)\) for observations generated from model (3), which can be written as follows:  
\[f(x_{i}|\theta)=\int f(x_{i},z_{i}|\theta)dz_{i}=\int f(x_{i}|z_{i},\theta) \phi(z_{i})dz_{i},\]  
where \(f(x_{i},z_{i}|\theta)\) is the joint probability distribution of (multivariate, observed) data \(x_{i}\) and (univariate, unobserved) random effects \(z_{i}\). In principle, the distribution \(Z\) of the \(z_{i}\) can be Gaussian or non-Gaussian. If it is non-Gaussian, then one will generally not be able to find an analytical form for the marginal density function. Following the approach for the linear latent variable in [16], we do not make any explicit assumptions regarding the distribution \(Z\) of the \(z_{i}\). Instead, we use the nonparametric maximum likelihood approach for mixture models, intensely discussed by [PERSON] et al. [27] and based on theoretical results leading back to [28], in order to deal with the integration over \(z_{i}\). Specifically, we replace the integral with a set of mass points \(z_{1},z_{2},\ldots,z_{k}\) and their corresponding masses \(\pi_{1},\pi_{2},\ldots,\pi_{k}\), where \(k=1,2,\ldots,K\), so that the approximated marginal probability density function can be written as follows:  
\[f(x_{i}|\theta)\approx\sum_{k=1}^{K}f(x_{i}|z_{k},\theta)\pi_{k},\]  
which is a mixture of Gaussian components, i.e.,  
\[x_{i}|z_{k},\theta\sim N_{m}(g(z_{k},\theta),\Sigma), \tag{4}\]  
and mixture probabilities \(\pi_{k},k=1,\ldots,K\).