Centered Partition Process: Informative Priors for Clustering

TL;DR

This paper introduces a Centered Partition process that incorporates prior expert knowledge into Bayesian clustering, allowing for more informed and flexible partitioning, demonstrated through simulations and an epidemiological case study.

Contribution

It proposes a novel Centered Partition prior that modifies EPPF to include prior knowledge, enhancing Bayesian clustering methods.

Findings

The CP prior effectively incorporates prior knowledge into clustering.

The methodology performs well in simulations and real epidemiological data.

The approach offers a flexible way to include expert input in Bayesian models.

Abstract

There is a very rich literature proposing Bayesian approaches for clustering starting with a prior probability distribution on partitions. Most approaches assume exchangeability, leading to simple representations in terms of Exchangeable Partition Probability Functions (EPPF). Gibbs-type priors encompass a broad class of such cases, including Dirichlet and Pitman-Yor processes. Even though there have been some proposals to relax the exchangeability assumption, allowing covariate-dependence and partial exchangeability, limited consideration has been given on how to include concrete prior knowledge on the partition. For example, we are motivated by an epidemiological application, in which we wish to cluster birth defects into groups and we have prior knowledge of an initial clustering provided by experts. As a general approach for including such prior knowledge, we propose a Centered Partition (CP) process that modifies the EPPF to favor partitions close to an initial one. Some properties of the CP prior are described, a general algorithm for posterior computation is developed, and we illustrate the methodology through simulation examples and an application to the motivating epidemiology study of birth defects.

Tables

Table 1. Table 3: Summary statistics of the distribution of congenital heart defects among cases. Defects are divided according the grouping provided from investigators.

Congenital Heart Defect	Abbreviation	Frequency	Percentage of cases
Septal
Atrial septal defect	ASD	765	0.15
Perimembranous ventricular septal defect	VSDPM	552	0.11
Atrial septal defect, type not specified	ASDNOS	225	0.04
Muscular ventricular septal defect	VSDMUSC	68	0.02
Ventricular septal defect, otherwise specified	VSDOS	12	0.00
Ventricular septal defect, type not specified	VSDNOS	8	0.00
Atrial septal defect, otherwise specified	ASDOS	4	0.00
Conotruncal
Tetralogy of Fallot	FALLOT	639	0.12
D-transposition of the great arteries	DTGA	406	0.08
Truncus arteriosus	COMMONTRUNCUS	61	0.01
Double outlet right ventricle	DORVTGA	35	0.01
Ventricular septal defect reported as conoventricular	VSDCONOV	32	0.01
D-transposition of the great arteries, other type	DORVOTHER	22	0.00
Interrupted aortic arch type B	IAATYPEB	13	0.00
Interrupted aortic arch, not otherwise specified	IAANOS	5	0.00
Left ventricular outflow
Hypoplastic left heart syndrome	HLHS	389	0.08
Coarctation of the aorta	COARCT	358	0.07
Aortic stenosis	AORTICSTENOSIS	224	0.04
Interrupted aortic arch type A	IAATYPEA	12	0.00
Right ventricular outflow
Pulmonary valve stenosis	PVS	678	0.13
Pulmonary atresia	PULMATRESIA	100	0.02
Ebstein anomaly	EBSTEIN	66	0.01
Tricuspid atresia	TRIATRESIA	46	0.01
Anomalous pulmonary venous return
Total anomalous pulmonary venous return	TAPVR	163	0.03
Partial anomalous pulmonary venous return	PAPVR	21	0.01
Atrioventricular septal defect
Atrioventricular septal defect	AVSD	112	0.02

Equations

$${\mathcal{S}}_{N,K}=\frac{1}{K!}\sum _{j=0}^K(-1{)}^j(\genfrac{}{}{0ex}{}{K}{j})(K-j{)}^N,$$

$$p(\mathbf{y}|\mathbf{c},\mathbf{\theta })=\prod _{k=1}^K\prod _{i\in B_k}p(y_i|{\theta }_k)=\prod _{k=1}^{K}p({\mathbf{y}}_k|{\theta }_k),$$

$$p(\mathbf{c}|\mathbf{y},\mathbf{\theta })\propto p(\mathbf{c})\prod _{k=1}^{K}p({\mathbf{y}}_k|{\theta }_k).$$

$$p(\mathbf{c}=\{B_1,\dots ,B_K\})\propto \prod _{j=1}^K\rho (B_j),$$

$$p(\mathbf{c}|{\mathbf{c}}_0,\psi )\propto p_0(\mathbf{c})e^{-\psi d(\mathbf{c},{\mathbf{c}}_0)},$$

$$s_l({\mathbf{c}}_0)=\{\mathbf{c}\in {\Pi }_N:d(\mathbf{c},{\mathbf{c}}_0)={\delta }_l\},l=0,1,\dots ,L.$$

$$\text{VI}(\mathbf{c},{\mathbf{c}}^{\prime })$$

$$=\sum _{j=1}^K\frac{{\lambda }_j}{N}\log \left(\frac{{\lambda }_j}{N}\right)+\sum _{l=1}^{K^{\prime }}\frac{{\lambda }_l^{\prime }}{N}\log \left(\frac{{\lambda }_l^{\prime }}{N}\right)-2\sum _{j=1}^K\sum _{l=1}^{K^{\prime }}\frac{{\lambda }_{jl}^{\wedge }}{N}\log \left(\frac{{\lambda }_{jl}^{\wedge }}{N}\right),$$

$$p(\mathbf{c}|{\mathbf{c}}_0,\psi )=\frac{e^{-\psi {\delta }_l}}{{\sum }_{u=0}^L|s_u({\mathbf{c}}_0)|e^{-\psi {\delta }_u}},\text{for }\mathbf{c}\in s_l({\mathbf{c}}_0),l=0,1,\dots ,L,$$

$$\frac{N!}{{\prod }_{j=1}^K{\lambda }_j!{\prod }_{i=1}^N{f}_i!}.$$

$$p(\mathbf{c}|{\mathbf{c}}_0,\psi )=\frac{g({\mathbf{\lambda }}_m)e^{-\psi {\delta }_l}}{{\sum }_{u=0}^L{\sum }_{v=1}^M{n}_{uv}g({\mathbf{\lambda }}_v)e^{-\psi {\delta }_u}},\text{for }\mathbf{c}\in s_{lm}({\mathbf{c}}_0),$$

$$p(c_i=k|{\mathbf{c}}^{-i},\mathbf{\theta },y_i)\propto \{\begin{array}{ll}& p(c_i=k|{\mathbf{c}}^{-i})p(y_i|{\theta }_k)k=1,\dots ,K^{-}\\ & p(c_i=k|{\mathbf{c}}^{-i})\int p(y_i|\theta )d{G}_0(\theta )k=K^{-}+1,\end{array}$$

$$p(c_i=k|{\mathbf{c}}^{-i},{\mathbf{c}}_0,\psi )\propto p_0(c_i=k|{\mathbf{c}}^{-i})\exp \{-\psi d(\mathbf{c},{\mathbf{c}}_0)\}k=1,\dots ,K^{-},K^{-}+1$$

$$p(\delta ={\delta }_l)=\sum _{\mathbf{c}\in {\Pi }_N}p(\mathbf{c})\mathcal{I}\left\{d(\mathbf{c},{\mathbf{c}}_0)={\delta }_l)\right\}=\sum _{\mathbf{c}\in s_l({\mathbf{c}}_0)}p(\mathbf{c})l=0,\dots ,L,$$

$$p(\delta ={\delta }_l)=\sum _{\mathbf{c}\in s_l({\mathbf{c}}_0)}\frac{p_0(\mathbf{c})e^{-\psi {\delta }_l}}{{\sum }_{u=0}^L{\sum }_{{\mathbf{c}}^{\ast }\in s_u({\mathbf{c}}_0)}{p}_0({\mathbf{c}}^{\ast })e^{-\psi {\delta }_u}}l=0,\dots ,L,$$

$$p(\delta ={\delta }_l)=\frac{{\sum }_{m=1}^M{n}_{lm}g({\mathbf{\lambda }}_m)e^{-\psi {\delta }_l}}{{\sum }_{u=0}^L{\sum }_{v=1}^M{n}_{uv}g({\mathbf{\lambda }}_v)e^{-\psi {\delta }_u}},l=0,\dots ,L.$$

$$p(\delta ={\delta }_l)=\frac{n_l{e}^{-\psi {\delta }_l}}{{\sum }_{u=0}^L{n}_u{e}^{-\psi {\delta }_u}}l=0,\dots ,L.$$

$${\mathcal{B}}_N=e^{-1}\sum _{k=1}^{\infty }\frac{k^N}{k!},$$

$$P(K=k)=e^{-1}\frac{k^N}{{\mathcal{B}}_{N}k!},$$

$${\hat{n}}_l={\mathcal{B}}^{\ast }\frac{1}{R^{\ast }}\sum _{r=1}^{R^{\ast }}\mathcal{I}\left\{d({\mathbf{c}}^{(r)},{\mathbf{c}}_0)=\widehat{{\delta }_l}\right\},$$

$${\hat{n}}_{lm}={\mathcal{B}}^{\ast }\frac{1}{R^{\ast }}\sum _{r=1}^{R^{\ast }}\mathcal{I}\left\{d({\mathbf{c}}^{(r)},{\mathbf{c}}_0)=\widehat{{\delta }_l}\right\}\mathcal{I}\left\{\mathbf{\Lambda }({\mathbf{c}}^{(r)})={\mathbf{\lambda }}_m\right\}.$$

