Bayesian Nonparametric Classification for Incomplete Data
With a High Missing Rate: an Application to Semiconductor Manufacturing Data

Algorithms for Gaussian mixture models on missing data have been studied in the last few decades. Ghahramani and Jordan (1994) used the Expectation-Maximization (EM) algorithm to find parameter values and missing components maximizing the likelihood of Gaussian mixture models. Zhang and Everson (2004) developed a Bayesian approach of mixture models using Gibbs sampler. This method utilizes full conditional distributions to obtain the joint posterior distribution of parameters and missing values. Williams et al. (2007) introduced variational inference based on the mean-field approximation for Bayesian mixture models. Both missing values and parameters are iteratively updated until the evidence lower bound (ELBO) converges.

We focus on the estimation of a Gaussian mixture model, permitting simultaneous inference of missing values through Gibbs sampling is a Markov Chain Monte Carlo (MCMC). Let $\textbf{x}_{i},i=1,\ldots,n$ be $n$ independent $p$-dimensional observations from the mixture distribution consists of $H$ Gaussian components. We partition an observation into two components $\textbf{x}_{i}=\{\textbf{x}_{i}^{o_{i}},\textbf{x}_{i}^{m_{i}}\}$, where $o_{i}\subset\{1,2,\ldots,p\}$ is an index set of observed variables and $m_{i}\subset\{1,2,\ldots,p\}$ is an index set of missing variables . That is, $\textbf{x}_{i}^{o_{i}}$ and $\textbf{x}_{i}^{m_{i}}$ indicate the observed components and missing components from the $i$th observation $\textbf{x}_{i}$, respectively. We can express the mixture distribution as follows:

$$f(\textbf{x}_{i})=\sum_{h=1}^{H}w_{h}\mathcal{N}_{p}(\textbf{x}_{i}|\boldsymbol{\mu}_{h},\boldsymbol{\Sigma}_{h})=\sum_{h=1}^{H}w_{h}\,\mathcal{N}_{p}\left(\begin{bmatrix}\textbf{x}_{i}^{o_{i}}\\ \textbf{x}_{i}^{m_{i}}\end{bmatrix}\Bigg{|}\begin{bmatrix}\boldsymbol{\mu}_{h}^{o_{i}}\\ \boldsymbol{\mu}_{h}^{m_{i}}\end{bmatrix},\begin{bmatrix}\Sigma_{h}^{o_{i},o_{i}}&\Sigma_{h}^{o_{i},m_{i}}\\ \Sigma_{h}^{m_{i},o_{i}}&\Sigma_{h}^{m_{i},m_{i}}\end{bmatrix}\right)$$

(1)

Here, $w_{h}$ is the non-negative mixing proportion and sum up to one, i.e., $\sum_{h=1}^{H}w_{h}=1$.

To implement the Gibbs sampling for mixture models with incomplete data, we need the full conditional posterior distributions of model parameters and missing values. In this paper we only cover the full conditional posterior for missing values and skip other parameters of mixture models; see the paper by Franzén (2006) for more details. If the data come from a multivariate normal distribution with mean $\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$, a full conditional density for missing values based on the observed data can derive easily the following equations (2) by using a special property of multivariate normal distribution.

	$\displaystyle f(\textbf{x}_{i}^{m_{i}}|\textbf{x}_{i}^{o_{i}},\text{others})\sim\mathcal{N}_{|m_{i}|}(\textbf{x}_{i}^{m_{i}};\boldsymbol{\mu}^{m_{i}|o_{i}},\boldsymbol{\Sigma}^{m_{i}|o_{i}}),$		(2)
	$\displaystyle\boldsymbol{\mu}^{m_{i}|o_{i}}=\boldsymbol{\mu}^{m_{i}}+\boldsymbol{\Sigma}^{m_{i},o_{i}}(\boldsymbol{\Sigma}^{o_{i},o_{i}})^{-1}(\textbf{x}_{i}^{o_{i}}-\boldsymbol{\mu}^{o_{i}}),$
	$\displaystyle\boldsymbol{\Sigma}^{m_{i}|o_{i}}=\boldsymbol{\Sigma}^{m_{i},m_{i}}-\boldsymbol{\Sigma}^{m_{i},o_{i}}(\boldsymbol{\Sigma}^{o_{i},o_{i}})^{-1}\boldsymbol{\Sigma}^{o_{i},m_{i}}.$

