A Bayesian Finite Mixture Model with Variable Selection for Data with Mixed-type Variables

TL;DR

This paper introduces a Bayesian finite mixture model that performs variable selection, handles censored biomarker data due to detection limits, and improves clustering robustness over traditional EM-based methods.

Contribution

It develops a Bayesian approach with Gibbs sampling for variable importance, censored data handling, and clustering, addressing key limitations of existing mixture models.

Findings

The model effectively identifies important variables in mixed-type data.

It accurately handles censored biomarker data in clinical datasets.

Simulation and real data analyses demonstrate improved clustering performance.

Abstract

Finite mixture model is an important branch of clustering methods and can be applied on data sets with mixed types of variables. However, challenges exist in its applications. First, it typically relies on the EM algorithm which could be sensitive to the choice of initial values. Second, biomarkers subject to limits of detection (LOD) are common to encounter in clinical data, which brings censored variables into finite mixture model. Additionally, researchers are recently getting more interest in variable importance due to the increasing number of variables that become available for clustering. To address these challenges, we propose a Bayesian finite mixture model to simultaneously conduct variable selection, account for biomarker LOD and obtain clustering results. We took a Bayesian approach to obtain parameter estimates and the cluster membership to bypass the limitation of the EM algorithm. To account for LOD, we added one more step in Gibbs sampling to iteratively fill in biomarker values below or above LODs. In addition, we put a spike-and-slab type of prior on each variable to obtain variable importance. Simulations across various scenarios were conducted to examine the performance of this method. Real data application on electronic health records was also conducted.

Figures

Equations