$$\log \left(\frac{\text{Pr}(y_{ij}=1|{\mathbf{x}}_{ij})}{\text{Pr}(y_{ij}=0|{\mathbf{x}}_{ij})}\right)=\text{logit}({\mathbf{\pi }}_{ij})={\alpha }_i+{\mathbf{x}}_{ij}^T{\mathbf{\beta }}_i,$$

$$y_{ij}$$

$${\alpha }_i$$

$$p(\mathbf{c})$$

$$p(K=k)=e^{-1}{k}^N/(k!{\mathcal{B}}_N),k\in \{1,\dots ,N\}.$$

$$-H(\mathbf{c})=$$

$$H(\mathbf{c}\wedge {\mathbf{c}}_0)=$$

$$\displaystyle+\frac{\lambda_{kc_{i,0}}^{-i}+1}{N}\log\left(\frac{\lambda_{kc_{i,0}}^{-i}+1}{N}\right)\Bigg{\}}$$

$$({\omega }_{ij}|-)\sim PG(1,{\alpha }_i+{\mathbf{x}}_{ij}^T{\mathbf{\beta }}_{c_i})$$

$$({\alpha }_i|-)\sim \mathcal{N}(a^{\ast },{\tau }^{\ast })$$

$$({\mathbf{\beta }}^{(k)}|-)\sim {\mathcal{N}}_p({\mathbf{b}}^{(k)},{\mathbf{Q}}^{(k)})$$

$$\text{Pr}(c_i=k|{\mathbf{c}}_{-i},-)\propto \text{Pr}(c_i=k|{\mathbf{c}}_{-i})\text{Pr}({\mathbf{y}}_i|{\mathbf{X}}_i,{\alpha }_i,c_i=k,{\mathbf{\beta }}^{(k)})$$

Full text

Bayesian clustering,

Bayesian nonparametrics,

centered process,

Dirichlet Process,

exchangeable probability partition function,

mixture model,

product partition model,

, , ,

and

1 Introduction

Clustering is one of the canonical data analysis goals in statistics. There are two main strategies that have been used for clustering; namely, distance and model-based clustering. Distance-based methods leverage upon a distance metric between data points and do not in general require a generative probability model of the data, while model-based methods rely on discrete mixture models, which model the data in different clusters as arising from kernels having different parameter values. The majority of the model-based literature leans on maximum likelihood estimation, commonly relying on the EM algorithm. Bayesian approaches that aim to approximate a full posterior distribution on the clusters have advantages in terms of uncertainty quantification, while also having the ability to incorporate prior information.

Although this article is motivated by providing a broad new class of methods for improving clustering performance in practice, we were initially motivated by a particular application involving birth defects epidemiology. In this context, there are $N$ different birth defects, which we can index using $i\in\{1,\ldots,N\}$, and we are interested in clustering these birth defects into mechanistic groups. This may be useful, for example, in that birth defects in the same group may have similar coefficients in logistic regression analysis relating different exposures to risk of developing the defect. Investigators have provided us with an initial partition $\boldsymbol{c}_{0}$ of the defects $\{1,\ldots,N\}$ into groups. It is appealing to combine this prior knowledge with information in the data from a grouped logistic regression to produce a posterior distribution on clusters, which characterizes uncertainty. The motivating question of this article is how to do this, with the resulting method ideally having broader impact to other types of centering of priors for clustering; for example, we may want to center the prior based on information on the number of clusters or cluster sizes.

With these goals in mind, we start by reviewing the relevant literature on clustering priors. Most of these methods assume exchangeability, which means that the prior probability of a partition $\boldsymbol{c}$ of $\{1,\ldots,N\}$ into clusters depends only on the number of clusters and the cluster sizes; the indices on the clusters play no role. Under the exchangeability assumption, one can define what is referred to in the literature as an Exchangeable Partition Probability Function (EPPF) (Pitman, 1995). This EPPF provides a prior distribution on the random partition $\boldsymbol{c}$. One direction to obtain a specific form for the EPPF is to start with a nonparametric Bayesian discrete mixture model with a prior for the mixing measure $P$, and then marginalize over this prior to obtain an induced prior on partitions. Standard choices for $P$, such as the Dirichlet (Ferguson, 1973) and Pitman-Yor process (Pitman and Yor, 1997), lead to relatively simple analytic forms for the EPPF. There has been some recent literature studying extensions to broad classes of Gibbs-type processes (Gnedin and Pitman, 2006; De Blasi et al., 2015), mostly focused on improving flexibility while maintaining the ability to predict the number of new clusters in a future sample of data.

There is also a rich literature on relaxing exchangeability in various ways. Most of the emphasis has been on the case in which a vector of features $\mathbf{x}_{i}$ is available for index $i$, motivating feature-dependent random partitions models. Building on the stick-breaking representation of the DP (Sethuraman, 1994), MacEachern (1999, 2000) proposed a class of fixed weight dependent DP (DDP) priors. Applications of this DDP framework have been employed in ANOVA modeling (De Iorio et al., 2004), spatial data analysis (Gelfand et al., 2005), time series (Caron et al., 2006) and functional data analysis applications (Petrone et al., 2009; Scarpa and Dunson, 2009) among many others, with some theoretical properties highlighted in Barrientos et al. (2012).

However such fixed weight DDPs lack flexibility in feature-dependent clustering, as noted in MacEachern (2000). This has motivated alternative formulations which allow the mixing weights to change with the features, with some examples including the order-based dependent Dirichlet process (Griffin and Steel, 2006), kernel- (Dunson and Park, 2008), and probit- (Rodriguez and Dunson, 2011) stick breaking processes.

Alternative approaches build on random partition models (RPMs), working directly with the probability distribution $p(\boldsymbol{c})$ on the partition $\boldsymbol{c}$ of indices $\{1,\ldots,N\}$ into clusters. Particular attention has been given to the class of product partition models (PPMs) (Barry and Hartigan, 1992; Hartigan, 1990) in which $p(\boldsymbol{c})$ can be factorized into a product of cluster-dependent functions, known as cohesion functions. A common strategy modifies the cohesion to allow feature-dependence; refer, for examples, to Park and Dunson (2010), Müller et al. (2011), Blei and Frazier (2011) and Dahl et al. (2017).

Our focus is fundamentally different. In particular, we do not have features $\mathbf{x}_{i}$ on indices $i$ but have access to an informed prior guess $\boldsymbol{c}_{0}$ for the partition $\boldsymbol{c}$; other than this information it is plausible to rely on exchangeable priors. To address this problem, we propose a general strategy to modify a baseline EPPF to include centering on $\boldsymbol{c}_{0}$. In particular, our proposed Centered Partition (CP) process defines the partition prior as proportional to an EPPF multiplied by an exponential factor which depends on a distance function $d(\boldsymbol{c},\boldsymbol{c}_{0})$ measuring how far $\boldsymbol{c}$ is from $\boldsymbol{c}_{0}$. The proposed framework should be broadly useful in including extra information into EPPFs, which tend to face issues in lacking incorporation of real prior information from applications.

The paper is organized as follows. Section 2 introduces concepts and notation related to Bayesian nonparametric clustering. In Section 3 we illustrate the general CP process formulation and describe an approach to posterior computation relying on Markov chain Monte Carlo (MCMC). Section 4 proposes a general strategy for prior calibration building on a targeted Monte Carlo procedure. Simulation studies and application to the motivating birth defects epidemiology study are provided in Section 5, with technical details included in an Appendix.

2 Clustering and Bayesian models

This section introduces some concepts related to the representation of the clustering space from a combinatorial perspective, which will be useful to define the Centered Partition process, along with an introduction to Bayesian nonparametric clustering models.

2.1 Set Partitions

Let $\boldsymbol{c}$ be a generic clustering of indices $[N]=\{1,\ldots,N\}$. It can be either represented as a vector of indices $\{c_{1},\ldots,c_{N}\}$, with $c_{i}\in\{1,\ldots,K\}$ for $i=1,\ldots,N$ and $c_{i}=c_{j}$ when $i$ and $j$ belong to the same cluster, or as a collection of disjoint subsets (blocks) $\{B_{1},B_{2},\ldots,B_{K}\}$ where $B_{k}$ contains all the indices of data points in the $k$-th cluster and $K$ is the number of clusters in the sample of size $N$. From a mathematical perspective $\boldsymbol{c}=\{B_{1},\ldots,B_{K}\}$ is a combinatorial object known as set partition of $[N]$. The collection of all possible set partitions of $[N]$, denoted with $\Pi_{N}$, is known as the partition lattice. We refer to Stanley (1997) and Davey and Priestley (2002) for an introduction to lattice theory, and to Meilă (2007) and Wade and Ghahramani (2018) for a review of the concepts from a more statistical perspective.

According to Knuth in Wilson and Watkins (2013), set partitions seem to have been studied first in Japan around A.D. 1500, due to a popular game in the upper class society known as genji-ko; five unknown incense sticks were burned and players were asked to identify which of the scents were the same, and which were different. Soon diagrams were developed to model all the $52$ outcomes, which corresponds to all the possible set partitions of $N=5$ elements. First results focused on enumerating the elements of the space. For example, for a fixed number of blocks $K$, the number of ways to assign $N$ elements to $K$ groups is described by the Stirling number of the second kind

[TABLE]

while the Bell number $\mathcal{B}_{N}=\sum_{K=1}^{N}\mathcal{S}_{N,K}$ describes the number of all possible set partitions of $N$ elements.

