Bayesian nonparametric mixture inconsistency for the number  of components: How worried should we be in practice?

Yannis ChaumenyYannis Chaumeny and Johan van der Molen Moris contributed equally<EMAIL_ADDRESS>
Department of Statistics

Ecole Polytechnique Fédérale de Lausanne

Lausanne, Switzerland
Johan van der Molen Moris^∗<EMAIL_ADDRESS>
MRC Biostatistics Unit

University of Cambridge

Cambridge, UK
Anthony C. Davison<EMAIL_ADDRESS>
Department of Statistics

Ecole Polytechnique Fédérale de Lausanne

Lausanne, Switzerland
Paul D. W. Kirk<EMAIL_ADDRESS>
MRC Biostatistics Unit

University of Cambridge

Cambridge, UK

Kevin Murphy and Bernhard Schölkopf

We consider the Bayesian mixture of finite mixtures (MFMs) and Dirichlet process mixture (DPM) models for clustering. Recent asymptotic theory has established that DPMs overestimate the number of clusters for large samples and that estimators from both classes of models are inconsistent for the number of clusters under misspecification, but  the implications for finite sample analyses are unclear. The final reported estimate after fitting these models is often a single representative clustering obtained using an MCMC summarisation technique, but it is unknown how well such a summary estimates the number of clusters. Here we investigate these practical considerations through simulations and an application to gene expression data, and find that (i) DPMs overestimate the number of clusters even in finite samples, but only to a limited degree that may be correctable using appropriate summaries, and (ii) misspecification can lead to considerable overestimation of the number of clusters in both DPMs and MFMs, but results are nevertheless often still interpretable. We provide recommendations on MCMC summarisation and suggest that although the more appealing asymptotic properties of MFMs provide strong motivation to prefer them, results obtained using MFMs and DPMs are often very similar in practice.

Bayesian nonparametric analysis, Cluster analysis, Dirichlet process, Number of clusters, Mixture model.

§ INTRODUCTION

Many areas of applied statistics use mixture models for cluster analysis.
One advantage of this approach compared to heuristic counterparts is the possibility of inference on both individual cluster allocations and the entire clustering structure. Mixture models have been developed for a wide variety of data types, with applications to gene expression analysis [Medvedovic and Sivaganesan, 2002, McLachlan et al., 2002, de Souto et al., 2008, Onogi et al., 2011, Crook et al., 2018], economics [Alfò et al., 2008, Wedel and Desarbo, 2002, Frühwirth-Schnatter et al., 2004], cell type classification [Chan et al., 2008, Prabhakaran et al., 2016], community detection [Geng et al., 2019, Legramanti et al., 2021], climate science [Gaffney et al., 2007], astrophysics [Mukherjee et al., 1998], health surveys [Ni et al., 2020], demographics [Carmona et al., 2019] and  finance [Dias et al., 2015].
Reliable estimation of the number of clusters is essential for accurate prediction and decision-making—for example, precision medicine relies on identifying disease subtypes [Sørlie et al., 2001, Lanza and Rhoades, 2013]—and various approaches to this have been developed [Dudoit and Fridlyand, 2002, James et al., 2001, Nobile, 2004, Woo and Sriram, 2006, Henna, 2008, Celeux et al., 2019].

Frequentist estimation of the number of
clusters is commonly treated as a model selection problem, approaches to which include likelihood ratio tests [McLachlan, 1987], information criteria [Biernacki et al., 2000, Celeux et al., 2006, Drton and Plummer, 2017] and overfitting with penalised likelihood mixture models [Manole and Khalili, 2021]. Identifiability issues are a major difficulty in both Bayesian and frequentist formulations [Liu and Shao, 2003]; see [Celeux et al., 2019] for an overview of the latter.

Bayesian mixture models treat the number of clusters as an unknown parameter for which principled uncertainty estimates may be obtained. We consider the mixture of finite mixtures [Nobile, 2004, Miller and Harrison, 2018] and Dirichlet process mixture models [MacEachern, 1994, Escobar and West, 1995, Neal, 2000], two common extensions of the finite mixture model that allow inference for the number of clusters. The first benefits from consistency results [Nobile, 1994], but estimation methods were inefficient until recently [Miller and Harrison, 2018]. The second is a nonparametric model that has become popular due to the existence of  efficient estimation algorithms [Neal, 2000, Jain and Neal, 2004].

Despite their attractions, recent theory has established that estimators of the numbers of clusters based on these models have undesirable asymptotic properties. Dirichlet process mixtures have been shown to give inconsistent estimators of the number of clusters, overestimation of which has been observed for finite samples due to the appearance of small extra clusters [Miller and Harrison, 2014]. Mixtures of finite mixtures have been suggested as a consistent alternative, but model misspecification can still result in inconsistency and poor estimation of the number of clusters in practice [Cai et al., 2021]. Moreover, high-dimensional data are particularly challenging for all clustering methods, which often incorrectly estimate the number of clusters [Chandra et al., 2021].

Although unsettling, poor asymptotic properties need not entail major problems in practical situations with finite samples, and our first aim in this work is to investigate to what extent the  asymptotic results sabotage the estimation of the number of clusters with Bayesian mixtures models with sample sizes that are typical for applications.

In practice inference is generally performed using Markov chain Monte Carlo (MCMC) algorithms [Neal, 2000]. This results in a large number of MCMC samples of cluster allocations, which cannot be interpreted directly owing to the presence of multiple modes over all the partitions of the data, changes in the number of components, and label switching, so  summarisation methods have been introduced to find a representative clustering [Binder, 1978, Medvedovic et al., 2004, Fritsch and Ickstadt, 2009, Wade and Ghahramani, 2018, Rastelli and Friel, 2018, Dahl et al., 2021]. To summarise MCMC samples despite label switching [Jasra et al., 2005], these post-processing methods are often based on Bayesian decision theory, that is, they minimise the posterior expectation of a suitable loss function.
The different choices among all possible combinations of loss function and optimisation method have led to the use of various summaries [Huelsenbeck and Andolfatto, 2007, Savage et al., 2010, Kirk et al., 2012, Lock and Dunson, 2013, Roth et al., 2014, Yuan et al., 2015, Carmona et al., 2019, Zafar et al., 2019, Fuentes-García et al., 2019, Ni et al., 2020], but the effect of summarisation methods has not yet been studied, especially with regard to the summarising number of clusters, which need not be representative of the corresponding marginal posterior distribution.
Hence a second goal of our work is to investigate how the summarisation method affects the number of clusters in light of the theoretical results. We focus on the low-dimensional case to avoid complications that can arise with clustering methods and MCMC samplers in high dimensions.

Another popular approach to estimating the number of components  is overfitting with finite mixture models. By deliberately fitting a model with an upper bound $K$ on the number of components we can estimate the true number of clusters $K_0$ by discarding empty components.
[Rousseau and Mengersen, 2011] show that overfitting can give  consistent estimators
, and that additional components vanish a posteriori at rate $N^{-1/2} \log(N)^\beta$ for some $\beta>0$.
We do not consider this approach, which requires
setting parameter $K$ which should be higher than the unknown $K_0$, and for $K$ high enough, this model is comparable to a Dirichlet process mixture [Frühwirth-Schnatter and Malsiner-Walli, 2019].

Non-decision-theoretic approaches intended to correct Dirichlet process mixture overestimation of the number of components also exist. For instance, [Guha et al., 2019] introduce an algorithm that merges small additional clusters of each MCMC sample in post-processing. They prove that their method is consistent for the number of components, but it involves a model-specific theoretical contraction rate to adjust the truncation of small clusters; they provide this rate for Gaussian mixture models. We do not consider such algorithms, because decision theory-based summaries are more common in practice.

In Section <ref>, we introduce the types of Bayesian mixture models that we consider, summarise known asymptotic (in)consistency results and present different strategies for summarising simulation output for these models. We perform a study on synthetic data in Section <ref> and illustrate the results using gene expression data. We conclude by discussing our results and give an overview of alternative approaches and potential future work.

§ BACKGROUND, THEORY AND SUMMARISATION METHODS

§.§ Bayesian mixture models

§.§.§ Finite mixture model

Mixture models presuppose that the data are generated from a mixture of component distributions belonging to a given parametric family, and clusters are determined by observations generated by the same component.
A Bayesian finite mixture model for observations $X_1, \ldots,  X_N$ with $K$ components can be defined hierarchically as

\begin{equation}
\begin{aligned}
\label{def:fmm}
X_i \mid \bm \theta, z_i &\sim F( \theta_{z_i} ), \quad i = 1, \ldots, N,\\
\theta_k  &\sim H, \quad k = 1,\ldots,K, \\
z_i \mid  \bm{\pi} &\sim \text{Categorical}_K(\bm{\pi}), \quad i = 1, \ldots, N, \\
\bm{\pi}  &\sim \text{Dirichlet}_K(\gamma, \ldots , \gamma),
\end{aligned}
\end{equation}

where the component distribution $F$ is parameterised by $\theta \in \Theta$ and $H$ is a prior distribution on the component parameters
. The component allocations $z_i$ indicate the component responsible for each observation, while the vector $\bm \pi$ represents the component proportions (also called the mixing distribution), i.e., $\pi_k$ is interpreted as the prior probability that an observation belongs to component $k$. The prior on $\bm \pi$ is here taken to be a symmetric Dirichlet distribution with a parameter $\gamma > 0$ commonly set to $1$ or $1/K$.
By marginalising over $\bm z$, the probability density function of an observation $X_i$ can be written as

\begin{equation}
\label{eq:mix}
p \left( X_i \mid \bm \theta, \bm \pi \right) = \sum_{k=1}^K \pi_k f( X_i \mid \theta_k),
\end{equation}

where $f$ is the component probability density function associated with $F$. Thus, we recover the standard representation of a mixture.

The most common choice of component distribution for continuous data  is multivariate Gaussian. For variables $ X_i \in \mathbbm R^p$ this corresponds to setting $F(\theta_k) = \mathcal N ( \mu_k,  \Sigma_k)$ such that $\theta_k = ( \mu_k, \Sigma_k)$ for $k=1, \ldots, K$, where $ \mu_k$ is a $p \times 1$ vector and $ \Sigma_k$ is a positive semi-definite $p \times p$ matrix.

§.§.§ Mixture of finite mixtures model

A natural extension of the finite mixture model (<ref>) is to treat the number of components $K$ as an unknown parameter with a prior $p_K$. The resulting mixture of finite mixtures (MFM) model can be defined hierarchically as

\begin{equation}
\begin{aligned}
\label{def:mfm}
K &\sim p_K, \\
\bm{\pi} \mid K &\sim \text{Dirichlet}_K(\gamma,\ldots,\gamma), \\
z_i \mid K, \bm{\pi} &\sim \text{Categorical}_K(\bm{\pi}), \quad i = 1, \ldots, N, \\
\theta_k \mid K &\sim H, \quad  k = 1,\ldots,K, \\
X_i \mid K, \bm{\theta}, z_i &\sim F( \theta_{z_i} ), \quad  i = 1, \ldots, N.
\end{aligned}
\end{equation}

If we set $p_K(k) = \mathbbm 1_{ \{ k = K_0 \} }$ we recover a finite mixture model with $K_0$ components. The mixture of finite mixtures model is sometimes also called a finite mixture model with a prior on the number of components [Nobile, 1994, Richardson and Green, 1997].