Missing values are filled with samples drawn from its full conditional distribution and combined with the fixed observed values to update subsequently other parameters from full conditionals. This updating process is called the multivariate normal imputation, first implemented by Schafer (1997), which is one of the multiple imputation methods. In the case of the mixture models, we assume all data points are generated from the mixture of multivariate normal distributions. A observation with the missing values belong to $h$th mixure component is imputed by sampling full conditional distribution for $\textbf{x}_{i}^{m_{i}}$ given $\boldsymbol{\mu}_{h}$ and $\boldsymbol{\Sigma}_{h}$ associated with $h$th mixture component. Mixture components that indicate clusters are determined by latent auxiliary variables. See Zhang and Everson (2004).

Dirichlet process (DP) introduced by Ferguson (1973) is the most popular Bayesian nonparametric model and used in many applications for clustering in the last two decades. Let $\mathcal{X}$ be a measurable space and $\mathcal{B}$ the Borel $\sigma$-field of subsets of $\mathcal{X}$. Then we say that the random probability measure $G$ on $(\mathcal{X},\mathcal{B})$ follows a Dirichlet process with a concentration parameter $\alpha>0$ and a baseline probability measure $G_{0}$, denoted by $G\sim DP(\alpha,G_{0})$, if for every finite disjoint partition $A_{1},\ldots A_{k}$ of $\mathcal{X}$,

$$(G(A_{1}),\ldots,G(A_{k}))\sim\text{Dir}(\alpha G_{0}(A_{1}),\ldots,\alpha G_{0}(A_{k})),$$

where Dir denotes the Dirichlet distribution. Dirichlet process can be represented in three different ways: (1) pólya urn scheme (Blackwell et al., 1973) (2) Chinese restaurant process (Aldous, 1985) (3) stick-breaking process (Sethuraman, 1994).

Sethuraman (1994) also showed that its realizations are discrete almost surely, even if $G_{0}$ is a continuous distribution. The discreteness of the DP implies that it is unsuitable prior for data generated from continuous distributions. To eliminate this drawback, Antoniak (1974) adopted Dirichlet process mixture models (DPMM) having the following hierarchical model formulations: for $i=1,\ldots,n,$

	$\displaystyle\textbf{x}_{i}|\theta_{i}$	$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}f(\textbf{x}_{i}|\theta_{i}),$
	$\displaystyle\theta_{i}$	$\displaystyle\stackrel{{\scriptstyle iid}}{{\sim}}G,$
	$\displaystyle G$	$\displaystyle\sim DP(\alpha,G_{0}),$

where $f$ is a parametric density function. The DPMM can be expressed as a limit of finite mixture models, where the number of mixture components is taken to infinity (Teh et al., 2006). For example, if a base measure $G_{0}$ is a multivariate normal Inverse-Wishart conjugate prior and $f(\textbf{x}_{i}|\theta_{i})$ is multivariate normal, the DPMM can be an infinite normal mixture model. The advantage of the DPMM is the number of mixture components is not fixed in advance of fitting the data and automatically infer the number of clusters from the data unlike finite normal mixture models. Various algorithms have been developed for posterior inference of the DPMM such as marginal Gibbs sampling (MacEachern, 1994; Escobar, 1994; Escobar and West, 1995; Bush and MacEachern, 1996; Neal, 2000), conditional Gibbs sampling (Ishwaran and James, 2001; Walker, 2007; Kalli et al., 2011; Ge et al., 2015), split-merge MCMC sampling (Jain and Neal, 2004), sequential updating and greedy search (SUGS) algorithms (Wang and Dunson, 2011) and variational approximation (Blei et al., 2006).

Generative models employ the Bayes theorem to build a classifier. Let the input or feature vector be X and the class labels be $Y$. We need the density of X conditioned on the class $k$, $P(\textbf{X}=\textbf{x}|Y=k)$ and the prior probability, $P(Y=k)$ to compute the posterior probability:

$$P(Y=k|\textbf{X}=\textbf{x})=\dfrac{P(Y=k)\cdot P(\textbf{X}=\textbf{x}|Y=k)}{\sum_{l=1}^{K}P(Y=l)\cdot P(\textbf{X}=\textbf{x}|Y=l)}.$$

(3)

Then, input x is assigned to the class having the highest posterior probability. Typical examples of generative models include linear discriminant analysis (LDA), quadratic discriminant, analysis (QDA), and Gaussian naive Bayes (GNB). The three models assume multivariate normal densities for $P(\textbf{X}=\textbf{x}|Y=k)$. Another example is to assume that the class-conditional density of X is a finite mixture of normals. This method is called mixture discriminant analysis (MDA) (Hastie and Tibshirani, 1996; Fraley and Raftery, 2002). Instead of using a finite mixture of normals, we propose to use the Dirichlet process mixture model, as the number of clusters is automatically determined and not limited.

We consider a binary classification problem and assume class labels $Y\in\{0,1\}$ has a binomial distribution with parameters $n$ and $p$. If we have $K$ $(>2)$ classes, $Y\in\{1,\ldots,K\}$ is assumed to have a multinomial distribution with $n$ and $\textbf{p}=(p_{1},\ldots,p_{K})$. The distribution of X conditioned on the class $k$ is modeled by the following hierarchical formulation for DPMM:

	$\displaystyle\textbf{x}_{i}|y_{i}=k$	$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}\mathcal{N}(\textbf{x}_{i}^{k};\boldsymbol{\mu}_{i}^{k},\boldsymbol{\Sigma}_{i}^{k}),$		(4)
	$\displaystyle(\boldsymbol{\mu}_{i}^{k},\boldsymbol{\Sigma}_{i}^{k})$	$\displaystyle\stackrel{{\scriptstyle iid}}{{\sim}}G_{k},\quad i=1,2,\ldots,n_{k},$
	$\displaystyle G_{k}$	$\displaystyle\sim DP(\alpha,M_{k}),\quad k=0,1,$

where $\textbf{x}_{i}^{k}$ is the $i$th feature vector and $\boldsymbol{\mu}^{k}$ and $\boldsymbol{\Sigma}^{k}$ are parameters of a normal distribution which belongs to only class $k$. Here, the base measure for class $k$, $M_{k}$ is the conjugate multivariate normal–inverse Wishart distribution, i.e. $M_{k}:=\mathcal{N}(\boldsymbol{\mu}^{k};\mathbf{m}_{0}^{k},\boldsymbol{\Sigma}^{k}/\tau_{0}^{k})\cdot\mathcal{I}\mathcal{W}(\boldsymbol{\Sigma}^{k};\mathbf{B}_{0}^{k},\nu_{0}^{k})$. Then the class-conditional density of X is given by

$$P(\textbf{X}=\textbf{x}|Y=k)=\sum_{h=1}^{\infty}w_{h}^{k}\mathcal{N}(\textbf{x};\boldsymbol{\mu}_{h}^{k},\boldsymbol{\Sigma}_{h}^{k}),\quad k\in\{0,1\}$$

(5)

using posterior samples for mixing proportion and cluster specific parameters. We use the improved slice sampler suggested by Ge et al. (2015) to generate samples from the posterior distribution. We describe the detailed MCMC algorithm for the DPMM on incomplete data in Appendix A. Since $Y$ has a binomial distribution, marginal probabilities over classes are estimated by

$$P(Y=k)=\frac{n_{k}}{n},\quad k\in\{0,1\},$$

(6)

where $n_{k}$ denotes the number of observations belonging to the class $k$. Putting (5) and (6) together in the equation (3), we can produce the following posterior probabilities for classes:

$$P(Y=k|\textbf{X}=\textbf{x})=\frac{n_{k}\times\sum_{h=1}^{\infty}w_{h}\mathcal{N}(\textbf{x};\boldsymbol{\mu}_{h}^{k},\boldsymbol{\Sigma}_{h}^{k})}{\sum_{l\in\{0,1\}}\left(n_{l}\times\sum_{h=1}^{\infty}w_{h}\mathcal{N}(\textbf{x};\boldsymbol{\mu}_{h}^{l},\boldsymbol{\Sigma}_{h}^{l})\right)}.$$

We compute posterior probabilities for all $k$ and choose only one class associated with the highest probability.

Both LDA and QDA assume multivariate normal densities for $P(\textbf{X}=\textbf{x}|Y=k)$ and GNB assumes that the features, $\textbf{x}_{1},\textbf{x}_{2},\ldots,\textbf{x}_{p}$ are conditional independent on $Y$ and $P(\textbf{x}_{j}|Y=k)$ is univariate normal density:

$$P(\textbf{X}=\textbf{x}|Y=k)=\prod_{j=1}^{p}P(\textbf{x}_{j}|Y=k).$$

Since three generative models are designed to different covariance structures of multivariate normal distributions, they have their respective classifiers. The decision boundaries for LDA are linear, while those for QDA are quadratic. GNB has not necessarily linear classification according to the data. However, MDA allows close approximation of not only linear but also nonlinear decision boundaries since mixture models can approximate arbitrary continuous distributions (Fraley and Raftery, 2002; Wang et al., 2010). Figure 1 shows that DPNB more accurately approximates most decision boundaries than do LDA, QDA, and GNB. In other words, DPNB is a much more flexible classifier than them. Furthermore, the classifier of DPNB provides comparable results with that of Support Vector Machine (SVM) and Random Forest (RF).

The DPNB can make inferences simultaneously on both unknown quantities and missing values. The same imputation approach described in the previous subsection 2.1 is capable of being applied in the mixture models with infinite components. For estimating (5), DPNB should divide the data into subsets where each subset belongs to only one class according to the inference process. Then it substitutes separately missing values with plausible values generated from samplers (2) based on respective subsets. We expect that imputation based on subsets is more accurate than that based on the full data. Because it is easy to find easily probable values estimated from homogeneous input data which belongs to only one class. If we spotlight imputation problems in the classification data, complete subsets are put together into a complete full dataset. In practice, the empirical results on 4 UCI datasets provide that this split and merge imputer has better performance than state-of-the-art imputation algorithms as the missing rate increases.

In reality, missing data may occur both training and test dataset. DPNB can make predictions for classes even if all elements of some features in the test set are absent. Let a new observed input vector denote $\textbf{x}_{\star}^{o_{\star}}$ and a new predicted class label $Y^{\star}$. The prediction rule for DPNB is given by

$$P(Y^{\star}=k|\textbf{x}_{\star}^{o_{\star}})=\dfrac{P(Y^{\star}=k)\cdot P(\textbf{x}_{\star}^{o_{\star}}|Y^{\star}=k)}{\sum_{l\in\{0,1\}}P(Y^{\star}=l)\cdot P(\textbf{x}_{\star}^{o_{\star}}|Y^{\star}=l)},\quad k\in\{0,1\}.$$

(7)

