# Directed Clustering of Multivariate Data Based on Linear or Quadratic Latent Variable Models
## 3 Methodology
### The M-Step for the Quadratic Model  
For model (2), taking partial derivatives of (6) w.r.t. \(\alpha\), \(\beta\), \(\eta\), and \(z_{k}\), leads to a complex set of expressions that we provide in Appendix A. However, these expressions involve many matrix inversions and are computationally not very stable to evaluate. Therefore, we follow the approach in [16], where, only within these expressions, we assume that \(\hat{\Sigma}\equiv\text{diag}(\sigma^{2})\), for some constant \(\sigma^{2}\), which does not need to be specified since it cancels out from the resulting simplified update equations for \(\alpha\), \(\beta\), and \(\eta\). The resulting simplified estimators are then given by the following:  
\[\hat{\alpha}=\frac{1}{n}\left(\sum_{i=1}^{n}x_{i}-\hat{\beta}\sum_{i=1}^{n} \sum_{k=1}^{K}w_{ik}\hat{z}_{k}-\eta\sum_{i=1}^{n}\sum_{k=1}^{K}w_{ik}\hat{z}_{ k}^{2}\right), \tag{7}\]  
\[\hat{\beta}=\frac{\sum_{i=1}^{n}\sum_{k=1}^{K}w_{ik}(x_{i}-\hat{\alpha}-\hat{ \eta}\hat{z}_{k}^{2})\hat{z}_{k}}{\sum_{i=1}^{n}\sum_{k=1}^{K}w_{ik}\hat{z}_{ k}^{2}}, \tag{8}\]  
\[\hat{\eta}=\frac{\sum_{i=1}^{n}\sum_{k=1}^{K}w_{ik}x_{i}\hat{z}_{k}^{2}-\hat{ \alpha}\sum_{i=1}^{n}\sum_{k=1}^{K}w_{ik}\hat{z}_{k}^{2}-\hat{\beta}\sum_{i=1 }^{n}\sum_{k=1}^{K}w_{ik}\hat{z}_{k}^{3}}{\sum_{i=1}^{n}\sum_{k=1}^{K}w_{ik} \hat{z}_{k}^{4}}. \tag{9}\]  
The estimates of \(z_{k}\) will be obtained by solving the cubic equation as follows:  
\[(\sum_{i=1}^{n}\sum_{j=1}^{m}2\cdot w_{ik}\hat{\eta}_{j}^{2}){z_{ k}}^{3}+3\cdot(\sum_{i=1}^{n}\sum_{j=1}^{m}w_{ik}\hat{\beta}_{j}\hat{\eta}_{j}){z_ {k}}^{2} \tag{10}\] \[-(\sum_{i=1}^{n}\sum_{j=1}^{m}2\cdot w_{ik}x_{ij}\hat{\eta}_{j}- \sum_{i=1}^{n}\sum_{j=1}^{m}2\cdot w_{ik}\hat{x}_{j}\hat{\eta}_{j}-\sum_{i=1 }^{n}\sum_{j=1}^{m}w_{ik}\hat{z}_{j}^{2})z_{k}\] \[-(\sum_{i=1}^{n}\sum_{j=1}^{m}w_{ik}x_{ij}\hat{\beta}_{j}-\sum_{ i=1}^{n}\sum_{j=1}^{m}w_{ik}\hat{x}_{i}\hat{\beta}_{j})=0.\]  
Unfortunately, there is still an added complication since, even when conditional on \(\hat{\alpha}\), \(\hat{\beta}\), and \(\hat{\eta}\), obtaining the analytical form of \(\hat{z}_{k}\) is not possible. However, the score equation (10) is a cubic equation of \(z_{k}\), for which solvers exist in standard software packages. We use the polyroot() function (available in the **base**[31] R package) to solve this cubic equation. It is important to note that the roots of a cubic equation fall into two scenarios: (i) one real root and two complex roots, in which case we consider the real root as the estimate of \(z_{k}\), and (ii) three real roots, where we always choose the smallest absolute root in the implementation.  
The estimates in Equations (7)-(10) still depend on each other, but they do not depend on \(\hat{\Sigma}\). Hence, within each M-step of the EM algorithm, one can iterate between them a small number of times, then estimate \(\pi_{k}\) via \(\hat{\pi}_{k}=\frac{1}{n}\sum_{i=1}^{n}w_{ik}\), and finally estimate the entries of \(\Sigma\) via the following:  
\[\hat{\sigma_{j}}^{2}=\frac{\sum_{i=1}^{n}\sum_{k=1}^{K}w_{ik}(x_{ij}-\hat{ \alpha}_{j}-\hat{\beta}_{j}\hat{z}_{k}-\hat{\eta}\hat{z}_{k}^{2})^{2}}{n}. \tag{11}\]  
The simulation study in the following section demonstrates that this procedure accurately estimates the model parameters. It is also noted that, as in [16], one can work with cluster-dependent variance matrices \(\Sigma_{k}\), \(k=1,\dots,K\), which vary in size and shape. However, in the context of directed clustering, the benefits of doing so are limited. Therefore, we relegate a brief consideration of cluster-dependent variance matrices to Appendix B but continue with the homoscedastic case in the main flow of this paper.