Performance Evaluation of Missing-Value Imputation Clustering Based on a Multivariate Gaussian Mixture Model.

BACKGROUND: It is challenging to deal with mixture models when missing values occur in clustering datasets. METHODS AND
RESULTS: We propose a dynamic clustering algorithm based on a multivariate Gaussian mixture model that efficiently imputes missing values to generate a "pseudo-complete" dataset. Parameters from different clusters and missing values are estimated according to the maximum likelihood implemented with an expectation-maximization algorithm, and multivariate individuals are clustered with Bayesian posterior probability. A simulation showed that our proposed method has a fast convergence speed and it accurately estimates missing values. Our proposed algorithm was further validated with Fisher's Iris dataset, the Yeast Cell-cycle Gene-expression dataset, and the CIFAR-10 images dataset. The results indicate that our algorithm offers highly accurate clustering, comparable to that using a complete dataset without missing values. Furthermore, our algorithm resulted in a lower misjudgment rate than both clustering algorithms with missing data deleted and with missing-value imputation by mean replacement.
CONCLUSION: We demonstrate that our missing-value imputation clustering algorithm is feasible and superior to both of these other clustering algorithms in certain situations.

Introduction

Clustering analysis, as a multivariate statistical method, refers to the process of classifying a set of observations into subsets, called clusters, such that observations in the same cluster are similar in certain respects [1-3]. Clustering is widely used in medical sciences, for instance when clustering diseases or gene-expression profiles. Clustering methods usually fall into two categories: hierarchical clustering methods [4], which are used for clustering datasets with small size [5,6]; and dynamic clustering methods, such as K-means [7,8] and self-organizing maps [9], which begin with an initial partitioning of the individuals and iteratively move individuals from one cluster to another until the criterion of convergence is met. With dynamic clustering, the number of clusters must be specified in advance [10]. Dynamic algorithms are mostly heuristically motivated, and they do not require an underlying statistical model. Nevertheless, selection of the “correct” number of clusters and of the best clustering method remains a topic for discussion. Model-based clustering methods [11-13] are a type of dynamic clustering based on the hypothesis that the whole dataset is a finite mixture of the same type of distribution with different sets of parameters, such as a finite mixture of a multivariate Gaussian distribution. Compared with “heuristic” algorithms, one obvious advantage to model-based clustering is that objective statistical criteria—such as the Akaike information criterion (AIC) or the Bayesian information criterion (BIC) [14]—are used to determine the number of clusters.

Missing data present a problem for medical research, given that the data are often supplied retrospectively and from various sources [15]. In particular, missing data are a frequent occurrence in microarray experiments and in medical research. Missing values are especially common in large-scale studies involving dozens of variables and hundreds of individuals. Indeed, many clustering methods require a full set of data. Individuals with missing values are either rejected or been estimated prior to the analysis. Consequently, several missing-value imputation methods have been developed [16-20], such as mean substitution, regression imputation, fuzzy c-means (FCM) clustering of incomplete data [21], and Gaussian mixture model-based missing-value imputation classification [22]. In this study, we propose a dynamic method for a model-based missing-value imputation clustering algorithm. With our proposed method, missing values are estimated iteratively until arithmetic convergence is reached during the process of clustering individuals. We used 12 simulated datasets. Each of these datasets had two versions: a complete version and a version where 10% of the individuals had at least one variable missing. We used these datasets to evaluate the clustering accuracy and the accuracy and precision of cluster-parameter estimators. We compared our proposed algorithm based on the missing-value imputation according to the maximum likelihood estimation (a “pseudo-complete” dataset) with a model-based clustering algorithm using the complete dataset and another two model-based clustering algorithms, one using a dataset with missing values deleted and the other imputing the missing values by mean replacement. In addition, we compared our algorithm with the FCM clustering method using real datasets—viz., Fisher’s Iris dataset, the Yeast Cell-cycle Gene-expression dataset, and the CIFAR-10 images dataset.

Methods

Multivariate Gaussian mixture model

The dataset is arranged as an n × k matrix denoted by Y, where n is the number of individuals and k is the number of variables. Let y be the observed value of the i-th individual in the j-th variable, for i = 1,2,…, n and j = 1,2,…, k. Let Y = (y,y,…,y)T be the i-th column in matrix YT—i.e., a k × 1 vector of the data for individual i under all variables. The value of y across all the k variables represents the expression level of the i-th individual. Under a finite multivariate Gaussian mixture model, each Y is assumed to follow a k-dimensional Gaussian mixture distribution. Mathematically, the mixture distribution for c clusters is as follows: where is the probability density function for the l-th k-dimensional Gaussian distribution with mean vector μ = (μ, μ,…, μ) and variance-covariance matrix , for l = 1,2,…,c. Moreover, P with , is the mixing proportion of cluster l. The mixing proportion p is defined as the proportion of individuals that belong to the l-th cluster. The joint log-likelihood function for n independent individual vectors is defined as follows:

Clustering algorithm for missing values

Most reports [23-25] consider only complete datasets (without missing values) for Gaussian mixture clustering. Indeed, it is challenging to estimate missing values accurately. An urgent problem concerns how to provide the optimal number of clusters c and how to cluster n individuals into the c clusters precisely. Therefore, we propose a missing-value imputation algorithm as follows.