$$\{\begin{array}{ll}{L}_m|\mathbf{X},\mathbf{Z},G,\mathbf{\beta }={\prod }_{i=1}^n{\prod }_{g=1}^G[f(x_{im}|{\mathbf{\beta }}_{mg}){]}^{z_{ig}}& \text{if variable }m\text{ is relevant}\\ L_m|\mathbf{X},\mathbf{\beta }={\prod }_{i=1}^{n}f(x_{im}|{\mathbf{\beta }}_m)& \text{if variable }m\text{ is irrelevant}\end{array}$$

$$l=\sum _{i=1}^{n}log\{\sum _{g=1}^G{\tau }_{g}f({\mathbf{x}}_i|{\mathbf{\beta }}_g)\}-Q_{\lambda }(\mathbf{\beta })$$

$$f({\mathbf{x}}_i|\mathbf{\tau },\mathbf{\theta },{\mathbf{A}}_1,\mathbf{\mu },{\mathbf{\sigma }}^2)=\sum _{g=1}^G{\tau }_{g}f({\mathbf{x}}_i|{\mathbf{\theta }}_g,{\mathbf{A}}_1,{\mathbf{\mu }}_g,{\mathbf{\sigma }}^2),$$

$$f({\mathbf{x}}_i|{\mathbf{z}}_i,\mathbf{\theta },{\mathbf{A}}_1,\mathbf{\mu },{\mathbf{\sigma }}^2)=\prod _{g=1}^G[f({\mathbf{x}}_i|{\mathbf{\theta }}_g,{\mathbf{A}}_1,{\mathbf{\mu }}_g,{\mathbf{\sigma }}^2){]}^{z_{ig}},$$

$$\mathbf{\tau }\sim Dir\left({\delta }_1,{\delta }_2,\dots {\delta }_G\right),$$

$$\{\begin{array}{ll}{\mu }_{mg}\sim N\left(0,{\sigma }_{{\Delta }_0}^2\right)& \text{if }{\Delta }_{mg}=0\\ {\mu }_{mg}\sim N\left(0,{\sigma }_{{\Delta }_1}^2\right)& \text{if }{\Delta }_{mg}=1,\end{array}$$

$${\sigma }_{{\Delta }_0}^2\sim Inv\Gamma \left(a_{{\Delta }_0},b_{{\Delta }_0}\right),$$

$${\gamma }_m\sim \Gamma \left(\tilde{a},\tilde{b}\right),$$

$$\{\begin{array}{ll}{\mathbf{\theta }}_{mg}\sim Dir\left({\mathbf{\alpha }}_{m{\Delta }_0}\right)& \text{if }{\Delta }_{mg}=0\\ {\mathbf{\theta }}_{mg}\sim Dir\left({\mathbf{\alpha }}_{{\Delta }_1}\right)& \text{if }{\Delta }_{mg}=1\end{array}$$

$$\{\begin{array}{ll}{\Delta }_{mg}\sim Bern\left(p_{1m}\right)& \text{if }m\in \{1,2,\dots ,q\}\\ {\Delta }_{mg}\sim Bern\left(p_{2m}\right)& \text{if }m\in \{q+1,q+2,\dots ,M\}\end{array}$$

$$p_{1m}\sim Beta(a_{p_1},b_{p_1})$$

$$p_{2m}\sim Beta(a_{p_2},b_{p_2})$$

$$f(Parameters|\mathbf{X})$$

$$\times \prod _{i=1}^{n}f({\mathbf{z}}_i|\mathbf{\tau })f(\mathbf{\tau })\prod _{m=1}^{q}N(A_{1m};{\mu }_A,{\sigma }_A^2)\prod _{m=1}^q\Gamma ({\gamma }_m;\tilde{a},\tilde{b})$$

$$\times \prod _{m=1}^q\prod _{g=2}^G[N({\mu }_{mg};0,{\sigma }_{{\Delta }_1}^2){\Delta }_{mg}+N({\mu }_{mg};0,{\sigma }_{{\Delta }_0}^2)(1-{\Delta }_{mg})]$$

$$\times Inv\Gamma ({\sigma }_{{\Delta }_0}^2;a_{{\Delta }_0},b_{{\Delta }_0})\prod _{m=1}^q\prod _{g=2}^{G}Bern({\Delta }_{mg};p_{1m})\prod _{m=1}^{q}Beta(p_{1m};a_{p_1},b_{p_1})$$

$$\times \prod _{m=q+1}^M\prod _{g=1}^G[Dir({\theta }_{mg};{\mathbf{\alpha }}_{{\Delta }_1}){\Delta }_{mg}+Dir({\theta }_{mg};{\mathbf{\alpha }}_{m{\Delta }_0})(1-{\Delta }_{mg})]$$

$$\times \prod _{m=q+1}^M\prod _{g=1}^{G}Bern({\Delta }_{mg};p_{2m})\prod _{m=q+1}^{M}Beta(p_{2m};a_{p_2},b_{p_2})$$

$${\mathbf{z}}_i|{\mathbf{x}}_i,\mathbf{\tau },{\mathbf{A}}_1,\mathbf{\mu },\mathbf{\gamma },\mathbf{\theta }\sim Multinomial\left(\frac{{\tau }_{1}f({\mathbf{x}}_i|{\mathbf{A}}_1,{\mathbf{\mu }}_1,{\mathbf{\gamma }}_1,{\mathbf{\theta }}_1)}{{\sum }_{g=1}^G{\tau }_{g}f({\mathbf{x}}_i|{\mathbf{A}}_1,{\mathbf{\mu }}_g,{\mathbf{\gamma }}_g,{\mathbf{\theta }}_g)},\dots ,\frac{{\tau }_{G}f({\mathbf{x}}_i|{\mathbf{A}}_1,{\mathbf{\mu }}_G,{\mathbf{\gamma }}_G,{\mathbf{\theta }}_G)}{{\sum }_{g=1}^G{\tau }_{g}f({\mathbf{x}}_i|{\mathbf{A}}_1,{\mathbf{\mu }}_g,{\mathbf{\gamma }}_g,{\mathbf{\theta }}_g)}\right).$$

$$\mathbf{\tau }|\mathbf{Z},{\delta }_1,\dots ,{\delta }_G\sim Dir\left({\delta }_1+\sum _{i=1}^n{z}_{i1},\dots ,{\delta }_G+\sum _{i=1}^n{z}_{iG}\right).$$

$$A_{1m}|\mathbf{X},\mathbf{Z},{\mathbf{\mu }}_m,{\gamma }_m,{\mu }_A,{\sigma }_A^2\sim N\left(\frac{{\sigma }_A^2({\sum }_{i=1}^n{x}_{im}-{\sum }_{i=1}^n{\sum }_{g=1}^G{z}_{ig}{\mu }_{mg})+({\mu }_A/{\gamma }_m)}{n{\sigma }_A^2+(1/{\gamma }_m)},\frac{{\sigma }_A^2/{\gamma }_m}{n{\sigma }_A^2+(1/{\gamma }_m)}\right).$$

$${\mu }_{mg}|\mathbf{X},\mathbf{Z},A_{1m},{\gamma }_m,{\sigma }_{{\Delta }_1}^2\sim N\left(\frac{{\sigma }_{{\Delta }_1}^2({\sum }_{i=1}^n{z}_{ig}{x}_{im}-{\sum }_{i=1}^n{z}_{ig}{A}_{1m})}{{\sigma }_{{\Delta }_1}^2{\sum }_{i=1}^n{z}_{ig}+(1/{\gamma }_m)},\frac{{\sigma }_{{\Delta }_1}^2/{\gamma }_m}{{\sigma }_{{\Delta }_1}^2{\sum }_{i=1}^n{z}_{ig}+(1/{\gamma }_m)}\right),$$

$${\mu }_{mg}|\mathbf{X},\mathbf{Z},A_{1m},{\gamma }_m,{\sigma }_{{\Delta }_0}^2\sim N\left(\frac{{\sigma }_{{\Delta }_0}^2({\sum }_{i=1}^n{z}_{ig}{x}_{im}-{\sum }_{i=1}^n{z}_{ig}{A}_{1m})}{{\sigma }_{{\Delta }_0}^2{\sum }_{i=1}^n{z}_{ig}+(1/{\gamma }_m)},\frac{{\sigma }_{{\Delta }_0}^2/{\gamma }_m}{{\sigma }_{{\Delta }_0}^2{\sum }_{i=1}^n{z}_{ig}+(1/{\gamma }_m)}\right).$$

$${\gamma }_m|\mathbf{X},\mathbf{Z},A_{1m},{\mathbf{\mu }}_m,\tilde{a},\tilde{b}\sim \Gamma \left(\tilde{a}+\frac{1}{2}n,\tilde{b}+\frac{1}{2}\sum _{i=1}^n\sum _{g=1}^G{z}_{ig}(x_{im}-A_{1m}-{\mu }_{mg}{)}^2\right).$$

$${\sigma }_{{\Delta }_0}^2|\mathbf{\Delta },\mathbf{\mu },a_{{\Delta }_0},b_{{\Delta }_0}\sim Inv\Gamma \left(a_{{\Delta }_0}+\frac{1}{2}\sum _{m=1}^q\sum _{g=2}^G(1-{\Delta }_{mg}),b_{{\Delta }_0}+\frac{1}{2}\sum _{m=1}^q\sum _{g=2}^G(1-{\Delta }_{mg}){\mu }_{mg}^2\right).$$

$${\mathbf{\theta }}_{mg}|\mathbf{X},\mathbf{Z},{\mathbf{\alpha }}_{{\Delta }_1}\sim Dir\left({\alpha }_{{\Delta }_{1}1}+\sum _{i=1}^n{x}_{im1}{z}_{ig},\dots ,{\alpha }_{{\Delta }_1{L}_m}+\sum _{i=1}^n{x}_{im{L}_m}{z}_{ig}\right)$$

$${\mathbf{\theta }}_{mg}|\mathbf{X},\mathbf{Z},{\mathbf{\alpha }}_{m{\Delta }_0}\sim Dir\left({\alpha }_{m{\Delta }_{0}1}+\sum _{i=1}^n{x}_{im1}{z}_{ig},\dots ,{\alpha }_{m{\Delta }_0{L}_m}+\sum _{i=1}^n{x}_{im{L}_m}{z}_{ig}\right).$$

$${\Delta }_{mg}|{\mu }_{mg},p_{1m},{\sigma }_{{\Delta }_0}^2,{\sigma }_{{\Delta }_1}^2\sim Bernoulli\left(\frac{p_{1m}N({\mu }_{mg};0,{\sigma }_{{\Delta }_1}^2)}{p_{1m}N({\mu }_{mg};0,{\sigma }_{{\Delta }_1}^2)+(1-p_{1m})N({\mu }_{mg};0,{\sigma }_{{\Delta }_0}^2)}\right)$$

$${\Delta }_{mg}|{\mu }_{mg},p_{2m},{\mathbf{\alpha }}_{m{\Delta }_0},{\mathbf{\alpha }}_{{\Delta }_1}\sim Bernoulli\left(\frac{p_{2m}Dir({\theta }_{mg};{\mathbf{\alpha }}_{{\Delta }_1})}{p_{2m}Dir({\theta }_{mg};{\mathbf{\alpha }}_{{\Delta }_1})+(1-p_{2m})Dir({\theta }_{mg};{\mathbf{\alpha }}_{m{\Delta }_0})}\right).$$

$$p_{1m}|{\mathbf{\Delta }}_m,a_{p_1},b_{p_1}\sim Beta\left(a_{p_1}+\sum _{g=2}^G{\Delta }_{mg},b_{p_1}+\sum _{g=2}^G(1-{\Delta }_{mg})\right).$$

$$p_{2m}|{\mathbf{\Delta }}_m,a_{p_2},b_{p_2}\sim Beta\left(a_{p_2}+\sum _{g=1}^G{\Delta }_{mg},b_{p_2}+\sum _{g=1}^G(1-{\Delta }_{mg})\right).$$

$$f(x_{i^{\prime }m}|A_{1m},{\mu }_{mg},{\gamma }_m)={\left(\frac{\sqrt{{\gamma }_m}h(x_{i^{\prime }m})}{\Phi [\sqrt{{\gamma }_m}(C_{Lm}-{\mu }_{mg}-A_{1m})]}\right)}^{z_{ig}}$$

$$f(x_{i^{\ast }m}|A_{1m},{\mu }_{mg},{\gamma }_m)={\left(\frac{\sqrt{{\gamma }_m}h(x_{i^{\ast }m})}{1-\Phi [\sqrt{{\gamma }_m}(C_{Um}-{\mu }_{mg}-A_{1m})]}\right)}^{z_{ig}},$$

Full text

A Bayesian Finite Mixture Model with Variable Selection for Data with Mixed-type Variables

Shu Wang

Department of Biostatistics, College of Public Health and Health Professions, University of Florida

University of Florida Health Cancer Center

Jonathan G. Yabes

Department of Biostatistics, Graduate School of Public Health, University of Pittsburgh

Department of Medicine, School of Medicine, University of Pittsburgh

Department of Clinical and Translational Science, School of Medicine, University of Pittsburgh

Chung-Chou H. Chang

Department of Biostatistics, Graduate School of Public Health, University of Pittsburgh

Department of Medicine, School of Medicine, University of Pittsburgh

Department of Clinical and Translational Science, School of Medicine, University of Pittsburgh

Finite mixture model is an important branch of clustering methods and can be applied on data sets with mixed types of variables. However, challenges exist in its applications. First, it typically relies on the EM algorithm which could be sensitive to the choice of initial values. Second, biomarkers subject to limits of detection (LOD) are common to encounter in clinical data, which brings censored variables into finite mixture model. Additionally, researchers are recently getting more interest in variable importance due to the increasing number of variables that become available for clustering.

To address these challenges, we propose a Bayesian finite mixture model to simultaneously conduct variable selection, account for biomarker LOD and obtain clustering results. We took a Bayesian approach to obtain parameter estimates and the cluster membership to bypass the limitation of the EM algorithm. To account for LOD, we added one more step in Gibbs sampling to iteratively fill in biomarker values below or above LODs. In addition, we put a spike-and-slab type of prior on each variable to obtain variable importance. Simulations across various scenarios were conducted to examine the performance of this method. Real data application on electronic health records was also conducted.

KEY WORDS: Censored biomarker, Clustering, Finite mixture model, Mixed data, Variable selection

1 Introduction

The finite mixture model (FMM) (McCutcheon,, 1987; Moustaki,, 1996; Nylund et al.,, 2007) has been used to uncover the latent mixture probability distributions in a combined statistical distribution of a population (Deb et al.,, 2008) when the population is heterogeneous in characteristics and consisting with a combination of several more homogeneous subgroups. This method has a natural interpretation of heterogeneity through the mixture of finite components with a distributional assumption conditioning on components for each variable.

Unlike the commonly used nonparametric clustering algorithms (e.g., distance-based, density-based), the FMM technique is a model-based method that can easily handle variables with mixed variable types and can do variable selection while performing clustering. This is especially important when the number of variables involved is large.

The current method of FMM has to deal with several challenges. The first challenge is related to the use of the EM algorithm (Dempster et al.,, 1977), which is the mostly commonly used estimation procedure for the FMMs. In the EM algorithm the convergence rate could be very slow and the solution could be highly dependent on the choices of initial values, especially in the multivariate settings (Biernacki et al.,, 2003; Karlis and Xekalaki,, 2003; McLachlan and Krishnan,, 2007).

The second challenge of using FMM is how to handle censored probability distributions. Oftentimes, covariates collected from medical setting include biomarker data of patients. Biomarkers are characteristics that can be accurately and reproducibly measured and accessed as an indicator of various biological processes (Group et al.,, 2001; Strimbu and Tavel,, 2010). Different biomarkers serve for different purposes. For example, temperature can be seen as a biomarker for fever; C-reactive protein (CRP) and Interleukin 6 (IL-6) are commonly used as biomarkers for sepsis (Pierrakos and Vincent,, 2010). Therefore, biomarkers usually contain important diagnosis information about subjects. However, many biomarkers are subject to a limit of detection (LOD). LOD could be lower detection limit (i.e., values below this limit could not be measured), higher detection limit (i.e., values above this limit could not be measured), or both. When biomarkers are outcome variables in data analysis, biomarker values more extreme than the detection limit are usually represented by the detection limit value and an additional binary variable is included to indicating whether the corresponding value is actually undetected or measured. A semiparametric censored regression model can then be used (Tobin,, 1958; Powell,, 1984; Honoré,, 1992). When biomarkers are predictors in analysis, multiple imputations are often used (Lubin et al.,, 2004; Lee et al.,, 2012; Bernhardt et al.,, 2015). If the objective of the analysis is to cluster the feature space when data containing censored biomarkers, the two aforementioned approaches cannot be applied. Therefore, conducting clustering for data with LODs is still not well-addressed.

The objective of this paper is threefold. First, to overcome the estimation limitation of the EM algorithm, we adopt Gibbs sampling, which has been shown to be a valid and practical way in estimation of the FMMs (Geman and Geman,, 1984), (Diebolt and Robert,, 1994). Second, we propose the use of a spike-and-slab type prior for categorical variables for the purpose of variable selection. Together with the traditional spike-and-slab prior for continuous variables, we incorporate variable selection into our Bayesian framework to provide quantitative information about variable importance. Third, we introduce an additional sampling step into our framework so that it is able to handle censored biomarker variables that are often encountered in clinical data.

Section 2 contains our review of currently used variable selection methods for the finite mixture models. We describe our proposed method in detail in Section 3. Section 4 includes simulation studies that are used to assess the performance of our methods and compare the performance to existing methods. In Section 5, our proposed method is applied to identify sepsis phenotypes using demographic, clinical, and biomarker data collected from electronic health records. Our conclusions and the summary of this study are in Section 6.

2 Existing methods review

Overall, there are two categories of variable selection methods for supervised or unsupervised machine learning: filter methods and wrapper methods (Blum and Langley,, 1997; Guyon and Elisseeff,, 2003; Fop et al.,, 2018). Filter methods refer to those whose feature selection procedures are conducted separately with clustering procedures. On the contrary, wrapper methods refer to those whose feature selection is conducted simultaneously with clustering procedures, like “wrapped” around clustering procedures. For example, the step 1 of the HyDaP algorithm (Wang et al.,, 2019) can be viewed as a filter method as it is conducted separately with the actual clustering algorithm. Wrapper methods are relatively more popular since it is naturally incorporated in clustering algorithms. In this section we focus on wrapper methods for FMMs.

Liu et al. proposed to conduct principle component analysis (PCA) before fitting Gaussian finite mixture models (Liu et al.,, 2003). This method assumes that only the first $K$ factors are relevant to clustering, where $K$ is a random variable that has a prior distribution. However, factors having larger eigen values in PCA do not necessarily contain more important information for clustering (Chang,, 1983).

Law et al. defined a binary indicator called feature saliency for each variable to reflect whether this variable is relevant to clustering or not. EM algorithm was used for estimation (Law et al.,, 2004). Let $\phi_{m}$ denote saliency for variable $m$. If $\phi_{m}=1$ then variable $m$ is relevant, otherwise variable $m$ is not, namely its distribution is independent of cluster labels. Let $\bm{X}$ denote a data matrix with $n$ subjects and $M$ variables. Let $x_{im}$ denote variable $m$ of subject $i$. Let $G$ denote number of clusters. Let $\bm{Z}$ denote cluster indicator matrix; $z_{ig}$ denote indicator variable of subject $i$ belonging to cluster $g$. Let $\bm{\beta}$ denote distributional parameters, $\bm{\beta}_{mg}$ denote parameters of variable $m$ in cluster $g$, $\bm{\beta}_{m}$ denote marginal parameter of variable $m$. Then $L_{m}$, the likelihood due to variable $m$, is defined as:

[TABLE]

Tadesse et al. and Li et al. used similar definitions of likelihood function and feature saliency. (Tadesse et al.,, 2005; Li et al.,, 2009). Different with the one proposed by Law et al., the method proposed by Tadesse et al. detects discriminating variables through reversible-jump MCMC instead of EM algorithm for high-dimensional Gaussian finite mixture models. Later, White et al. applied this idea on latent class analysis (White et al.,, 2016) for Bayesian variable selection. While the method proposed by Li et al. adopted Variational Learning of Bayesian approximation (VB) for inference. They claimed that although using EM algorithm and VB usually lead to identical results, VB could avoid the situation of getting infinite likelihood when there is singular cluster, which may encounter using EM. Later Sun et.al. used the same framework and extended Gaussian mixture model to Student’s t mixture model to better handle outliers (Sun et al.,, 2018).

Another category of methods is penalization approach (Fop et al.,, 2018). The general idea is to maximize a penalized log-likelihood which is defined as:

[TABLE]

where $\bm{x}_{i}$ is vector $(x_{i1},x_{i2},\dots,x_{iM})^{T}$; $\bm{\beta}$ is vector $(\bm{\beta}_{1},\bm{\beta}_{2},\dots,\bm{\beta}_{G})^{T}$ and each $\bm{\beta}_{g}$ represents distributional parameters for cluster $g$; $Q$ is a function of distributional parameters; $\lambda$ is penalty parameter. Different penalty functions were proposed including $L1$ penalty (Pan and Shen,, 2007), sample size weighted $L1$ penalty (Bhattacharya and McNicholas,, 2014), $L_{\infty}$ penalty (Wang and Zhu,, 2008) and other variations to achieve the goal of variable selection. These methods all assume that the differences of mixture components lie in mean parameters, which assumes similar variance within all clusters. This assumption was later relaxed by Xie et al. by adding two penalty terms so that the variance covariance structure is cluster-specific diagonal matrices (Xie et al.,, 2008). Later this method was further extended to unconstrained variance covariance matrices (Zhou et al.,, 2009).

Another category of methods that is commonly applied is spike and slab prior (Mitchell and Beauchamp,, 1988; Madigan and Raftery,, 1994; George and McCulloch,, 1997; Ishwaran et al.,, 2005). It is demonstrated in Figure 1. Spike and slab prior is originally proposed for variable selection in linear regressions, but it can be naturally applied on Gaussian mixture models for unsupervised clustering. This method first chooses one cluster as reference, and obtain mean difference between other clusters with this reference cluster. Similar to feature saliency, a binary indicator representing variable importance is also defined. For other methods defined feature saliency, each variable has only one feature saliency indicating whether it is relevant to clustering or not. While under spike and slab prior, each variable first has a cluster-specific importance, for cluster $2,3,\dots,G$ and then we aggregate these cluster-specific importance values into one overall importance or weight value. Let $\Delta_{mg}$ denote importance of variable $m$ within cluster $g$, $\mu_{mg}$ the mean difference between cluster $g$ and cluster $1$ of variable $m$, where $g=2,3,\dots,G$. $\Delta_{mg}=0$ indicates that for variable $m$, cluster $g$ is not different with reference cluster, namely variable $m$ within cluster $g$ is not important. Then in next iteration, we assign $\mu_{mg}$ a “spike” prior (the grey density in Figure 1): a Gaussian prior centered at 0 with very small variance. In this way, a value close to 0 is very likely to be sampled as the updated value of $\mu_{mg}$. Otherwise if $\Delta_{mg}=1$, it indicates that for variable $m$, cluster $g$ is different with reference cluster, namely variable $m$ within cluster $g$ is important. Then we assign $\mu_{mg}$ a “slab” prior (the red density in Figure 1): a Gaussian prior centered at 0 with very large variance. Therefore, any value is possible to be sampled as the updated value of $\mu_{mg}$. Let $\Delta_{m}$ denote final weight of variable $m$. Within each iteration, as long as $\Delta_{mg}=1$ for at least one cluster, we assign 1 to $\Delta_{m}$ for that iteration to represent that variable $m$ is important. After all iterations, we calculate $p(\Delta_{m}=1)$ to be the weight of variable $m$.

Although researchers have explored a lot of variable selection methods for finite mixture models, all above methods only focus on single type of variables. Not many methods incorporate both types of variables. The method proposed by Raftery and Dean is able to handle various types of variables through stepwise regression like procedures (Raftery and Dean,, 2006). They used approximate Bayes factor as variable selection criteria and obtained final optimal variable sets through a greedy search algorithm. This method naturally takes advantage of the parametric form of finite mixture models, but it could be extremely slow if a data set contains a large number of variables. Besides, this method only identifies that whether a variable is important or not instead of providing quantitative values for importance. Later it was extended through combining LASSO-like procedures and model selection procedures (Celeux et al.,, 2018) to be more efficient. But similarly, this method does not provide real-valued variable weights.

3 Bayesian finite mixture model with variable selection

3.1 Notation and proposed model

We will first define some notations that will be used to define our proposed Bayesian finite mixture model.

Let $\bm{X}$, $\mathbf{x}_{i}=(x_{i1},x_{i2},\dots,x_{iM})^{T}$ denote the vector of covariates for subject $i$ and $\bm{X}=(\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{n})^{T}$ be an $n\times M$ data matrix, where $n$ is the number of subjects and $M$ is the number of variables. Let $\mathbf{z}_{i}=(z_{i1},z_{i2},\dots,z_{iG})^{T}$ denote the vector of cluster-membership indicators for subject $i$, where $z_{ig}=1$ if the subject $i$ belongs to cluster $g$ and [math] otherwise. $\sum_{g=1}^{G}z_{ig}=1$, where $G$ is the number of clusters. Without loss of generality, we let the first $q$ variables to be continuous normally distributed and the rest $q+1^{th}$ to $M^{th}$ variables to be categorical.

Let $\bm{A}_{1}=(A_{11},A_{12},\dots,A_{1q})^{T}$ denote the vector of mean values of all continuous variables in cluster $1$, where $A_{1m}$ is the mean of variable $m$ in cluster $1$, $m=1,2,\dots,q$. We also define $\bm{\sigma}^{2}=(\sigma_{1}^{2},\sigma_{2}^{2},\dots,\sigma_{q}^{2})^{T}$ as the vector of variances for all continuous variables, where $\sigma^{2}_{m}$ is the variance of variable $m$, $m=1,2,\dots,q$.

We further define $\bm{\mu}$ as a $q\times G$ matrix with its $(m,g)$ element $\mu_{mg}$ representing the mean difference of variable $m$ between cluster $g$ and cluster $1$, where $m=1,2,\dots,q$. Note that cluster $1$ is the reference cluster. Then, vector $\bm{\mu}_{g}=(\mu_{1g},\mu_{2g},\dots,\mu_{qg})^{T}$ represents the mean difference between cluster $g$ and cluster $1$ of all normally distributed variables and vector $\bm{\mu}_{m}=(\mu_{m1},\mu_{m2},\dots,\mu_{mG})^{T}$ represents the mean differences of variable $m$ for all $G$ clusters. For identifiability, $\bm{\mu}_{1}$ is set to be a vector with all elements being $1$.

For categorical variables, let $\bm{\theta}$ be a $(M-q)\times G$ matrix with its $(m,g)$ element $\bm{\theta}_{mg}$ representing the distributional parameters of variable $m$ within cluster $g$, where $m=q+1,q+2,\dots,M$. Vector $\bm{\theta}_{g}=(\bm{\theta}_{(q+1)g},\bm{\theta}_{(q+2)g},\dots,\bm{\theta}_{Mg})^{T}$ represents the parameters of all categorical variables within cluster $g$ and vector $\bm{\theta}_{m}=(\bm{\theta}_{m1},\bm{\theta}_{m2},\dots,\bm{\theta}_{mG})^{T}$ represents the parameters of variable $m$ for $m=q+1,q+2,\dots,M$.

We define density function of the finite mixture model (FMM) with $G$ clusters as

[TABLE]

where $\sum_{g=1}^{G}\tau_{g}=1$ and $\tau_{g}$ is the probability that a subject belongs to cluster $g$. We let $\bm{\tau}=(\tau_{1},\tau_{2},\dots,\tau_{G})^{T}$ be the vector of cluster mixture probabilities.

For subject $i$, the corresponding density function for data $(\mathbf{x}_{i}$ given its cluster membership $\mathbf{z}_{i})$ could be written as:

[TABLE]

where $f(\mathbf{x}_{i}|\bm{\theta}_{g},\bm{A}_{1},\bm{\mu}_{g},\bm{\sigma}^{2})=\prod_{m=1}^{M}f(x_{im}|\theta_{mg},A_{1m},\mu_{mg},\sigma_{m}^{2})$ based on conditional independence, which assumes variables are independent with each other conditional on the cluster labels.

3.2 Priors

We define conjugate priors for all parameters specified in the above-defined FMM. Details of the prior distribution for each parameter is specified below.

For the cluster indicator matrix $\mathbf{Z}$, we specify a multinomial prior distribution with the form: $\mathbf{z}_{i}\sim Multinomial\left(G,\tau_{1},\tau_{2},\dots\tau_{G}\right)$. We also let the probability of cluster membership $\bm{\tau}$ follow a Dirichlet distribution

[TABLE]

where $\delta_{1},\delta_{2},\dots\delta_{G}$ are the hyper-parameters of $\bm{\tau}$.

Let parameter $A_{1m}$ follows $N\left(\mu_{A},\sigma_{A}^{2}\right)$, where $\mu_{A}$ and $\sigma_{A}^{2}$ are both hyper-parameters.

Let $\bm{\Delta}$ be a $G\times M$ matrix representing the collection of importance indicator for each variable, where $\Delta_{mg}$ is the importance indicator of variable $m$ within cluster $g$. Vector $\bm{\Delta}_{g}$ is $(\Delta_{1g},\Delta_{2g},\dots,\Delta_{Mg})^{T}$ representing the importance indicator of all variables within cluster $g$; and vector $\bm{\Delta}_{m}$ is $(\Delta_{m1},\Delta_{m2},\dots,\Delta_{mG})^{T}$ representing the importance indicator of variable $m$. For the purpose of variable selection, we apply a spike-and-slab priors for parameter $\mu_{mg}$ with the form

[TABLE]

where $\sigma_{\Delta_{0}}^{2}$ and $\sigma_{\Delta_{1}}^{2}$ are both hyper-parameters. For identifiability, we define $\mu_{mg}=0$ when $g=1$. When $\Delta_{mg}=1$, namely cluster $g$ is different from cluster $1$, $\mu_{mg}$ should be away from 0. Therefore, $\sigma_{\Delta_{1}}^{2}$ should be a very large number, for increasing the variability of cluster mean difference, so that it is likely to be different with [math]. When $\Delta_{mg}=0$, which occurs if cluster $g$ is not different from cluster $1$, $\mu_{mg}$ should be close to 0. Therefore, $\sigma_{\Delta_{0}}^{2}$ should be a very small number so that the mean differences among cluster means would be close to zero. We used an Inverse-Gamma prior for parameter $\sigma_{\Delta_{0}}^{2}$:

[TABLE]

where $a_{\Delta_{0}}$ and $b_{\Delta_{0}}$ are both hyper-parameters.

The precision parameter $\gamma_{m}=\sigma_{m}^{-2}$ is assigned a prior following a Gamma distribution,

[TABLE]

where $\tilde{a}$ and $\tilde{b}$ are both hyper-parameters.

Inspired by the usual spike-and-slab priors for continuous variables, we propose similar priors for categorical variables. The challenge of building a spike-and-slab prior for a categorical variable include (1) we do not have an appropriate distribution for the probability difference; (2) a categorical variable may contain multiple levels so it is hard to compare all levels altogether between two clusters. Therefore, for a categorical variable, instead of comparing to a reference clustering group, we compare the distribution of the levels within each of the $G$ clusters to the overall marginal distribution of this categorical variable. If $\Delta_{mg}=0$, this indicates that the distribution of the categories within cluster $g$ is identical to the marginal distribution. If $\Delta_{mg}=1$, it implies that the distribution within cluster $g$ is different from the marginal distribution.

We define a Dirichlet distribution for categorical variable parameter $\bm{\theta}_{mg}$:

[TABLE]

where $\bm{\alpha}_{m\Delta_{0}}$ and $\bm{\alpha}_{\Delta_{1}}$ are hyper-parameters. Vector $\bm{\alpha}_{m\Delta_{0}}=(\alpha_{m\Delta_{0}1},\alpha_{m\Delta_{0}2},\dots,\alpha_{m\Delta_{0}L_{m}})^{T}$ is proportional to $\bm{\theta}_{mg}$ to make the corresponding prior center at marginal parameters of variable $m$, where $L_{m}$ is the number of categories of variable $m$; its elements are relatively larger numbers so that the prior is “spike” at marginal parameters of variable $m$ when $\Delta_{mg}=0$. Elements of vector $\bm{\alpha}_{\Delta_{1}}$ are all $1$ thus the pdf of this Dirichlet distribution is a constant. Therefore, this prior becomes a “slab” one when $\Delta_{mg}=1$.

We used Bernoulli distribution for the importance indicator variable $\Delta_{mg}$:

[TABLE]

where $p_{1m}$ and $p_{2m}$ are hyper-parameters of $\Delta_{mg}$.

Beta distribution for $p_{1m}$ and $p_{2m}$:

[TABLE]

[TABLE]

where $a_{p_{1}}$, $b_{p_{1}}$, $a_{p_{2}}$, $b_{p_{2}}$ are hyper-parameters.

A graphical depiction of our model is shown in Figure 2.

3.3 Posteriors

With specified prior distributions of all parameters, the corresponding posterior distribution then can be obtained with the form:

[TABLE]

We can derive posterior distribution for each parameter as below:

For cluster indicator matrix $\mathbf{Z}$:

[TABLE]

For cluster mixture proportion $\bm{\tau}$:

[TABLE]

For $A_{1m}$ in non-censored normal distributed variable $m$ within cluster $g$, where $m=1,2,\dots,q$:

[TABLE]

For mean difference $\mu_{mg}$, where $m=1,2,\dots,q$:

when $\Delta_{mg}=1$,

[TABLE]

when $\Delta_{mg}=0$,

[TABLE]

For precision parameter $\gamma_{m}$, where $m=1,2,\dots,q$:

[TABLE]

For hyper-parameter $\sigma_{\Delta_{0}}^{2}$ of $\mu_{mg}$:

[TABLE]

For parameter $\bm{\theta}_{mg}$, where $m=q+1,q+2,\dots,M$:

When $\Delta_{mg}=1$,

[TABLE]

When $\Delta_{mg}=0$,

[TABLE]

For importance indicator variable $\bm{\Delta}$:

When $m\in{1,2,\dots,q}$:

[TABLE]

When $m\in{q+1,q+2,\dots,M}$:

[TABLE]

For $p_{1m}$ for continuous variables, where $m=1,2,\dots,q$:

[TABLE]

For $p_{2m}$ for categorical variables, where $m=q+1,q+2,\dots,M$:

[TABLE]

We set different hyper-parameters for $p_{1m}$ and $p_{2m}$ to make variable selection more flexible. In this way, we can control the extent of shrinkage for continuous and categorical variables separately. In the next section, we will introduce how different choices of hyper-parameters could affect variable selection with more details.

For continuous variables with detection limit, we add one more sampling step: sample values below the lower detection limit and above the upper detection limit from the remainder of truncated normal distributions, respectively.

[TABLE]

[TABLE]

where $C_{Um}$ is the upper limit and $C_{Lm}$ is the lower limit for variable $m$; $i^{\prime}$ represents subjects whose real values of variable $m$ are lower than $C_{Lm}$; $i^{*}$ represents subjects whose real values of variable $m$ are higher than $C_{Um}$; $h$ is standard normal distribution; $\Phi$ is CDF of standard normal distribution.

3.4 Hyper-parameters

From the Bayesian framework of our proposed model, we will need to specify the values of the following hyper-parameters before starting the sampling procedures: $\bm{\delta}$ for $\bm{\tau}$; $\mu_{A}$ and $\sigma_{A}^{2}$ for $\bm{A_{1}}$; $a_{\Delta_{0}}$ and $b_{\Delta_{0}}$ for $\sigma_{\Delta_{0}}^{2}$; $\sigma_{\Delta_{1}}^{2}$ for $\bm{\mu}$; $\tilde{a}$ and $\tilde{b}$ for $\bm{\gamma}$; $\bm{\alpha}_{\Delta_{1}}$ and $\bm{\alpha}_{m\Delta_{0}}$ for $\bm{\theta}$; $a_{p_{1}}$ and $b_{p_{1}}$ for $p_{1m}$; and $a_{p_{2}}$ and $b_{p_{2}}$ for $p_{2m}$.

We can control the extent of shrinkage through controlling hyper-parameters $a_{p_{1}}$ and $b_{p_{1}}$ for continuous variables, and $a_{p_{2}}$ and $b_{p_{2}}$ for categorical variables. We can start with $a_{p_{1}}=b_{p_{1}}=a_{p_{2}}=b_{p_{2}}=1$, namely $Beta(1,1)$ to be the prior distribution of $p_{1m}$ and $p_{2m}$. If we prefer more shrinkage, then we could increase $b_{p_{1}}$ and $b_{p_{2}}$; otherwise we could decrease $b_{p_{1}}$ and $b_{p_{2}}$. In Figure 3 we show probability density functions of $Beta(1,1)$, $Beta(1,3)$ and $Beta(3,1)$ as an example. We can find that $Beta(1,1)$ (the grey line in Figure 3) is the same as $Unif(0,1)$, namely it is a non-informative flat prior. Under prior $Beta(1,3)$ (the red curve in Figure 3), it is more likely to sample smaller numbers for corresponding $\Delta_{mg}$ so that we can achieve more shrinkage. While under prior $Beta(3,1)$ (the blue curve in Figure 3), it is more likely to sample larger numbers for the corresponding $\Delta_{mg}$ so that we put less shrinkage. Assigning different hyper-parameters to $p_{1m}$ and $p_{2m}$ makes the shrinkage more flexible since we can control continuous and categorical variable separately. If we prefer the same shrinkage for all variables, we can simply set $a_{p_{1}}=a_{p_{2}}$ and $b_{p_{1}}=b_{p_{2}}$.

3.5 Estimation procedure

Since all parameters in the model have conjugate priors, we can apply Gibbs sampling to estimate unknown parameters and the procedure is summarized as follows.

Step 1: Set initial values for all parameters.

Step 2: For each variable in the dataset (e.g., age), update its distributional parameters (e.g., mean and standard deviation), namely $\bm{A_{1}}$, $\bm{\mu}$, and $\bm{\gamma}$ for a continuous variable or $\bm{\theta}$ for a categorical variable by sampling from the corresponding posterior distributions of parameters. If censored biomarker variables are encountered, update the censored variables before updating the distributional parameters.

Step 3: For each variable, update $\bm{\Delta}$ given the values of all other parameters.

Step 4: For each variable, update $p_{1m}$ or $p_{2m}$ given the values of all other parameters.

Step 5: Update $\sigma_{\Delta_{0}}^{2}$, $\mathbf{Z}$, and $\bm{\tau}$ given the values of all other parameters.

Step 6: Repeat Steps 2 to 5 for many iterations.

Step 7: Estimate each parameter using the mean value of its posterior distribution. Also, calculate the posterior probabilities of cluster membership for each subject and assign the cluster with the highest probability to that subject.

3.6 Label switching in Gibbs sampling

Label switching is a common problem in Gibbs sampling. It has been extensively studied and discussed in the literatures (Diebolt and Robert,, 1994; Stephens,, 2000; Frühwirth-Schnatter,, 2001; Marin et al.,, 2005; Papastamoulis and Iliopoulos,, 2010). Scrambling the cluster labels will not affect the likelihood function, the priors, and the posterior distributions. In Gibbs sampling of the cluster-specific parameters, there is a possibility of assigning wrong cluster labels to these parameters. This is called the label switching issue, which could cause biased estimation of the cluster-specific parameters. We adopt Stephen’s method (Stephens,, 2000; Papastamoulis,, 2015) to resolve this issue.

Let $\bm{P}^{(t)}$ be a $n\times G$ matrix with element $p^{(t)}_{ig}$ representing posterior probabilities of a subject $i$ belonging to cluster $g$ at MCMC iteration $t$; $t=1,2,\dots,T$, where $n$ is total number of subjects, $T$ is the total MCMC iteration times and $G$ is the total number of clusters. The basic idea of Stephen’s method is trying to permute the estimates from the MCMC iterations if needed so that those posterior probability matrix $\bm{P}^{(t)}$ agree each other for $t=1,2,\dots,T$. Specifically, Stephen’s Method minimizes the Kullback-Leibler divergence (also called relative entropy) between $\frac{1}{T}\sum_{t=1}^{T}\bm{P}^{(t)}$ and $\bm{P}^{(t)}$ for each $t$.

 

The relabeling algorithm is summarized as follows:

1. Assign an initial cluster labeling for all the MCMC iterations. We let $\bm{U}$ be a $T\times G$ matrix whose $t$th row represents the cluster labels for the $t$th iteration. For example, if $t=1$ and $G=3$, the cluster labels for the first MCMC iteration could be $(1,2,3)$, $(1,3,2)$, $(2,1,3)$, $(2,3,1)$, $(3,1,2)$ or $(3,2,1)$. To assign an initial clustering labeling is equivalent to assigning an initial values of the matrix $\bm{U}$.

2. Based on the current cluster labels, we calculate the mean posterior probability of cluster membership for all subjects across all iterations by averaging $\bm{P}^{(1)},\bm{P}^{(2)},\dots,\bm{P}^{(T)}$. That is, the mean posterior probability of the $g$th cluster membership for subject $i$ has the form:

[TABLE]

for $i=1,2,...,n$ and $g=1,2,...,G,$ where $p_{ig}^{(t)}$ is the $(i,g)$th element of the matrix $\bm{P}^{(t)}$.

3. Update each row of the matrix $\bm{U}$, the cluster labels of an iteration, by minimizing

[TABLE]

for $t=1,2,\dots,T$.

4. Repeat Steps 2 and 3 until matrix $\bm{U}$ is unchanged. Let $\bm{\Psi}$ be the $T\times G$ matrix of a parameter that needs to be estimated, where $T$ is the total number of MCMC iterations and $G$ is the total number of clusters. We will obtain the final parameter estimates for each iteration (each row of $\bm{\Psi}$) by permuting the cluster labels based on the values in the matrix $\bm{U}$.

4 Simulation studies

In this section we use simulations to evaluate the performance of the proposed Bayesian FMM relative to the other existing approaches. As reviewed in Section 2, most existing FMM frameworks with the ability of variable selection can only be applied on data consists of single variable type. Therefore, we chose methods that can handle mixed data: the HyDaP algorithm (Wang et al.,, 2019), Partition Around Medoids (PAM) with Gower distance (Gower,, 1971), K-prototypes (Huang,, 1998), regular FMM that uses the EM algorithm for estimation, and PAM with distance defined in factorial analysis of mixed data (FAMD) (Pagès,, 2014) as comparison methods. Clustering performance is evaluated by Adjusted Rand Index (ARI).

Assuming that there are 3 underlying true clusters with cluster sizes of 100, 100, and 100. In terms of variable importance, we designed 4 scenarios covering (1) both variable types contribute to clustering, (2) only continuous variables contribute to clustering, and (3) only categorical variables contribute to clustering. They are similar to the settings used in the paper that proposed the HyDaP algorithm. Details of these simulation settings can be found in Table 1. To evaluate clustering performance of different methods in the presence of censored biomarkers, we focused on the first simulation scenario and examined combinations of varying censoring proportions and varying numbers of censored variables. For each simulation scenario, 500 datasets were generated.

Our hyper-parameters used in simulations are as follows: $(0.33,0.33,0.33)^{T}$ was used for $\bm{\delta}$ assuming we do not have much informative on $\bm{\tau}$; 0 and 100 were used for $\mu_{A}$ and $\sigma_{A}^{2}$ assuming we do not have much informative on $A_{1m}$; 2 and 0.0001 were used for $a_{\Delta_{0}}$ and $b_{\Delta_{0}}$ as we believe $\sigma_{\Delta_{0}}^{2}$ should be a small number so that the corresponding prior for $\mu_{mg}$ is a “spike” one; 1000 was used for $\sigma_{\Delta_{1}}^{2}$ as we believe $\sigma_{\Delta_{1}}^{2}$ should be a large number so that the corresponding prior for $\mu_{mg}$ is a “slab” one; 2 and 1 were used for $\tilde{a}$ and $\tilde{b}$ assuming we do not have much informative on $\gamma_{m}$; 1 for $\bm{\alpha}_{\Delta_{1}}$ so that the corresponding prior for $\theta_{mg}$ is a “slab” one. We used $(10,10,10)^{T}$ for $\bm{\alpha}_{m\Delta_{0}}$, where $m=1,2,\dots,M$, since every categorical variable has 3 categories with true marginal probabilities $(0.33,0.33,0.33)^{T}$. We used 1 and 2 for $a_{p_{1}}$ and $b_{p_{1}}$, 1 and 2 for $a_{p_{2}}$ and $b_{p_{2}}$ as we would like some extent of shrinkage on both continuous and categorical variables.

4.1 Without existence of censored biomarker variables

Clustering performances in terms of ARI are shown in Table 2. The first row, Bayesian FMM, is our proposed Bayesian FMM with variable selection. In addition, the median and $2.5^{th}$ to $97.5^{th}$ percentile interval of variable weights obtained from our proposed Bayesian FMM with variable selection across all simulated datasets are shown in Table 3.

In simulation 1(a), we simulated a total of 5 variables: 4 continuous and 1 categorical. All except one continuous variable truly contribute to clustering. The sole categorical variable also contributes to clustering. Our Bayesian FMM with variable selection performed the best (ARI 1.00 [0.98, 1.00]) compared with other methods. Table 3 shows that the median weights of variables $x_{1}$, $x_{2}$, $x_{3}$ and $x_{5}$ are $1.00$ with $2.5^{th}$ and $97.5^{th}$ percentile interval $(1.00,1.00)$, indicating that these variables were found to be important for clustering. Meanwhile, variable $x_{4}$ had low weights (0.01 [0.01, 0.06]). These weights correspond to the true setting.

In simulation 1(b), we simulated a total of 10 variables: 7 continuous and 3 categorical. Six out of seven continuous variables truly contribute to clustering; two out of three categorical variables contribute to clustering. The difference between this setting and 1(a) is that in simulation 1(a), continuous variables with dominant influences exist while such variables do not exist in simulation 1(b). Our Bayesian FMM with variable selection performed the best (ARI 0.99 [0.96, 1.00]) compared with other methods. In Table 3, variables $x_{1}$, $x_{4}$ and $x_{7}$ had weights 0.05 (0.04, 0.07), 0.04 (0.03, 0.04), 0.05 (0.04, 0.07) indicating none of them has contribution to clustering; variables $x_{3}$, $x_{5}$, $x_{6}$ and $x_{9}$ had weights 0.10 (0.05, 0.19), 0.12 (0.06, 0.16), 0.11 (0.06, 0.17) and 0.12 (0.08, 0.22) respectively indicating that they have small contributions to clustering; variable $x_{2}$ had weight 0.30 (0.11, 0.41) indicating that it had slightly larger importance than $x_{3}$, $x_{5}$, $x_{6}$ and $x_{9}$. Variables $x_{8}$ and $x_{10}$ had weights 1.00 (1.00, 1.00) indicating they’re dominant variables that are highly relevant to clustering. These weights correctly reflect the true setting.

In simulation 2, we simulated a total of 5 variables: 4 continuous and 1 categorical. This setting is the same as setting 1(a) except that the sole categorical variable does not contribute to clustering. Our Bayesian FMM with variable selection performed the best (ARI 1.00 [0.98, 1.00]) compared with other methods. Variables $x_{1}$ to $x_{4}$ exhibited similar weights as those in simulation 1(a), but variable $x_{5}$ now had a low weight (0.12 [0.08, 0.23]), reflecting the underlying setting.

In simulation 3, we simulated a total of 7 variables: 4 continuous and 3 categorical. None of the continuous variables truly contribute to clustering. Two out of three categorical variables contribute to clustering. Our Bayesian FMM with variable selection performed similarly (ARI 0.74 [0.65, 0.83]) with the best performer: the HyDaP algorithm (ARI 0.75 [0.63, 0.85]). Table 3 shows that variables $x_{1}$ to $x_{4}$ all had weights close to 0 as in the true setting none of them are distinguishable across clusters. While categorical variable $x_{5}$ had very low weight as it has very small differences across clusters in true setting. Meanwhile, variables $x_{6}$ and $x_{7}$ which truly are associated with clustering had high weights (1.00 [1.00, 1.00]).

In summary, the clustering performance of our Bayesian FMM is always the top performer across all simulation scenarios. In addition, our Bayesian FMM with variable selection is able to provide quantitative variable importance which is more informative than the dichotomous information we obtained using the HyDaP algorithm, especially for simulation 1(b).

4.2 With existence of censored biomarker variables

We evaluated the clustering performance of our proposed Bayesian FMM with variable selection under simulation setting 1(a), where existing methods performed well. We compared the performances of proposed method with (1) our Bayesian FMM together with the naive method which uses half of lower detection limit to fill in undetected biomarker values; (2) the HyDaP algorithm together with the naive method; (3) PAM with Gower distance together with the naive method; and (4) K-prototypes together with the naive method. We designed 4 scenarios based on the original simulation setting 1(a): (1) only variable $x_{3}$ and $x_{4}$ are censored and both with 20% censoring; (2) only variable $x_{3}$ and $x_{4}$ are censored and both with 50% censoring; (1) all continuous variables are censored and each with 20% censoring; (2) all continuous variables are censored and each with 50% censoring. Clustering results are shown in Table 4 and obtained variable weights using proposed Bayesian FMM with variable selection are shown in Table 5.

In Table 4 we can observe that except the last column, our Bayesian FMM with variable selection performed the best. When all the continuous variables have very high censoring proportions, our method performs typically well but with higher variability from simulation to simulation (ARI 0.92 [0.42, 0.98]). Using Bayesian FMM with naive method to fill in censored values instead of using our embedded sampling approach, the performance is worse. Meanwhile, performance of the HyDaP algorithm was consistently good, though not the best, across all scenarios. When all the continuous variables have large censoring proportions, it can still yield high ARI with narrower percentile interval (ARI 0.85 [0.75, 0.94]). PAM with Gower distance performed poorly in all scenarios. Performance of K-prototypes is satisfactory when only a few continuous variables are censored, but becomes worse when all continuous variables are censored. In terms of variable weights, Table 5 indicates that our Bayesian FMM with embedded sampling approach yields weights that reflect the true setting for all scenarios.

In summary, when most continuous variables are censored with high probability, we suggest using the HyDaP algorithm since it would provide robust performance under this scenario. In all other scenarios, our proposed Bayesian FMM with variable selection is preferred as it would yield the best results.

5 Application

We analyzed the SENECA data and assumed that the data contains three latent clusters, as shown in  (Wang et al.,, 2019). Variables used in SENECA are listed in the appendix.

Due to the assumption of conditional independence, namely variables are independent with each other conditional on cluster labels, we dropped variables heart rate, blood sodium level (Na), aspartate aminotransferase (AST), and blood urea nitrogen (BUN) because of their high correlations to other variables in the dataset. The goals of our analysis are to obtain quantitative weights for all variables in the analytic dataset and to assess whether the identified clusters are related to those identified by the HyDaP algorithm.

Our hyper-parameters used in the analysis are as follows: $\bm{\delta}=(0.33,0.33,0.33)^{T}$ was used assuming that we do not have much information on $\bm{\tau}$. We set $\mu_{A}=0$ and $\sigma_{A}^{2}=1,000$ assuming that we do not have much information on $A_{1m}$. To set $a_{\Delta_{0}}=2$ and $b_{\Delta_{0}}=0.005$ as we believe that $\sigma_{\Delta_{0}}^{2}$ should be a small number so that the corresponding prior for $\mu_{mg}$ is a spike one. We set $\sigma_{\Delta_{1}}^{2}=100$ as we believe that $\sigma_{\Delta_{1}}^{2}$ should be a large number so that the corresponding prior for $\mu_{mg}$ is a slab one. We also let $\tilde{a}=2$ and $\tilde{b}=1$ assuming that we do not have much information on $\gamma_{m}$. $\bm{\alpha}_{\Delta_{1}}=1$ so that the corresponding prior for $\theta_{mg}$ is a slab one. We set $\bm{\alpha}_{\Delta_{0}}=(10.2,10)^{T}$ for gender and $\bm{\alpha}_{\Delta_{0}}=(15.5,2.4,2.1)^{T}$ for race based on their marginal distributions. We would like some shrinkage on both continuous and categorical variables, so we set $a_{p_{1}}=1$, $b_{p_{1}}=2$, $a_{p_{2}}=1$, and $b_{p_{2}}=2$.

By applying our proposed Bayesian FMM, we obtained 3 clusters with sample sizes $9,983$, $7,447$, and $2,759$, respectively. Table 6 shows these results cross tabulating with the 4 clusters that were identified using the HyDaP algorithm. Cluster 1 and Cluster 3 identified by the Bayesian FMM are similar to Cluster 1 and Cluster 4 identified by the HyDaP algorithm, respectively. Cluster 2 of the Bayesian FMM is similar to Clusters 2 and 3 altogether from the HyDaP algorithm. The results show that the 3 clusters identified by using the proposed Bayesian FMM are consistent with those identified from the HyDaP algorithm.

Table 7 summarizes the clustering weight of each variable obtained from the Bayesian FMM. None of the variables had driving influences on the final clustering results as their weights were low. This finding is consistent with what we found using the HyDaP algorithm. An advantage of using our proposed Bayesian FMM in ranking variable importance is its ability to provide a real-valued weight for each variable, not just dichotomizing certain variables as important or not. In SENECA, we found that ESR, troponin, lactate, sex, albumin, bicarbonate, GCS, and INR had relatively high weights, which mean that these variables contributed more in forming the 3 clusters than other variables did. Such information provides more details than the general importance obtained from the HyDaP algorithm. This is especially useful for clinical data which usually do not contain any dominant variables in terms of clustering.

Table 8 and Figure 4 show the distributions of clinical endpoints across the 3 clusters obtained by Bayesian FMM with variable selection. We can observe that Cluster 1 has the lowest proportion for all clinical endpoints while Cluster 2 has the second lowest ones. Cluster 3 has the highest proportions. We observe that these also correspond to the results obtained in the HyDap algorithm but our method is able to provide quantitative weights of all involved variables so that we have better knowledge about their importance and better interpretation of clustering results.

6 Discussion

Clustering has received a lot of attention and been applied in various areas these days. However, clustering methods that can handle mixed types of variables (both continuous and categorical) are still limited. Finite mixture model is a branch of clustering methods that is able cluster mixed data. But most existing FMM frameworks with variable selection are limited to single data type. Therefore, the goal of this paper is to develop FMM that can cluster data with mixed types of variables and perform variable selection in clustering.

We proposed a Bayesian FMM that can simultaneously cluster variables with mixed types, calculate variable weights, as well as handle censored biomarker variables. In this method we apply a Bayesian framework in order to bypass the limitations in the EM algorithm which is the standard estimation method for a FMM. In addition to identifying clusters, our model can provide real-valued variable weights which are more informative than a dichotomy of a variable being important vs. not in clustering. In addition, our proposed method is able to handle censored biomarker variables through recovering their underlying distributions. For a dataset of variables with mixed types but without censoring, our proposed Bayesian FMM performs better than the HyDaP algorithm and other existing methods across various simulation settings. If censored variables exist in the data, proposed Bayesian FMM with embedded sampling approach performs better than other clustering algorithms with naive fill-in methods. However, when variables have high censoring proportions, the proposed Bayesian FMM with embedded sampling approach may not consistently outperform other approaches. In this case, the HyDaP algorithm with ad-hoc imputations provides better and robust results. Under all scenarios, our proposed Bayesian FMM with variable selection is able to provide reasonable weights of all variables.

Users of the proposed Bayesian FMM model need to be aware of certain limitations. Same as finite mixture models, the proposed model also has the disadvantage of unverifiable distributional assumptions. Besides, the model assumes conditional independence (i.e., variables are independent conditional on cluster membership) so it may not perform well when variables are subject to within-cluster correlations. In computations, because that we need to specify the values for all hyper-parameters before running the algorithm, the computation time for the proposed model is usually longer than that for the EM algorithm.

Future work can be done on the top of the development of this paper. For the Bayesian FMM with variable selection, the framework can be extended to account for more distributions. Approaches of handling censored variables can also be further developed.

Appendix A

Variables used in SENECA data analysis

Age

Gender: categorical variable; 2 levels (male/female)

Race: categorical variables; 3 levels (white/black/hispanic)

Maximum temperature within 6 hours of ER presentation

Maximum heart rate within 6 hours of ER presentation

Minimum systolic blood pressure within 6 hours of ER presentation

Maximum respiration rate within 6 hours of ER presentation

Maximum albumin within 6 hours of ER presentation

Maximum Cl within 6 hours of ER presentation

Maximum erythrocyte sedimentation rate (ESR) within 6 hours of ER presentation

Maximum hemoglobin within 6 hours of ER presentation

Maximum bicarbonate within 6 hours of ER presentation

Maximum Sodium within 6 hours of ER presentation

Minimum Glasgow Coma Scale (GCS) within 6 hours of ER presentation

Elixhauser Score

Maximum white blood cell within 6 hours of ER presentation

Maximum bands within 6 hours of ER presentation

Maximum creatinine within 6 hours of ER presentation

Maximum bilirubin within 6 hours of ER presentation

Maximum troponin within 6 hours of ER presentation

Maximum lactate within 6 hours of ER presentation

Maximum alanine aminotransferase (ALT) within 6 hours of ER presentation

Maximum aspartate aminotransferase (AST) within 6 hours of ER presentation

Maximum C-reactive protein within 6 hours of ER presentation

Maximum international normalized ratio (INR) within 6 hours of ER presentation

Maximum glucose within 6 hours of ER presentation

Maximum Platelets within 6 hours of ER presentation

Maximum blood urea nitrogen (BUN) within 6 hours of ER presentation

Oxygen saturation (SaO2)

Minimum partial pressure of oxygen (PaO2) within 6 hours of ER presentation