Interest progressively shifted towards characterizing the structure of the space of partitions using the notion of partial order. Consider $\Pi_{N}$ endowed with the set containment relation $\leq$, meaning that for $\boldsymbol{c}=\{B_{1},\ldots,B_{K}\},\boldsymbol{c}^{\prime}=\{B^{\prime}_{1},\ldots,B^{\prime}_{K^{\prime}}\}$ belonging to $\Pi_{N}$, $\boldsymbol{c}\leq\boldsymbol{c}^{\prime}$ if for all $i=1,\ldots,K,B_{i}\subseteq B^{\prime}_{j}$ for some $j\in\{1,\ldots,K^{\prime}\}$. Then the space $(\Pi_{N},\leq)$ is a partially ordered set (poset), which satisfies the following properties:

Reflexivity: for every $\boldsymbol{c},\boldsymbol{c}^{\prime}\in\Pi_{N}$, $\boldsymbol{c}\leq\boldsymbol{c}$, 2. 2.

Antisymmetry: if $\boldsymbol{c}\leq\boldsymbol{c}^{\prime}$ and $\boldsymbol{c}^{\prime}\leq\boldsymbol{c}$ then $\boldsymbol{c}=\boldsymbol{c}^{\prime}$, 3. 3.

Transitivity: if $\boldsymbol{c}\leq\boldsymbol{c}^{\prime}$ and $\boldsymbol{c}^{\prime}\leq\boldsymbol{c}^{\prime\prime}$, then $\boldsymbol{c}\leq\boldsymbol{c}^{\prime\prime}$.

Moreover, for any $\boldsymbol{c},\boldsymbol{c}^{\prime}\in\Pi_{N}$, it is said that $\boldsymbol{c}$ is covered (or refined) by $\boldsymbol{c}^{\prime}$ if $\boldsymbol{c}\leq\boldsymbol{c}^{\prime}$ and there is no $\boldsymbol{c}^{\prime\prime}$ such that $\boldsymbol{c}<\boldsymbol{c}^{\prime\prime}<\boldsymbol{c}^{\prime}$. Such a relation is indicated by $\boldsymbol{c}\prec\boldsymbol{c}^{\prime}$. This covering relation allows one to represent the space of partitions using a Hasse diagram, in which the elements of $\Pi_{N}$ correspond to nodes in a graph and a line is drawn from $\boldsymbol{c}$ to $\boldsymbol{c}^{\prime}$ when $\boldsymbol{c}\prec\boldsymbol{c}^{\prime}$; there is a connection from a partition $\boldsymbol{c}$ to another one when the second can be obtained by splitting or merging one of the blocks in $\boldsymbol{c}$. See Fig. 1 for an example of the Hasse diagram of $\Pi_{4}$. Conventionally, the partition with just one cluster is represented at the top of the diagram and denoted as $\boldsymbol{1}$, while the partition having every observation in its own cluster is at the bottom and indicated with $\boldsymbol{0}$.

This representation of the set partitions space $\Pi_{N}$ as a partially ordered set provides a useful framework to characterize its elements. As already mentioned, two partitions connected in the Hasse diagram can be obtained from one another by means of a single operation of split or merge; a sequence of connections is a path, linking the two extreme partitions $\boldsymbol{0}$ and $\boldsymbol{1}$. A path starting from $\boldsymbol{0}$ connects partitions with an increasing number of blocks, from $1$ to $N$, which is referred as the rank of the partition. Set partitions with the same rank may differ in terms of their configuration $\boldsymbol{\Lambda}(\boldsymbol{c})$, the sequence of block cardinalities $\{|B_{1}|,\ldots,|B_{K}|\}$, which corresponds to another combinatorial object known as an integer partition of $N$. In combinatorics, an integer partition is defined as the multiset of positive integers $\{\lambda_{1}\ldots\lambda_{K}\}$, listed in decreasing order by convention, such that $\sum_{i=1}^{K}\lambda_{i}=N$. Also the associated space of all possible integer partitions $I_{N}$ is a partially ordered set, making the definition of configuration a poset mapping $\boldsymbol{\Lambda}(\cdot):\boldsymbol{c}\in\Pi_{N}\rightarrow\boldsymbol{\lambda}\in I_{N}$.

Finally, the space $\Pi_{N}$ is a lattice, based on the fact that every pair of elements has a greatest lower bound (g.l.b.) and a least upper bound (l.u.b.) indicated with the “meet” $\land$ and the “join” $\lor$ operators, i.e. $\boldsymbol{c}\land\boldsymbol{c}^{\prime}=\text{g.l.b.}(\boldsymbol{c},\boldsymbol{c}^{\prime})$ and $\boldsymbol{c}\lor\boldsymbol{c}^{\prime}=\text{l.u.b.}(\boldsymbol{c},\boldsymbol{c}^{\prime})$ and equality holds under a permutation of the cluster labels. An element $\boldsymbol{c}\in\Pi_{N}$ is an upper bound for a subset $\boldsymbol{S}\subseteq\Pi_{N}$ if $\boldsymbol{s}\leq\boldsymbol{c}$ for all $\boldsymbol{s}\in\boldsymbol{S}$, and it is the least upper bound for a subset $\boldsymbol{S}\subseteq\Pi_{N}$ if $\boldsymbol{c}$ is an upper bound for $\boldsymbol{S}$ and $\boldsymbol{c}\leq\boldsymbol{c}^{\prime}$ for all upper bounds $\boldsymbol{c}^{\prime}$ of $\boldsymbol{S}$. The lower bound and the greatest lower bound are defined similarly, and the definition applies also to the elements of the space $I_{N}$. Consider, as an example, $\boldsymbol{c}=\{1\}\{2,3,4\}$ and $\boldsymbol{c}^{\prime}=\{3\}\{1,2,4\}$; their greatest lower bound is $\boldsymbol{c}\land\boldsymbol{c}^{\prime}=\{1\}\{3\}\{2,4\}$ while the lowest upper bound is $\boldsymbol{c}\lor\boldsymbol{c}^{\prime}=\{1,2,3,4\}$. Considering the Hasse diagram in Fig 1 the g.l.b. and l.u.b. are the two partitions which reach both $\boldsymbol{c}$ and $\boldsymbol{c}^{\prime}$ through the shortest path, respectively from below and from above.

2.2 Bayesian mixture models

From a statistical perspective, set partitions are key elements in a Bayesian mixture model framework. The main underlying assumption is that observations $y_{1},\ldots,y_{N}$ are independent conditional on the partition $\boldsymbol{c}$, and their joint probability density can be expressed as

[TABLE]

with $\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{K})$ a vector of unknown parameters indexing the distribution of observations $\mathbf{y}_{k}=\{y_{i}\}_{i\in B_{k}}$ for each cluster $k=1,\ldots,K$. In a Bayesian formulation, a prior distribution is assigned to each possible partition $\boldsymbol{c}$, leading to a posterior of the form

[TABLE]

Hence the set partition $\boldsymbol{c}$ is conceived as a random object and elicitation of its prior distribution is a critical issue in Bayesian modeling.

The first distribution one may use, in the absence of prior information, is the uniform distribution, which gives the same probability to every partition with $p(\boldsymbol{c})=1/\mathcal{B}_{N}$; however, even for small values of $N$ the Bell number $\mathcal{B}_{N}$ is very large, making computation of the posterior intractable even for simple choices of the likelihood. This motivated the definition of alternative prior distributions based on different concepts of uniformity, with the Jensen and Liu (2008) prior favoring uniform placement of new observations in one of the existing clusters, and Casella et al. (2014) proposing a hierarchical uniform prior, which gives equal probability to set partitions having the same configuration.

Usual Bayesian nonparametric procedures build instead on discrete nonparametric priors, i.e. priors that have discrete realization almost surely. Dirichlet and Pitman-Yor processes are well known to have this property, as does the broader class of Gibbs-type priors. Any discrete random probability measure $\tilde{p}$ can induce an exchangeable random partition. Due to the discreteness of the process, $\tilde{p}$ induces a partition of the observations $y_{1},\ldots,y_{N}$ which can be characterized via an Exchangeable Probability Partition Function. For both Dirichlet and Pitman-Yor processes, the EPPF is available in closed form as reported in Table 1 along with the case of the finite mixture model with $\kappa$ components and a symmetric Dirichlet prior with parameters $(\gamma/\kappa,\ldots,\gamma/\kappa)$. Notice that $\lambda_{j}=|B_{j}|$ is the cardinality of the clusters composing the partition, while notation $(x)_{r}$ is for the rising factorial $x(x+1)\cdots(x+r-1)$.

There is a strong connection with the exchangeable random partitions induced by Gibbs-type priors and product partition models. A product partition model assumes that the prior probability for the partition $\boldsymbol{c}$ has the following form

[TABLE]

with $\rho(\cdot)$ known as the cohesion function. The underlying assumption is that the prior distribution for the set partition $\boldsymbol{c}$ can be factorized as the product of functions that depend only on the blocks composing it. Such a definition, in conjunction with formulation (1) for the data likelihood, guarantees the property that the posterior distribution for $\boldsymbol{c}$ is still in the class of product partition models.

Distributions in Table 1 are all characterized by a cohesion function that depends on the blocks through their cardinality. Although the parameters can control the expected number of clusters, this assumption is too strict in many applied contexts in which prior information is available about the grouping. In particular, the same probability is given to partitions with the same configuration but having a totally different composition.

3 Centered Partition Processes