Suppose that there are r missing values in Y. Thus, y, y,…, y are missing. The Y, μ, and S are then divided into where

Suppose further that Y is from the l-th k-dimensional Gaussian distribution and that Y is known. The conditional expectation of Y, which belongs to the l-th cluster, is derived as follows:

We provide the initial value for μ, S, and S. Then, based on the criterion of the minimum mean-square deviations from Eq (3), the conditional expectation E(Y|Y) can be calculated, representing the best predicted function for Y [26]. Under the mixture distribution of the individual vector Y, the conditional expectation for Y based on Y is where E(Y|Y) denotes the estimators for the missing values Y. Subsequently, a “pseudo-complete” dataset can be constructed.

Iterative clustering algorithm

The proposed model-based clustering algorithm for datasets with missing values assigns each individual to one of the c clusters with a certain probability. We define this probability as p, which is that of the i-th individual belonging to the l-th cluster. An individual is assigned to the l-th cluster if p is greater than a certain pre-determined threshold. The probability can be calculated using the expectation-maximization (EM) algorithm [27] to achieve the maximum likelihood for the objective function derived with Eq (2). The EM algorithm begins with some initial parameters that are set in advance. Then, each parameter is iteratively updated in the algorithm until convergence to a minimizer is reached. An EM iteration includes the following steps:

Initialize the prior probabilities of the cluster assignment and the cluster parameters:

Calculate the missing values using Eqs (3) and (4) to generate the “pseudo-complete” dataset for .

Update the posterior probabilities of the cluster assignment:

Update the cluster proportions, mean vectors, and variance-covariance matrices:

Repeat steps (1)–(4) until convergence is reached.

The number of clusters c can also be treated as an unknown parameter and inferred by the BIC or AIC tests. The BIC is derived as follows: where q = c(k+1)(k+2)/2−1 is the number of independent parameters to be estimated in the model, is the likelihood of —i.e., the vector for the maximum likelihood estimation of the parameters—and n is the size of the dataset. The number of clusters c is determined by the maximum BIC value.

Results

Simulation analysis

A simulation was designed to evaluate the feasibility and accuracy of the proposed missing-values imputation algorithm. The individuals belonging to each cluster were known exactly and without subjective errors. Indeed, when analyzing real data, the misjudgment rate (MR) may reflect confounding errors in the experiment and the subjective errors. Without loss of generality, we simulated datasets with a two-dimensional Gaussian distribution and some missing values. All simulations were performed using statistical SAS (version 9.3; SAS/IML) software.

Design of the simulation

Suppose that there are 500 individuals derived from one of two (two-dimensional) Gaussian populations, denoted by MVN (μ1,S1) and MVN (μ2,S2), from which 10% of the individuals have at least one variable randomly removed. The entire dataset contains 1000 individuals. Each individual has two variables, and 100 individuals have at least one variable missing. The cluster’s mean vectors are denoted by μ1 and μ2, which were simulated at three different levels: A1: μ1 = (0, 0), μ2 = (2.5, 2.5), A2: μ1 = (0, 0), μ2 = (2.0, 2.0), and A3: μ1 = (0, 0), μ2 = (1.5, 1.5). The variance-covariance matrices are denoted by S1 and S2 (without loss of generality, supposing that S1 = S2 = S), and these were set at four levels: , , , and . The total number of treatment combinations (datasets) is therefore 3 × 4 = 12, that is, A1B1, A1B2, A1B3, A1B4, A2B1, A2B2, A2B3, A2B4, A3B1, A3B2, A3B3, and A3B4. Twenty replicated simulations were conducted for each of the twelve scenarios.

Our missing-value imputation cluster algorithm (denoted by M-3, generating a “pseudo-complete” dataset) was compared with three other cluster algorithms: a cluster algorithm using the complete dataset (denoted by M-1, i.e., the above simulation dataset without any missing values); a cluster algorithm using a dataset from which all individuals with missing variables are deleted (denoted by M-2, i.e., a cropped dataset); and a cluster algorithm using missing-value imputation by mean replacement in which the missing values are replaced by the mean value of the whole dataset (denoted by M-4, another “pseudo-complete” dataset). Thus, we used these four algorithms to cluster simulated datasets. These algorithms are all based on the multivariate Gaussian mixture model, using the maximum likelihood estimation with the EM algorithm. As such, the clustering results could be fairly compared. We evaluated M-1, M-2, M-3, and M-4 using the following metrics: (1) the average convergence rate; (2) the accuracy and precision of their respective parameter estimates; and (3) the total misjudgment rate (MR), where MR is the ratio of all misjudged individuals to the total number of individuals in the 20 replicates and MR = 1 - accuracy. A chi-square test and multiple comparisons of Scheffé’s Confidence Interval were used to test the differences in MR among the four algorithms.

Simulation results

The average parameter estimates include the mean vectors (± se), variance-covariance matrices (± se), the likelihood value of the maximum-likelihood function, and the total MR for the four clustering algorithms. The results are presented in Tables 1–3. Our proposed missing-value imputation clustering algorithm (M-3) outperformed the other two algorithms (M-2 and M-4) in terms of the convergence rate for most of the simulation datasets. Indeed, its superiority was most obvious when the mean vectors of two clusters were near each other and the variance-covariance of each cluster was more disperse, such as with datasets A2B4, A3B3, and A3B4. Moreover, the proposed algorithm resulted in the most accurate and precise parameter estimations, closely approximating the true values.