Some examples of discrete priors $p_K$ are a Poisson distribution for $K-1$ [Nobile, 2007] or a Geometric distribution [Miller and Harrison, 2018].

§.§.§ Dirichlet process mixture model

Another extension of the finite mixture model is the Dirichlet process mixture model, under which

each observation $X_i$ has an associated parameter, $\theta_i$, with $\theta_1, \ldots, \theta_N$ assumed to be independent and identically distributed according to a random distribution $G$.
The prior on $G$ is a Dirichlet process, $G \sim \text{DP} (\alpha, H)$, where $\alpha$ is the concentration parameter and $H$ is the base measure.

The Dirichlet process mixture (DPM) model can therefore be written hierarchically as

\begin{equation}
\begin{aligned}
\label{def:dpm}
G &\sim \text{DP}(\alpha, H),\\
\theta _i \mid G &\sim G, \quad i = 1,\ldots,N, \\
X_i \mid \theta_i &\sim F( \theta_i ), \quad  i = 1, \ldots, N.
\end{aligned}
\end{equation}

A practical definition of the resulting distribution of $\theta$ is provided by the Pólya urn scheme [Blackwell and MacQueen, 1973].

In contrast to the MFM, under which the prior density of the number of components is fixed as $p_K$, the Dirichlet process prior allows an increasing number of components as the number of observations $N$ grows. Indeed, the prior expectation of $T$, the number of components for $X_1, \ldots, X_N$ or equivalently the number of unique values of $\bm z$, satisfies

\begin{equation}
\mathbbm E (T) = \sum_{i = 1}^N \frac{\alpha}{\alpha + i -1} = O\left(\alpha \log \frac{N}{\alpha} \right).
\end{equation}

For more properties of  the Dirichlet process see [Ghosal, 2010], for example.

§.§ Known asymptotic properties

There are various theoretical results for the estimation of the number of components with Bayesian mixture models. In general we present asymptotic results as the number of observations $N$ tends to infinity.
Due to the generality of mixture models, useful results require further assumptions, but these hold in most practical settings and in particular for those we study in Section <ref>.

Although we focus on component allocation and the number of components, an underlying issue is density estimation: although the posterior number of components may be inconsistent, the mixture density $\sum_{k=1} \pi_k f(\bm X \mid  \theta_k)$ can still converge to the true generating density.   [Nguyen, 2013], for example, shows that the posterior distribution of a DPM model converges to the true distribution in the Wasserstein metric, and [Guha et al., 2019] show that the posterior distribution of a MFM model contracts to the true generating parameters at an optimal rate under correct specification, up to relabeling of the components.  These authors also show that even under misspecification both types of models can asymptotically recover the true parameters under certain conditions.
Under an identifiability condition on the true generating distribution, they show in particular that a multivariate normal mixture distribution converges in the Wasserstein metric, though the convergence is slower than for the well-specified case. However, as we describe in Section <ref>, under misspecification the estimated number of components is  inconsistent and diverges. This is a recurring trade-off for mixture models: it is possible to obtain an arbitrarily good density estimate, but at the cost of overestimating the number of components.

§.§.§ Consistency of mixture of finite mixtures

For the mixture of finite mixtures model, [Nobile, 1994] shows that the posterior converges to the true number of components, if the mixture is identifiable and the component distribution is continuous. Mixture-identifiability means that a mixture of $K$ components is strictly more expressive than any mixture with $K' < K$ components and it is uniquely identified up to relabeling of the components. Many models satisfy this condition, including Gaussian mixtures and most common continuous component distributions.

§.§.§ Inconsistency of Dirichlet process mixture

[Miller and Harrison, 2014] show that the posterior distribution of the number of components of a Dirichlet process mixture, or more generally a Pitman–Yor mixture model does not contract to the true number of components when the observations arise from a finite mixture.

However it is still unclear what the posterior distribution actually does, for instance, it might contract around a wrong value, or diverge entirely. [Yang et al., 2020] show that the posterior diverges when a uniform or Gaussian prior is taken for $\bm \theta$.
Since we do not know the convergence speed it is difficult to assess the finite-sample repercussions of this, and other aspects such as the choice of $\alpha$ may be more important when $N$ is finite.
In fact, [Miller and Harrison, 2014] do not cover the case where $\alpha$ is also a parameter with a hyperprior, though they suspect their result to hold regardless. However, more recent work has shown that consistency can be achieved by putting a prior on $\alpha$ [Ascolani et al., 2022], or allowing it to depend on the sample size [Ohn and Lin, 2022], although it is unclear whether these results hold for mixtures of multivariate Gaussian densities with unknown covariance matrix, which are most commonly used for clustering.

§.§.§ Inconsistency under misspecification

The consistency result for the mixture of finite mixtures model in Section <ref> holds when the generating and modelling component distribution families are the same, but we may ask whether it remains true under misspecification, which is frequent, perhaps even invariable, in practice.
We saw in Section <ref> that misspecification of the allocation prior can lead to inconsistency. [Cai et al., 2021] show that mixtures of finite mixtures with misspecified component distributions are also inconsistent for the number of components. Indeed, they show that the posterior diverges for any misspecified location-scale family of component distributions, if the corresponding mixture is identifiable and absolutely continuous with respect to the parameter $ \bm \theta$, and the data generating mixture lies in the Kullback–Liebler support of the prior. Under an assumption of degenerate limits they also prove a more general statement for any continuous density/mass function.

§.§ Summarisation of MCMC samples

Analytic computations for Bayesian mixture models are intractable, so the usual approach is to use Markov chain Monte Carlo (MCMC) sampling, which results in a time series $\bm z^{(1)},\ldots,\bm z^{(M)}$ from the posterior $p(\bm z,\bm  \theta \mid X_1, \ldots, X_N)$, where $M$ is chosen adequately large in practice.
Even if the $\bm z^{(m)}$ are sampled from the exact posterior distribution, the results cannot be interpreted directly  using sample averages or other standard summary statistics owing to two major  identifiability issues for the components, namely label-switching and variation in the number of components.