Our focus is on incorporating structured knowledge about clustering of the finite set of indices $[N]=\{1,\ldots,N\}$ in the prior distribution within a Bayesian mixture model framework. We consider as a first source of information a given potential clustering, but our approach can also accommodate prior information on summary statistics such as the number of clusters and cluster sizes.

3.1 General formulation

Assume that a potential clustering $\boldsymbol{c}_{0}$ is given and we wish to include this information in the prior distribution. To address this problem, we propose a general strategy to modify a baseline EPPF to shrink towards $\boldsymbol{c}_{0}$. In particular, our proposed CP process defines the prior on set partitions as proportional to a baseline EPPF multiplied by a penalization term of the type

[TABLE]

with $\psi>0$ a penalization parameter, $d(\boldsymbol{c},\boldsymbol{c}_{0})$ a suitable distance measuring how far $\boldsymbol{c}$ is from $\boldsymbol{c}_{0}$ and $p_{0}(\boldsymbol{c})$ indicates a baseline EPPF, that may depend on some parameters that are not of interest at the moment. For $\psi\rightarrow 0$, $p(\boldsymbol{c}|\boldsymbol{c}_{0},\psi)$ corresponds to the baseline EPPF $p(\boldsymbol{c}_{0})$, while for $\psi\rightarrow\infty$, $p(\boldsymbol{c}=\boldsymbol{c}_{0})\rightarrow 1$.

Note that $d(\boldsymbol{c},\boldsymbol{c}_{0})$ takes a finite number of discrete values $\Delta=\{\delta_{0},\ldots,\delta_{L}\}$, with $L$ depending on $\boldsymbol{c}_{0}$ and on the distance $d(\cdot,\cdot)$. We can define sets of partitions having the same fixed distance from $\boldsymbol{c}_{0}$ as

[TABLE]

Hence, for $\delta_{0}=0$, $s_{0}(\boldsymbol{c}_{0})$ denotes the set of partitions equal to the base one, meaning that they differ from $\boldsymbol{c}_{0}$ only by a permutation of the cluster labels. Then $s_{1}(\boldsymbol{c}_{0})$ denotes the set of partitions with minimum distance $\delta_{1}$ from $\boldsymbol{c}_{0}$, $s_{2}(\boldsymbol{c}_{0})$ the set of partitions with the second minimum distance $\delta_{2}$ from $\boldsymbol{c}_{0}$ and so on. The introduced exponential term penalizes equally partitions in the same set $s_{l}(\boldsymbol{c}_{0})$ for a given $\delta_{l}$, but the resulting probabilities may differ depending on the chosen baseline EPPF.

3.2 Choices of distance function

The proposed CP process modifies a baseline EPPF to include a distance-based penalization term, which aims to shrink the prior distribution towards a prior partition guess. The choice of distance plays a key role in determining the behavior of the prior distribution. A variety of different distances and indices have been employed in clustering procedures and comparisons. We consider in this paper the Variation of Information (VI), obtained axiomatically in Meilă (2007) using information theory, and shown to nicely characterize neighborhoods of a given partition by Wade and Ghahramani (2018). The Variation of Information is based on the Shannon entropy $H(\cdot)$, and can be computed as

[TABLE]

where $\log$ denotes $\log$ base 2, and $\lambda^{\land}_{jl}$ the size of blocks of the intersection $\boldsymbol{c}\land\boldsymbol{c}^{\prime}$ and hence the number of indices in block $j$ under partition $\boldsymbol{c}$ and block $l$ under $\boldsymbol{c}^{\prime}$. Notice that VI ranges from [math] to $\log_{2}(N)$. Although normalized versions have been proposed (Vinh et al., 2010), some desirable properties are lost under normalization. We refer to Meilă (2007) and Wade and Ghahramani (2018) for additional properties and empirical evaluations.

An alternative definition of the VI can be derived from lattice theory, exploiting the concepts provided in Section 2.1. We refer to Monjardet (1981) for general theory about metrics on lattices and ordered sets, and Rossi (2015) for a more recent review focused on set partitions. In general, a distance between two different partitions $\boldsymbol{c},\boldsymbol{c}^{\prime}\in\Pi_{N}$ can be defined by means of the Hasse diagram via the minimum weighted path, which corresponds to the shortest path length when edges are equally weighted. Instead, when edges depend on the entropy function through $w(\boldsymbol{c},\boldsymbol{c}^{\prime})=|H(\boldsymbol{c})-H(\boldsymbol{c}^{\prime})|$, the minimum weighted path between two partitions is the Variation of Information. Notice that two partitions are connected when in a covering relation, then $\boldsymbol{c}\land\boldsymbol{c}^{\prime}$ is either equal to $\boldsymbol{c}$ or $\boldsymbol{c}^{\prime}$ and $VI(\boldsymbol{c},\boldsymbol{c}^{\prime})=w(\boldsymbol{c},\boldsymbol{c}^{\prime})$. The minimum weight $w(\boldsymbol{c},\boldsymbol{c}^{\prime})$ corresponds to $2/N$ which is attained when two singleton clusters are merged, or conversely, a cluster consisting of two points is split (see Meilă, 2007).

3.3 Effect of the prior penalization

We first consider the important special case in which the baseline EPPF is $p_{0}(\boldsymbol{c})=1/\mathcal{B}_{N}$ and the CP process reduces to $p(\boldsymbol{c}|\boldsymbol{c}_{0},\psi)\propto\exp\{-\psi d(\boldsymbol{c},\boldsymbol{c}_{0})\}$ with equation (4) simplifying to

[TABLE]

where $|\cdot|$ indicates the cardinality and $s_{l}(\boldsymbol{c}_{0})$ is defined in (5).

Considering $N=5$, there are $52$ possible set partitions; Figure 2 shows the prior probabilities assigned to partitions under the CP process for different values of $\psi\in(0,3)$ with $\psi=0$ corresponding to the uniform prior. Notice that base partitions with the same configuration (e.g. for $\boldsymbol{c}_{0}=\{1,2\}\{3,4,5\}$ all the partitions with blocks sizes $\{3,2\}$), will behave in the same way, with the same probabilities assigned to partitions different in composition.

Non-zero values of $\psi$ increase the prior probability of partitions $\boldsymbol{c}$ that are relatively close to the chosen $\boldsymbol{c}_{0}$. However, the effect is not uniform but depends on the structure of both $\boldsymbol{c}$ and $\boldsymbol{c}_{0}$. For example, consider the inflation that occurs in the blue region as $\psi$ increases from [math] to $3$. When $\boldsymbol{c}_{0}$ has $2$ blocks (Figure 2a) versus $4$ (Figure 2d) there is a bigger increase. This is because the space of set partitions $\Pi_{N}$ is not “uniform”, since given a fixed configuration there is a heterogeneous number of partitions. Rewriting $\boldsymbol{\lambda}=(\lambda_{1},\ldots,\lambda_{K})$ as $\boldsymbol{\lambda}=(1^{f_{1}},2^{f_{2}},\ldots,K^{f_{K}})$, with the notation indicating that there are $f_{i}$ elements of $\boldsymbol{\lambda}$ equal to $i$, the number of set partitions with configuration $\boldsymbol{\lambda}$ is

[TABLE]

For example, for $\{221\}=1^{1}2^{2}3^{0}4^{0}5^{0}$, the number of corresponding set partitions is $15$, while there are $10$ set partitions of type $\{311\}$.

While the uniform distribution gives the same probability to each partition in the space, the EPPF induced by Gibbs-type priors distinguishes between different configurations, but not among partitions with the same configuration. We focus on the Dirichlet process case, being the most popular process employed in applications. Under the DP the induced EPPF $p_{0}(\boldsymbol{c})\propto\alpha^{K}\prod_{j=1}^{K}\Gamma(\lambda_{j})$ is a function of the configuration $\boldsymbol{\Lambda}(\boldsymbol{c})$, which is one of $\{\boldsymbol{\lambda}_{1},\ldots,\boldsymbol{\lambda}_{M}\}$ since the possible configurations are finite and correspond to the number of integer partitions. Letting $g(\boldsymbol{\Lambda}(\boldsymbol{c}))=\alpha^{K}\prod_{j=1}^{K}\Gamma(\lambda_{j})$, the formulation in (4) can be written as

[TABLE]

where $s_{lm}(\boldsymbol{c}_{0})=\{\boldsymbol{c}\in\Pi_{N}:d(\boldsymbol{c},\boldsymbol{c}_{0})=\delta_{l},\boldsymbol{\Lambda}(\boldsymbol{c})=\boldsymbol{\lambda}_{m}\}$, the set of partitions with distance $\delta_{l}$ from $\boldsymbol{c}_{0}$ and configuration $\boldsymbol{\lambda}_{m}$ for $l=0,1,\ldots,L$ and $m=1,\ldots,M$, with $n_{lm}$ indicating the cardinality. The factorization (7) applies for the family of Gibbs-type priors in general, with different expressions of $g(\boldsymbol{\Lambda}(\boldsymbol{c}))$.

In Figure 3 we consider the prior distribution induced by the CP process when the baseline EPPF $p_{0}(\boldsymbol{c})$ comes from a Dirichlet process with concentration parameter $\alpha=1$, considering the same base partitions and values for $\psi$ as in Figure 2. For the same values of the parameter $\psi$, the behavior of the CP process changes significantly due to the effect of the base prior. In particular, in the top left panel the CP process is centered on $\boldsymbol{c}_{0}=\{1,2,3,4,5\}$, the partition with only one cluster, which is a priori the most likely one for $\psi=0$. In general, for small values of $\psi$ the clustering process will most closely resemble that for a DP, and as $\psi$ increases the DP prior probabilities are decreased for partitions relatively far from $\boldsymbol{c}_{0}$ and increased for $\boldsymbol{c}_{0}$ relatively close.