We need to compute the predictive density of $\textbf{x}_{\star}^{o_{\star}}$ conditioned on the class $k$ to complete the equation (7), the posterior probability of the new class label. It can be approximated by using the Monte Carlo Integration from the posterior samples in the training process as shown in (8). Therefore, DPNB can build classifiers regardless of missingness and missing rate on both training set and test set. Experiments also support that it predicts classes more accurately than do competitive models.

	$\displaystyle P(\textbf{x}_{\star}^{o_{\star}}|Y^{\star}=k)$	$\displaystyle=\int P(\textbf{x}_{\star}^{o_{\star}},\textbf{x}_{\star}^{m_{\star}}\,|Y^{\star}=k)\,\,d\textbf{x}_{\star}^{m_{\star}}$		(8)
		$\displaystyle=\int P(\textbf{x}_{\star}^{o_{\star}},\textbf{x}_{\star}^{m_{\star}}\,|\,\boldsymbol{\mu},\boldsymbol{\Sigma},Y^{\star}=k)\pi(\boldsymbol{\mu},\boldsymbol{\Sigma}\,|Y^{\star}=k)\,d\boldsymbol{\mu}\,d\boldsymbol{\Sigma}d\textbf{x}_{\star}^{m_{\star}}$
		$\displaystyle=\int\sum_{h}\pi_{h}^{k}\,P(\textbf{x}_{\star}^{o_{\star}},\textbf{x}_{\star}^{m_{\star}}\,|\boldsymbol{\mu}_{h},\boldsymbol{\Sigma}_{h},Y^{\star}=k)\pi(\boldsymbol{\mu}_{h},\boldsymbol{\Sigma}_{h}\,|Y^{\star}=k)\,d\boldsymbol{\mu}_{h}\,d\boldsymbol{\Sigma}_{h}\,d\textbf{x}_{\star}^{m_{\star}}$
		$\displaystyle\approx\frac{1}{M}\sum_{j=1}^{M}\sum_{h_{j}}\pi_{h_{j}}^{k}\int P(\textbf{x}_{\star}^{o_{\star}},\textbf{x}_{\star}^{m_{\star}}\,|\boldsymbol{\mu}_{h_{j}},\boldsymbol{\Sigma}_{h_{j}},Y^{\star}=k)\,d\textbf{x}_{\star}^{m_{\star}}\,\,(\because\text{MC integration})$
		$\displaystyle=\frac{1}{M}\sum_{j=1}^{M}\sum_{h_{j}}\pi_{h_{j}}^{k}P(\textbf{x}_{\star}^{o_{\star}}\,|\boldsymbol{\mu}_{h_{j}}^{o_{\star}},\boldsymbol{\Sigma}_{h_{j}}^{o_{\star},o_{\star}},Y^{\star}=k)$
		$\displaystyle=\frac{1}{M}\sum_{j=1}^{M}\sum_{h_{j}}\pi_{h_{j}}^{k}\cdot\mathcal{N}_{|o_{\star}|}\left(\textbf{x}_{\star}^{o_{\star}}\,|\,[\boldsymbol{\mu}_{h_{j}}^{k}]^{o_{\star}},[\boldsymbol{\Sigma}_{h_{j}}^{k}]^{o_{\star},o_{\star}}\right),$

where $M$ is the number of posterior samples generated by MCMC.

For each UCI dataset, we generate 100 different incomplete datasets removing 10% to 60% of the complete values based on the MCAR or MAR missing assumptions. We make use of the normalized root mean squared errors (NRMSE) as the imputation accuracy measure along with its standard deviation across the 100 replicated datasets. The NRMSE is defined as

$$\text{NRMSE}=\sqrt{\dfrac{\text{mean}((X^{\text{true}}-X^{\text{imp}})^{2})}{\text{Var}(X^{\text{true}})}},$$

where $X^{\text{true}}$ is the original data and $X^{\text{imp}}$ the imputed data.

We compare our proposed methods with several popular imputation methods such as multivariate imputation by chained equations using the predictive mean matching; Buuren and Groothuis-Oudshoorn (2010) (denoted by MICE), random forest based imputation; Stekhoven and Bühlmann (2012) (denoted by RF), KNN based imputation with the optimal number of neighbors which minimizes cross-validation errors; Troyanskaya et al. (2001) (denoted by KNN), imputer utilizing deep denoising autoencoders; Gondara and Wang (2018) (denoted by MIDA), and imputation technique by adapting generative adversarial nets; Yoon et al. (2018) (denoted by GAIN).

Figure 2 shows that the imputer of the DPNB model performs well in most datasets with either the lowest or the second lowest average NRMSE values across 100 replicates irrespective of the missingness mechanism (see Table 9 of Appendix B for details). In particular, we find out from both Table 2 and Table 3 that the proposed model is more accurate than other imputation algorithms as missing rates increase. The DPNB model can fill missing values via different covariance structures constructed by all available observations. However, methods to form conditional distributions using RF and MICE are less accurate than the DPNB because they utilize partial observations and variables instead of full information. Deep learning-based imputation methods such as GAIN and MIDA have poor performance due to model complexity relative to the number of observed data points. Among all models, the worst-performing method is the KNN based imputation.