Table 1

Average parameter estimates under 4 different simulation datasets A1B1-A1B4 in 20 replicates.

Treatment		Iterative time	Likelihood value	Probability estimate		Mean vector estimate				Covariance matrix estimate				MR (%)
				P^_1	P^_2	μ^_1 ± se		μ^_2 ± se		S^_1 ± se		S^_2 ± se
A1B1	M-1	51	-1876.32	0.50	0.50	0.00±0.01	0.00±0.01	2.49±0.07	2.50±0.06	0.25±0.04	0.15±0.03	0.25±0.04	0.15±0.02	0.20a
										0.15±0.03	0.26±0.04	0.15±0.02	0.25±0.04
	M-2	58	-1864.94	0.50	0.50	0.00±0.01	-0.01±0.02	2.47±0.08	2.53±0.07	0.24±0.05	0.14±0.03	0.25±0.03	0.15±0.02	0.30a
										0.14±0.03	0.26±0.04	0.15±0.02	0.24±0.04
	M-3	64	-1897.05	0.50	0.50	0.00±0.02	0.00±0.01	2.49±0.07	2.50±0.06	0.24±0.04	0.14±0.02	0.24±0.04	0.14±0.04	0.20a
										0.14±0.02	0.26±0.03	0.14±0.04	0.25±0.04
	M-4	55	-1800.38	0.49	0.49	-0.01±0.03	0.00±0.02	2.46±0.08	2.46±0.09	0.22±0.07	0.13±0.04	0.23±0.06	0.14±0.05	1.50b
										0.13±0.04	0.23±0.06	0.14±0.05	0.21±0.07
A1B2	M-1	70	-2578.47	0.50	0.50	0.01±0.03	0.00±0.04	2.47±0.14	2.47±0.14	0.48±0.06	0.28±0.04	0.47±0.07	0.28±0.04	2.40a
										0.28±0.04	0.53±0.08	0.28±0.04	0.48±0.06
	M-2	80	-2362.69	0.50	0.50	-0.01±0.04	0.01±0.04	2.52±0.15	2.48±0.16	0.48±0.07	0.32±0.04	0.47±0.07	0.26±0.06	2.80a
										0.32±0.04	0.54±0.08	0.26±0.06	0.46±0.08
	M-3	87	-2591.61	0.50	0.50	-0.01±0.04	0.01±0.03	2.47±0.14	2.46±0.16	0.47±0.06	0.27±0.05	0.46±0.07	0.28±0.04	2.60a
										0.27±0.05	0.53±0.07	0.28±0.04	0.48±0.06
	M-4	93	-2653.93	0.48	0.48	0.02±0.09	0.01±0.08	2.43±0.26	2.44±0.23	0.42±0.14	0.25±0.11	0.44±0.12	0.25±0.11	5.64b
										0.25±0.11	0.43±0.13	0.25±0.11	0.44±0.12
A1B3	M-1	107	-2837.52	0.49	0.51	0.02±0.07	-0.01±0.08	2.47±0.21	2.47±0.21	0.68±0.16	0.42±0.10	0.71±0.14	0.43±0.09	5.90 ab
										0.42±0.10	0.78±0.15	0.43±0.09	0.75±0.14
	M-2	142	-2569.69	0.50	0.50	-0.02±0.06	0.01±0.06	2.48±0.20	2.52±0.19	0.73±0.16	0.47±0.08	0.69±0.16	0.40±0.09	5.45a
										0.47±0.08	0.81±0.15	0.40±0.09	0.67±0.17
	M-3	134	-2822.26	0.49	0.51	-0.01±0.07	-0.02±0.06	2.45±0.21	2.46±0.20	0.65±0.17	0.39±0.08	0.73±0.14	0.42±0.07	6.12ab
										0.39±0.08	0.76±0.12	0.42±0.07	0.68±0.15
	M-4	138	-2859.62	0.45	0.55	-0.05±0.15	0.04±0.12	2.40±0.37	2.43±0.32	0.63±0.21	0.40±0.16	0.70±0.19	0.39±0.11	8.94b
										0.40±0.16	0.60±0.20	0.39±0.11	0.71±0.19
A1B4	M-1	180	-3211.91	0.49	0.52	-0.05±0.10	-0.03±0.12	2.40±0.25	2.42±0.28	0.86±0.24	0.51±0.17	0.94±0.24	0.59±0.15	8.50a
										0.51±0.17	1.05±0.26	0.59±0.15	1.02±0.23
	M-2	214	-2852.78	0.48	0.52	-0.04±0.10	-0.06±0.14	2.38±0.30	2.37±0.32	0.84±0.25	0.55±0.16	0.97±0.19	0.56±0.16	9.20a
										0.55±0.16	1.05±0.26	0.56±0.16	0.94±0.25
	M-3	215	-3100.00	0.48	0.52	-0.05±0.12	-0.06±0.13	2.39±0.26	2.38±0.30	0.84±0.23	0.48±0.19	0.90±0.25	0.57±0.14	8.70a
										0.48±0.19	1.02±0.21	0.57±0.14	0.99±0.20
	M-4	267	-3289.32	0.43	0.57	-0.07±0.20	-0.08±0.25	2.17±0.48	2.22±0.52	0.77±0.32	0.43±0.28	0.82±0.27	0.48±0.15	18.90b
										0.43±0.28	1.09±0.23	0.48±0.15	1.07±0.25