3.4 Posterior computation under Gibbs-type priors

Certain MCMC algorithms for Bayesian nonparametric mixture models can be easily modified for posterior computation in CP process models. In particular, we adapt the so-called “marginal algorithms” developed for Dirichlet and Pitman-Yor processes. These methods are called marginal since the mixing measure $P$ is integrated out of the model and the predictive distribution is used within a MCMC sampler. In the following, we recall Algorithm 2 in Neal (2000) and illustrate how it can be adapted to sample from the CP process posterior. We refer to Neal (2000) and references therein for an overview and discussion of methods for both conjugate and nonconjugate cases, and to Fall and Barat (2014) for adaptation to Pitman-Yor processes.

Let $\boldsymbol{c}$ be represented as an $N$-dimensional vector of indices $\{c_{1},\ldots,c_{N}\}$ encoding cluster allocation and let $\theta_{k}$ be the set of parameters currently associated to cluster $k$. The prior predictive distribution for a single $c_{i}$ conditionally on $\boldsymbol{c}^{-i}=\{c_{1},\ldots,c_{i-1},c_{i+1},\ldots,c_{N}\}$ is exploited to perform the Gibbs sampling step allocating observations to either a new cluster or one of the existing ones. Algorithm $2$ in Neal (2000) updates each $c_{i}$ sequentially for $i=1,\ldots,N$ via a reseating procedure, according to the conditional posterior distribution

[TABLE]

with $K^{-}$ the number of clusters after removing observation $i$. The conditional distribution $p(c_{i}=k|\boldsymbol{c}^{-i})$ is reported in Table 2 for different choices of the prior EPPF. Notice that, for the case of finite Dirichlet prior, the update consists only in the first line of equation (8), since the number of classes is fixed. For Dirichlet and Pitman-Yor processes, when observation $i$ is associated to a new cluster, a new value for $\theta$ is sampled from its posterior distribution based on the base measure $G_{0}$ and the observation $y_{i}$. This approach is straightforward when we can compute the integral $\int p(y_{i}|\theta)dG_{0}(\theta)$, as will generally be the case when $G_{0}$ is a conjugate prior.

Considering the proposed CP process, the conditional distribution for $c_{i}$ given $\boldsymbol{c}^{-i}$ can still be computed, but it depends both on the base prior and the penalization term accounting for the distance between the base partition $\boldsymbol{c}_{0}$ and the one obtained by assigning the observation $i$ to either one of the existing classes $k\in\{1,\ldots,K^{-}\}$ or a new one. Hence, the step in equation (8) can be easily adapted by substituting the conditional distribution for $p(c_{i}=k|\boldsymbol{c}^{-i})$ with

[TABLE]

with $\boldsymbol{c}=\{\boldsymbol{c}^{-i}\cup\{c_{i}=k\}\}$ the current state of the clustering and $p_{0}(c_{i}=k|\boldsymbol{c}^{-i})$ one of the conditional distributions in Table 2. Additional steps on the implementation using the variation of information as a distance are given in the Appendix (Algorithm Marginal sampling using variation of information).

Extension to the non-conjugate context can be similarly handled exploiting Algorithm $8$ in Neal (2000) based on auxiliary parameters, which avoids the computation of the integral $\int p(y_{i}|\theta)dG_{0}(\theta)$. The only difference is that, when $c_{i}$ is updated, $m$ temporary auxiliary variables are introduced to represent possible values of components parameters that are not associated with any other observations. Such variables are simply sampled from the base measure $G_{0}$, with the probabilities of a new cluster in Table 2 changing into $(\alpha/m)/(\alpha+N-1)$ for the Dirichlet process and to $[(\alpha+\sigma K^{-})/m]/(\alpha+N-1)$ for the Pitman-Yor process, for $k=K^{-}+1,\ldots,K^{-}+m$.

4 Prior calibration

As the number of observations $N$ increases, the number of partitions explodes, and higher values of $\psi$ are needed to place non-negligible prior probability in small to moderate neighborhoods around $\boldsymbol{c}_{0}$. The prior concentration around $\boldsymbol{c}_{0}$ depends on three main factors: i) $N$ through $\mathcal{B}_{N}$, i.e. the cardinality of the space of set partitions, ii) the baseline EPPF $p_{0}(\boldsymbol{c}_{0})$ and iii) where $\boldsymbol{c}_{0}$ is located in the space. We hence propose a general method to evaluate the prior behavior under different settings, while suggesting how to choose the parameter $\psi$.

One may evaluate the prior distribution for different values of $\psi$ and check its behavior using graphs such as those in Section 3.3, however they become difficult to interpret as the space of partitions grows. We propose to evaluate the probability distribution of the distances $\delta=d(\boldsymbol{c},\boldsymbol{c}_{0})$ from the known partition $\boldsymbol{c}_{0}$. The probability assigned to different distances by the prior is

[TABLE]

with $\mathcal{I}(\cdot)$ the indicator function and $s_{l}(\boldsymbol{c}_{0})$ denoting the set of partitions distance $\delta_{l}$ from $\boldsymbol{c}_{0}$, as defined in (5). Consider the uniform distribution on set partitions, $p(\delta=\delta_{l})=|s_{l}(\boldsymbol{c}_{0})|/\mathcal{B}_{N}$, the proportion of partitions distance $\delta_{l}$ from $\boldsymbol{c}_{0}$. Under the general definition of the CP process, the resulting distribution becomes

[TABLE]

with the case of Gibbs-type EPPF corresponding to

[TABLE]

Notice that the uniform EPPF case is recovered when $g(\boldsymbol{\lambda}_{m})=1$ for $m=0,\ldots,M$, so that $\sum_{m=1}^{M}n_{lm}=n_{l}$. Hence the probability in (9) simplifies to

[TABLE]

In general, since distances are naturally ordered, the corresponding cumulative distribution function can be simply defined as $F(\delta)=\sum_{\delta_{l}\leq\delta}p(\delta_{l})$ for $\delta\in\{\delta_{0},\ldots,\delta_{L}\}$ and used to assess how much mass is placed in different size neighborhoods around $\boldsymbol{c}_{0}$ under different values of $\psi$. Hence we can choose $\psi$ to place a specified probability $q$ (e.g. $q=0.9$) on partitions within a specified distance $\delta^{*}$ from $\boldsymbol{c}_{0}$. This would correspond to calibrating $\psi$ so that $F(\delta^{*})\approx q$, with $F(\delta^{*})\geq q$. In other words, partitions generated from the prior would have at least probability $q$ of being within distance $\delta^{*}$ from $\boldsymbol{c}_{0}$.

The main problem is in computing the probabilities in equations (10)-(11), which depend on all the set partitions in the space. In fact, one needs to count all the partitions having distance $\delta_{l}$ for $l=0,\ldots,L$ when the base EPPF is uniform, while taking account of configurations in the case of the Gibbs-type priors. Even if there are quite efficient algorithms to list all the possible set partitions of $N$ (see Knuth, 2005; Nijenhuis and Wilf, 2014), it becomes computationally infeasible due to the extremely rapid growth of the space; for example from $N=12$ to $13$, the number of set partitions grows from $\mathcal{B}_{12}=4,213,597$ to $\mathcal{B}_{13}=27,644,437$.

We propose a general strategy to approximate prior probabilities assigned to different distances from $\boldsymbol{c}_{0}$ focused on obtaining estimates of distance values and related counts, which represent the sufficient quantities to compute (10)-(11) under different values of $\psi$. We consider a targeted Monte Carlo procedure which augments uniform sampling on the space of set partitions with a deterministic local search using the Hasse diagram to estimate the counts for small values of the distance.

4.1 Deterministic local search

Poset theory provides a nice representation of the space of set partitions by means of the Hasse diagram illustrated in Section 2.1, along with suitable definition of metrics. A known partition $\boldsymbol{c}_{0}$ can be characterized in terms of number of blocks $K_{0}$ and configuration $\boldsymbol{\Lambda}(\boldsymbol{c}_{0})$. These elements allows one to locate $\boldsymbol{c}_{0}$ in the Hasse diagram and then explore connected partitions by means of split and merge operations on the clusters in $\boldsymbol{c}_{0}$.

As an illustrative example, consider the Hasse diagram of $\Pi_{4}$ in Figure 4 and $\boldsymbol{c}_{0}=\{1\}\{2,3,4\}$, having $2$ clusters and configuration $\boldsymbol{\Lambda}({\boldsymbol{c}_{0}})=\{31\}$. Let $\mathcal{N}_{1}(\boldsymbol{c}_{0})$ denote the sets of partitions directly connected with $\boldsymbol{c}_{0}$, i.e. partitions covering $\boldsymbol{c}_{0}$ and those covered by $\boldsymbol{c}_{0}$. In general, a partition $\boldsymbol{c}_{0}$ with $K_{0}$ clusters is covered by ${K_{0}}\choose{2}$ partitions and covers $\sum_{j=1}^{K_{0}}2^{\lambda_{j}-1}-1$. In the example, $\mathcal{N}_{1}(\boldsymbol{c}_{0})$ contains $\{1,2,3,4\}$ obtained from $\boldsymbol{c}_{0}$ with a merge operation on the two clusters, and all the partitions obtained by splitting the cluster $\{2,3,4\}$ in any possible way. The base idea underlying the proposed local search, consists in exploiting the Hasse diagram representation to find all the partitions in increasing distance neighborhoods of $\boldsymbol{c}_{0}$. One can list partitions at $T$ connections from $\boldsymbol{c}_{0}$ starting from $\mathcal{N}_{1}(\boldsymbol{c}_{0})$ by recursively applying split and merge operations on the set of partitions explored at each step. Potentially, with enough operations one can reach all the set partitions, since the space is finite with lower and upper bounds.