Label-switching creates an incoherence in the values of the posterior samples of $\bm z$. Any permutation of the component labels leaves the posterior density unchanged,
so taking a summary clustering, such as the average of the samples, is not helpful, because  two samples $\bm z^{(m)}$ and $\bm z^{(m')}$ may represent the same component allocation with permuted labels: with $N=3$ and $K=2$, for example, the labelings $1,1,2$ and $2,2,1$ lead to the same clustering.

Moreover, the number of components vary during sampling, further complicating the challenge of summarisation.

Formally, a summarisation method is a function from $\bm z^{(1)},\ldots,\bm z^{(M)}$ to a single summary clustering $\bm z^*$ of the $N$ observations, i.e.,

\begin{equation}
\label{def:summ_clust}
\begin{aligned}
\mathbbm N^{N\times M} &\rightarrow \mathbbm N ^N\\
\big\{\bm  z^{(m)}  \big\}_{m = 1, \ldots, M } &\rightarrow \bm z^*.
\end{aligned}
\end{equation}

Figure <ref> illustrates the summarisation task and shows the difference between the posterior number of components and the number of clusters of a summary clustering.  Post-processing methods are intended to provide a “representative” $\bm z^*$.  The concept of representativity is not uniquely defined, but [Binder, 1978] introduced a formal definition within a loss minimisation framework. For most summarisation methods we can define a loss function $L : \mathbbm N^N \times  \mathbbm N^N \rightarrow \R_{\geq 0}$ such that

\begin{equation}
\label{def:loss}
\bm z^* = \argmin_{\hat{\bm z}} \mathbbm E \left[ L(\bm z,\hat{\bm z}) \mid \bm X \right],
\end{equation}

where the expectation is taken over $\bm z$.
In order that the loss $L$ be invariant to relabeling, clusterings should be treated as partitions of the set $\{1, \ldots,N \}$, and then a direct MCMC approximation of (<ref>) is

\begin{equation}
\label{eq:loss_approx}
\bm z^* = \argmin_{\hat{\bm z}} \frac{1}{M} \sum_{m = 1}^M L \left( \bm z^{(m)},\hat{\bm z} \right). %\tilde p \left[ \bm z^{(m)} \mid \bm X \right],
\end{equation}

A central quantity for post-processing samples is the posterior similarity matrix (PSM) $\bm P$ with elements

\begin{equation}
\label{def:psm}
P_{ij} = p( z_i = z_j \mid X )
\approx \frac{1}{M} \sum_{m=1}^M \mathbbm 1 \left[ z_i^{(m)} = z_j^{(m)} \right], \quad i,j = 1, \ldots , N,
\end{equation}

where $\mathbbm 1$ is the indicator function, which measures the frequency of co-clustering of pairs of observations and is invariant to label-switching or to changes in the numbers of components during sampling.
To simplify notation in this section we write $K(\bm z)$ for the number of clusters in $\bm z$ and denote the contingency counts of two clusterings $\bm z$ and $\hat{\bm  z}\bm $ as

\begin{equation}
n_{kk'} = \sum_{i,j = 1}^N \mathbbm 1(z_i = k) \mathbbm 1(\hat z_j = k'), \quad
n_k = \sum_{i=1}^N \mathbbm 1(z_i = k),\quad k = 1, \ldots, K(\bm z), k ' = 1, \ldots, K(\hat{\bm z}).
\end{equation}

These form a $K(\bm z) \times K (\hat{\bm z})$ matrix representation of the information relevant to cluster comparisons.

Illustration of MCMC samples and an MCMC summary clustering. Cluster labels are represented by different colours.

In the following section we present various loss functions including the most used in practice. Given a loss function, we then need to optimise (<ref>) to obtain a summary clustering. We present algorithms for approximate optimisation in Section <ref>.

§.§ Loss functions

§.§.§ Maximum a posteriori

A straightforward approach to obtaining a summary clustering $\bm z^*$ is to take the sample $\bm z^{(m)}$ that maximizes the unnormalised posterior density. This corresponds to taking the 0-1 loss in (<ref>), i.e., $L_{\text{0-1}}(\bm z,\hat{ \bm z}) = \mathbbm 1 ( \hat{\bm  z} = \bm z)$ where $\hat{\bm z} = \bm z$ is defined up to relabeling. This entirely circumvents the label switching problem by only considering a single sample, rather than comparing clusterings with possibly different labelings.

§.§.§ Binder's loss

The first loss function for summarising MCMC samples is due to [Binder, 1978], which uses pairwise co-clustering disagreement to measure differences between two clusterings via

\begin{equation}
L_{\text{B}} (\bm z, \hat{\bm z}) = \mathlarger{\sum}_{i < j}
l_1 \mathbbm 1 (z_i = z_j) \mathbbm 1 ( \hat z_i \neq \hat z_j)
+ l_2 \mathbbm 1 (z_i \neq z_j) \mathbbm 1 (\hat z_i = \hat z_j),
\end{equation}

where $l_1, l_2 > 0 $ are commonly set to $l_1=l_2=1$, though other values could be used in case of asymmetric classification error. We set $l_1 = l_2 = 1$, but [Dahl et al., 2021] show that changing the ratio $l_1/l_2$ can significantly affect the number of clusters.

The expected posterior loss can be written in terms of the posterior similarity matrix as

\begin{equation}
\mathbbm E \left[ L_{\text{B}} (\bm z, \hat{\bm z} ) \mid X \right] = \sum_{i < j} \lvert \mathbbm 1 ( \hat z_i = \hat z_j) - P_{ij} \rvert .
\end{equation}

Binder's loss has a direct relationship with the commonly-used Rand index $\text R(\cdot,\cdot)$  [Rand, 1971]. Indeed, for $l_1 = l_2 = 1$,

\begin{equation}
\label{eq:binder_rand}
L_{\text{B}} (\bm z, \hat{\bm z}) = {N \choose 2} \left[ 1 - \text R(\bm z, \hat{\bm z}) \right],
\end{equation}

so minimising the expected posterior Binder loss is equivalent to maximising the expected posterior Rand index.

§.§.§ Posterior Expected Adjusted Rand (PEAR)

[Fritsch and Ickstadt, 2009] propose using the  summary clustering that maximises the posterior adjusted Rand index [Hubert and Arabie, 1985]. While they do not explicitly introduce a loss function, we define the adjusted Rand loss as

\begin{equation}
L_{\text{AR}} (\bm z, \hat{\bm z}) = 1 - \text{AR}(\bm z, \hat{\bm z}),
\end{equation}

where the adjusted Rand index AR$(\cdot,\cdot)$  corrects the Rand index for chance under permutations within clusters.
For computational reasons [Fritsch and Ickstadt, 2009] suggest taking the approximation $\mathbbm 1 (z_i = z_j) \approx P_{ij}$ in AR$(\bm z, \hat{\bm z})$, but [Dahl et al., 2021] argue that this  is unnecessary and may lead to inaccuracies. Nevertheless, we will use
the approximate optimisation, as it is typically used in practice.

§.§.§ Variation of Information

A loss function can also be based on a distance metric on the space of clusterings. [Wade and Ghahramani, 2018] suggest using as a loss function the variation of information [Meila, 2007], which has the advantage of inheriting theoretical properties from the distance space. The variation of information distance between two clusterings is

\begin{equation}
\hspace{-3mm}
\begin{aligned}
L_\text{VI}(\bm z,\hat{\bm z}) &= -H(\bm z) - H( \hat{\bm  z}) + 2 H(\bm z, \hat{\bm z}) \\
&=  \sum_{k=1}^{K(\bm z)} \frac{n_k}{N} \log_2\left(\frac{n_k}{N}\right) +
\sum_{k'=1}^{K(\hat{\bm z}) } \frac{\hat n_{k'}}{N} \log_2\left(\frac{\hat n_{k'}}{N}\right)
- 2 \sum_{k=1}^{K(\bm z)} \sum_{k'=1}^{K (\hat{\bm z})} \frac{n_{kk'}}{N}
\log_2\left(\frac{n_{kk'}}{N}\right),
\end{aligned}
\end{equation}

where $H$ is called the entropy  and $H(\bm z, \hat{\bm z} )$ is the joint entropy for cluster allocations. Intuitively, this distance is small if the amount of information two clusterings share, captured by the joint entropy, is close to the individual information of each clustering, captured by their entropy.

The posterior expectation of $L_\text{VI}$ by (<ref>) involves $\mathbbm E \left[ \log_2  \left( n_{kk'}\right) \mid \bm X \right] $, computation of which requires the contingency counts for each sample. Instead, [Wade and Ghahramani, 2018] find a lower bound on the posterior expected loss, which only requires the computation of the posterior similarity matrix and leads to the loss function

\begin{equation}
L_{\text{VILB}}(\bm z,\hat{\bm z})    = \sum_{k'=1}^{K(\hat{\bm z})} n_{k'} \log_2 n_{k'} - 2 \sum_{i =1}^N \log_2 \left[ \sum_{j=1}^N \mathbbm 1\left(\hat z_i = \hat z_j \right) P_{ij} \right],
\end{equation}

where we omit terms that do not depend on $\hat{\bm  z}$. We use the lower bound approximation in our implementation.

§.§ Optimisation

Loss optimisation over all possible clusterings is infeasible for computational reasons. Indeed, if $\mathcal C (N)$ is the set of all clusterings of $\{1, \ldots, N\}$, then $|\mathcal C (N)|$ grows very rapidly as a function of $N$. The number of partitions into two clusters is $2^{N-1}-1$, so a simple lower bound is $2^{N-1}-1 \leq |\mathcal C (N)|$, indicating that optimisation of the loss function is NP-hard  and it is necessary to consider a reduced set of clusterings. We approximate the optimisation as

\begin{equation}
\label{def:opt}
\bm z^* = \argmin_{\hat{\bm z} \in \mathcal C (N)} \mathbbm E \left[ L(\bm z,\hat{\bm z}) \mid \bm X \right]
\approx \argmin_{\hat{\bm z} \in S} \mathbbm E \left[ L(\bm z,\hat{\bm z}) \mid \bm X \right],
\end{equation}

where $S \subset \mathcal C(N)$. In the following subsections we consider various choices of $S$.

We also mention more recent optimisation methods that we do not consider for brevity. [Rastelli and Friel, 2018] use a greedy search algorithm and avoid computing the posterior similarity matrix, which is a costly intermediate step when $N$ is large, but nevertheless require certain types of loss functions which can be updated efficiently when a single $z_i$ is changed. [Dahl et al., 2021] extend this optimisation algorithm with parallel programming and further search steps. However, these two approaches have not yet found wide use in practice, so we restrict ourselves to popular methods.

§.§.§ Hierarchical clustering

A common post-processing approach is to convert the posterior similarity matrix $\bm P$ into a metric between any two observations and then apply standard hierarchical clustering methods [Medvedovic et al., 2004, Fritsch and Ickstadt, 2009].
For this we define the posterior dissimilarity matrix $\bm D$ as $D_{ij} = 1 - P_{ij}$ for $i,j = 1, \ldots,N$, and perform agglomerative clustering on $\bm D$ with either complete or average linkage.

Denoting the cluster allocation with $k$ clusters by $\hat{\bm z}(k)$, [Fritsch and Ickstadt, 2009] suggest optimising the loss function in (<ref>) with $S = \left\{ \hat{ \bm z}(k) : k = 1, \ldots, k_{\text{max}} \right\}$ for some upper bound $k_{\text{max}}$ such as $\lfloor N / 8 \rfloor$. If the posterior similarity matrix is already computed, this approach is very efficient, as we can update the linkage distances in constant time at every merging step, for a total of $O(N)$ operations. However, this yields only a single element of $S$ per value of $k$, whereas a desirable property of $S$ would be concentration around an appropriate numbers of clusters.

[Medvedovic et al., 2004] introduced hierarchical clustering without loss function optimisation, instead obtaining $\bm z^*$ directly as the clustering that exceeds a linkage distance $1-\varepsilon$ across any pair of clusters. [Fritsch and Ickstadt, 2009] suggest  $1-\varepsilon = 0.99$, which we consider as a default. Unlike other methods we present, Medvedovic clustering does not depend on a loss function and does not seem to be expressible in decision-theoretic terms.

§.§.§ Partitioning around medoids

[Liverani et al., 2015] suggest optimisation using the popular $k$-medoids algorithm, which generalises $k$-means for any metric. More precisely, they use $D_{ij} = 1 - P_{ij}$ as a distance metric and choose partitioning around medoids as an efficient implementation of $k$-medoids.   [Kaufman and Rousseeuw, 1990] improve this approach by instead successively swapping medoids to minimize a global cost function. Nonetheless, when $N$ is large, partitioning around medoids is also inefficient as it requires $O\left[ k(N-k)^2 \right]$ iterations. We may therefore consider more efficient approaches such as that of [Schubert and Rousseeuw, 2019].

By denoting the output clustering of k-medoids with $k$ clusters by $\hat{\bm z}(k)$, [Liverani et al., 2015] suggest optimising the loss function in (<ref>) with
$S = \left\{ \hat{\bm  z}(k) : k = 1, \ldots, k_{\text{max}} \right\}$ for some upper bound $k_{\text{max}}$. We choose $k_{\text{max}} = \lfloor N/8 \rfloor$ in our implementation.

§.§.§ Sample search

A straightforward choice for $S$ is the set of MCMC samples $S = \left\{ \bm z^{(1)}, \ldots , \bm z^{(M)} \right\}$ [Dahl, 2006]. The resulting summary clustering is a posterior sample $\bm z^{(m)}$, which also provides associated sample component parameters $\bm \theta^{(m)}$. However, if all samples inherit an undesired property from the sampling method or the model itself,  such as overestimation of the number of components, the summary clustering will also have this property.

§ EVALUATION

We now summarise the results of our analysis applying MFMs and DPMs, with and without conjugate hyperpriors on the component parameters, and with different summarisation methods.  We use three scenarios with simulated data and one with real data.

The first scenario corresponds to a sample with a moderate sample size for which MFMs and DPMs are similar. Since the marginal posteriors for the number of clusters are highly concentrated at the truth, this easy case sets a baseline for comparing the different summarisation methods.

The second scenario corresponds to a larger sample size at which the DPM and MFM start to differ and which is more challenging for the summarisation methods. We focus on the DPM with hyperpriors since the corresponding posterior for the number of clusters is the flattest, with a mode that is higher than the truth in some of the simulated datasets.

The third scenario is used to explore how a common type of misspecification affects summarisation methods, by modelling correlated data using multivariate Gaussian components with diagonal covariance matrices.

Finally, we show results for the gene expression dataset used in Cai et al., 2021, where we expect to observe both the effects of DPM inconsistency and of misspecification. Again, we compare the outcomes of the different summarisation methods and use cell types as a proxy for true labels.

§.§ Modelling specification

We consider continuous data and multivariate normal mixture models. In the multivariate models we consider either a model with full covariance matrix or one with a diagonal covariance matrix. Both assumptions are common in practice, and in particular the restriction to a diagonal covariance has valued computational advantages.
The component distributions with a full covariance specification are

\begin{equation}
\begin{aligned}
\label{def:mvn}
\mu_k  &\sim \mathcal{N}_p\left(\widehat{ \mu}, \widehat{ C} \right), \\
\Sigma_k  &\sim \text{InverseWishart}\left(p,\widehat{C}\right), \\
X_i \mid  z_i,  \bm \mu, \bm \Sigma &\sim \mathcal{N}_p\left(  \mu_{z_i}, \Sigma_{z_i} \right), \quad  i = 1, \ldots, N,
\end{aligned}
\end{equation}

where $\widehat{ \mu} = N^{-1}\sum_{i=1}^N X_i$ is the sample mean, $\widehat{ C} = N^{-1}\sum_{i=1}^N (X_i - \widehat{ \mu}) (X_i - \widehat{ \mu})^T $ is the sample covariance and $k$ is the component index. As it is common to include hyperpriors to give models more flexibility or to incorporate prior assumptions, we also consider an alternative model with conjugate hyperpriors on the component parameters, i.e.,

\begin{equation}
\begin{aligned}
\label{def:mvnh}
m_k &\sim \mathcal N_p\left(\widehat{\mu}, \widehat{C} \right),\\
C_k &\sim  \text{InverseWishart}_p\left(p,\widehat{C}\right), \\
\nu_k - p + 1 &\sim \text{Gamma}\left(\alpha = 2, \beta = 2\right),\\
W_k &\sim \text{InverseWishart}_p\left(p,\widehat{ C}\right), \\
\mu_k \mid m_k, C_k &\sim \mathcal{N}_p\left(m_k, C_k\right), \\
\Sigma_k \mid \nu_k, W_k  &\sim \text{InverseWishart}_p\left(\nu_k,W_k\right), \\
X_i \mid  z_i, \bm \mu, \bm \Sigma &\sim \mathcal{N}_p\left( \bm \mu_{z_i}, \bm \Sigma_{z_i} \right), \quad  i = 1, \ldots, N.
\end{aligned}
\end{equation}

For the diagonal covariance model we standardize the observations beforehand, so our modelling assumptions can rely on standardised data. Hence, we have

\begin{equation}
\begin{aligned}
\label{def:mvnaac}
\lambda_{kd} &\sim \text{Gamma}(1,1),\\
\mu_{kd} \mid \lambda_{z_id}^{-1}  &\sim \mathcal{N}\left(0, \lambda_{z_id}^{-1}\right), \\
X_{id} \mid  z_i, \bm \mu, \bm \lambda &\sim \mathcal{N}\left( \mu_{z_id}, \lambda_{z_id}^{-1}\right),  \quad i = 1, \ldots, N,
\end{aligned}
\end{equation}

for each dimension $d=1,\ldots,p$ independently.

We consider two types of allocation prior, either a mixture of finite mixtures (<ref>) or a Dirichlet process mixture (<ref>). For the first, $\bm{\pi} \mid K \sim \text{Dirichlet}_K(1,\ldots,1)$ and $K \sim \text{Geometric}(0.1)$, following [Miller and Harrison, 2018]. For the second, $\alpha \sim \text{Exponential}(1)$.

For both models we obtain MCMC samples with the split-merge sample with four independent chains, each with $2000$ iterations, and the $100$ first iterations removed as burn-in. We thin the MCMC samples by keeping every second iteration, giving $M = 3800$ samples in total. [Miller and Harrison, 2018] show fast convergence of the split-merge sampler, so we choose the sample size sufficiently large to be representative while limiting the running time. We assess convergence using the Geweke and the Gelman–Rubin $\hat R$ diagnostics  across multiple chains [Geweke, 1991, Gelman and Rubin, 1992, Vehtari et al., 2020].

§.§ Simulation study

In order that the true number of components is known, we generate synthetic data using a finite mixture model from [Miller and Harrison, 2018] with four clusters. Our synthetic bivariate observations are distributed as

\begin{equation}\label{eq:sim-data}
\begin{aligned}
z_i &\sim \text{Categorical}( \bm \pi ), \quad i=1, \ldots,N, \\
X_i &\sim \mathcal N_2 \left(\mu_{z_i}, \Sigma_{z_i}\right), \quad i=1, \ldots,N,
\end{aligned}
\end{equation}

\begin{equation}
\begin{gathered}
\bm \pi = (0.44, 0.3, 0.25,0.01),\\[1pt]
\mu_1 = (4,4)^T, \, \mu_2 = (7,4)^T, \, \mu_3 = (6,2)^T, \, \mu_4 = (8,11)^T, \\[1pt]
\Sigma_1 = \begin{pmatrix}
2 & 0\\
0 & 2
\end{pmatrix},\, \Sigma_2 = R \begin{pmatrix}
2 & 0\\
0 & 2
\end{pmatrix}R^T, \, R = \begin{pmatrix}
\cos(\pi/4) & -\sin(\pi/4)\\
\sin(\pi/4) & \cos(\pi/4)
\end{pmatrix},\\[1pt]
\Sigma_3 = \begin{pmatrix}
3 & 0\\
0 & 0.1
\end{pmatrix}, \,
\Sigma_4 = \begin{pmatrix}
0.1 & 0\\
0 & 0.1
\end{pmatrix}.
\end{gathered}
\end{equation}

§.§.§ Moderate sample size

To obtain representative results we generated $50$ independent datasets with $N= 500$ observations from model (<ref>). Figure <ref> shows one of these datasets, which display overlap between the component densities. Uncertainty between clusters is often seen in applications, so by fixing $\bm \pi$, $\bm \mu$ and $\bm \Sigma$ we control the difficulty of the clustering task. If instead we generate data from the model with hyperpriors, then we do not control distances between component means, which may confound our interpretation of results.

Synthetic Gaussian data with four clusters and $N = 500$ observations.

We perform estimation with the full covariance models (<ref>) and (<ref>) and four model types: the mixture of finite mixtures with and without a hyperprior on the parameters of the mixture components and the Dirichlet process mixture with and without hyperpriors, which we respectively abbreviate to MFMH, MFM, DPMH and DPM. At each iteration we obtain a sample of component allocations $\bm z$ and the number of active components $T$.

Figure <ref> shows every resulting posterior on the number of components. For clarity, we maintain a consistent colouring of the four model types in all following figures. To obtain the posterior of $K$ for the mixture of finite mixtures model, we use the posterior on $K$ given the number of active components $T$, provided by [Miller and Harrison, 2018]. Although the dataset in Figure <ref> is well-specified by the full covariance model and $N = 500$ is a moderate sample size, there are differences between the model types: the hyperpriors seem to flatten the posterior number of components towards higher values, and  the Dirichlet process mixture allows exploration of higher numbers of components during MCMC sampling.

The identifiability issues make diagnostics for MCMC samples of $\bm \theta$, $\bm z$ and $\bm \pi$ challenging, and we can only consider the number of components or $\alpha$. Figure <ref> in the Appendix gives an example of the number of components for a single chain before thinning and including the burn-in samples when fitting each of our models to the data shown in Figure <ref>.

Posteriors on the number of components for a mixture of finite mixtures (a) without and (b) with a hyperprior,  and for a Dirichlet process mixture (c) without and (d) with a hyperprior, for $50$ synthetic datasets with $N = 500$. There are in fact four components.

For every dataset and model type we derive the summary clusterings described in Section <ref>, each of which gives a number of clusters.
As an example, Figure <ref> summarises the results for the mixture of finite mixtures model without a hyperprior. The posteriors at the right are those in Figure <ref>. The summarisation methods with a loss function are Binder's loss (Binder), the posterior expected adjusted Rand index (PEAR), the lower bound of the variation of information (VI-LB), average and complete linkage hierarchical clustering (Average and Complete), sample search (Samples) and partitioning around medoids (PAM). Although this case is rather simple, some methods vastly overestimate the number of clusters.

Number of clusters of summary clustering methods and posterior number of components for the mixture of finite mixtures (MFM) model without a hyperprior, for $50$ independent synthetic datasets with $N=500$. Summary names are composed of the loss function and the optimisation method, except for Medvedovic clustering.  The true number of components, four, is shown by the horizontal dashed line.

Our main focus is the number of clusters, but  this is linked to the underlying cluster allocations.
Nonetheless, an accurate clustering can coexist with an inaccurate number of clusters, at least when accuracy is measured with popular scores such as the adjusted Rand index.
The addition of small clusters has a minor impact on these scores provided the large clusters are correctly identified, as illustrated in the Appendix; see Figure <ref>.

§.§.§ Dirichlet process mixture inconsistency

To investigate the Dirichlet process mixture and asymptotic results of Section <ref> we set $N = 10^4$. We consider the Dirichlet process mixture with hyperpriors, as it shows the highest overestimation of the number of components even when $N = 500$, and in Figure <ref> show the posterior number of components for $50$ independent datasets.  The Dirichlet process mixture model tends to overestimate the number of components, in particular compared to the mixture of finite mixtures model.

Posteriors on the number of components for a mixture of finite mixtures (a) without and (b) with a hyperprior, and for a Dirichlet process mixture (c) without and (d) with a hyperprior, for $50$ synthetic datasets with $N = 10^4$. There are in fact four components.

The Dirichlet process mixture model tends to introduce small extra clusters, as noticed by [Miller and Harrison, 2014] and illustrated in Figure <ref> in the Appendix.

Figure <ref> gives the number of clusters for a selection of summary clusterings. We exclude results for Binder's loss and PEAR loss with average linkage, which can both massively overestimate the number of clusters, and we do not use partitioning around medoids because of its much longer running time when $N$ is large.
As expected, optimisation with sample search tends to follow the posterior distribution and leads to overestimation of the number of clusters, whereas Medvedovic clustering and variation of information both give the correct number of clusters for all $50$ datasets.  This confirms that summary clustering methods may eliminate any small extra clusters introduced by the DPM model. Indeed, this holds even for data whose posterior mode is appreciably larger than four and the posterior credibility of four components is low.

Figure <ref> of the Appendix gives two examples of summary clusterings; PEAR loss with average linkage is representative of methods not included in Figure <ref>, as it gives $161$ clusters, while variation of information leads to four clusters only.
Both Binder's loss and PEAR loss with complete linkage optimisation tend to create small clusters in the overlaps between components.
This occurs independently of the model type, so it is a different phenomenon from Dirichlet process mixture inconsistency; indeed, [Wade and Ghahramani, 2018] already notice this for Binder's loss.
Our results in Section <ref> show that Binder's loss leads to poor results for the number of clusters that for $N = 500$.

For PEAR loss we could investigate whether the approximation by [Fritsch and Ickstadt, 2009] is the reason for these additional clusters.  Variation of information successfully adjusts for small extra clusters introduced during sampling, but the resulting clustering differs from the generating component allocation, as the clusters do not have Gaussian shapes.
We conjecture that the borders between clusters approximate the contour lines of equal component densities, e.g., $\{ x: f(x \mid \theta_1) = f(x \mid \theta_2) \}$, which would explain the straight boundaries between clusters.

Number of clusters of summary clustering methods and posterior number of components for the Dirichlet process mixture model with a hyperprior (DPMH), for $50$ independent synthetic datasets with $N=10^4$. Summary names are composed of the loss function and the optimisation method, except for Medvedovic clustering.

§.§.§ Misspecification inconsistency

We now model the previous correlated data with $N = 10^4$ using the diagonal covariance model (<ref>), which is popular for computational reasons, to investigate asymptotic results from Section <ref>. The model (<ref>) is misspecified, because it cannot express the correlation from $\Sigma_2$. Indeed, with this choice we can specify a conjugate prior on the component distributions and sample directly from the component prior, thereby avoiding the conditional Gibbs update during the split-merge algorithm. Moreover the number of covariance parameters reduces from $p(p+1)/2$ to $p$ for every component, so sampling should be more stable in high-dimensional cases.

We again generate $50$ datasets of correlated Gaussian mixture observations and sample from the diagonal covariance model (<ref>). Figure <ref> shows   the posterior numbesr of components for two model types we consider. The MFM model gives concentrated posterior distributions for the number of components that contract onto a single value, though this varies across datasets. It is unclear why the model is sensitive to data generation, as we would expect almost identical data when $N$ is so large. MCMC samples, not shown here, show the creation of nearly identical overlapping components sharing observations of the well-specified components. One reason for this could be that, as $K$ gets large, the symmetric Dirichlet prior attributes more density to allocations with similar proportions, thereby giving potentially lower posterior credibility to uneven proportions. To investigate this we could choose a smaller $\gamma$, such as $\gamma = 1/K$. The DPM model posterior distributions also vary across datasets, but are less concentrated than those for the MFM model. This difference may arise because of the overestimation described in Section <ref>.

Posterior on the number of components under misspecification for (a) the mixture of finite mixtures; and (b) Dirichlet process mixture, with diagonal covariance on $50$ independent synthetic datasets with $N= 10^4$.  There are in fact four components.

Figure <ref> shows the numbers of clusters for summarisation methods for the MFM and DPM models.
As in the previous section, we omit methods with significant overestimation.
Overall summary clusterings seem to minimise the number of components compared to the posterior but they do not recover the true number of components.
Among these methods, variation of information leads to the most consistent reduction, giving seven clusters in most cases independently of the model type.
Sample search optimisation seems to follow the posterior distribution. The adjusted Rand index of summary clusterings, not shown here, is similar for all methods except for the slightly less accurate sample search, similarly to Figure <ref>.

Dirichlet process mixture model.

Mixture of finite mixtures model.

Number of clusters of summary clustering methods and posterior number of components for (a) the misspecified mixture of finite mixtures (MFM) model and (b) the misspecified Dirichlet process mixture (DPM) model, for $50$ synthetic datasets with $N=10^4$. Summary names are composed of the loss function and the optimisation method, except for Medvedovic clustering.  The true number of components, four, is shown by the horizontal dashed line.

The only component that cannot be expressed using the diagonal covariance model is the second component.  As $\Sigma_2$ is non-diagonal, the misspecified model introduces additional components. Figure <ref> shows an MCMC sample selected by the variation of information loss in which four components approximate the second component. The summary clustering with hierarchical clustering optimisation inherits the three additional components from the MCMC samples. It seems that summary clusterings do not adjust for misspecification in this example, but only avoid small additional clusters, in particular for the DPM model.

Summary clusterings with variation of information optimised with (a) sample search and (b) complete linkage hierarchical clustering for the misspecified mixture of finite mixtures (MFM) MCMC samples on synthetic data with $N=10^4$.

§.§ Gene expression data

We now consider gene expression data from single-cell RNA sequencing of somatosensory cells [Zeisel et al., 2015], which have $N = 3005$ to address asymptotic theory, and which were also used by [Cai et al., 2021] to illustrate the inconsistency result presented in Section <ref>.
The data are 19,972 gene measurements from 3005 cells with seven cell types. We pre-process the data following [Prabhakaran et al., 2016] and use the $p=10$ genes with the largest standard deviations because we aim to use our full covariance model in addition to the the diagonal covariance model used by [Cai et al., 2021]; [Prabhakaran et al., 2016] took $p=558$.
Modelling using the full covariance model is much slower for large $p$ due to the additional $p(p-1)/2$ correlation parameters, but sampling from both allows us to compare the diagonal covariance and more expressive full covariance models, though even the latter can be expected to be misspecified.   According to [Zeisel et al., 2015], there are seven known cell types, which we will consider as ground truth.

Figure <ref> gives the posterior distributions of the number of components for models (<ref>) and (<ref>).
The MFM and DPM models give similar posteriors for the full covariance model, with a mode at eight components, suggesting that the DPM overestimation of the number of components tends to be moderate for sample sizes in the order of $10^3$ and confirming our results in Section <ref>.
The MFM posterior gives slightly more credibility to higher numbers of components than does that of the DPM.  The diagonal covariance specification leads to a significant overestimation of the number of components, confirming the theoretical results from Section <ref> and our simulation results in Section <ref>.
Furthermore, the DPM model contributes to the overestimation of the number of components, with a mode at 18 components compared to 16 for the MFM model.
Both model types show bimodal posterior distributions with a second local mode and have similar shapes, though the DPM posterior is shifted towards higher numbers of components.

MFM (full covariance model).

DPM (full covariance model).

MFM (diagonal covariance model).

DPM (diagonal covariance model).

Posterior numbers of components for the mixture of finite mixtures (MFM) and Dirichlet process mixture (DPM) with full or diagonal covariance on mouse single-cell RNA-sequence data. There are seven true clusters according to [Zeisel et al., 2015].

Table <ref> gives the number of clusters of summarisation methods. We do not consider methods based on partitioning around medoids due to computational limitations. Overall we confirm our simulation results, with overestimation of the number of clusters for Binder and PEAR losses with hierarchical clustering, especially with average linkage.  Medvedovic clustering and variation of information correct the overestimation of the Dirichlet process mixture model with diagonal covariance, giving numbers of clusters similar to the MFM model. However, as in Section <ref>, summarisation methods do not correct the overestimation due to misspecification.  The adjusted Rand indices show that the full covariance model is more accurate than the misspecified diagonal covariance model.  We cannot expect perfect recovery of the true seven cell types, even with the full covariance model, because we use a very general model without consideration of the data type and without further prior assumptions.

1c     4cNumber of clusters      4cAdjusted Rand index

1cV1.5     2cFull      2cV1.5Diagonal
2cFull      2cV1.5Diagonal

1cV1.5     1cMFM     DPM     MFM     DPM     MFM     DPM     MFM     DPM

Medvedovic     8     8     15     16     0.38     0.35     0.34     0.35

Binder+Average     30     36     59     55     0.38     0.38     0.33     0.33

Binder+Complete     11     12     19     17     0.38     0.38     0.34     0.35

Binder+Samples     9     9     15     18     0.37     0.38     0.34     0.34

PEAR+Average     25     29     33     32     0.38     0.38     0.33     0.33

PEAR+Complete     11     12     18     17     0.38     0.38     0.34     0.35

PEAR+Samples     9     9     15      18     0.37     0.38     0.34      0.34

VI-LB+Average     9     8     14     15     0.38     0.38     0.33     0.35

VI-LB+Complete     9     9     15     15     0.38     0.38     0.34     0.35

VI-LB+Samples     9     8     15     19     0.37     0.37     0.34     0.33

Posterior mode     8     8     16     18     4c

gray!15 Truth      gray!157     gray!157     gray!157     gray!157     4c

Number of clusters and adjusted Rand indices of summary clusterings for the MFM and DPM models with full or diagonal covariance on mouse single-cell RNA-sequence data. The adjusted Rand index compares the summary clustering with the true clustering with seven clusters given by [Zeisel et al., 2015].

We use the Uniform Manifold Approximation and Projection (UMAP) [McInnes et al., 2020] to visualise our clusterings. Figure <ref> illustrates summary clusterings for the MFM model on a two-dimensional projection of our data, which we compare to the true clustering according to [Zeisel et al., 2015].
The misspecified model seems to create subclusters of the full covariance, i.e., the clusters of the diagonal covariance model seem to split up the larger clusters of the full covariance model. For synthetic data, we showed in Section <ref> that the diagonal covariance model splits up clusters with correlated data. We confirm our conjecture in Figure <ref> by visualising the contingency counts for summary clusterings of the MFM model. By reordering the clusters of the diagonal covariance model, we roughly split up clusters of the full covariance model (rows) into clusters of the diagonal covariance (columns). For example, cluster 1 of the full diagonal covariance model can be split into clusters 4, 1 and 6 of the diagonal covariance model.

Full covariance model.

Diagonal covariance model.

True clusters.

Summary clusterings of variation of information average linkage hierarchical clustering for the MFM model with a full (a) or diagonal (b) covariance specification on the mouse single-cell RNA-sequence data. We visualise the 10-dimensional data with the Uniform Manifold Approximation and Projection (UMAP) on two dimensions, denoted by $( Y_1, Y_2)$. The true clusters (c) are given by [Zeisel et al., 2015].

Normalised contingency counts between summary clusterings of the MFM model with a full and diagonal covariance on the mouse single-cell RNA-sequence data. We reorder clusters of the diagonal covariance by their largest count and normalise columns for readability, so that entries correspond to proportions relative to the diagonal covariance cluster sizes. We use variation of information average linkage hierarchical clustering for both models.

In conclusion, our analysis of the mouse single-cell RNA-sequence data confirms some of our simulated results. Misspecification of the component distributions seems to have a bigger effect on the overestimation of the number of components than does the Dirichlet process prior. Summarisation methods can correct the overestimation of the DPM model compared to the MFM model but they do not correct for misspecification of the component distributions. However, misspecification of Gaussian component distributions seems to split components of their well-specified counterpart, giving reasonably interpretable summary clusterings.

§ DISCUSSION

Our simulation results confirm theoretical inconsistencies for the number of components of Bayesian mixture models from recent literature for finite samples. MCMC summarisation methods may correct estimation errors, but the summarisation method matters, as some popular methods can lead to severe overestimation. Our gene expression analysis supports these findings and shows the importance of appropriate model specification.

Dirichlet process mixture models overestimate the number of components asymptotically in the sample size [Miller and Harrison, 2014], but  our results show limited overestimation for practical sample sizes. Moreover, some summarisation methods consistently correct overestimation and give the correct number of clusters. In our simulations, for example, variation of information [Wade and Ghahramani, 2018] consistently led to accurate summary clusterings. We believe that summarisation methods can find the true number of clusters because additional clusters introduced by the DPM model tend to be small and diffused across MCMC samples.
We therefore expect that past studies using this model are mostly unaffected by these theoretical drawbacks, especially if they use summary clusterings for final interpretations.
Nevertheless, we recommend the mixture of finite mixtures model in general, due to its theoretical consistency and equally efficient sampling methods.

Bayesian mixture models with misspecified component distributions overestimate the number of components asymptotically in the sample size [Cai et al., 2021], and our results show that this overestimation can be large for Gaussian mixture models, and cannot be corrected by summary clusterings.
Nonetheless, we conjecture that misspecified mixture models tend to only split true clusters, leading to moderate errors in interpretation.  However, our results are limited to misspecification by constraining component distributions to be less expressive, rather than specifying entirely different families of distributions.

Some summary clustering methods lead to overestimation of the number of clusters, even for well-specified models and simple examples, where MCMC samples concentrate around the true number of clusters. We have shown that summary clusterings obtained with Binder's loss [Binder, 1978] or the posterior expected average adjusted Rand index [Fritsch and Ickstadt, 2009] overestimate the true number of clusters on synthetic and real data, when using common optimisation methods. These results agree with those from simulations in [Rastelli and Friel, 2018] and [Wade and Ghahramani, 2018], where they also observed overestimation for Binder's loss (using equal penalties of misclassification) and correct estimation by the lower bound of the variation of information loss. Moreover, we have shown that the choice of optimisation methods for summary clustering can greatly influence the number of clusters. Here we have compared some of the most popular methods, but acknowledge that there are many such methods <cit.> and hence our analyses are not exhaustive. By exploring a bigger part of the space of all possible partitions, these more recent methods are able to further minimise the loss function, but interestingly this seems to have a negative impact on the estimation of the number of clusters.

Our results only consider Gaussian mixture models, so we could extend our work to other continuous models, such as Laplace location mixtures or regression models, or to discrete mixture models, such as multinomial mixture models. Such extensions, however, would require other Gibbs updates for split-merge sampling. An alternative would be to generate data from a different family of distributions and keep using Gaussian mixture models.

Our results show that the summary number of clusters can greatly depend on the summarisation method, so theoretical work might help explain under what conditions different methods perform best.

Finally, when misspecification is chosen over efficiency, we could investigate guided summarisation methods that account for misspecification by merging additional clusters. Introducing intuitive adjustable settings to summarisation methods would accommodate further assumptions more generally, and allow the user to choose  appropriate representative clusterings  according to the purpose of subsequent analysis.

PDWK and JvdMM acknowledge MRC grant MC_UU_00002_13.  This work was supported by the Swiss National Science Foundation and the National Institute for Health Research [Cambridge Biomedical Research Centre at the Cambridge University Hospitals NHS Foundation Trust]. The views expressed are those of the authors and not necessarily those of the NHS, the NIHR or the Department of Health and Social Care.

urlstyle[Alfò et al., 2008]
M. Alfò, G. Trovato, and R. J. Waldmann.
Testing for country heterogeneity in growth models using a finite
mixture approach.
Journal of Applied Econometrics, 230 (4):0
487–514, 2008.
URL <https://onlinelibrary.wiley.com/doi/abs/10.1002/jae.1008>.

[Ascolani et al., 2022]
F. Ascolani, A. Lijoi, G. Rebaudo, and G. Zanella.
Clustering consistency with Dirichlet process mixtures, May 2022.
URL <http://arxiv.org/abs/2205.12924>.
Number: arXiv:2205.12924 arXiv:2205.12924 [math, stat].

[Biernacki et al., 2000]
C. Biernacki, G. Celeux, and G. Govaert.
Assessing a mixture model for clustering with the integrated
completed likelihood.
IEEE Transactions on Pattern Analysis and Machine
Intelligence, 220 (7):0 719–725, 2000.

[Binder, 1978]
D. A. Binder.
Bayesian cluster analysis.
Biometrika, 650 (1):0 31–38, 04 1978.
ISSN 0006-3444.
URL <https://doi.org/10.1093/biomet/65.1.31>.

[Blackwell and MacQueen, 1973]
D. Blackwell and J. B. MacQueen.
Ferguson distributions via Pólya urn schemes.
The Annals of Statistics, 10 (2):0 353–355,
URL <https://doi.org/10.1214/aos/1176342372>.

[Cai et al., 2021]
D. Cai, T. Campbell, and T. Broderick.
Finite mixture models do not reliably learn the number of components.
In M. Meila and T. Zhang, editors, Proceedings of the 38th
International Conference on Machine Learning, volume 139 of
Proceedings of Machine Learning Research, pages 1158–1169, 18–24 Jul
URL <https://proceedings.mlr.press/v139/cai21a.html>.

[Carmona et al., 2019]
C. Carmona, L. Nieto-Barajas, and A. Canale.
Model-based approach for household clustering with mixed scale
Advances in Data Analysis and Classification, 130
(2):0 559–583, June 2019.

[Celeux et al., 2006]
G. Celeux, F. Forbes, C. P. Robert, and D. M. Titterington.
Deviance information criteria for missing data models.
Bayesian Analysis, 10 (4):0 651–673, 2006.
URL <https://doi.org/10.1214/06-BA122>.

[Celeux et al., 2019]
G. Celeux, S. Frühwirth-Schnatter, and C. P. Robert.
Model selection for mixture models – perspectives and strategies.
In G. Celeux, S. Frühwirth-Schnatter, and C. P. Robert, editors,
Handbook of Mixture Analysis, pages 117–154. Chapman and Hall/CRC,

[Chan et al., 2008]
C. Chan, F. Feng, J. Ottinger, D. Foster, M. West, and T. B. Kepler.
Statistical mixture modeling for cell subtype identification in flow
Cytometry Part A, 73A0 (8):0 693–701, 2008.

[Chandra et al., 2021]
N. K. Chandra, A. Canale, and D. B. Dunson.
Escaping the curse of dimensionality in Bayesian model based
clustering, 2021.
arXiv:2006.02700 [stat.ME].

[Crook et al., 2018]
O. M. Crook, C. M. Mulvey, P. D. W. Kirk, K. S. Lilley, and L. Gatto.
A Bayesian mixture modelling approach for spatial proteomics.
PLOS Computational Biology, 140 (11):0
e1006516, Nov. 2018.
ISSN 1553-7358.
URL <https://dx.plos.org/10.1371/journal.pcbi.1006516>.

[Dahl, 2006]
D. B. Dahl.
Model-based clustering for expression data via a Dirichlet process
mixture model.
In K.-A. Do, P. Müller, and M. Vannucci, editors, Bayesian
Inference for Gene Expression and Proteomics, pages 201–218. Cambridge
University Press, 2006.

[Dahl et al., 2021]
D. B. Dahl, D. J. Johnson, and P. Mueller.
Search algorithms and loss functions for Bayesian clustering, 2021.
arXiv:2105.04451 [stat.ME].

[de Souto et al., 2008]
M. C. de Souto, I. G. Costa, D. S. de Araujo, T. B. Ludermir, and A. Schliep.
Clustering cancer gene expression data: a comparative study.
BMC Bioinformatics, 90 (1):0 497, Dec. 2008.
ISSN 1471-2105.

[Dias et al., 2015]
J. G. Dias, J. K. Vermunt, and S. Ramos.
Clustering financial time series: New insights from an extended
hidden Markov model.
European Journal of Operational Research, 2430
(3):0 852–864, 2015.
ISSN 0377-2217.

[Drton and Plummer, 2017]
M. Drton and M. Plummer.
A Bayesian information criterion for singular models.
Journal of the Royal Statistical Society. Series B (Statistical
Methodology), 790 (2):0 323–380, 2017.
ISSN 13697412, 14679868.
URL <http://www.jstor.org/stable/44682516>.

[Dudoit and Fridlyand, 2002]
S. Dudoit and J. Fridlyand.
A prediction-based resampling method for estimating the number of
clusters in a dataset.
Genome Biology, 30 (7):0 1–21, 2002.
URL <https://doi.org/10.1186/gb-2002-3-7-research0036>.

[Escobar and West, 1995]
M. D. Escobar and M. West.
Bayesian density estimation and inference using mixtures.
Journal of the American Statistical Association, 900
(430):0 577–588, 1995.
ISSN 01621459.
URL <http://www.jstor.org/stable/2291069>.

[Fritsch and Ickstadt, 2009]
A. Fritsch and K. Ickstadt.
Improved criteria for clustering based on the posterior similarity
Bayesian Analysis, 40 (2):0 367–391, 2009.
URL <https://doi.org/10.1214/09-BA414>.

[Frühwirth-Schnatter and Malsiner-Walli, 2019]
S. Frühwirth-Schnatter and G. Malsiner-Walli.
From here to infinity: sparse finite versus Dirichlet process
mixtures in model-based clustering.
Advances in Data Analysis and Classification, 130
(1):0 33–64, 2019.
ISSN 1862-5347.
URL <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6448299/>.

[Frühwirth-Schnatter et al., 2004]
S. Frühwirth-Schnatter, R. Tüchler, and T. Otter.
Bayesian analysis of the heterogeneity model.
Journal of Business & Economic Statistics, 220
(1):0 2–15, 2004.
URL <https://doi.org/10.1198/073500103288619331>.

[Fuentes-García et al., 2019]
R. Fuentes-García, R. H. Mena, and S. G. Walker.
Modal posterior clustering motivated by Hopfield’s network.
Computational Statistics & Data Analysis, 137:0
92–100, 2019.
ISSN 0167-9473.

[Gaffney et al., 2007]
S. J. Gaffney, A. W. Robertson, P. Smyth, S. J. Camargo, and M. Ghil.
Probabilistic clustering of extratropical cyclones using regression
mixture models.
Climate dynamics, 290 (4):0 423–440, 2007.

[Gelman and Rubin, 1992]
A. Gelman and D. B. Rubin.
Inference from iterative simulation using multiple sequences.
Statistical Science, 70 (4):0 457–472, 1992.
URL <https://doi.org/10.1214/ss/1177011136>.

[Geng et al., 2019]
J. Geng, A. Bhattacharya, and D. Pati.
Probabilistic community detection with unknown number of communities.
Journal of the American Statistical Association, 1140
(526):0 893–905, 2019.
URL <https://doi.org/10.1080/01621459.2018.1458618>.

[Geweke, 1991]
J. F. Geweke.
Evaluating the accuracy of sampling-based approaches to the
calculation of posterior moments.
Staff Report 148, Federal Reserve Bank of Minneapolis, 1991.
URL <https://ideas.repec.org/p/fip/fedmsr/148.html>.

[Ghosal, 2010]
S. Ghosal.
The Dirichlet process, related priors and posterior asymptotics.
In N. L. Hjort, C. Holmes, P. Müller, and S. G. Walker, editors,
Bayesian Nonparametrics, Cambridge Series in Statistical and
Probabilistic Mathematics, pages 35–79. Cambridge University Press,
Cambridge, 2010.

[Guha et al., 2019]
A. Guha, N. Ho, and X. Nguyen.
On posterior contraction of parameters and interpretability in
Bayesian mixture modeling, 2019.
arXiv:1901.05078 [math.ST].

[Henna, 2008]
J. Henna.
An estimator of the number of components of a finite mixture of
multivariate distributions.
Journal of the Japan Statistical Society, 380
(3):0 379–389, 2008.
URL <https://doi.org/10.14490/jjss.38.379>.

[Hubert and Arabie, 1985]
L. Hubert and P. Arabie.
Comparing partitions.
Journal of Classification, 20 (1):0 193–218,

[Huelsenbeck and Andolfatto, 2007]
J. P. Huelsenbeck and P. Andolfatto.
Inference of population structure under a Dirichlet process model.
Genetics, 1750 (4):0 1787–1802, 04 2007.
ISSN 1943-2631.
URL <https://doi.org/10.1534/genetics.106.061317>.

[Jain and Neal, 2004]
S. Jain and R. M. Neal.
A split-merge Markov chain Monte Carlo procedure for the Dirichlet
process mixture model.
Journal of Computational and Graphical Statistics, 130
(1):0 158–182, 2004.
URL <https://doi.org/10.1198/1061860043001>.

[James et al., 2001]
L. F. James, D. J. Marchette, and C. E. Priebe.
Consistent estimation of mixture complexity.
The Annals of Statistics, 290 (5):0
1281–1296, 2001.
URL <https://doi.org/10.1214/aos/1013203454>.

[Jasra et al., 2005]
A. Jasra, C. C. Holmes, and D. A. Stephens.
Markov chain Monte Carlo methods and the label switching problem in
Bayesian mixture modeling.
Statistical Science, 200 (1):0 50–67, 2005.
URL <https://doi.org/10.1214/088342305000000016>.

[Kaufman and Rousseeuw, 1990]
L. Kaufman and P. J. Rousseeuw.
Partitioning Around Medoids (Program PAM), chapter 2, pages
John Wiley & Sons, Ltd, 1990.

[Kirk et al., 2012]
P. Kirk, J. E. Griffin, R. S. Savage, Z. Ghahramani, and D. L. Wild.
Bayesian correlated clustering to integrate multiple datasets.
Bioinformatics, 280 (24):0 3290–3297, 10
ISSN 1367-4803.
URL <https://doi.org/10.1093/bioinformatics/bts595>.

[Lanza and Rhoades, 2013]
S. T. Lanza and B. L. Rhoades.
Latent Class Analysis: An alternative perspective on subgroup
analysis in prevention and treatment.
Prevention Science, 140 (2):0 157–168, Apr.
ISSN 1389-4986, 1573-6695.
URL <http://link.springer.com/10.1007/s11121-011-0201-1>.

[Legramanti et al., 2021]
S. Legramanti, T. Rigon, D. Durante, and D. B. Dunson.
Extended stochastic block models with application to criminal
networks, 2021.
arXiv:2007.08569 [stat.ME].

[Liu and Shao, 2003]
X. Liu and Y. Shao.
Asymptotics for likelihood ratio tests under loss of
The Annals of Statistics, 310 (3):0 807–832,
URL <https://doi.org/10.1214/aos/1056562463>.

[Liverani et al., 2015]
S. Liverani, D. I. Hastie, L. Azizi, M. Papathomas, and S. Richardson.
PReMiuM: An R package for profile regression mixture models using
Dirichlet processes.
Journal of Statistical Software, 640 (7):0
1–30, 2015.

[Lock and Dunson, 2013]
E. F. Lock and D. B. Dunson.
Bayesian consensus clustering.
Bioinformatics, 290 (20):0 2610–2616, 08
ISSN 1367-4803.
URL <https://doi.org/10.1093/bioinformatics/btt425>.

[MacEachern, 1994]
S. MacEachern.
Estimating normal means with a conjugate style Dirichlet process
Communications in Statistics - Simulation and Computation,
23:0 727–741, 01 1994.

[Manole and Khalili, 2021]
T. Manole and A. Khalili.
Estimating the number of components in finite mixture models via the
Group-Sort-Fuse procedure.
The Annals of Statistics, 490 (6):0
3043–3069, 2021.
URL <https://doi.org/10.1214/21-AOS2072>.

[McInnes et al., 2020]
L. McInnes, J. Healy, and J. Melville.
UMAP: Uniform Manifold Approximation and Projection for dimension
reduction, 2020.
arXiv:1802.03426 [stat.ML].

[McLachlan, 1987]
G. J. McLachlan.
On bootstrapping the likelihood ratio test stastistic for the number
of components in a normal mixture.
Journal of the Royal Statistical Society. Series C (Applied
Statistics), 360 (3):0 318–324, 1987.
ISSN 00359254, 14679876.
URL <http://www.jstor.org/stable/2347790>.

[McLachlan et al., 2002]
G. J. McLachlan, R. W. Bean, and D. Peel.
A mixture model-based approach to the clustering of microarray
expression data.
Bioinformatics, 180 (3):0 413–422, 03 2002.
ISSN 1367-4803.
URL <https://doi.org/10.1093/bioinformatics/18.3.413>.

[Medvedovic and Sivaganesan, 2002]
M. Medvedovic and S. Sivaganesan.
Bayesian infinite mixture model based clustering of gene expression
Bioinformatics, 180 (9):0 1194–1206, 09 2002.
ISSN 1367-4803.
URL <https://doi.org/10.1093/bioinformatics/18.9.1194>.

[Medvedovic et al., 2004]
M. Medvedovic, K. Yeung, and R. Bumgarner.
Bayesian mixture model based clustering of replicated microarray
Bioinformatics, 200 (8):0 1222–1232, Feb.
URL <https://doi.org/10.1093/bioinformatics/bth068>.

[Meila, 2007]
M. Meila.
Comparing clusterings—an information based distance.
Journal of Multivariate Analysis, 980 (5):0
873–895, 2007.
ISSN 0047-259X.

[Miller and Harrison, 2014]
J. W. Miller and M. T. Harrison.
Inconsistency of Pitman–Yor process mixtures for the number of
Journal of Machine Learning Research, 150
(96):0 3333–3370, 2014.
URL <http://jmlr.org/papers/v15/miller14a.html>.

[Miller and Harrison, 2018]
J. W. Miller and M. T. Harrison.
Mixture models with a prior on the number of components.
Journal of the American Statistical Association, 1130
(521):0 340–356, 2018.
URL <https://doi.org/10.1080/01621459.2016.1255636>.

[Mukherjee et al., 1998]
S. Mukherjee, E. D. Feigelson, G. J. Babu, F. Murtagh, C. Fraley, and
A. Raftery.
Three types of gamma-ray bursts.
The Astrophysical Journal, 5080 (1):0 314,

[Neal, 2000]
R. M. Neal.
Markov chain sampling methods for Dirichlet process mixture models.
Journal of Computational and Graphical Statistics, 90
(2):0 249–265, 2000.
ISSN 10618600.
URL <http://www.jstor.org/stable/1390653>.

[Nguyen, 2013]
X. Nguyen.
Convergence of latent mixing measures in finite and infinite mixture
The Annals of Statistics, 410 (1):0 370–400,
URL <https://doi.org/10.1214/12-AOS1065>.

[Ni et al., 2020]
Y. Ni, P. Müller, M. Diesendruck, S. Williamson, Y. Zhu, and Y. Ji.
Scalable Bayesian nonparametric clustering and classification.
Journal of Computational and Graphical Statistics, 290
(1):0 53–65, 2020.
URL <https://doi.org/10.1080/10618600.2019.1624366>.

[Nobile, 1994]
A. Nobile.
Bayesian analysis of finite mixture distributions.
Ph.D. dissertation, Carnegie Mellon University, United States –
Pennsylvania, 1994.
ISBN: 9798607381097.

[Nobile, 2004]
A. Nobile.
On the posterior distribution of the number of components in a
finite mixture.
The Annals of Statistics, 320 (5):0
2044–2073, 2004.
URL <https://doi.org/10.1214/009053604000000788>.

[Nobile, 2007]
A. Nobile.
Bayesian finite mixtures: a note on prior specification and posterior
arXiv:0711.0458 [stat], Nov. 2007.
URL <http://arxiv.org/abs/0711.0458>.

[Ohn and Lin, 2022]
I. Ohn and L. Lin.
Optimal Bayesian estimation of Gaussian mixtures with growing
number of components, Mar. 2022.
URL <http://arxiv.org/abs/2007.09284>.
Number: arXiv:2007.09284 arXiv:2007.09284 [math, stat].

[Onogi et al., 2011]
A. Onogi, M. Nurimoto, and M. Morita.
Characterization of a Bayesian genetic clustering algorithm based
on a Dirichlet process prior and comparison among Bayesian clustering
BMC Bioinformatics, 12:0 263, 06 2011.

[Prabhakaran et al., 2016]
S. Prabhakaran, E. Azizi, A. Carr, and D. Pe'er.
Dirichlet process mixture model for correcting technical variation in
single-cell gene expression data.
In M. F. Balcan and K. Q. Weinberger, editors, Proceedings of
the 33rd International Conference on International Conference on Machine
Learning - Volume 48, pages 1070––1079. JMLR, 2016.

[Rand, 1971]
W. M. Rand.
Objective criteria for the evaluation of clustering methods.
Journal of the American Statistical Association, 660
(336):0 846–850, 1971.

[Rastelli and Friel, 2018]
R. Rastelli and N. Friel.
Optimal Bayesian estimators for latent variable cluster models.
Statistics and Computing, 280 (6):0
1169–1186, nov 2018.
ISSN 0960-3174.
URL <https://doi.org/10.1007/s11222-017-9786-y>.

[Richardson and Green, 1997]
S. Richardson and P. J. Green.
On Bayesian analysis of mixtures with an unknown number of
components (with discussion).
Journal of the Royal Statistical Society. Series B (Statistical
Methodology), 590 (4):0 731–792, 1997.

[Roth et al., 2014]
A. Roth, J. Khattra, D. Yap, A. Wan, E. Laks, J. Biele, G. Ha, S. Aparicio,
A. Bouchard-Côté, and S. P. Shah.
PyClone: statistical inference of clonal population structure in
Nature Methods, 110 (4):0 396–398, Apr. 2014.
ISSN 1548-7091, 1548-7105.
URL <http://www.nature.com/articles/nmeth.2883>.

[Rousseau and Mengersen, 2011]
J. Rousseau and K. Mengersen.
Asymptotic behaviour of the posterior distribution in overfitted
mixture models.
Journal of the Royal Statistical Society. Series B (Statistical
Methodology), 730 (5):0 689–710, 2011.
ISSN 13697412, 14679868.
URL <http://www.jstor.org/stable/41262270>.

[Savage et al., 2010]
R. S. Savage, Z. Ghahramani, J. E. Griffin, B. J. de la Cruz, and D. L. Wild.
Discovering transcriptional modules by Bayesian data integration.
Bioinformatics, 260 (12):0 158–167, 06 2010.
ISSN 1367-4803.
URL <https://doi.org/10.1093/bioinformatics/btq210>.

[Schubert and Rousseeuw, 2019]
E. Schubert and P. J. Rousseeuw.
Faster $k$-medoids clustering: Improving the PAM, CLARA, and
CLARANS algorithms.
In G. Amato, C. Gennaro, V. Oria, and M. Radovanović, editors,
Similarity Search and Applications, pages 171–187. Springer
International Publishing, 2019.

[Sørlie et al., 2001]
T. Sørlie, C. M. Perou, R. Tibshirani, T. Aas, S. Geisler, H. Johnsen,
T. Hastie, M. B. Eisen, M. van de Rijn, S. S. Jeffrey, T. Thorsen, H. Quist,
J. C. Matese, P. O. Brown, D. Botstein, P. E. Lønning, and A.-L.
Gene expression patterns of breast carcinomas distinguish tumor
subclasses with clinical implications.
Proceedings of the National Academy of Sciences, 980
(19):0 10869–10874, 2001.
ISSN 0027-8424.
URL <https://www.pnas.org/content/98/19/10869>.

[Vehtari et al., 2020]
A. Vehtari, A. Gelman, D. P. Simpson, B. Carpenter, and P.-C. Burkner.
Rank-normalization, folding, and localization: An improved $\hat R$
for assessing convergence of MCMC (with discussion).
Bayesian Analysis, 160 (2):0 667–718, 2020.

[Wade and Ghahramani, 2018]
S. Wade and Z. Ghahramani.
Bayesian cluster analysis: Point estimation and credible balls
(with discussion).
Bayesian Analysis, 130 (2):0 559–626, 2018.
URL <https://doi.org/10.1214/17-BA1073>.

[Wedel and Desarbo, 2002]
M. Wedel and W. S. Desarbo.
Market segment derivation and profiling via a finite mixture model
Marketing Letters, 130 (1):0 17–25, 2002.
ISSN 09230645, 1573059X.
URL <http://www.jstor.org/stable/40216618>.

[Woo and Sriram, 2006]
M.-J. Woo and T. N. Sriram.
Robust estimation of mixture complexity.
Journal of the American Statistical Association, 1010
(476):0 1475–1486, 2006.
ISSN 01621459.
URL <http://www.jstor.org/stable/27639766>.

[Yang et al., 2020]
C.-Y. Yang, E. Xia, N. Ho, and M. I. Jordan.
Posterior Distribution for the Number of Clusters in
Dirichlet Process Mixture Models.
arXiv:1905.09959 [cs, math, stat], Oct. 2020.
URL <http://arxiv.org/abs/1905.09959>.

[Yuan et al., 2015]
K. Yuan, T. Sakoparnig, F. Markowetz, and N. Beerenwinkel.
BitPhylogeny: a probabilistic framework for reconstructing
intra-tumor phylogenies.
Genome Biology, 160 (1):0 36, Dec. 2015.
ISSN 1474-760X.

[Zafar et al., 2019]
H. Zafar, N. Navin, K. Chen, and L. Nakhleh.
SiCloneFit: Bayesian inference of population structure, genotype,
and phylogeny of tumor clones from single-cell genome sequencing data.
Genome Research, 290 (11):0 1847–1859, Nov.
ISSN 1088-9051, 1549-5469.
URL <http://genome.cshlp.org/lookup/doi/10.1101/gr.243121.118>.

[Zeisel et al., 2015]
A. Zeisel, A. B. Muñoz-Manchado, S. Codeluppi, P. Lönnerberg, G. L. Manno,
A. Juréus, S. Marques, H. Munguba, L. He, C. Betsholtz, C. Rolny,
G. Castelo-Branco, J. Hjerling-Leffler, and S. Linnarsson.
Cell types in the mouse cortex and hippocampus revealed by
single-cell RNA-seq.
Science, 3470 (6226):0 1138–1142, 2015.
URL <https://www.science.org/doi/abs/10.1126/science.aaa1934>.

§ APPENDIX A.

Number of components of MCMC sample allocation for a mixture of finite mixtures with (MFMH) and without (MFM) a hyperprior and for a Dirichlet process mixture with (DPMH) and without (DPM) a hyperprior, jittered for visibility, on synthetic data with $N = 500$. The 100 samples in red are burn-in and not included in posterior analysis.

Adjusted Rand index of summary clusterings with the true generating cluster allocation for the MFM model without a hyperprior on $50$ independent synthetic datasets with $N = 500$. Summary names are composed of the loss function and the optimisation method, except for Medvedovic clustering.

Example of a MCMC allocation sample for the Dirichlet process mixture model with a hyperprior on synthetic data with $N= 10^4$. Observations are coloured by their cluster allocation. One can discern six clusters, including a green one on the right and a blue one roughly aligned on $X_1 = 7.5$.

Summary clusterings of PEAR loss (a) and variation of information (b) both optimised with average linkage hierarchical clustering on MCMC samples of the DPM model with a hyperprior on synthetic data with $N=10^4$.


11institutetext: Ping An Health Technology, Beijing, China 22institutetext:
Ping An Property & Casualty Insurance Company, Shenzhen, China
33institutetext: Central University of Finance and Economics, Beijing, China
44institutetext: Ping An Health Cloud Company Limited., Shenzhen, China
55institutetext: Ping An International Smart City Technology Co., Ltd.,
Shenzhen, China

# PASH at TREC 2021 Deep Learning Track: Generative Enhanced Model for Multi-
stage Ranking

Yixuan Qiao 11** Hao Chen 11 Jun Wang 11 Yongquan Lai 22 Tuozhen Liu 33
Xianbin Ye 11 Xin Tang 22 Rui Fang 22 Wenfeng Xie 22 Guotong Xie 114455**

###### Abstract

This paper describes the PASH participation in TREC 2021 Deep Learning Track.
In the recall stage, we adopt a scheme combining sparse and dense retrieval
method. In the multi-stage ranking phase, point-wise and pair-wise ranking
strategies are used one after another based on model continual pre-trained on
general knowledge and document-level data. Compared to TREC 2020 Deep Learning
Track, we have additionally introduced the generative model T5 to further
enhance the performance.

###### Keywords:

expansion method dense retrieval transfer learning multi-stage ranking model
parallel.

††∗Corresponding Author. Email: {qiaoyixuan528<EMAIL_ADDRESS>
##  1 Introduction

TREC 2021 is the thirtieth edition of the Text REtrieval Conference (TREC).
The main goal of TREC is to create the evaluation infrastructure required for
large-scale testing of information retrieval (IR) technology. This includes
research on best methods for evaluation as well as development of the
evaluation materials themselves. “Retrieval technology” is broadly interpreted
to include a variety of techniques that enable and/or facilitate access to
information that is not specifically structured for machine use. Specifically,
the Deep Learning track focuses on IR tasks where a large training set is
available, allowing us to compare a variety of retrieval approaches including
deep neural networks and strong non-neural approaches, to see what works best
in a large-data regime.

The Deep Learning track had two tasks, Document Ranking and Passage Ranking.
For each task, participants could either do their own retrieval from the full
collection or re-rank an initial retrieved set provided by the track
organizers. Both tasks used the same set of 477 test questions.

##  2 Methodology

In this section, we briefly describe the components in the multi-stage ranking
pipeline. Details of some methods can be found in TREC 2020 notebook paper.
The passage ranking task contains 285,328 queries on a collection of
138,364,198 passages, totally 296,419 query passage pairs annotated as
positive for relevance. Few queries are matched with multiple relevant
passages. The document ranking task contains 331,748 queries on a collection
of 11,959,635 documents, totally 341,836 query document pairs annotated as
positive for relevance. Few queries are matched with multiple relevant
documents. NDCG (Normalized Discounted Cumulative Gain) is the main official
evaluation index.

###  2.1 Multi-way Matching

####  2.1.1 Sparse retrieval

We use the docTTTTTquery[1] (also known as docT5query or doc2query-T5) to
generate queries for which the passage might be relevant. For passage, We
sample 40 queries per passage with T5-base. For document, we first segment
each document into passages by applying a sliding window of ten sentences with
a stride of five. Each passage is then prepended with the url and title of the
document. Finally, we generate ten queries per passage. Each document
prepended with the url and title info was appended with all its predicted
queries are indexed by BM25[2] as before.

####  2.1.2 Dense retrieval

Dense retrieval methods have shown great promise over sparse retreval methods
in a range of NLP (Natural Language Processing) problems. We choose ColBERT[3]
as our Dense Retriever due to the effectiveness/efficiency trade-offs.

###  2.2 Multi-stage Ranking

We inherit the model from TREC 2020 Deep Learning Track including BERT-
Large[4], ALBERT-XXLarge[5], ELECTRA-Large[6] and XLNet-Large[7] model. In
particular, BERT-Large uses the same pre-training and fine-tuning strategies
as last year but using this year’s data set. Compared with last year, we have
made the following further improvements:

  * •

We use pair-wise loss as the goal of the second stage of ranking. So we can
use more data instead of TREC labeled 4-class data.

  * •

Generative model T5[8] is introduced to further enhance the performance

    * a)

Implementation based on Megatron-LM[9] repository.

    * b)

Same training strategy as [10].

    * c)

Training the T5-3B model takes approximately 12 days on 32 V100 GPUs consuming
about 12.8M training examples.

    * d)

The T5-11B model takes about 30 days under the same situation. An amazing fact
is that the zero-shot NDCG@5 of the T5-11B is already 0.3651.

###  2.3 Ensemble

Benefited from multi-way matching, we just sequentially train different re-
rankers with different random seeds for ensemble learning.

##  3 Results

We submitted six official runs in each passage and document ranking sub-task.
Table 1 and 2 present the results of our runs in passage ranking task. Table 3
and 4 present the results of our runs in document ranking task. We use
$\\_$r${\ast}$ to denote the runs for re-ranking, $\\_$f${\ast}$ for full
ranking.

Table 1: Results on the TREC 2021 passage re-rank task. RUN | MAP@100 | NDCG@5 | NDCG@10 | R@100
---|---|---|---|---
pash_r1 | 0.2362 | 0.7190 | 0.6951 | 0.3261
pash_r2 | 0.2389 | 0.7390 | 0.7076 | 0.3261
pash_r3 | 0.2385 | 0.7390 | 0.7072 | 0.3261
Table 2: Results on the TREC 2021 passage full rank task. RUN | MAP@100 | NDCG@5 | NDCG@10 | R@100
---|---|---|---|---
pash_f1 | 0.3193 | 0.7596 | 0.7494 | 0.4868
pash_f2 | 0.3318 | 0.7596 | 0.7494 | 0.5300
pash_f3 | 0.3378 | 0.7596 | 0.7494 | 0.5500
Table 3: Results on the TREC 2021 document re-rank task. RUN | MAP@100 | NDCG@5 | NDCG@10 | R@100
---|---|---|---|---
pash_doc_r1 | 0.2665 | 0.7452 | 0.7150 | 0.3195
pash_doc_r2 | 0.2640 | 0.7308 | 0.7076 | 0.3195
pash_doc_r3 | 0.2672 | 0.7383 | 0.7164 | 0.3195
Table 4: Results on the TREC 2021 document full rank task. RUN | MAP@100 | NDCG@5 | NDCG@10 | R@100
---|---|---|---|---
pash_doc_f1 | 0.3111 | 0.6745 | 0.5832 | 0.3683
pash_doc_f4 | 0.3498 | 0.7516 | 0.7404 | 0.4412
pash_doc_f5 | 0.3521 | 0.7508 | 0.7368 | 0.4419

##  4 Conclusion

The experiments demonstrate that the expressiveness of ColBERT is well
integrated with BM25+docTTTTTquery. Compared with the transformer-based model,
the generative model T5 has a higher NDCG@5 and a lower NDCG@10, which means
that relevant passages are already in top5. This phenomenon largely hints at
the thirst for multi-stage ranking in generative models, which can be expected
to show powerful effects with the pair-wise or list-wise training. We believe
that the interweaving sequence of different methods would be a focus priority
in the future work.

## References

  * [1] Nogueira, R., Yang, W., Lin, J. & Cho, K.  Document expansion by query prediction.  _arXiv preprint arXiv:1904.08375_ (2019).
  * [2] Zeng, Z. & Sakai, T.  Bm25 pseudo relevance feedback using anserini at waseda university.  In _OSIRRC@ SIGIR_ , 62–63 (2019).
  * [3] Khattab, O. & Zaharia, M.  Colbert: Efficient and effective passage search via contextualized late interaction over bert.  In _Proceedings of the 43rd International ACM SIGIR conference on research and development in Information Retrieval_ , 39–48 (2020).
  * [4] Devlin, J., Chang, M.-W., Lee, K. & Toutanova, K.  Bert: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , 4171–4186 (2019).
  * [5] Lan, Z. _et al._ Albert: A lite bert for self-supervised learning of language representations.  In _International Conference on Learning Representations_ (2019).
  * [6] Clark, K., Luong, M.-T., Le, Q. V. & Manning, C. D.  Electra: Pre-training text encoders as discriminators rather than generators.  _arXiv preprint arXiv:2003.10555_ (2020).
  * [7] Yang, Z. _et al._ Xlnet: Generalized autoregressive pretraining for language understanding.  _Advances in neural information processing systems_ 32 (2019).
  * [8] Raffel, C. _et al._ Exploring the limits of transfer learning with a unified text-to-text transformer.  _Journal of Machine Learning Research_ 21, 1–67 (2020).
  * [9] Narayanan, D. _et al._ Efficient large-scale language model training on gpu clusters using megatron-lm.  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , 1–15 (2021).
  * [10] Pradeep, R., Nogueira, R. & Lin, J.  The expando-mono-duo design pattern for text ranking with pretrained sequence-to-sequence models.  _arXiv preprint arXiv:2101.05667_ (2021).