Response to ‘Letter to the Editor: on the stability and internal consistency of component-wise sparse mixture regression based clustering’, Zhang et al.

We have recently published a clustering method, namely CSMR, based on mixture regression in the context of high-dimensional predictors [1]. The motivation of our work was to provide a scalable method to deal with high-dimensional molecular features when there exist heterogeneous relationships between the molecular features and a disease phenotype of interest. In our original work, we simulated different data environments, including sample size , number of cluster-differentiating and cluster-specific predictors , number of clusters , and noise level . We evaluated CSMR on the simulation datasets, based on its accuracy of clustering using Rand index, and feature selection using true positive/negative rate, as well as the consistency of predicted and observed response values using Pearson correlation. In a real-world Cancer Cell Line Encyclopedia (CCLE) dataset [2], we evaluated CSMR based on the consistency of predicted and observed response values using cross-validation.

The letter by Zhang et al. raised the concerns that the performance of CSMR drops significantly when the number of cluster-differentiating predictors increases from 5 to 20, as demonstrated by evaluation metrics including the adjusted Rand index (ARI) and the clustering internal consistency (IC), both of which are advocated by Zhang et al. to be used for evaluating clustering methods. We appreciate Zhang et al.’s interests in our method and agree that their observations of the performance drop in the large case are indeed true. We would like to mention that mixture regression models, with the presence of both low- and high-dimensional features, rely on Expectation–Maximization (EM) algorithm or its variants for a solution. However, EM algorithm often converges to the maximum likelihood estimate of the mixture parameters locally [3], and using EM algorithm to solve mixture regression problem might yield poor clustering performance if the parameters are not initialized properly, and finding a good initial value is always a challenge. In other words, EM algorithm is highly sensitive to initialization, and different initial values may lead to different solutions and hence the instability. Many works have been published on establishing the convergence and stability criterion for EM algorithm [4–7]. For example, a guarantee for the EM algorithm to converge to the unique global optimum is when the likelihood is unimodal with certain regularity conditions [4]. However, the likelihood function could be multi-modal in many cases, for which only local optimum could be guaranteed; and what’s worse, sometimes the local optimum is a poor one that is far away from any global optimum of the likelihood [7]. In all, while EM algorithm has enjoyed popular use in solving mixture model, there is still space for improved theoretical understanding.

Here, we must clarify that in our setup, the total number of features, , is equal to . Here is the number of features that have no power to differentiate the different clusters, called the noise features. In our simulation data, is usually much larger than , to mimic the high-dimensional feature space setting. Obviously, there are only features that truly define the clusters. When there is no noise features, i.e. , under the case of , which is the case the authors pointed out to be least stable, we have , which is relatively large. While there is no direct empirical or theoretical study on how the number of features might impact the stability of the EM algorithm, a recent study investigated the effect of sample size on the estimation efficiency of mixture regression model using EM algorithm [8]. Their simulation studies suggest that the stability of the results will be reduced with low class separation and small sample size. We argue that for a fixed sample size 400, and cluster number , is a relatively large feature space, making a sample size of 400 a relatively small sample size. Hence, according to Jaki et al. [8], the stability of the mixture regression results will be impacted. We believe this is the root reason for the instability of CSMR.

At the initialized clustering assignment, as well as at the end of each of the iteration of CSMR, after a noise feature filtering step, CSMR will call an external mixture regression solver to fit the ‘low-dimensional’ model to update the posterior probability. We observe that for , N = 400, when the number of predictors increases to 60, the low-dimensional version of mixture regression fails to stably produce clustering results, even for the same dataset, under the same parameter setting, as evidenced by the low pair-wise ARI for repeated runs of (low-d) mixture regression on the same dataset using the same parameter setting (Figure 1). Here, we used flexmix function from the ‘flexmix’ R package as the low-d mixture regression solver, as it is currently the most popular and widely used implementation of mixture regression. Such instability is caused by the existence of many local optima, and the clustering results differ due to different starting values. Figure 1 suggests that the clustering performance deteriorates drastically with the increase of the size of the feature space, . The increase of in some sense makes the sample size relatively smaller. On the CCLE real dataset, the number of molecular features is large, and they are strongly collinear as many genes perform similar or related functions. In addition, it is unclear how well-separated the hidden clusters are, which depends on the level of heterogeneity of the different cell lines. In this case, it is very likely that the feature filtering step in CSMR would retain a large number of features, making the low-d mixture regression model fail to perform stably, and hence resulting in the unstable clustering results as demonstrated by Zhang et al. And the multi-collinearity problem among the features makes the coefficient estimation efficiency even worse. We are working to improve the algorithm of CSMR by limiting the retained molecular features that is fed into the low-d mixture regression solver, such that the solver will be stable even with relatively larger number of predictor space. Currently, we have taken a simple but effective approach, by setting a limit on the maximum number of features entering into the low-d mixture model to be 40 or , whichever is smaller. The retained features are selected to be the ones with top Pearson correlation strengths with the response variable. While the Pearson correlation seems to go against the mixture relations, it is only needed at the initialization step. The users could certainly set up different criterion for feature filtering. Overall, with this simple step, we will not allow a large number of features to enter into the low-d mixture model, avoiding the possible instability of the low-d mixture regression solver. As shown in Figure 2, the clustering performance significantly improved with such a feature number limiting step. However, we must argue that the instability issue may persist with increased (cluster-differentiating) predictor size, which violates the sparsity assumption of CSMR.

In summary, the instability issue of CSMR is caused by the instability of the EM-algorithm-based low-d mixture regression solver when facing with relatively large number of predictors, which is common to all EM algorithm-based mixture regression solver due to the existence of many local optima. In application of the CSMR algorithm, one underlying assumption is the sparsity of the differentiating predictors. Clearly, 60 true predictors with respect to sample size of 400 are not very sparse. However, in real-world biomedical data, no one could foresee how many differentiating genomic predictors there is, that is hidden in the large amount of feature space. It warrants further investigation using more comprehensive simulation data settings to explore the impact of multicollinearity, cluster separation and predictor dimensions on the performance of EM algorithm in solving mixture regression model.