In practice, the space is too huge to be explored entirely, and a truncation is needed. From the example in Figure 4, $\mathcal{N}_{1}(\boldsymbol{c}_{0})$ contains $3$ partitions with distance $0.69$ from $\boldsymbol{c}_{0}$ and one with distance $1.19$. Although $\mathcal{N}_{2}(\boldsymbol{c}_{0})$ may contain partitions closer to $\boldsymbol{c}_{0}$ than this last, the definition of distance in Section 3.2 guarantees that there are no other partitions with distance from $\boldsymbol{c}_{0}$ less than $0.69$. Since the VI is the minimum weighted path between two partitions, all the partitions reached at the second exploration step add a nonzero weight to distance computation. This consideration extends to an arbitrary number of explorations $T$, with $\delta_{L^{*}}=\text{min}\{d(\boldsymbol{c}^{*},\boldsymbol{c}_{0})\}_{\boldsymbol{c}^{*}\in\mathcal{N}_{T}(\boldsymbol{c}_{0})}$ being the upper bound on the distance value. By discarding all partitions with distance greater that $\delta_{L^{*}}$, one can compute exactly the counts in equations (10)-(11) related to distances $\delta_{0},\ldots,\delta_{L^{*}}$. Notice that $2/N$ is the minimum distance between two different partitions, and $2T/N$ is a general lower bound on the distances from $\boldsymbol{c}_{0}$ that can be reached in $T$ iterations.

4.2 Monte Carlo approximation

We pair the local exploration with a Monte Carlo procedure to estimate the counts and distances greater that $\delta_{L^{*}}$, in order to obtain a more refined representation of the prior distance probabilities. Sampling uniformly from the space of partitions is not in general a trivial problem, but a nice strategy has been proposed in Stam (1983), in which the probability of a partition with $K$ clusters is used to sample partitions via an urn model. Derivation of the algorithm starts from the Dobiński formula (Dobiński, 1877) for the Bell numbers

[TABLE]

which from a probabilistic perspective corresponds to the $k$-th moment of the Poisson distribution with expected value equal to $1$. Then a probability distribution for the number of clusters $K\in\{1,\ldots,N\}$ of a set partition can be defined as

[TABLE]

which is a well defined law thanks to (12). To simulate a uniform law over $\Pi_{N}$, Stam (1983)’s algorithm first generates the number of clusters $K$ according to (13) and, conditionally on the sampled value, it allocates observations to the clusters according a discrete uniform distribution over $\{1,\ldots,K\}$. We refer to Stam (1983) and Pitman (1997) for derivations and proof of the validity of the algorithm.

We adapt the uniform sampling to account for the values already computed by rejecting all the partitions with distance less that $\delta_{L^{*}}$, restricting the space to $\left\{\Pi_{N}\setminus\{\mathcal{N}_{t}(\boldsymbol{c}_{0})\}_{t=0}^{T}\right\}$. In practice, few samples are discarded since the probability to sample one such partition corresponds to $|\{\mathcal{N}_{t}(\boldsymbol{c}_{0})\}_{t=0}^{T}|/\mathcal{B}_{N}$, which is negligible for small values of exploration steps $T$ that are generally used in the local search. A sample of partitions $\boldsymbol{c}^{(1)},\ldots,\boldsymbol{c}^{(R)}$, can be used to provide an estimate of the counts. Let $R^{*}$ denote the number of accepted partitions and $\mathcal{B}^{*}=\mathcal{B}_{N}-|\{\mathcal{N}_{t}(\boldsymbol{c}_{0})\}_{t=0}^{T}|$ be the number of partitions in the restricted space. Conditionally on the observed values of distances in the sample, $\hat{\delta}_{(L^{*}+1)},\ldots,\hat{\delta}_{L}$, an estimate of the number of partitions with distance $\hat{\delta}_{l}$ to use in the uniform EPPF case is

[TABLE]

obtained by multiplying the proportions of partitions in the sample by the total known number of partitions. For the Gibbs-type EPPF case one needs also to account for the configurations $\boldsymbol{\lambda}_{1},\ldots,\boldsymbol{\lambda}_{M}$ in a given orbital of the distance; hence, the estimates are

[TABLE]

Pairing these estimates with the counts obtained via the local search, one can evaluate the distributions in equations (10)-(11) for different values of $\psi$. The entire procedure is summarized in Algorithm Prior calibration in the Appendix. Although it requires a considerable number of steps, the procedure can be performed one single time providing information for different choices of $\psi$ and EPPFs. Moreover the local search can be implemented in parallel to reduce computational costs.

We consider an example for $N=12$ and $\boldsymbol{c}_{0}$ with configuration $\{3,3,3,3\}$. Figure 5 shows the resulting cumulative probability estimates of the CP process under uniform and DP($\alpha=1$) base distributions, estimated with $m=4$ iterations of the local search and $20,000$ samples. Dots represent values of the cumulative probabilities, with different colors in correspondence to different values of the parameter $\psi$. Using these estimates one can assess how much probability is placed in different distance neighborhoods of $\boldsymbol{c}_{0}$; tables in Figure 5 show the distance values defining neighborhoods around $\boldsymbol{c}_{0}$ with $90$% prior probability. If one wishes to place such probability mass on partitions within distance $1$ from $\boldsymbol{c}_{0}$, a value of $\psi$ around $10$ and $15$ is needed, respectively, under uniform and DP base prior. We suggest, when performing the analysis, to consider also neighborhood values of the chosen $\psi$, in order to assess the sensitivity of the results.

5 The National Birth Defects Prevention Study

The National Birth Defects Prevention Study (NBDPS) is a multi-state population-based, case-control study of birth defects in the United States (Yoon et al., 2001). Infants were identified using birth defects surveillance systems in recruitment areas within ten US states (Arkansas, California, Georgia, Iowa, Massachusetts, New Jersey, New York, North Carolina, Texas, and Utah), which cover roughly 10% of US births. Diagnostic case information was obtained from medical records and verified by a standardized clinician review specific to the study (Rasmussen et al., 2003). Participants in the study included mothers with expected dates of delivery from 1997-2009. Controls were identified from birth certificates or hospital records and were live-born infants without any known birth defects. Each state site attempted to recruit $300$ cases and $100$ (unmatched) controls annually. A telephone interview was conducted with case and control mothers to solicit a wide range of demographic, lifestyle, medical, nutrition, occupational and environmental exposure history information.

Because birth defects are highly heterogeneous, a relatively large number of defects of unknown etiology are included in the NBDPS. We are particularly interested in congenital heart defects (CHD), the most common type of birth defect and the leading cause of infant death due to birth defects. Because some of these defects are relatively rare, in many cases we lack precision for investigating associations between potential risk factors and individual birth defects. For this reason, researchers typically lump embryologically distinct and potentially etiologically heterogeneous defects in order to increase power (e.g., grouping all heart defects together), even knowing the underlying mechanisms may differ substantially. In fact, how best to group defects is subject to uncertainty, despite a variety of proposed groupings available in the literature (Lin et al., 1999).

In this particular application, we consider $26$ individual heart defects, which have been previously grouped into $6$ categories by investigators (Botto et al., 2007). The prior grouping is shown in Table 3, along with basic summary statistics of the distribution of defects in the analyzed data. We are interested in evaluating the association between heart defects and about $90$ potential risk factors related to mothers’ health status, pregnancy experience, lifestyle and family history. We considered a subset of data from NBDPS, excluding observations with missing covariates, obtaining a dataset with $8,125$ controls, while all heart defects together comprise $4,947$ cases.

5.1 Modeling birth defects

Standard approaches assessing the impact of exposure factors on the risk to develop a birth defect often rely on logistic regression analysis. Let $i=1,\ldots,N$ index birth defects, while $j=1,\ldots,n_{i}$ indicates observations related to birth defect $i$, with $y_{ij}=1$ if observation $j$ has birth defect $i$ and $y_{ij}=0$ if observation $j$ is a control, i.e. does not have any birth defect. Let $\mathbf{X}_{i}$ denote the data matrix associated to defect $i$, with each row $\mathbf{x}^{T}_{ij}=(x_{ij1},\ldots,x_{ijp})$ being the vector of the observed values of $p$ categorical variables for the $j$th observation. At first one may consider $N$ separate logistic regressions of the type

[TABLE]

with $\alpha_{i}$ denoting the defect-specific intercept, and $\boldsymbol{\beta}_{i}$ the $p\times 1$ vector of regression coefficients. However, Table 3 highlights the heterogeneity of heart defect prevalences, with some of them being so few as to preclude separate analyses.

A first step in introducing uncertainty about clustering of the defects may rely on a standard Bayesian nonparametric approach, placing a Dirichlet process prior on the distribution of regression coefficient vector $\boldsymbol{\beta}_{i}$ in order to borrow information across multiple defects while letting the data inform on the number and composition of the clusters. A similar approach has been previously proposed in MacLehose and Dunson (2010), with the aim being to shrink the coefficient estimates towards multiple unknown means. In our setting, an informed guess on the group structure is available through $\boldsymbol{c}_{0}$, reported in Table 3.