M-1 indicates a complete data clustering algorithm; M-2 indicates a missing-data-deleted clustering algorithm; M-3 indicates our missing-value imputation clustering algorithm

M-4 indicates the clustering algorithm for missing-value imputation by mean replacement. These four algorithms are based on the multivariate Gaussian mixture model.

indicates the multiple comparisons of the differences in MR among the algorithms M-1, M-2, M-3, and M-4: having the different letters indicate that there is statistical significance between these two groups (P<0.05), and vice versa.

Table 3

Average parameter estimates under 12 different simulation datasets A3B1-A3B4 in 20 replicates.

Treatment		Iterative time	Likelihood value	Probability estimate		Mean vector estimate				Covariance matrix estimate				MR (%)
				P^_1	P^_2	μ^_1 ± se		μ^_2 ± se		S^_1 ± se		S^_2 ± se
A3B1	M-1	129	-1816.34	0.50	0.50	0.00±0.06	0.01±0.06	1.47±0.15	1.48±0.14	0.23±0.08	0.13±0.05	0.23±0.09	0.14±0.06	4.90a
										0.13±0.05	0.26±0.09	0.14±0.06	0.25±0.10
	M-2	337	-1637.81	0.51	0.49	0.01±0.07	0.01±0.06	1.53±0.15	1.47±0.16	0.23±0.08	0.14±0.06	0.23±0.10	0.12±0.06	5.35a
										0.14±0.06	0.27±0.10	0.12±0.06	0.22±0.09
	M-3	134	-1762.64	0.50	0.50	-0.01±0.07	0.00±0.06	1.46±0.16	1.46±0.16	0.21±0.08	0.13±0.05	0.21±0.09	0.14±0.05	5.00a
										0.13±0.05	0.26±0.08	0.14±0.05	0.25±0.08
	M-4	359	-1800.32	0.52	0.48	-0.02±0.09	0.03±0.08	1.44±0.25	1.45±0.26	0.20±0.14	0.13±0.07	0.21±0.13	0.14±0.06	9.20b
										0.13±0.07	0.24±0.13	0.14±0.06	0.25±0.13
A3B2	M-1	312	-2322.43	0.45	0.55	-0.06±0.10	-0.03±0.07	1.43±0.26	1.45±0.25	0.45±0.23	0.22±0.11	0.45±0.23	0.32±0.09	11.50a
										0.22±0.11	0.55±0.23	0.32±0.09	0.53±0.22
	M-2	444	-2102.31	0.45	0.55	-0.08±0.11	0.07±0.10	1.41±0.27	1.42±0.26	0.41±0.25	0.23±0.13	0.45±0.24	0.30±0.10	13.45ab
										0.23±0.13	0.47±0.25	0.30±0.10	0.54±0.23
	M-3	300	-2296.72	0.46	0.54	-0.06±0.10	-0.03±0.08	1.42±0.26	1.43±0.25	0.44±0.23	0.25±0.10	0.46±0.22	0.30±0.09	12.00a
										0.25±0.10	0.54±0.22	0.30±0.09	0.54±0.21
	M-4	489	-2358.99	0.43	0.57	-0.09±0.25	0.07±0.26	1.39±0.33	1.40±0.35	0.41±0.28	0.25±0.14	0.43±0.29	0.30±0.13	16.50 b
										0.25±0.14	0.54±0.29	0.30±0.13	0.52±0.28
A3B3	M-1	584	-2586.53	0.42	0.60	-0.15±0.28	-0.14±0.26	1.38±0.46	1.40±0.47	0.68±0.32	0.49±0.18	0.80±0.31	0.49±0.16	18.00a
										0.49±0.18	0.72±0.31	0.49±0.16	0.78±0.30
	M-2	603	-2364.52	0.39	0.63	-0.19±0.35	-0.18±0.33	1.38±0.52	1.37±0.50	0.64±0.35	0.48±0.17	0.81±0.32	0.46±0.18	23.30 b
										0.48±0.17	0.65±0.34	0.46±0.18	0.75±0.32
	M-3	310	-2578.77	0.44	0.58	-0.10±0.25	-0.10±0.26	1.40±0.45	1.39±0.45	0.65±0.31	0.44±0.17	0.74±0.29	0.47±0.15	18.00a
										0.44±0.17	0.78±0.32	0.47±0.15	0.78±0.29
	M-4	590	-2600.62	0.37	0.63	-0.18±0.38	-0.17±0.41	1.38±0.55	1.37±0.54	0.63±0.39	0.42±0.25	0.82±0.39	0.48±0.21	25.00b
										0.42±0.25	0.68±0.40	0.48±0.21	0.75±0.39
