{'fraction_non_alphanumeric': 0.03819975612459816, 'fraction_numerical': 0.012127258618778405, 'mean_word_length': 4.299729761485137, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hidden Markov Models on Variable Blocks with a Modal Clustering Algorithm and Applications

Motivated by high-throughput single-cell cytometry data with applications to vaccine development and immunological research, we consider statistical clustering in large-scale data that contain multiple rare clusters. We propose a new hierarchical mixture model, namely Hidden Markov Model on Variable Blocks (HMM-VB), and a new mode search algorithm called Modal Baum-Welch (MBW) for efficient clustering. Exploiting the widely accepted chain-like dependence among groups of variables in the cytometry data, we propose to treat the hierarchy of variable groups as a figurative time line and employ a HMM-type model, namely HMM-VB. We also propose to use mode-based clustering, aka modal clustering, and overcome the exponential computational complexity by MBW. In a series of experiments on simulated data HMM-VB and MBW have better performance than existing methods. We also apply our method to identify rare cell subsets in cytometry data and examine its strengths and limitations.


Introduction
Multivariate mixture model is arguably the most widely used statistical method for clusters identification. In this paper, we develop a novel hierarchical mixture modeling approach aimed at large-scale moderately high dimensional data sets that contain multiple rare clusters. The work is motivated by single-cell studies with applications to vaccine development and immunological research. The conventional mixture model approach is not up to the challenge of capturing clusters with highly unbalanced sizes even though the dimensionality of such data is in the tens. In order to model the data accurately without missing the rare clusters, a large number of mixture components is needed. On the other hand, the curse of dimensionality prevents the use of many components.
Furthermore, the increase of components results in higher computational complexity. To resolve this dilemma, we leverage a natural chain-like dependence among groups of variables in such data, a key fact that biologists and clinicians have been utilizing in exploratory studies. Existing statistical modeling techniques do not exploit this important characteristic. Our method is an effort to bridge the gap.

Background on Single-cell Cytometry Data Analysis
Interrogation of cell population heterogeneity has been made possible by the recent advancements in single-cell cytometry technologies (Perfetto et al. 2004;Bandura et al. 2009;Maecker et al. 2012;Chattopadhyay et al. 2014;Spitzer and Nolan 2016). For example, current high-throughput flow cytometry experiments measure 10 ∼ 20 parameters (cell markers/variables) routinely, including both phenotypic and functional markers, on a large number of single cells (hundreds of thousands to several millions). The current mass cytometry (CyTOF) can now measure up to 50 parameters at a single cell level.
A key first step to analyze this wealth of data is to partition the data (cells) from a sample (typically blood or tissue) into clusters based on the measured cell markers. The identified clusters are usually referred to as (cell) subsets. In most studies, the sample sizes are large, reaching to several millions. In addition, cell subsets of interest are typically in low frequencies (e.g., ∼ 0.01% of total cells). There is a need for detecting cell heterogeneity, especially very low frequency cell subsets for subsequent analysis. For example, association studies which help understand the association between cellular heterogeneity and disease progression (Ciuffreda et al. 2008;Lin et al. 2015a;Seshadri et al. 2015;Corey et al. 2015).
Traditionally, the cell subsets are identified by a manual gating strategy. Specifically, a manual sequential process that visually demarcates cells in bounded regions (called gates) on histogram or 2-D scatter plot projections. Figure 1 provides a simple illustration for the manual gating analysis on one flow cytometry data. Red lines are the gates, and cells within the region defined by the gates are identified as a specific cell subset. For example, to discriminate CD4+ T cells, which is one cell subset, a sequence of subsetting procedures is performed. Two physical markers, Forward (FSC-A) and side (SSC-A) light scatter, are first used to distinguish lymphocytes from all the live cells.
Lymphocytes can then be further partitioned based on 3 fluorescence parameters: CD3, CD4 and CD8 cell-surface markers. CD4+ T cells are the subclass of lymphocytes having high values of CD3 and CD4 but low value of CD8. Within CD4+ T-cell populations, additional functional markers such as intracellular makers (IL2 and IFNg) can further distinguish many functionally different CD4+ T-cell subsets. The sequence of markers to use is called gating hierarchy, which is determined by expert knowledge. The exact shape and location of the gates are manually drawn. The gates are typically used to dichotomize the continuous marker expressions into binary value: positive and negative. Therefore, the manual gating analysis is subjective and hard to reproduce. The red lines on each subplots are called gates. Cells within the red lines are the subset of interest, which is subsetted and projected on the next subplot in the sequence. The percentages are the frequencies of the identified cell subsets relative to the total number of cells.
It is evident that the manual gating analysis is essentially a hierarchical clustering procedure.
The hierarchy is formed on a chain of variable groups. More specifically, consider the sequential visualization process, any move one step further means a new group of variables are examined and used to select a more refined subset of data from an existing set chosen based on the previous groups of variables. By using the well-established gating hierarchy that projects cells on lower dimensions, each subplot provides a finer resolution of the cellular heterogeneity. In other words, the variables are grouped and ordered based on expert knowledge, and then the data are zoomed in by progressively putting constraints on more and more variable groups.

Automated Clustering Using Mixture Models
To overcome the limitations of the manual gating strategy, probabilistic clustering based on statistical mixture models (Boedigheimer and Ferbas 2008;Lo et al. 2008;Chan et al. 2008;Pyne et al. 2009;Aghaeepour et al. 2013;Lin et al. 2016) have been established to identify cell subsets objectively and automatically. One advantage of the mixture model is that the goodness of fit to any data sample can be improved by increasing the number of mixture components. The simplest approach to identify cell subsets based on a mixture model is to assign each mixture component to an individual cell subset. However, cell subsets can have any arbitrary shape. The parametric distribution of each mixture component is often inadequate to capture different shapes of the subsets. Various strategies have been proposed to merge multiple mixture components so that an individual cell subset can be better modeled (Pyne et al. 2009;Finak et al. 2009;Chan et al. 2010;Aghaeepour et al. 2011;Lin et al. 2016). In the statistical learning literature, the methodology of merging multiple mixture components into one cluster based on the modes of mixture density is called modal clustering (Li et al. 2007).
Although the typical mixture modeling approach has achieved some success, it is not powerful enough to detect clusters of very low probabilities. This is because such low probability mixture components are "consumed" by large background clusters in the data. One may naturally ask why not to increase the number of mixture components dramatically. However, to robustly estimate a mixture model, the number of components is limited by some maximum tolerable value depending on the data size (again curse of dimensionality). Our observation indicates that even the maximum tolerable value is not adequate for detecting those low frequency clusters. Moreover, a practical issue with using more components is the computational burden. It is thus imperative to design and implement effective and computationally accessible novel statistical models that can flexibly and robustly fit data on both high and low probability regions, in addition to the capability of handling non-Gaussian shape subsets.
Comparing the existing statistical mixture modeling approaches with the manual gating strategy, it is not hard to notice that the chain-like dependence among groups of variable (aka gating hierarchy), vital for the feasibility of the manual gating strategy, is not exploited by the aforemen-tioned mixture models. Clearly, if the dependence can be considered, more effective modeling of high dimensional data is possible. The recent work of Lin et al. (2013), which develops a hierarchical nonparametric Bayesian mixture model, explores along this direction. We see it as a prelude of our current work. The key idea of that paper is to partition cell features into two groups, also called subvectors. The first subvector is modeled by a Gaussian mixture model. Given the latent component identity of the first subvector, the second subvector is conditionally independent from the first subvector, and follows another Gaussian mixture distribution. In the conditional mixture distribution of the second subvector, the prior probabilities of the components depend on the latent component identity of the first subvector, while the parameters of individual normal components do not. This modeling feature is achieved via the use of the hierarchical Dirichlet process prior. The automatic clustering mechanism resembles the zoom-in on data in a hierarchical way by manual gating. The latent component identity of the first subvector determines the first level clustering. Clusters obtained from the first level are then further divided based on the latent component identity of the second subvector. It is found that this two-subvector model is more efficient and accurate in rare subset identification than the traditional generative mixture models.

Our Contributions
The next logical extension of Lin et al. (2013) is to accommodate more than two groups of variables, similarly as in the manual gating where a sequence of variable groups are visited. Albeit obvious in concept, there are technical hurdles seemingly impossible to conquer by the approach of hierarchical nonparametric Bayesian. First, when the number of variable groups grows, the computational complexity for estimating the latent component identities grows exponentially. The Markov chain Monte Carlo (MCMC) simulation for the two-subvector model is already highly intensive. There is also the classic issue of label switching when MCMC is applied to estimate the mixture models (e.g., Richardson and Green (1997); Celeux et al. (2000); Stephens (2000)). Practically, the existing method is already hitting a wall just to go beyond two subvectors. Second, a more subtle issue is the exponential growth of combinations of the latent component identities. Suppose there are k variable groups and each group is modeled by a mixture of M components. Then there are M k combinations of components over all the variables, each can correspond to one cluster. This exponential increase in the number of clusters is unrealistic, raising doubt on the legitimacy of the result. Even if we perform modal clustering, which could dramatically reduce the number of clusters, such algorithm to find the modes again has an exponential complexity.
In this paper, we develop a general statistical model, namely Hidden Markov Model on Variable Blocks (HMM-VB), and overcome both of the aforementioned difficulties by deriving new algorithms for model estimation and mode identification. Although our new model is not a HMM in the conventional sense, a point elaborated in Section 3, the latent graphical dependence structure for the variables is essentially the same as that of the sequential vector data in HMM. We thus put our model in the light of HMM and leverage the powerful computational tools for HMM. To the best of our knowledge, our work is the first to exploit the hierarchical dependence structure among variable groups by HMM.
The estimation algorithm for HMM can resolve the estimation difficulty encountered by the MCMC approach. One natural idea to use HMM-VB to cluster data is to find the underlying sequence of states with maximum posterior. However, the number of possible state sequences is exponential in the sequence length, in practice, often much larger than the data size. To tackle this issue, we cast the HMM-VB back as a Gaussian mixture model and use modal clustering to merge the components in Gaussian mixture model. However, the number of components in the equivalent Gaussian mixture model equals the number of possible state sequences. As a result, a straightforward application of a mode-seeking algorithm, such as Modal EM (MEM) (Li et al. 2007), still faces exponential computational complexity. We discover a new numerical algorithm with linear complexity for finding modes of a HMM-VB, which is proven to be equivalent to MEM. This is one key for the success of clustering based on HMM-VB.
Our experiments on both simulated and real data sets show that the approach of HMM-VB is effective in discovering rare clusters and is robust to non-Gaussian shaped clusters. It naturally incorporates domain knowledge on the dependence structure of the variables. We also explore the potential of HMM-VB as a general tool for parsimonious mixture modeling. In particular, we test HMM-VB on data not complying to the chain-like dependence of variable blocks and comparisons with the general mixture models are performed.
The rest of the paper is organized as follows. In Section 2, we introduce notations and overview existing techniques most relevant to our proposed methods. In Section 3, we describe the proposed model and develop efficient algorithms for model fitting and for performing modal clustering. Section 4 presents the results of both simulation studies and the application to large-scale cytometry data. Comparisons are made with some competing models. Finally, we summarize and discuss our contributions in Section 5.

Preliminaries
The finite Gaussian mixture model (GMM) is commonly used for clustering. For random variable X ∈ R d , where d is the number of dimensions, a GMM with M components has density function: where π k is the mixture component prior probability and φ(· | θ k ) is the multivariate normal density parameterized by θ k = (µ k , Σ k ), µ k being the d-dimensional mean vector and Σ k the d × d covariance matrix. The collection of parameters θ = (π 1 , ..., π M , θ 1 , ..., θ M ).
Specifically, Z = k if, and only if, X comes from component k. Z is also called the component identity of X. To perform clustering, the usual approach is to compute the posterior probability P (Z = k | X = x), and assign X to the cluster k with the maximum posterior. However, this approach is inadequate to model clusters with arbitrary shapes and cannot ensure that the clusters are reasonably separated. One major idea explored in the literature is to merge multiple mixture components for a better and more flexible representation of an individual cluster. We refer to (Melnykov and Maitra 2010) for a thorough review on clustering based on finite mixture models.

MEM for Gaussian Mixtures and Modal Clustering
The Modal EM (MEM) algorithm developed by Li et al. (2007) performs efficient merging of mixture components. It resembles the expectation-maximization (EM) algorithm (Dempster et al. 1977), as reflected by the name "modal EM". However, the optimization objective of MEM is the local maxima over x for f (x | θ) under a given θ, while EM is to find local maxima over θ for is the density of the kth mixture component. Staring from any initial value x [0] , MEM generates a local maximum of the mixture density by performing the following two iterative steps: (1) At iteration r, 2. M-step: solve The computational efficiency of MEM enabled the development of a new clustering approach by Li et al. (2007), referred to as modal clustering. In Li et al. (2007), a nonparametric Gaussian kernel density estimate is used; and initialized by each data point, MEM is applied to find the mode associated with every point. Data points associated to the same mode are classified to be in the same cluster. In Lee and Li (2012), a natural extension is developed for the general GMM.
To save computation, instead of applying MEM to every data point, it is applied to the means of the mixture components. Components (hence data points) with mean vectors associated to the same mode are merged into one cluster. Whether a point-wise mode association or a componentwise mode association is preferred depends on the nature of the application and the computational resource. In practice, the difference in the clustering results we have observed is quite small. For the finite mixture model, the number of components is usually much smaller than the data size. Hence component-wise mode association saves computation. In some situation, the number of mixture components can be enormous (much more than the data size), as is the case for our hierarchical mixture model to be introduced in Section 3. Clearly, a point-wise mode association would be preferred. However, finding a mode starting from every data point can still be costly for big data.
As a result, we take the scheme of component-wise mode search, but we do not find the associated mode for every component mean. Instead, we first identify all the components that have been chosen by at least one data point according to the maximum a posteriori criterion. The number of such components is upper bounded by the data size and is often much smaller. The modes are then computed with initialization set by each of these chosen component means.
Under the framework of modal clustering, the role of a mixture component is primarily for achieving good density estimation. We no longer rely on a component to represent a cluster. This provides flexibility in modeling. We can allow components to overlap substantially for the sake of more accurate modeling of the density. When it comes to partitioning data, mode association ensures that different groups of data are sufficiently separated. On a practical side, the choice of the number of components is not as crucial as the conventional mixture model based clustering, where that number determines how many clusters are generated. In fact, in the modal clustering algorithm in Li et al. (2007), the Gaussian kernel density is used, where the number of components equals the data size. More elaborated discussions on the advantage of modal clustering and comparisons with the usual mixture model clustering are provided in that paper. Interestingly, the modal clustering technique has been applied to cytometric data analysis by Ray and Pyne (2012).
A computational framework called flowScape was developed to mimic the analytical actions of human experts in the process of manual gating (Ray and Pyne 2012). The purpose of that work is to enhance the manual gating process, while our work belongs to the school of automated clustering.
The flexibility provided by modal clustering for modeling data is precisely what we need for the applications we consider. As will be explained in Section 3, we exploit a HMM-type model which can be cast as a mixture model with an enormous number of components, even exceeding the data size. This causes no difficulty in clustering via mode association.

Hidden Markov Model
k, l ∈ S. Let the initial probabilities of states be denoted by π k = P (s 1 = k), k ∈ S.
2. Given the hidden state s t , the observation x t is conditionally independent from s t and x t for any t = t; and the distribution of x t given s t depends on s t , but not t. Denote the conditional In summary: The parameters to be estimated in a HMM are the transition probabilities: a k,l , k, l = 1, ..., M , the initial probabilities: π k , k = 1, ..., M , and µ k , Σ k for each state k = 1, ..., M . HMM is usually estimated by the EM algorithm. However, because the cardinality ofS grows exponentially with the sequence length, the computational complexity of a direct application of EM is of exponential complexity. This technical hurdle was overcome by the Baum-Welch (BW) algorithm that achieves without suffering optimality a complexity linear in the sequence length and quadratic in the number of states. The BW algorithm, a special instance of EM, was developed in the 1960's before the general EM algorithm was developed in the 1970's. As a result, we still call the estimation algorithm Baum-Welch, following the convention of the literature on HMM.
Under a set of parameters, let L k (t) be the conditional probability of being in state k at position t given the entire observed sequence x = {x 1 , x 2 , ..., x T }. Let I(·) be the indicator function that equals 1 when the argument is true and 0 otherwise. Then Let H k,l (t) be the conditional probability of being in state k at position t and being in state l at position t + 1, i.e., seeing a transition from k to l at t, given the entire observed sequence x.
Since s ∈S and the |S| = M T , it is infeasible to compute L k (t) and H k,l (t) by the above equations directly. As part of the BW algorithm, the forwardbackward algorithm is used to compute L k (t) and H k,l (t) efficiently. The amount of computation needed is at the order of M 2 T ; and memory required is at the order of M T . We present the forwardbackward algorithm in Appendix A. For the current discussion, suppose the quantities have been solved.

M-step: Update parameters by
The initial probabilities of states π k are often manually determined. We can also estimate them by π k ∝ T t=1 L k (t), subject to M k=1 π k = 1, or π k ∝ L k (1). So far we present the BW algorithm in the case of estimation based on a single sequence. The extension to estimation from multiple sequences is relatively easy, which we describe in Appendix A. In our work, as will be explained in the next section, we always handle multiple sequences. For a detailed exposure to HMM, we refer to (Young et al. 1997).

Hidden Markov Model on Variable Blocks
Consider X ∈ R d . We propose a dependence structure among groups of variables that in essence matches the statistical dependence among observations at different time spots assumed by HMM. We call the model Hidden Markov Model on Variable Blocks (HMM-VB). Suppose the variables are divided into blocks t = 1, 2, ..., T , where T is the total number of blocks. Let the number of In a nutshell, HMM-VB treats the block index t as "time". The variables in the same block are treated as multivariate data and their joint density is assumed to be a GMM. The Gaussian component identity of a variable block acts as the underlying state for this block. To capture the dependence between the variable blocks, the states are assumed to follow a Markov chain. As discussed in Section 1.3, such a dependence structure is motivated by the gating hierarchy among groups of variables, which has been well accepted by researchers on cytometry data analysis.
Without loss of generality, suppose we have ordered the variables in a way such that the first d 1 variables are in block 1, the next d 2 variables are in block 2, and so on. Let variable blocks where M t is the number of mixture components for variable block x (t) , t = 1, ..., T . Let the set of all possible sequences beŜ = S 1 × S 2 · · · × S T . |Ŝ| = T t=1 M t . We assume: 1. {s 1 , s 2 , ..., s T } follow a Markov chain. Let π k = P (s 1 = k), k ∈ S 1 . Let the transition probability matrix A t = (a (t) k,l ) between s t and s t+1 be defined by a (t) k,l = P (s t+1 = l|s t = k), k ∈ S t , l ∈ S t+1 .
2. Given s t , x (t) is conditionally independent from the other s t and x (t ) , t = t. We also assume that given s t = k, the conditional density of Denote s = {s 1 , ..., s T }. To summarize, the density of HMM-VB is given by Remarks on comparison with the conventional HMM: 1. The variable blocks x (t) 's are not from the same vector space. Hence, the parameters of the distribution of x (t) given s t = k depend not only on k but also on t.
2. The underlying Markov chain for {s 1 , ..., s T } is not time invariant. In fact, the the state space S t varies with t.

Maximum Likelihood Estimation
We derive the corresponding BW algorithm for HMM-VB. It is an instance of the EM algorithm, similarly as the BW algorithm for the usual HMM. The proof of the BW algorithm as an equivalence of the EM algorithm for HMM-VB follows closely that for the usual HMM. We thus omit it here.
Clearly, it is not meaningful to estimate HMM-VB using a single sequence, aka, a single data point.
HMM-VB is after all a model for X ∈ R d . The imposed sequence structure is for the sake of modeling dependence among the variables, motivated by the nature of the data encountered in applications.
Let x = (x (1) , x (2) , ..., x (T ) ) ∈ R d (assume column-wise concatenation) be the full-dimensional data. Consider estimation of HMM-VB based on a data set {x 1 , x 2 , ..., x n }, x i ∈ R d , i = 1, ..., n. Use i ∈ R dt to denote the tth variable block of x i , and the notation x Making the dependence on x explicit in the notation, we define L k (x, t) and H k,l (x, t) similarly as in Eq. (4) and Eq. (5).
The above equations easily extend to the case of weighted sample points. It can occur in practice that each sample point is assigned with a weight. For instance, quantization is often used to reduce the data size significantly. Instead of the original data, one may use the quantized points, each of which can represent a different number of original points and hence is assigned with a weight proportional to that number. Suppose weight w i is assigned to sample x i . The E-step is not affected.
In the update of parameters, we can simply multiply w i in front of each summand appeared in the equations above.
The forward-backward algorithm for computing L k (x i , t) and H k,l (x i , t) is essentially the same as the forward-backward algorithm for the usual HMM. The fact that the variable blocks are not from the same vectors space and the state spaces vary with t does not cause any intrinsic difference.
Define the forward probability α k (x, t) as the joint probability of observing the first t variable blocks x (τ ) , τ = 1, ..., t, and being in state k at time t: This probability can be evaluated by the following recursive formula: Define the backward probability β k (x, t) as the conditional probability of observing the variable blocks after time t, x (τ ) , τ = t + 1, ..., T , given the state at block t is k: Set β k (x, T ) = 1 , for all k ∈ S T .
The backward probability can be evaluated using the following recursion: The probabilities L k (x, t) and H k,l (x, t) are solved by The normalizing factor P (x) = k∈St α k (x, t)β k (x, t) for any t (equation holds for any t).
To initialize the model, we design several schemes. In our experiments, models from different initializations are estimated and the one is chosen with maximum likelihood. In our baseline initialization scheme, k-means clustering is applied individually to each variable block using all the data instances. Based on the clustering result of k-means, we take every cluster as one mixture component and compute the sample mean and sample covariance matrix of data in that cluster. To reduce the sensitivity to the initial clustering result, we also compute the pooled common sample covariance matrix for the clusters. The initial covariance matrix of a component is then set to be a convex combination of the cluster-specific sample covariance and the common sample covariance.
The transition probabilities are always initialized to be uniform. Under the second initialization scheme, we randomly sample a subset from the whole data and apply the baseline initialization to the subset. Under the third initialization scheme, we randomly pick a subset from the data and treat points in this subset as the cluster centroids of k-means. These centroids will induce a cluster partition of the whole data, based on which we initialize the component means and covariance matrices likewise, as in the baseline method. Both the second and the third initialization schemes are repeated several times with different random start.

Modal Baum-Welch Algorithm
HMM-VB can be viewed as a special case of a GMM where each component of the GMM corresponds to a particular sequence of states s = {s 1 , ..., s T }, that is, a combination of states for all the variable blocks. We call this equivalent GMM the GMM mapped from HMM-VB. Each component is a Gaussian distribution with mean µ s = (µ (1) s 2 , ..., µ s T ) (column-wise stack of vectors) and a covariance matrix, denoted by Σ s , that is block diagonal. The tth diagonal block in Σ s is Σ We can thus readily apply the modal clustering framework for GMM to data modeled by HMM-VB. However, as the number of possible sequences grows exponentially with the number of variable blocks T , a direct application is computationally infeasible. We discover that because of the block diagonal structure of the covariance matrix of the GMM mapped from HMM-VB, we can in fact avoid computing the posterior of x belonging to each component (exponentially many of them!). Instead, we only need L k (x, t) for all k and t when updating x in the M-step of MEM. Because the BW algorithm solves L k (x, t) at a complexity linear in T , we can achieve linear complexity for MEM on HMM-VB as well. We call this new algorithm Modal Baum-Welch (MBW) Algorithm.
For brevity of notation, we use x (t),r to denote the value of the tth variable block at iteration r.
And x [r] = (x (1),r , x (2),r , ..., x (T ),r ) is the concatenated full vector at iteration r. The equivalence of MBW and the Modal EM algorithm is ensured by Theorem 1 below, which is proved in Appendix B.

Theorem 1. For the mapped GMM from a HMM-VB, suppose the solution of the M-step in the MEM algorithm given by Eq. (3) is divided into blocks
. Then The MBW algorithm iterates the following two steps: 1. E-step: Compute L k (x [r] , t), for k ∈ S t , t = 1, ..., T .
2. M-step: The clustering method based on MBW is straightforward. We first find the state sequence s Since different x i 's may yield the same sequence, we then identify the collection of distinct s * i . For each distinct sequence, say s . Use µ s * i as an initialization for MBW to find the mode associated with it. If µ s * i and µ s * j are brought to the same mode by MBW, the data point x i and x j are put in the same cluster. When |Ŝ| is very large, the number of different s * i 's can become close to the data size. Hence, the computation we can save by seeking modes starting from µ s * i 's instead of the original data diminishes. As a result, in such cases, we recommend seeking modes directly from the original data. In fact, in the experiment described in Section 4.2, the modes are found in this way.

Experiments
In this section, we present experiments on several simulated data sets (4.1) and one CyTOF data (4.2). For each data set, the BW algorithm was run repeatedly starting from multiple initial models. In Section 3.1, the different ways of initialization are described. Among the final models, we choose the one yielding the maximum likelihood.

Simulation with Various Types of Variable Block Structures
We first conduct simulation studies to examine the effectiveness of HMM-VB at capturing small clusters and its robustness against different parameter settings.

Two variable blocks
Using the same set-up as in Lin et al. (2013), a sample of size 10, 000 with dimension d = 8 is drawn from a hierarchical mixture model. There are two variable blocks. Following the notations in the previous section, x i is divided into two variable blocks, also called subvectors, x i was generated from, as in the assumption of HMM-VB. The data is designed to have a distinct cluster after standardization (subtract mean and divided by the standard error). In particular, the standardized data have a well-separated region that the five dimensions x 2 , x 3 , x 6 , x 7 , x 8 are of high positive values, and the rest are negative. The subset of interest with size 140 is indicated in red in Figure 2. Lin et al. (2013) showed the standard GMM failed completely to identify the subset of interest.
Here we contrast the hierarchical Bayesian mixture model approach (Lin et al. 2013) with results from analysis using the new HMM-VB. To fit HMM-VB, we only need to specify the number of mixture components for each variable block, in this example, M 1 and M 2 . If casting as a GMM, the HMM-VB has M 1 · M 2 components for the full dimensional data. Model selection using BIC is conducted to select the optimal M 1 and M 2 . Summaries on various model specifications are listed in Table 1. The model with M 1 = 7 and M 2 = 10 has the lowest BIC, thus it is selected as the optimal model.
Recall that we search modes by MBW starting from the mean of every component that has been chosen by a data point according to the maximum a posteriori rule. We call the components that have not been chosen by any data point the empty mixture components, while the others the nonempty. Table 1 shows that even though the total number of nonempty mixture components increases with M 1 and M 2 , the number of clusters after modal clustering remains relatively stable. Figure 3, HMM-VB with M 1 = 7 and M 2 = 10 correctly identifies 135 observations out of the 140 target sample. It yields similar clustering performance as the hierarchical mixture model developed in Lin et al. (2013), while the hierarchical mixture model can only accommodate two variable blocks.

Multiple Variable Blocks
A sample of size 10, 000 with d = 13 dimensions is drawn such that the first 5 dimensions are generated from a single multivariate normal distribution with zero mean and diagonal covariance matrix 3I. The rest 8 dimensions, divided into 2 variable blocks, are generated according to Section 4.1.1. In summary, the simulated data here have 3 variable blocks with d 1 = 5, d 2 = 5, d 3 = 3. Similarly, the data contain a well-separated region such that after data standardization, the five dimensions Next, we compare the clustering result with the generic GMM (see Eq. (1)), which can be taken as a HMM-VB with one variable block. We also apply modal clustering to GMM. Similarly, BIC is computed for models with different numbers of mixture components M . Results are listed in Table 3. BIC suggests the 5−component GMM is the optimal model which identified in total 4 clusters. The subset of interest is completely masked as is shown in Figure 5 middle panel.
The 5−component GMM is unable to identify the correct subset region. We then dramatically increased the number of mixture components to 100. The 100−component GMM identified in total 69 clusters, significantly larger than the number of clusters determined by the optimal HMM-VB.

No Information on Variable Blocks
Lastly, we study the clustering performance of HMM-VB when there is no clear-cut variable block structures. A sample of size 10, 000 with d = 10 dimensions is drawn from a 50−component GMM.
The component priors for the first 6 components are π 1:6 = (0.0025, 0.005, 0.1, 0.1, 0.05, 0.03), and the other priors are equal, that is, π s = 0.7125/44 ≈ 0.016, s = 7, ..., 50. Furthermore, we make the first two rare components to represent two distinct clusters by setting the first mean vector µ 1 to a vector of 5's and µ 2 to a vector of −5's. The other mean vectors µ 3:50 are independently generated from a multivariate normal distribution with 0 mean and identity covariance matrix.
The first two covariance matrices Σ 1:2 are both identity matrix. The other covariance matrices Σ 3:50 are independently generated from an inverse Wishart distribution with 15 degrees of freedom and diagonal scale matrix 5I. Figure 6 shows the scatter plots of the generated data on selected dimensions. The two components are masked by the background data in some dimensions. For example, projected to the 5th and 9th dimensions, the first component (red) cannot be completely separated from the background data (grey).
We first fit the data with HMM-VB using different model specifications. Specifically, we let Model 1 define a HMM-VB with 10 variable blocks such that each variable forms one block. The data is fitted hierarchically from the first variable to the last one, with M 1:10 = 10. Model 2 is defined similarly to Model 1, but we reverse the ordering of the variables such that x = (x 10 , x 9 , ..., x 1 ).
Lastly, we let Model 3 define a HMM-VB with 10 variable blocks such that we randomly reorder x = (x 10 , x 9 , x 6 , x 5 , x 4 , x 3 , x 2 , x 1 , x 8 , x 7 ). Table 4 shows the confusion matrices for the three models respectively. In particular, we assign all the background data that do not belong to the first two components to cluster 0. Table 4 suggests that HMM-VB is quite robust to the ordering of the  variables. All the three HMM-VB models can consistently identify the two rare components with high accuracy.  We see that even in the lack of the variable block structure, HMM-VB can still outperform GMM in identifying multiple rare clusters. This result suggests that HMM-VB is a strategy to form parsimonious versions of GMM which can be more effective than the common way of varying the number of components.

Study of CyTOF Data
In this section, we study the performance of HMM-VB on a higher dimensional data set obtained from CyTOF experiment (Becher et al. 2014). The particular data set that we analyze here is from mouse lung sample obtained from three C57BI6 wild-type mice and three Csf2rb −/− mice, which in total contains 46, 204 single cells with 39 measured cell markers. According to the gating hierarchy provided in Becher et al. (2014), it defines roughly 11 variable blocks, with maximum block size 8 and minimum block size 1. Becher et al. (2014) performed automated clustering on a projected (latent) 2−dimensional space by first using nonlinear dimension reduction technique on the original 39−dimensional space. However, it has been studied, e.g. Lin et al. (2015b), that dimension reduction generates a "cluttered" display, and this can prevent density estimation from accurately representing low probability regions.
We first standardize the data. The reason is that the moderately high-dimensional data has a nearly singular covariance matrix (that is, when a single Gaussian is fit). This prevents the direct use of GMM. It should be noted that HMM-VB encounters no difficulty in fitting the data because of the smaller dimensions of the individual variable blocks. In order to compare GMM and HMM-VB, we work on the standardized data through the rest of the section.
We fit the data using HMM-VB with specification as follows. There are 11 variable blocks. For the variable blocks i with dimension lower than 4, we set the corresponding number of mixture components M i = 5. For variable blocks j with dimension between 5 and 7, we set the corresponding number of mixture components M j = 10. For the other variable blocks, we set the corresponding number of mixture components to 15. After modal clustering, this model results in 825 clusters.
The average CPU time per model fitting is 11.4 min on an iMac with Intel i7 Processor at 3.0 GHz with 8GB memory. We next fit the data using GMM with M = 500. On the same iMac, the model fitting takes 284.2 min, which is 25 times longer than did HMM-VB. After modal clustering, it results in 316 clusters.
We compare the performance of the two models for clustering a relatively low probability region (6% of total cells), which is shown in Figure 7. Specifically, Figure  indicates that GMM has difficulty in estimating the structure of this moderately high-dimensional data.

Discussion
We develop and explore a novel hierarchical mixture model (HMM-VB) with the goal of automated clustering of large-scale data sets that contain multiple rare clusters. The method exploits the natural structure among groups of variables for more effective modeling and clustering analysis. Both the simulated and real data examples demonstrate the effectiveness of HMM-VB for identifying multiple rare clusters in large-scale data sets. The practical motivation is from the automated analysis of single-cell cytometry data. The standard visual gating approach suffers severe limitations including non-reproducibility and labor intensiveness. The cytometry field is increasingly interested in using more consistent automated statistical methods. Standard mixture models, however, are not powerful enough to detect clusters of very low probabilities when applied to large data sets. One key feature of the new model is by design the ability to quantify the inherent chain-like dependence among groups of variables for more effective clustering, especially in identifying rare clusters that are almost undetectable by existing mixture modeling approaches.
Technically, our clustering method integrates two new algorithms, one for fitting HMM-VB and the other for performing modal clustering, both necessary to make HMM-VB a practical tool for large-scale data analysis. In theory, the number of components grows exponentially with the num-ber of variable groups. Existing methods for estimation and mode identification thus bear exponential complexity, rendering them infeasible on real-world data. We derive and implement algorithms with linear complexity in the number of variable blocks for both tasks.
The development of more structured mixture models for clustering is necessary for large-scale data, especially single-cell data as the technologies advance and the number of markers that can be measured increases. However, due to the huge number of mixture components, many tiny clusters can be generated. The total number of points in those clusters can be a small fraction of the entire data. For instance, nearly half of the clusters identified in the analysis of CyTOF data are of size one.
In this situation, more sophisticated methods that measure the separability of these tiny clusters and other larger ones can be used for modal clustering (e.g., Lee and Li (2012)). The development of HMM-VB can also be used to determine the optimal variable groupings and orderings based on some model selection criterion, such as BIC. This can be useful for developing new cytometry experiment for providing objective optimal gating hierarchy.

Appendix A
We present the forward-backward algorithm for computing L k (t) and H k,l (t) efficiently. Define the forward probability α k (t) as the joint probability of observing the first t vectors x τ , τ = 1, ..., t, and being in state k at time t: α k (t) = P (x 1 , x 2 , ..., x t , s t = k).
This probability can be evaluated by the following recursive formula: α l (t − 1)a l,k , 1 < t ≤ T, 1 ≤ k ≤ M.
The probabilities L k (t) and H k,l (t) are solved by H k,l (t) = P (s t = k, s t+1 = l | u) = P (x, s t = k, s t+1 = l) P (x) = 1 P (x) α k (t)a k,l b l (x t+1 )β l (t + 1), where P (x) = M k=1 α k (t)β k (t). We now provide the Baum-Welch algorithm for the case of multiple sequences. For brevity, assume all the sequences are of length T . Denote the ith sequence by x i = {x i,1 , x i,2 , ..., x i,T }, i = 1, ..., n. In each iteration, we compute the forward and backward probabilities for each sequence separately in the same way as previously described. We also compute L k (t) and H k,l (t) separately for each sequence. As a general pattern of notations, we put a superscript (i) to indicate the quantities for the ith sequence.
3. Update the parameters using L .

Appendix B
We prove the equivalence of MBM and MEM for HMM-VB. It is clear that the density of HMM-VB in Eq. (6) is a mixture model if we take the state sequence s ∈Ŝ as the index for the mixture component. Each component is a Gaussian distribution with mean µ s = (µ (1) s 2 , ..., µ s T ) (column-wise stack of vectors) and a covariance matrix, denoted by Σ s , that is block diagonal. The tth diagonal block in Σ s is Σ (2) s 2 0 · · · 0 · · · · · · · · · · · · · · · If we apply MEM directly to HMM-VB and keep in mind that s is the index for the mixture component, we need to compute the posterior P (s | x) in the E-step and in the M-step. The computational hurdle is that the number of possible sequences s, that is, |Ŝ|, grows exponentially with T (assuming similar |S t |'s).
Because Σ s is block diagonal, we have Hence the tth variable block of x is given by st   , t = 1, 2, ..., T.
Let I(·) be the indicator function that equals 1 when the argument is true. Note that s∈Ŝ P (s | x) Σ (t) according to the definition of L k (x, t) in Eq. (8). By the same technique, we can show that Thus we have proved the equivalence of MBW and MEM for HMM-VB.