We consider a simple approach building on the Bayesian version of the model in (16), and allowing the exposure coefficients $\boldsymbol{\beta}_{i}$ for $i=1,\ldots,N$ to be shared across regressions, while accounting for $\boldsymbol{c}_{0}$. The model written in a hierarchical form is

[TABLE]

where $CP(\boldsymbol{c}_{0},\psi,p_{0}(\boldsymbol{c}_{0}))$ indicates the Centered Partition process, with base partition $\boldsymbol{c}_{0}$, tuning parameter $\psi$ and baseline EPPF $p_{0}(\boldsymbol{c}_{0})$. We specify the baseline EPPF so that when $\psi=0$ the prior distribution reduces to a Dirichlet Process with concentration parameter $\alpha$. Instead, for $\psi\rightarrow\infty$ the model corresponds to $K$ separate logistic regressions, one for each group composing $\boldsymbol{c}_{0}$. The model estimation can be performed by leveraging a Pòlya-Gamma data-augmentation strategy for Bayesian logistic regression (Polson et al., 2013), combined with the procedure illustrated in Section 3.4 for the clustering update step. The Gibbs sampler is detailed in the Appendix (Algorithm Gibbs sampling for shared logistic regression).

5.2 Simulation study

We conduct a simulation study to evaluate the performance of our approach in accurately estimating the impact of the covariates across regressions with common effects, under different prior guesses. In simulating data, we choose a scenario mimicking the structure of our application, considering a number of defects $N=12$ equally partitioned in $4$ groups. We consider $p=10$ dichotomous explanatory variables and assume that defects in the same group have the same covariates effects. We take a different number of observations across defects, with $\{n_{1},n_{2},n_{3}\}=\{100,600,200\}$, $\{n_{4},n_{5},n_{6}\}=\{300,100,100\}$, $\{n_{7},n_{8},n_{9}\}=\{500,100,200\}$, $\{n_{10},n_{11},n_{12}\}=\{200,200,200\}$. For each defect $i$ with $i=1,\ldots,12$ we generate a data matrix $\mathbf{X}_{i}$ by sampling each of the variables from a Bernoulli distribution with probability of success equal to $0.5$. We set most of coefficients $\beta_{i1},\ldots,\beta_{i10}$ to [math], while defining a challenging scenario with small to moderate changes across different groups. In particular we fix $\{\beta_{1},\beta_{2},\beta_{3},\beta_{4}\}=\{0.7,-1.2,0.5,0.5\}$ for group $1$, $\{\beta_{4},\beta_{5},\beta_{6}\}=\{0.7,-0.7,0.7\}$ for group $2$, $\{\beta_{9},\beta_{10}\}=\{0.7,-1.2\}$ for group $3$ and $\{\beta_{1},\beta_{2},\beta_{9},\beta_{10}\}=\{0.7,-0.7,0.7,-0.7\}$ for group $4$. Finally response variables $\mathbf{y}_{i}$ for $i=1,\ldots,12$ are drawn from a Bernoulli distribution with probability of success $p_{i}=\text{logit}(\mathbf{X}_{i}^{T}\boldsymbol{\beta}_{i})$.

We compare coefficients and partition estimates from a grouped logistic regression using a DP prior with $\alpha=1$ and using a CP prior with DP base EPPF with $\alpha=1$. In evaluating the CP prior performances, we consider both the true known partition and a wrong guess. Posterior estimates are obtained using the Gibbs sampler described in the Appendix. We consider a multivariate normal distribution with zero mean vector and covariance matrix $\mathbf{Q}=\text{diag}_{p}(2)$ as base measure for the DP, while we assume the defect-specific intercepts $\alpha_{i}\sim N(0,2)$ for $i=1,\ldots,12$. We run the algorithm for $5,000$ iterations discarding the first $1,000$ as burn-in, with inspection of trace-plots suggesting convergence of the parameters.

In evaluating the resulting estimates under different settings, we take as baseline values for coefficients the maximum likelihood estimates obtained under the true grouping. Figure 6 shows the posterior similarity matrices obtained under the Dirichlet and Centered Partition processes, along with boxplots of the distribution of differences between the coefficients posterior mean estimates and their baseline values, for each of the $12$ simulated defects. We first centered the CP prior on the true known grouping and, according to the considerations made in Section 4.2, we fixed the value of $\psi$ to $15$ for the CP process prior, founding the maximum a posteriori estimate of the partition almost recovering the true underlying grouping expect for merging together the third and fourth group.

We also considered other values for $\psi$ close to $15$, and report the case for $\psi=17$ in Figure 6, for which the true grouping is recovered, with resulting mean posterior estimates of the coefficients almost identical to the baseline. When considering the Dirichlet process, although borrowing information across the defects, it does not distinguish between all the groups but individuate only the first one, while the CP process recovers the true grouping, with better performances in estimating the coefficients.

Finally, we evaluate the CP prior performances when centered on a wrong guess $\boldsymbol{c}_{0}^{\prime}$ of the base partition. In particular, we set $\boldsymbol{c}_{0}^{\prime}=\{1,5,9\}\{2,6,10\}\{3,7,11\}\{4,8,12\}$. Despite having the same configuration of $\boldsymbol{c}_{0}$, it has distance from $\boldsymbol{c}_{0}$ of approximately $3.16$, where the maximum possible distance is $\log_{2}(12)=4.70$. Under such setting we estimate the partition $\hat{\boldsymbol{c}}=\{1,2,3,5\}\{4,6,7,8,9,10,11,12\}$ via maximum at posteriori, obtaining two clusters. Although we center the prior in $\boldsymbol{c}_{0}^{\prime}$, the estimated partition results to be closer to the one induced by the DP ($0.65$) than $\boldsymbol{c}_{0}^{\prime}$ ($2.45$), with also similar performances in the coefficient estimation, which may be interpreted as a suggestion that the chosen base partition is not supported by the data.

5.3 Application to NBDPS data

We estimated the model in (5.1) on the NBDPS data, considering the controls as shared with the aim of grouping cases into informed groups on the basis of the available $\boldsymbol{c}_{0}$. In order to choose a value for the penalization parameter, we consider the prior calibration illustrated in Section 4, finding a value of $\psi=40$ assigning a $90$% probability to partitions within a distance around $0.8$, where the maximum possible distance is equal to $4.70$. In terms of moves on the Hasse diagram we are assigning $90\%$ prior probability to partitions at most at $11$ split/merge operations from $\boldsymbol{c}_{0}$, given that the minimum distance from $\boldsymbol{c}_{0}$ is $2/N\approx 0.07$. To assess sensitivity of the results, we performed the analysis under different values of $\psi\in\{0,40,80,120,\infty\}$. In particular, for $\psi=0$ the clustering behavior is governed by a Dirichlet process prior, while $\psi\rightarrow\infty$ corresponds to fixing the groups to $\boldsymbol{c}_{0}$.

In analyzing the data we run the Gibbs sampler for $10,000$ iterations and use a burn-in of $4,000$, under the same prior settings as in Section 5.2. Figure 8 summarizes the posterior estimates of the allocation matrices under different values of $\psi$, with colored dots emphasizing differences with the base partition $\boldsymbol{c}_{0}$. Under the DP process ($\psi=0$) the estimated partition differs substantially from the given prior clustering. Due to the immense space of the possible clusterings, this is likely reflective of limited information in the data, combined with the tendency of the DP to strongly favor certain types of partitions, typically characterized from few large clusters along many small ones. When increasing the value of the tuning parameter $\psi$ the estimated clustering is closer to $\boldsymbol{c}_{0}$, with a tendency in favoring a total number of three clusters. In particular, for $\psi=120$ one of the groups in $\boldsymbol{c}_{0}$ is recovered (left ventricular outflow), while the others are merged in two different groups. It is worth noticing that AVSD, which is placed in its own group under $\boldsymbol{c}_{0}$, is always grouped with other defects with a preference for ones in the septal class (blue color). Also two defects of this last class, ASD and ASDOS, happen to be lumped together across different values of $\psi$, and are in fact two closely related defects.

Details on the results for each of the estimated models are given in the Appendix (Figures 10-14) and summarized here. Figure 9 shows a heatmap of the mean posterior log odds-ratios for increasing values of the penalization parameter $\psi$, with dots indicating if they are significant according to a $95\%$ credibility interval. In general, the sign of the effects does not change for most of the exposure factors across the different clusterings. Figure 9 focuses on pharmaceutical use in the period from $1$ month before the pregnancy and $3$ months during, along with some exposures related to maternal behavior and health status.

We found consistent results for known risk factors for CHD in general, including for diabetes (Correa et al., 2008) and obesity (Waller et al., 2007). The finding that nausea is associated with positive outcomes is consistent with prior literature (Koren et al., 2014). The association between use of SSRIs and pulmonary atresia was also noted in Reefhuis et al. (2015). It is worth noticing that estimates obtained under the DP prior are less consistent with prior work. In particular, there apparent artifacts such as the protective effect of alcohol consumption related to defects in the bigger cluster, which is mitigated from an informed borrowing across the defects. On the other side, estimates under separate models for AVSD or PAPVR , which corresponds to $0.02\%$ and $0.01\%$ of cases respectively, show how a separate analysis of cases with low prevalence misses even widely assessed risk factors, as for example diabetes.

Discussion