We perform 10-fold cross-validation to obtain an estimate of the classification accuracy using simulated missing datasets generated from the previous subsection 4.1. The cross-validation process is repeated 10 times for fair comparisons. Since four UCI datasets have various class distributions and binary or multi-class classification problems, we use three types of classification performance metrics to measure the performance of classifiers:(a) accuracy rate, (b) area under the ROC curve (AUC), and (c) F1 score. The classification metrics corresponding to datasets are also shown in Table 4.

We compare the DPNB model with other prediction models on the incomplete dataset. Other approaches predict test cases after the imputation stage. Some details of the procedures regarding post-imputation prediction are as follows. First, we divide the dataset with missing values into training and test sets. Imputation algorithms utilized in subsection 4.1 fill missing values of training sets. Second, we build the support vector machines (SVM) with radial basis function (RBF) kernels as a benchmark classifier on the imputed training datasets. Third, missing values in test sets are replaced with the mean of the features of the imputed training set. Finally, we make predictions on the imputed test set using the trained model. We call all competing methods types of imputation algorithms used in both training and test phases.

Figure 3 provides that the DPNB model comes up with the best performance except for the Breast Cancer Wisconsin (Diagnostic) dataset, where it is competitive. Both Table 5 and Table 6 also show that the DPNB has a better classification performance than the others on a dataset for higher missing rates. As missing rates increase, the prediction accuracy of the DPNB is not greatly reduced unlike other competitive approaches. This supports that the DPNB is less affected by both missing data mechanisms and proportions missing proportions than other methods. See Table 10 of Appendix B for further details about these results.

To assess the imputation accuracy based on the real world dataset, we consider three scenarios by deleting artificially 50, 100, and 500 complete values in the semiconductor dataset under MCAR assumption. In every cases, we made 100 replicates, respectively. For performance comparisons, we then compute the average of normalized root mean squared errors obtained by each method in all cases.

The results of the DPNB and competing models are provided in Figure 5 and Table 7. They demonstrate that the DPNB model provides more accurate values and consistent performance through both minimum mean and the lowest standard deviation values over 100 replications of NRMSEs. As shown in the Breast Cancer example, both RF and MICE yield similar imputation performance with the DPNB model. Imputation methods based on deep architectures including MIDA and GAIN have difficulty selecting appropriate parameters to prevent the problem of overfitting. These methods have then lower imputation accuracy than the DPNB, RF and MICE.

In this experiment, we randomly partitioned the real dataset into training and test sets with different ratios. The training set was set from 50% to 90% of the overall dataset and the test set was accordingly set from 50% to 10% of it. For example, if 90% of the data is used as the training set, then the remaining 10% is used as the test set. Likewise, we created 100 replicate datasets in all cases. Furthermore, the bootstrapping-based oversampling technique that replicates observations from minority class was applied into training sets in order to balance the two classes. Note that the same indices of observations selected from the minority class were used in all methods given percentages of the training sets.

We add the SVM with a linear kernel as another base classifier. We will name competitive methods by combining imputation techniques and kernels of support vector machine, e.g., “RF+L” means that the RF imputes missing values of training sets and the SVM with a linear kernel is fitted on the imputed training sets and predicts the test set values. The results of this experiment are given in Table 8. The columns of competing methods represent their average F1-score over 100 replicate datasets, normalized by the average F1-score of the DPNB model. The values higher than one indicate that the methods provide better performance than the DPNB and values lower than one indicate worse performance in imbalanced settings.

As shown in Table 8, the DPNB method performs well in practice since it ranks high in almost all cases. In particular, it yields the best prediction performance on two out of the five ratios. We can say that the DPNB model provides stable performance by obtaining the highest average normalized F1-score and average rank and stable inference in the incomplete data with high missing rates. However, since the semiconductor dataset has very high missing rates, the DPNB does not provide remarkable performance against other methods as shown in the previous section. The methods using the KNN imputer have poor performance irrespective of the type of classifiers. We also present the actual average F1-score in Table 11 of Appendix B for more details.