A3B4	M-1	645	-3013.37	0.35	0.65	-0.13±0.46	-0.17±0.38	1.30±0.72	1.30±0.70	0.92±0.38	0.50±0.20	1.14±0.40	0.68±0.24	28.00b
										0.50±0.20	1.08±0.34	0.68±0.24	0.97±0.39
	M-2	690	-2633.70	0.20	0.80	-0.22±0.74	-0.15±0.62	1.00±0.99	1.20±0.98	0.70±0.57	0.38±0.40	1.26±0.62	1.01±0.38	48.32c
										0.38±0.40	0.79±0.51	1.01±0.38	1.23±0.58
	M-3	520	-2892.21	0.40	0.60	-0.11±0.34	0.15±0.31	1.34±0.66	1.37±0.63	0.93±0.36	0.62±0.18	1.07±0.37	0.69±0.20	21.30a
										0.62±0.18	1.07±0.39	0.69±0.20	0.96±0.41
	M-4	835	-2901.33	0.10	0.90	-0.58±0.74	-0.55±0.80	0.90±1.04	0.97±1.07	0.61±0.72	0.21±0.58	1.26±0.77	0.89±0.47	56.00d
										0.21±0.58	0.62±0.74	0.89±0.47	1.23±0.77

M-1 indicates a complete data clustering algorithm; M-2 indicates a missing-data-deleted clustering algorithm; M-3 indicates our missing-value imputation clustering algorithm

M-4 indicates the clustering algorithm for missing-value imputation by mean replacement. These four algorithms are based on the multivariate Gaussian mixture model.

indicates the multiple comparisons of the differences in MR among the algorithms M-1, M-2, M-3, and M-4: having the different letters indicate that there is statistical significance between these two groups (P<0.05), and vice versa.

Thus, our proposed missing-value imputation clustering algorithm correctly clusters individuals with missing variables. The total MR resulting from the proposed algorithm (M-3) demonstrates that its performance is superior to that of M-4 in all simulation trials, with significant statistical difference (p < 0.05), and superior to that of M-2 in simulation trials A2B4, A3B3, and A3B4, with significant statistical difference (p < 0.05). Moreover, the clustering accuracy of our proposed algorithm was comparable to that of the complete dataset (M-1) (p > 0.05). In general, the more similar the mean vectors of two clusters or the more disperse their variance-covariance, the higher the MR of all clustering methods and the better the performance of our proposed algorithm compared to M-2 and M-4. Interestingly, however, our algorithm performed well in such cases. For example, in scenario A3B4, the speed of convergence with M-2 and M-4 was lower in both cases than for other scenarios and had a total MR of 48.32% and 56.00%, respectively, whereas the total MR with our proposed algorithm M-3 was 21.30%. In fact, this result was slightly better than M-1, which resulted in an MR of 28.00%.

Real data analysis

Iris dataset