There is a very rich literature on priors for clustering, with almost all of the emphasis on exchangeable approaches, with a smaller literature focused on including dependence on known features (eg., temporal or spatial structure or covariates). The main contribution of this article is to propose what is seemingly a first attempt at including prior information on an informed guess at the clustering structure. We were particularly motivated by a concrete application to a birth defects study in proposing our method, which is based on shrinking an initial clustering prior towards the prior guess.

There are many immediate interesting directions for future research. One thread pertains to developing better theoretical insight and analytical tractability into the new class of priors. For existing approaches, such as product partition models and Gibbs-type partitions, there is a substantial literature providing simple forms of prediction rules and other properties. It is an open question whether such properties can be modified to our new class. This may yield additional insight into the relative roles of the base prior, centering value and hyperparameters in controlling the behavior of the prior and its impact on the posterior.

Another important thread relates to applications of the proposed framework beyond the setting in which we have an exact guess at the complete clustering structure. In many cases, we may have an informed guess or initial clustering in a subset of the objects under study, with the remaining objects (including future ones) completely unknown. Conceptually the proposed approach can be used directly in such cases, and also when one has different types of prior information on the clustering structure than simply which objects are clustered together.

acknowledgement

The authors gratefully acknowledge this work was supported in part through cooperative agreements from the Centers for Disease Control and Prevention to the centers participating in the National Birth Defects Prevention Study and by the National Institutes of Health (R01ES027498; U50CCU422096; 5U01DD001036; PA96043; PA 02081; FOA DD09-001).

Appendix

Prior calibration

Algorithm 1 : Estimation of counts statistics related to distances neighborhoods of $\boldsymbol{c}_{0}$

  Local search

  0. Start from the base partition $\boldsymbol{c}_{0}$ with $|K_{0}|$ clusters and configuration $\boldsymbol{\lambda}_{m_{0}}$ and set $\delta_{0}=0$ and $\mathcal{N}_{0}(\boldsymbol{c}_{0})=\boldsymbol{c}_{0}$.

  for $t=1,\ldots,T$ do

     Obtain $\mathcal{N}_{t}(\boldsymbol{c}_{0})$ from partitions in $\mathcal{N}_{t-1}(\boldsymbol{c}_{0})$ by exploring all directed connections, i.e. partitions obtained with one operation of split/merge on elements $\mathcal{N}_{t-1}(\boldsymbol{c}_{0})$.

  end for

  2. Compute the distance from $\boldsymbol{c}_{0}$ and all partitions in $\mathcal{N}_{T}(\boldsymbol{c}_{0})$ and take the minimum distance, $\delta_{L^{*}}$; discard all partitions having distances greater than $\delta_{L^{*}}$.

  3. Obtain counts $n_{l}$ and $n_{lm}$ relative to distances $\delta_{1},\ldots,\delta_{L^{*}}$ for $m=1,\ldots,M$.

  Monte Carlo approximation

  for $r=1,\ldots,R$ do

     4. Sample the number of clusters $K$ from the discrete probability distribution

[TABLE]

     5. Conditional on $K$ generate a partition $\boldsymbol{c}^{(r)}=\{c_{1}^{(r)},\ldots,c_{N}^{(r)}\}$ by sampling each $c^{(r)}_{i}$ from a discrete uniform distribution on $\{1,\ldots,K\}$.

     6. If $d(\boldsymbol{c}^{(r)},\boldsymbol{c}_{0})>\delta_{L^{*}}$ reject the partition.

  end for

  7. Let $R^{*}$ be the number of accepted partitions, and estimate counts $\hat{n}_{l}$ and $\hat{n}_{lm}$ for $m=1,\ldots,M$ and according to (14)-(15) conditional on the observed distance values $\hat{\delta}_{(L^{*}+1)},\ldots,\hat{\delta}_{L}$.

  8. Using $R^{*}$ be the number of accepted partitions, and estimate counts $\hat{n}_{l}$ and $\hat{n}_{lm}$ relative to distances $\hat{\delta}_{L^{*}+1},\ldots,\hat{\delta}_{L}$ for $m=1,\ldots,M$.

Marginal sampling using variation of information

We describe how to compute the penalization term in the marginal sampling step described in Section 3.4 using the Variation of Information as a distance, but the same procedure applies when using other distances based on blocks sizes. Let $K^{-}$ and $K_{0}^{-}$ denote respectively the number of clusters in $\boldsymbol{c}^{-i}$ and $\boldsymbol{c}_{0}^{-i}$, i.e. partitions $\boldsymbol{c}$ and $\boldsymbol{c}_{0}$ after removing the $i$ observation.

Algorithm 2 : Computation strategy for the penalization term in marginal sampling

  Let $K^{-}$ and $K_{0}^{-}$ denote respectively the number of clusters in $\boldsymbol{c}^{-i}$ and $\boldsymbol{c}_{0}^{-i}$, i.e. partitions $\boldsymbol{c}$ and $\boldsymbol{c}_{0}$ after removing the $i$ observation.

  for $i=1,\ldots,N$ do

     1. Compute cardinalities $\{\lambda_{1}^{-i},\ldots,\lambda_{K^{-}}^{-i}\}$ representing the number of observations in each cluster for $\boldsymbol{c}^{-i}$.

     2. Compute $\lambda_{lm}^{-i}$, the number of observations in cluster $l$ under $\boldsymbol{c}^{-i}$ and cluster $m$ under $\boldsymbol{c}_{0}^{-i}$ for $l=1,\ldots,K^{-}$ and $m=1,\ldots,K_{0}^{-}$.

     for $k=1,\ldots,K^{-},K^{-}+1$ do

        Let $c_{i,0}$ be the cluster of index $i$ under partition $\boldsymbol{c}_{0}$.Compute $d(\boldsymbol{c},\boldsymbol{c}_{0})\propto-H(\boldsymbol{c})+2H(\boldsymbol{c}\land\boldsymbol{c}_{0})$ for $\boldsymbol{c}=\{\boldsymbol{c}^{-i}\cup k\}$ using

[TABLE]

     end for

  end for

Gibbs sampling for shared logistic regression

In estimating the model, a Pólya-gamma data augmentation strategy is employed; for each $y_{ij}$ we introduce a latent variable $\omega_{ij}\sim PG(1,\alpha_{i}+\mathbf{x}_{ij}^{T}\boldsymbol{\beta}_{c_{i}})$ for each observation $j$ in defect-specific dataset $i$ for $i=1,\ldots,N$.

Algorithm 3 : Gibbs sampling for posterior computation

  Conditionally on the cluster allocation vector $\mathbf{c}=(c_{1},\ldots,c_{n})$ and data $\{\mathbf{y}_{i},\mathbf{X}_{i}\}$ for $i=1,\ldots,N$, update mixture related parameters and Pólya-gamma latent variables as follows.

  ——————————————————————————————————–

   [1] Sample Pólya-gamma latent variables for each observation in each dataset

  for $i=1,\ldots,N$ and $j=1,\ldots,n_{i}$ do

     

[TABLE]

  end for

  ——————————————————————————————————–

   [2] Update defect-specific intercept, exploiting Pólya-gamma conjugancy

  for $i=1,\ldots,N$ do

     

[TABLE]

with $\tau^{*}=\tau_{0}+\sum_{j=1}^{n_{i}}\omega_{ij}$ and $a^{*}=[a_{0}\tau_{0}+\sum_{j=1}^{n_{i}}(y_{ij}-1/2-\omega_{ij}\mathbf{x}_{ij}^{T}\boldsymbol{\beta}_{c_{i}})]/\tau^{*}$

  end for

  ——————————————————————————————————–

   [3] Defining $\kappa_{ij}:=y_{ij}-1/2-w_{ij}\alpha_{i}$, then the vector ($\kappa_{ij}/\omega_{ij}|c_{i}=h,\omega_{ij})\sim\mathcal{N}(\mathbf{x}_{i}^{T}\boldsymbol{\beta}^{(k)},1/\omega_{ij})$, and each cluster-specific coefficient vector $\boldsymbol{\beta}_{h}$ can be updated by aggregating all observations and augmented data relative to birth defects that are in the same cluster.

  for $k=1,\ldots,K$ do

     Let $\mathbf{X}^{(k)}$, $\mathbf{y}^{(k)}$, $\boldsymbol{\kappa}^{(k)}$ be the obtained quantities relative to cluster $h$, and $\mathbf{\Omega}^{(k)}$ a diagonal matrix with the corresponding Pólya-gamma augmented variables. Then update cluster-specific coefficients vector from

[TABLE]

with $\mathbf{Q}^{(k)}=(\mathbf{X}^{(h)T}\mathbf{\Omega}^{(k)}\mathbf{X}^{(k)}+\mathbf{Q}^{-1})^{-1}$ and $\mathbf{b}^{(k)}=\mathbf{Q}^{(k)}(\mathbf{X}^{(k)T}\boldsymbol{\kappa}^{(k)}+\mathbf{Q}^{-1}\mathbf{b})$.

  end for

  —————————————————————————————————–

   [4] Allocate each birth defect $i$ to one of the clusters

  for $i=1,\ldots,N$ do

     Sample the class indicator $c_{i}$ conditionally on $\boldsymbol{c}_{-i}=(c_{1},\ldots,c_{i-1},c_{i+1},\ldots,c_{n})$ from the discrete distribution with probabilities

[TABLE]

with

[TABLE]

being the model likelihood evaluated for cluster $h$ and $\text{Pr}(c_{i}=h|\boldsymbol{c}^{-(i)})$ computed as described in Section 3.4.

  end for

Results for NBDPS data application