Fisher’s Iris dataset (Fisher, 1936) [12] is perhaps the best-known database for comparing clustering algorithms. The dataset contains three clusters with 50 individuals each, where each cluster is a species of the genus Iris: namely, Iris setosa, Iris versicolor, and Iris virginica. These three species are thus segregated. The dataset contains four variables: sepal length, sepal width, petal length, and petal width. Fisher’s Iris dataset is available on the Internet (http://archive.ics.uci.edu/ml/datasets/Iris). We randomly removed one of the four variables in 10% of the individual observations of Iris data; the resulting dataset constituted a dataset with missing values. Table 4 lists the clustering results from five clustering algorithms: the above four clustering algorithms [the clustering algorithm with a complete dataset (M-1), the clustering algorithm in which individuals with missing variables are deleted (M-2), our proposed missing-value imputation algorithm (M-3), and the clustering algorithm for missing-value imputation by mean replacement (M-4)] and the FCM clustering method applied to the dataset with data missing from 10% of the individuals. In the M-4 and FCM methods, the missing values are replaced by the mean of the whole dataset. The superiority of our proposed method is obvious insofar as M-3 resulted in the lowest MR of the five algorithms. Because of Gaussian mixture distribution of the dataset, our algorithm makes full use of the known data based on the Gaussian mixture model to estimate missing values, the accuracy of our method was higher than those of the FCM method and the M-2 and M-4 algorithms, and it closely approximated the accuracy obtained using the complete-dataset clustering algorithm (M-1). Moreover, the model-based algorithms (M-1, M-2, M-3, and M-4) resulted in the highest BIC values in the three clusters, validating the effectiveness of using the BIC to determine the number of clusters.

Table 4

The comparison among four algorithms and FCM method for Fisher’s Iris dataset.

Method	n	Misjudgment rate (MR)	Iterative time	BIC value in three clusters
M-1	150	3 (2.0%)	37	-307.848
M-2	135	6 (4.4%)	43	-273.539
M-3	150	3 (2.0%)	36	-305.211
M-4	150	7 (4.7%)	49	-307.382
FCM	150	5 (3.3%)

M-1 indicates the complete data clustering algorithm; M-2 indicates missing data deleted clustering algorithm; M-3 indicates our missing-value imputation clustering algorithm; M-4 indicates clustering algorithm for missing-value imputation by mean replacement; FCM indicates Fuzzy C-Means clustering.

Yeast Cell-cycle dataset

The Yeast Cell-cycle dataset [28] shows the fluctuation in the expression levels of the whole genome (6220 genes) over 17 time points taken at 10-minute intervals, covering nearly two cell cycles. The entire raw dataset is available at http://genome-www.stanford.edu/cellcycle/. Among the effectively expressed genes, 416 genes were found to have consistent periodic changes in expression level. Yeung et al. picked out 384 genes over 15 time points that reached their peaks at five (C = 5) cell phases on the basis of the size of the buds, after which five non-intersecting clusters were classified [29]. These 384 gene expressions, a subset with a moderate sample size, are available at http://www.cs.washington.edu/homes/kayee/model. All 384 genes were assigned to one of the five clusters by the original researchers [28,29]. Therefore, we can use this dataset to test the performance of the clustering algorithms introduced here and compare them with the performance of existing methods. We randomly removed one of the 15 variable values in 10% of the genes to create a dataset with missing values. Table 5 lists the clustering results from five clustering algorithms: the above four model-based clustering algorithms (M-1, M-2, M-3, and M-4) and the FCM clustering method. Further validating the results from Fisher’s Iris dataset, we found that our missing-value imputation algorithm had the lowest MR. In addition, because our missing-value imputation based on the Gaussian mixture distribution is compatible with the Yeast Cell-cycle dataset, the accuracy of our method was significantly higher than that of the FCM method, and it closely approximated the accuracy from using the complete-dataset clustering algorithm (M-1).

Table 5

The comparison among four algorithms and FCM method for Yeast Cell-cycle dataset.

Method	n	Misjudgment rate (MR)	Iterative time	BIC value in five clusters
M-1	384	65 (16.9%)	1265	-6472.448
M-2	346	62 (17.9%)	1340	-5892.953
M-3	384	62 (16.1%)	1260	-6316.730
M-4	384	78 (20.3%)	1307	-6490.706
FCM	384	80 (20.8%)

M-1 indicates the complete data clustering algorithm; M-2 indicates missing data deleted clustering algorithm; M-3 indicates our missing-value imputation clustering algorithm; M-4 indicates the clustering algorithm for missing-value imputation by mean replacement; FCM indicates Fuzzy C-Means clustering.

CIFAR-10 images dataset

The CIFAR-10 dataset [30] contains 60,000 32 × 32 color images (3072 variables) divided into ten classes, which represent different kinds of pictures. The classes are designed to be completely mutually exclusive. For example, neither automobile pictures nor truck pictures contain images of pickup trucks. Having a large sample size, each class contains 5000 training images and 1000 test images. The entire dataset was used to test the performance of test images according to training images in supervised classification, and it is available at http://www.cs.toronto.edu/~kriz/cifar.html. We used the training images to test the performance of our unsupervised clustering algorithm and compared it to existing clustering algorithms. We randomly removed 10 of the 3072 variable values in 10% of the color images to create a dataset with missing values.

Table 6 lists the clustering results from five clustering algorithms: the above four model-based clustering algorithms (M-1, M-2, M-3, and M-4) and the FCM clustering method. We found that our missing-value imputation algorithm had a lower total MR than algorithms M-2 and M-4. The accuracy of our method closely approximated the accuracy obtained using the clustering algorithm with complete data (M-1) but was lower than that of the FCM method. (It should be noted that most clusters of the CIFAR-10 data have deviants from the normal distribution.) Thus, the performance of our model-based algorithm does not offer clustering that is more accurate than the FCM method. Moreover, the speed of computations was much lower in algorithms M-1, M-2, M-3, and M-4.

Table 6

The comparison among four algorithms and FCM method for CIFAR-10 dataset.

Method	MR
	total	Cluster 1	Cluster 2	Cluster 3	Cluster 4	Cluster 5	Cluster 6	Cluster 7	Cluster 8	Cluster 9	Cluster 10
M-1	39.9%	38.2%	41.4%	35.6%	43.5%	40.5%	37.9%	38.5%	42.1%	35.4%	45.6%
M-2	42.4%	42.6%	40.8%	35.2%	45.8%	42.0%	46.3%	40.4%	46.7%	38.4%	45.4%
M-3	39.8%	42.0%	41.0%	36.2%	39.6%	39.7%	38.9%	38.0%	43.2%	34.2%	45.0%
M-4	44.0%	45.3%	43.6%	37.4%	45.9%	43.7%	48.5%	42.8%	47.9%	41.3%	43.9%
FCM	35.9%	35.8%	35.6%	35.5%	36.7%	34.1%	37.6%	32.6%	37.6%	35.7%	37.8%

M-1 indicates the complete data clustering algorithm; M-2 indicates missing data deleted clustering algorithm; M-3 indicates our missing-value imputation clustering algorithm; M-4 indicates clustering algorithm for missing-value imputation by mean replacement; FCM indicates Fuzzy C-Means clustering.

Discussion

Missing values are unavoidable in experimental data and data from field investigations [31]. Many imputation methods [16-20] have been proposed, such as the mean-substitution method, the regression imputation method, the EM estimation method, and the FCM method. Clustering is a widely used statistical method, especially in gene expression profiles. Most classical clustering methods, such as K-means clustering and hierarchical clustering, are unable to deal with missing values. Therefore, we developed a model-based clustering algorithm for imputing missing values. The proposed method utilizes information from existing data, and it is based on the Gaussian mixture model. With our imputation algorithm, missing values are estimated by calculating the conditional expectation, from which a “pseudo-complete” dataset can be generated and used for clustering.

The performance of our proposed algorithm was tested and compared with other algorithms using both real and simulated data. We found that our proposed method outperformed other algorithms. In particular, our proposed method produced clustering results comparable to an algorithm that uses a complete dataset, and it outperformed an algorithm in which samples with missing data are removed and another algorithm in which missing values are imputed by mean replacement. The results are consistent with Hunt and Jorgensen’s reports [32]. Even with datasets in which the boundaries of the clusters overlap, such as datasets A3B3, A2B4, and A3B4, our algorithm performed well.

Furthermore, we found that the proposed algorithm outperformed the FCM clustering method on Fisher’s Iris dataset and on the Yeast Cell-cycle dataset. The FCM clustering method, a heuristic approach, has been verified as a successful missing-value clustering method [33]. FCM clusters data based on the characteristics of the samples. The criteria for FCM clustering are inconsistent, however, because the characteristics of dynamic algorithms lead to uncertainties about the samples, degrading FCM’s ability to precisely cluster data [34-36]. Our study shows that the proposed method is more likely to result in a low MR when the dataset is exactly a multivariable Gaussian mixture distribution. That is, it is more accurate than other methods with both simulated data and real data, e.g., Fisher’s Iris dataset and the Yeast Cell-cycle dataset. One possible reason for this is that our method estimates missing values with a conditional expectation based on the Gaussian mixture distribution and generates a “pseudo-complete” dataset. In addition, our method combines dynamic clustering with discriminant analysis and uses the EM algorithm to achieve a maximum-likelihood function. It can thus better utilize the information from known data to discover the intrinsic and implicit properties of a dataset. Usually, model-based methods are better than heuristic approaches, even though the former are always more time consuming. Model-based methods depend on a probability model, which is usually proposed based on the experience and feasibility of the model. The convergence and accuracy of the EM algorithm have already been established [25]. The Gaussian mixture model is both convenient and robust. Moreover, a model-based method provides substantial information—e.g., the posterior probabilities of individuals and the parameters of each cluster—and this information can be exploited to more accurately impute missing values with iterative updates to the algorithm until convergence to a minimizer is reach.

Our method uses the Bayesian posterior probability to cluster individuals; individuals are assigned to the cluster having the highest posterior probability. It is more reasonable to assign individuals to a cluster having a high posterior probability, such as 0.8 or 0.9, and this improves the clustering accuracy. Whereas traditional dynamic clustering methods cluster individuals directly to a certain category, this process is difficult when the boundaries of the clusters overlap. Another advantage to a model-based approach is that statistical criteria, such as AIC and BIC [14], can be used to discover the number of clusters that best suit the data. We confirmed the effectiveness of such criteria with results indicating three clusters for Fisher’s Iris dataset (a dataset with exactly three species of Iris) and five clusters for the Yeast Cell-cycle Gene-expression dataset. However, the reversible-jump Markov-chain Monte Carlo (MCMC) method was previously developed to infer the number of distributions in a mixture [37]. MCMC estimates the optimal number of clusters from a Bayesian perspective. Our future research shall compare MCMC with the BIC criteria used in this study. Finally, a prerequisite for raw data collected by experimenters is that they be normalized. Indeed, our method is more sensitive to data that have not been normalized. Therefore, either a log or some other form of data transformation should be applied in order to distribute the data in a normal fashion. Furthermore, because our proposal relies on the information from known individuals when imputing missing values, a large sample size is desirable for more accuracy when estimating missing values.

Conclusions

In this paper, we have presented a model-based missing-data imputation algorithm for clustering numerical data with missing values. We achieved a high convergence rate with our algorithm when it was tested on simulated and real datasets. In addition, the statistical power of our algorithm is high, and it can accurately estimate parameters and correctly cluster individuals with missing variables. Moreover, the results of our evaluation show a low MR, even when the mean vectors of two clusters are close to each other or when their variance-covariance is disperse. The algorithm was applied successfully to datasets (both simulated and real) having missing data. Its performance is superior to that of an algorithm in which samples with missing data were removed and another algorithm in which missing values were imputed by mean replacement. A comparative evaluation demonstrated that our proposed algorithm outperforms the FCM method when the data from each cluster fit a multivariate Gaussian distribution.

Table 2

Average parameter estimates under 4 different simulation datasets A2B1-A2B4 in 20 replicates.

Treatment		Iterative time	Likelihood value	Probability estimate		Mean vector estimate				Covariance matrix estimate				MR (%)
				P^_1	P^_2	μ^_1 ± se		μ^_2 ± se		S^_1 ± se		S^_2 ± se
A2B1	M-1	81	-1903.64	0.50	0.50	-0.01±0.04	0.00±0.03	2.00±0.12	2.00±0.13	0.25±0.05	0.15±0.04	0.24±0.04	0.14±0.04	1.64a
										0.15±0.04	0.25±0.05	0.14±0.04	0.24±0.05
	M-2	117	-1696.53	0.51	0.49	0.01±0.04	-0.01±0.04	2.00±0.14	1.99±0.13	0.24±0.05	0.15±0.03	0.25±0.05	0.14±0.04	1.86a
										0.15±0.03	0.26±0.04	0.14±0.04	0.24±0.06
	M-3	92	-1885.63	0.50	0.50	-0.01±0.04	0.01±0.04	1.99±0.14	1.99±0.12	0.23±0.05	0.14±0.04	0.23±0.06	0.15±0.04	1.60a
										0.14±0.04	0.26±0.04	0.15±0.04	0.24±0.05
	M-4	126	-2008.47	0.48	0.52	-0.02±0.07	0.01±0.06	1.95±0.19	1.96±0.20	0.22±0.10	0.12±0.07	0.21±0.10	0.11±0.05	2.26a
										0.12±0.07	0.27±0.12	0.11±0.05	0.27±0.07
A2B2	M-1	120	-2372.14	0.50	0.50	-0.02±0.07	-0.01±0.06	2.01±0.16	1.97±0.17	0.45±0.13	0.27±0.09	0.46±0.14	0.28±0.10	5.60a
										0.27±0.09	0.52±0.11	0.28±0.10	0.50±0.14
	M-2	179	-2128.49	0.51	0.49	-0.02±0.08	0.02±0.06	1.95±0.18	1.96±0.18	0.48±0.11	0.27±0.08	0.47±0.14	0.25±0.09	6.30ab
										0.27±0.08	0.53±0.13	0.25±0.09	0.45±0.13
	M-3	135	-2462.73	0.51	0.49	0.01±0.07	0.02±0.06	1.95±0.19	1.97±0.17	0.46±0.11	0.28±0.08	0.44±0.09	0.26±0.09	5.40a
										0.28±0.08	0.53±0.10	0.26±0.09	0.46±0.08
	M-4	266	-2410.50	0.45	0.55	-0.03±0.10	0.02±0.11	1.91±0.28	1.92±0.29	0.44±0.17	0.25±0.10	0.44±0.18	0.24±0.11	9.00 b
										0.25±0.10	0.54±0.16	0.24±0.11	0.43±0.19
A2B3	M-1	219	-2699.35	0.48	0.52	-0.05±0.08	-0.03±0.08	1.92±0.27	1.94±0.24	0.67±0.20	0.40±0.10	0.75±0.19	0.42±0.11	10.90a
										0.40±0.10	0.78±0.21	0.42±0.11	0.78±0.19
	M-2	387	-2489.38	0.47	0.53	-0.07±0.09	0.05±0.10	1.90±0.27	1.90±0.28	0.64±0.25	0.38±0.13	0.72±0.22	0.41±0.12	11.40a
										0.38±0.13	0.73±0.24	0.41±0.12	0.79±0.21
	M-3	242	-2709.34	0.47	0.53	-0.07±0.10	-0.05±0.10	1.90±0.28	1.90±0.30	0.66±0.18	0.38±0.09	0.76±0.20	0.47±0.11	11.00a
										0.42±0.09	0.77±0.18	0.47±0.11	0.75±0.20
	M-4	450	-2646.83	0.42	0.58	-0.14±0.20	0.12±0.18	1.89±0.41	1.88±0.45	0.64±0.29	0.38±0.17	0.72±0.27	0.40±0.15	16.96b
										0.38±0.17	0.78±0.27	0.40±0.15	0.77±0.23
A2B4	M-1	360	-3014.36	0.45	0.55	-0.11±0.16	0.09±0.17	1.85±0.34	1.80±0.32	0.86±0.26	0.52±0.18	1.01±0.23	0.63±0.17	12.80a
										0.52±0.18	1.04±0.25	0.63±0.17	1.05±0.23
	M-2	486	-2677.25	0.43	0.57	-0.15±0.25	-0.15±0.25	1.84±0.35	1.82±0.32	0.82±0.30	0.51±0.19	1.04±0.29	0.61±0.18	18.00b
										0.51±0.19	0.94±0.25	0.61±0.18	1.01±0.29
	M-3	342	-2907.18	0.44	0.56	-0.10±0.17	0.08±0.14	1.86±0.34	1.82±0.32	0.87±0.27	0.53±0.18	0.92±0.28	0.62±0.17	12.30a
										0.53±0.18	1.03±0.25	0.62±0.17	1.03±0.29
	M-4	440	-2977.46	0.40	0.60	-0.20±0.35	0.19±0.33	1.82±0.57	1.81±0.62	0.76±0.41	0.51±0.27	1.05±0.33	0.65±0.27	21.80c
										0.51±0.27	0.90±0.39	0.65±0.27	1.07±0.34

M-1 indicates a complete data clustering algorithm; M-2 indicates a missing-data-deleted clustering algorithm; M-3 indicates our missing-value imputation clustering algorithm

M-4 indicates the clustering algorithm for missing-value imputation by mean replacement. These four algorithms are based on the multivariate Gaussian mixture model.

indicates the multiple comparisons of the differences in MR among the algorithms M-1, M-2, M-3, and M-4: having the different letters indicate that there is statistical significance between these two groups (P<0.05), and vice versa.