# Graphical Dirichlet Process for Clustering Non-Exchangeable Grouped Data

Arhit Chakrabarti Department of Statistics, Texas A&M University, College
Station, TX Yang Ni Department of Statistics, Texas A&M University, College
Station, TX All correspondence should be addressed to<EMAIL_ADDRESS>Ellen
Ruth A. Morris Department of Nutrition, Texas A&M University, College Station,
TX  Program in Integrative Nutrition & Complex Diseases, Texas A&M University,
College Station, TX  Current address: Texas A&M Veterinary Medical Diagnostic
Laboratory, College Station, TX  Michael L. Salinas Department of Nutrition,
Texas A&M University, College Station, TX  Program in Integrative Nutrition &
Complex Diseases, Texas A&M University, College Station, TX  Robert S. Chapkin
Department of Nutrition, Texas A&M University, College Station, TX  Program in
Integrative Nutrition & Complex Diseases, Texas A&M University, College
Station, TX  Bani K. Mallick Department of Statistics, Texas A&M University,
College Station, TX

###### Abstract

We consider the problem of clustering grouped data with possibly non-
exchangeable groups whose dependencies can be characterized by a known
directed acyclic graph. To allow the sharing of clusters among the non-
exchangeable groups, we propose a Bayesian nonparametric approach, termed
graphical Dirichlet process, that jointly models the dependent group-specific
random measures by assuming each random measure to be distributed as a
Dirichlet process whose concentration parameter and base probability measure
depend on those of its parent groups. The resulting joint stochastic process
respects the Markov property of the directed acyclic graph that links the
groups. We characterize the graphical Dirichlet process using a novel
hypergraph representation as well as the stick-breaking representation, the
restaurant-type representation, and the representation as a limit of a finite
mixture model. We develop an efficient posterior inference algorithm and
illustrate our model with simulations and a real grouped single-cell dataset.

Keywords: Bayesian nonparametrics, clustering, directed acyclic graph, family-
owned restaurant process, non-exchangeable groups.

##  1 Introduction

This article considers clustering of grouped data where the groups are _non-
exchangeable_. We are interested in settings where the data are _partially
exchangeable_ (de Finetti, 1938), which entails the exchangeability of the
observations within each group but not across the groups. We consider
dependent group–specific random probability measures, thereby allowing the
borrowing of information across non-exchangeable groups. We represent the
dependencies among groups through a known _directed acyclic graph_ (DAG) with
nodes denoting groups and directed edges denoting the group dependencies. Such
data are abundant in many areas such as genomics. For example, our motivating
application is a single-cell RNA-sequencing (scRNA-seq) study that aimed to
investigate intestinal stem cell differentiation processes in mice with
colorectal cancer. The experiments started from a baseline group where the
mice were genetically wild-type, fed with a normal diet, and treated with no
cancer therapy (placebo). Then to understand the main effects of genotype,
diet, and cancer therapy on colonic crypt and tumor niche cell composition,
the experimenters introduced three new groups of mice, each differing from the
baseline group by exactly one factor (Apc knock-out, a high-fat diet, or a new
cancer treatment AdipoRon). To determine the two-way interaction effects,
three additional groups of mice were studied, each of which differed from the
baseline group by two factors (e.g., mice with Apc knock-out, a high-fat diet,
and no cancer treatment). Lastly, for a three-way interaction, they introduced
the eighth group of mice with Apc knock-out, a high-fat diet, and the new
treatment AdipoRon. The progression of these experiments from baseline to the
study of main effects, two-way interactions, and three-way interactions
manifests the non-exchangeability of the experimental groups (e.g., the
baseline group is expected to be more similar to the “main effect” groups than
the “three-way interaction” group). With this grouped scRNA-seq dataset, our
goal is to cluster cells based on gene expression at the single-cell level
within each experimental group while allowing information to be shared across
these non-exchangeable groups with a novel DAG-based Bayesian nonparametric
model.

The _Dirichlet process_ (DP, Ferguson, 1973) and its variations (De Blasi et
al., 2013; Barrios et al., 2013) have been the backbone of numerous model-
based Bayesian nonparametric clustering methods (Hjort et al., 2010; Müller et
al., 2015). The DP, $DP(\alpha_{0},G_{0})$, is a probability measure on
probability measures, where $\alpha_{0}>0$ is the concentration parameter and
$G_{0}$ is a base probability measure. There have been extensive studies on DP
mixture models (Antoniak, 1974; Escobar and West, 1995; MacEachern and Müller,
1998), which enable clustering without having to fix the number of clusters _a
priori_. When there are groups present in the data, naively, one could
consider either a separate DP mixture model for each group on one extreme or a
single DP mixture model ignoring the groups on the other extreme. However, it
is often desirable to identify group-specific clusters while allowing the
groups to be linked so that clusters are comparable across groups. Given the
goal of clustering the observations within each group, consider a set of
random probability measures, $G_{j}$, one for each group $j$, where each
$G_{j}$ is distributed as $DP(\alpha_{0j},G_{0j})$ with group-specific
concentration parameter $\alpha_{0j}$ and base probability measure $G_{0j}$.
Many methods have been proposed to link these group-specific DPs to induce
dependencies through the parameter $\alpha_{0j}$ and/or $G_{0j}$ (Cifarelli
and Regazzini, 1978; Mallick and Walker, 1997; Kleinman and Ibrahim, 1998;
Müller et al., 2004). Perhaps one of the most well-known methods is the
hierarchical Dirichlet process (HDP, Teh et al., 2006), which falls in the
general framework of dependent DP (MacEachern, 1999, 2000) and assumes each
group-specific $G_{j}$ is distributed as $DP(\alpha_{0},G_{0})$ where
$\alpha_{0}$ is the shared concentration parameter and $G_{0}$ is the shared
base probability measure for all groups. They further assume that $G_{0}$
follows another DP, $G_{0}\sim DP(\gamma,H)$. Since draws from a DP are
discrete with probability one (Sethuraman, 1994), the base measure $G_{0}$ is
almost surely discrete, which ensures that the group-specific probability
measure $G_{j}$ shares the same set of atoms. The corresponding HDP mixture
model is thus capable of identifying group-specific clusters while borrowing
strength across groups. By construction, HDP mixture model assumes that both
the observations within each group and the groups are exchangeable. A similar
approach with a different scope, the nested DP (Rodríguez et al., 2008),
assumes $G_{j}$ follows a DP-distributed random probability measure with
another DP as the base measure, $G_{j}\sim Q$ and $Q\sim
DP(\alpha_{0},DP(\gamma,H))$. The nested structure allows for the clustering
of groups but restricts the clusters of observations within each group to be
either identical or completely unrelated across groups. Similarly to HDP,
nested DP also assumes both the observations within each group and the groups
to be exchangeable. Several recent works (Camerlenghi et al., 2019; Lijoi et
al., 2022) have been proposed to take advantage of the cluster-sharing feature
of the HDP and the group-clustering feature of the nested DP. Dependent DP has
also been extensively used to model random distributions with various other
types of dependencies such as spatial and temporal dependencies (Iorio et al.,
2004, 2006; Dunson and Herring, 2005; Gelfand et al., 2005; Griffin and Steel,
2006; Nieto-Barajas and Contreras-Cristán, 2014; Dahl et al., 2017); see
Quintana et al., 2020 for a recent review of different dependent DPs.

In this paper, we are interested in modeling a set of group-specific random
distributions of which the (conditional) dependencies can be characterized by
a DAG whose nodes represent the groups. More precisely, we assume that the
joint distribution of the set of group-specific random distributions
factorizes with respect to a DAG and, therefore, respects its Markov property
(i.e., conditional independencies). We call such graph-dependent DP, the
_graphical Dirichlet process_ (GDP). Using GDP as a mixing distribution, the
GDP mixture model gives rise to group-specific clusters, which depend directly
on their Markov blanket. The known flexibility of DAG in representing
conditional dependencies renders the generality of the proposed GDP for
modeling dependent random distributions and group-specific clusters beyond
exchangeable groups. The use of DAGs in Bayesian nonparametrics has been
considered in recent literature. Dey et al., 2022 proposed a graphical
Gaussian process to parsimoniously model multivariate spatial data by
incorporating conditional independencies among variables encoded by a DAG. Gu
and Dunson, 2023 proposes a pyramid-shaped deep latent variable model for
categorical data using a DAG to represent the layer-wise latent conditional
dependency structure. These works showcased the usefulness of DAGs through
their factorization in Bayesian nonparametric models. We also exploit such
factorization in this paper but our model is significantly different from
theirs in both approaches and scopes. For example, their graphs link variables
whereas ours link groups, and they focus on the modeling of multivariate
spatial fields or generative models for categorical data whereas we focus on
clustering non-exchangeable grouped data. The proposed GDP is a general model.
The well-known HDP is a special case of GDP with a specific type of DAG – a
fork, i.e., one parent node and many children nodes (detailed in Section 3.1);
see Figure 1(a). Several existing works on time-evolving topic models can also
be reformulated using a DAG to capture the time-dependency structure (Srebro
and Roweis, 2005; Ren et al., 2008; Zhang et al., 2010).

In this paper, we will characterize the proposed GDP by a novel _hypergraph_
representation. We will also provide several other representations analogous
to those for the HDP, i.e., a stick-breaking representation, a restaurant-type
representation, and a representation as an infinite limit of a finite mixture
model. We develop efficient posterior sampling based on the SALTSampler
(Director et al., 2017) and a Blocked Gibbs sampler for DP/HDP (Ishwaran and
James, 2001; Das et al., 2023). Simulations and the motivating grouped single-
cell data are used to demonstrate our method. In summary, our main
contribution is three-fold:

  1. 1.

We propose a general Bayesian nonparametric approach, GDP, to incorporate non-
exchangeable group dependencies for clustering.

  2. 2.

We provide several characterizations of GDP, each providing a different
perspective.

  3. 3.

We develop a Metropolis-within-blocked-Gibbs sampler for posterior inference.
Since HDP is a special case of GDP, this also contributes to a new sampler for
HDP. The difficulty of sampling the global weights for HDP is mitigated by
using the specialized proposal of SALTSampler (Director et al., 2017).

$G_{1}^{(0)}$$G_{0}$$G_{1}^{(1)}$$G_{2}^{(1)}$$G_{l_{1}-1}^{(1)}$$G_{l_{1}}^{(1)}$…
(a) HDP

$G_{1}^{(0)}$$G_{0}$$G_{1}^{(1)}$$G_{2}^{(1)}$$G_{l_{1}-1}^{(1)}$$G_{l_{1}}^{(1)}$$G_{1}^{(2)}$$G_{2}^{(2)}$$G_{l_{2}-1}^{(2)}$$G_{l_{2}}^{(2)}$……$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$
(b) GDP for DAG with a unique root node.

Figure 1: Schematic illustration of HDP and GDP. HDP is a special case of GDP
when the DAG is a fork.

The remainder of the paper is organized as follows. Section 2 provides a brief
overview of some preliminaries needed for the remainder of the paper. Section
3 introduces the proposed GDP and the corresponding nonparametric mixture
model. We introduce the hyperpriors of our model and also present two lemmas,
which are the backbone of our main result in Theorem 1. In Section 4, we
present different representations of the proposed GDP. In Section 5, we
provide simulations to illustrate our method. Section 6 presents a real data
analysis using the proposed method on the motivating single-cell data. The
paper concludes with a brief discussion in Section 7.

##  2 Preliminaries

###  2.1 Directed acyclic graph

We first provide a brief background on DAG. Let $D=\left(V,E\right)$ be a DAG
consisting of a set of nodes $V=\left\\{1,2,\dots,p\right\\}$ and a set of
directed edges $E\subset V\times V$ that does not contain any directed cycles.
We denote a directed edge from the node $i$ to node $j$ by $j\leftarrow i$ and
call $i$ a _parent_ of $j$. A node without parents is called a _root_. For a
DAG, there exists at least one root. Let $\bm{Y}=\\{Y_{1},\dots,Y_{p}\\}$ be a
set of random variables. Every node $j\in V$ represents a random variable
$Y_{j}$; later in this paper, $Y_{j}$ will be a random probability measure. In
a DAG model, also known as a Bayesian network, the probability distribution
$\mathcal{P}(\bm{Y})$ is assumed to factorize over $D$,
$\mathcal{P}(\bm{Y})=\prod_{j=1}^{p}\mathcal{P}(Y_{j}\mid Y_{pa(j)})$, where
$pa(j)=\\{k\in V|j\leftarrow k\\}$ denotes the collection of parents of node
$j$. This DAG factorization implies that the distribution $\mathcal{P}$
respects the conditional independence relationships encoded by the graph $D$
via the notion of d-separation (Pearl, 2009); and vice versa. For instance,
any node is conditional independent of its non-descendants given its parents,
i.e., $Y_{j}\perp Y_{nd(j)}|Y_{pa(j)}$ for any $j\in V$ where $\perp$ denotes
independence, $nd(j)=V\backslash de(j)\backslash\\{j\\}$ denotes the non-
descendants of node $j$, and $de(j)=\\{k\in V|k\leftarrow\cdots\leftarrow
j\\}$ denotes the descendants of node $j$. A _Markov blanket_ of any node $j$
from $V$ is any subset $V_{1}$ of $V$ such that $Y_{j}\perp Y_{V\backslash
V_{1}}|Y_{V_{1}}$. In other words, $V_{1}$ contains all the information in $V$
about the node $j$. DAG models are convenient tools to parsimoniously specify
a multivariate distribution through its conditionals, which is especially
useful in this paper for specifying a multivariate distribution of a set of
random probability measures.

###  2.2 Infinite mixture model

Next, we present a brief overview of infinite mixture models for a single
population, the DP mixture model, and for multiple exchangeable populations,
the HDP mixture model.

####  2.2.1 Dirichlet process mixture model

For a single population, let $x_{i}$ denote the $i$th realization of a random
variable $X$. We consider a mixture model,

$\displaystyle\theta_{i}\mid G$ $\displaystyle\overset{iid}{\sim}G,$ (1)
$\displaystyle x_{i}\mid\theta_{i}$
$\displaystyle\overset{ind}{\sim}F(\theta_{i}),$

where $F(\theta_{i})$ denotes the distribution of $x_{i}$ parameterized by
$\theta_{i}$. The parameters $\theta_{i}$’s are conditionally independent
given the prior distribution $G$. In a DP mixture model, $G$ is assigned a DP
prior, $G\sim DP(\alpha_{0},G_{0})$ with concentration $\alpha_{0}$ and base
probability measure $G_{0}$.

Sethuraman, 1994 presented the _stick-breaking representation_ of the DP based
on independent sequences of i.i.d. random variables
$(\pi_{k}^{\prime})_{k=1}^{\infty}$ and $(\phi_{k})_{k=1}^{\infty}$, which is
given by,

$\displaystyle\pi_{k}^{\prime}$
$\displaystyle\mid\alpha_{0}\overset{iid}{\sim}Beta(1,\alpha_{0}),$
$\displaystyle\phi_{k}\mid G_{0}\overset{iid}{\sim}G_{0},$ (2)
$\displaystyle\pi_{k}$
$\displaystyle=\pi_{k}^{\prime}\prod_{l=1}^{k-1}(1-\pi_{l}^{\prime}),$
$\displaystyle G=\sum_{k=1}^{\infty}\pi_{k}\delta_{\phi_{k}},$ (3)

where $\delta_{\phi}$ is a point mass at $\phi$ and $\phi_{k}$’s are called
the _atoms_ of $G$. The sequence of random weights
$\bm{\pi}=(\pi_{k})_{k=1}^{\infty}$ constructed from Eq. 2 and Eq. 3 satisfies
$\sum_{k=1}^{\infty}\pi_{k}=1$ with probability one. The random probability
measure on the set of integers is denoted by
$\bm{\pi}\sim\text{GEM}(\alpha_{0})$ for convenience where GEM stands for
Griffiths, Engen and McCloskey (Pitman, 2002). It is clear from Eq. (1) and
Eq. (3) that $\theta_{i}$ takes the value $\phi_{k}$ with probability
$\pi_{k}$. Let $z_{i}$ be a categorical variable such that $z_{i}=k$ if
$\theta_{i}=\phi_{k}$. An equivalent representation of a Dirichlet process
mixture is given by,

$\displaystyle\bm{\pi}\mid\alpha_{0}$
$\displaystyle\sim\text{GEM}(\alpha_{0}),$ $\displaystyle z_{i}\mid\bm{\pi}$
$\displaystyle\overset{iid}{\sim}\bm{\pi},$ (4) $\displaystyle\phi_{k}\mid
G_{0}$ $\displaystyle\overset{iid}{\sim}G_{0},$ $\displaystyle x_{i}\mid
z_{i},(\phi_{k})_{k=1}^{\infty}$
$\displaystyle\overset{ind}{\sim}F(\phi_{z_{i}}).$

####  2.2.2 Hierarchical Dirichlet process mixture model

Suppose observations are now organized into multiple exchangeable groups. Let
$x_{ji}$ denote the observation $i$ from group $j$ and $\theta_{ji}$ denote
the parameter specifying the mixture component associated with the
corresponding observation. Let $F(\theta_{ji})$ denote the distribution of
$x_{ji}$ given $\theta_{ji}$ and $G_{j}$ denote a prior distribution for
$\theta_{ji}$. The group-specific mixture model is given by,

$\displaystyle\theta_{ji}\mid G_{j}$ $\displaystyle\overset{ind}{\sim}G_{j},$
(5) $\displaystyle x_{ji}\mid\theta_{ji}$
$\displaystyle\overset{ind}{\sim}F(\theta_{ji}).$

As with the DP mixture model, when the random measures $G_{j}$’s are assigned
an HDP prior,

$\displaystyle G_{0}\mid\gamma,H$ $\displaystyle\sim DP(\gamma,H),$ (6)
$\displaystyle G_{j}\mid\alpha_{0},G_{0}$ $\displaystyle\sim
DP(\alpha_{0},G_{0}),$

the corresponding mixture model is referred to as the HDP mixture model. The
global random probability measure $G_{0}$ is distributed as a DP with
concentration parameter $\gamma$ and base probability measure $H$. The group-
specific random measures $G_{j}$’s are conditionally independent given $G_{0}$
and hence are exchangeable (de Finetti, 1938). They are distributed as DP with
the base measure $G_{0}$ and some concentration parameter $\alpha_{0}$. The
probability model (5) along with (6) completes the specification of an HDP
mixture model. Because DP-distributed $G_{0}$ is almost surely discrete, the
atoms of $G_{j}$’s and hence the group-specific clusters are necessarily
shared across groups.

##  3 Graphical Dirichlet Process

When groups are non-exchangeable (e.g., due to study design), the joint
distribution of $G_{j}$’s specified by (6) may not be appropriate. Our
approach to the problem of sharing clusters among non-exchangeable groups is
through specifying a general joint distribution of $G_{j}$’s that respect the
Markov property of a DAG $D$ that links the groups. We assume that the
underlying DAG $D$ is known and we define the appropriate prior on the nodes
of the DAG and refer to the resulting stochastic process on the graph as the
graphical Dirichlet process (GDP). We show how this prior can be used in the
non-exchangeable grouped mixture model setting.

###  3.1 The Proposed GDP

Let the nodes $V$ of DAG $D=(V,E)$ now represent the group-specific random
probability measures $G_{j}$’s. The edges $E$ represent the conditional
dependence of $G_{j}$’s. Then the joint distribution of the random probability
measures follows the DAG factorization
$\mathcal{P}(G_{1},\dots,G_{p}|D)=\prod_{j=1}^{p}\mathcal{P}(G_{j}\mid
G_{pa(j)})$, where $G_{pa(j)}$ is the set of random probability measures
indexed by the parents $pa(j)$ of node $j$. For convenience, we assume $D$ has
a unique root; see Figure 1(b). This assumption does not diminish the
generality of our approach as a DAG with multiple roots can always be
converted, without losing any conditional dependencies, to a DAG with a unique
root by simply augmenting the DAG with a hidden common parent of the roots;
that hidden common parent becomes the unique root of the new DAG (Figure 2).
The augmentation only changes the Markov blanket of the original root nodes.
Specifically, the Markov blanket of any original root node is simply augmented
with the hidden parent node. As the Markov blanket of any other node remains
unchanged, the distributions of all other nodes remain the same, and hence
this augmentation does not alter the conditional dependencies of the original
DAG.

Let us introduce a few terms before describing the proposed GDP. We denote the
root node, which may be hidden, as the _layer 0_ of DAG $D$. The child nodes
of the root node are termed as the _layer-1_ nodes, and we assume that there
are $l_{1}$ of them. Similarly, we assume that there are a total of $l_{2}$
child nodes from the layer-1 nodes, which we refer to as the _layer-2_ nodes.
We assume that there are $K$ layers in the given DAG $D$ and at any layer $k$,
there are $l_{k}$ nodes. The total number of non-root nodes is
$\sum_{k=1}^{K}l_{k}=p$. We define the concentration parameters and random
measures of node $j$ in the layer $k$ of DAG $D$ as $\alpha_{j}^{(k)}$ and
$G_{j}^{(k)},\ \ j=1,\dots,l_{k}$. We denote by $an^{(k,l)}(j)$ the collection
of generation-$l$ ancestors of node $j$ in layer $k$ of the DAG. For example,
$an^{(k,1)}(j)$ denotes the parents (generation-$1$ ancestors) of the node $j$
in layer $k$, and $an^{(k,2)}(j)$ denotes the collection of the parents of the
nodes in $an^{(k,1)}(j)$ or in other words, $an^{(k,2)}(j)$ denotes the
collection of “grand-parents” (generation-$2$ ancestors) of node $j$ in layer
$k$ of the DAG.

$G_{1}^{(0)}$$G_{1}^{(1)}$$G_{2}^{(1)}$$G_{l_{1}-1}^{(1)}$$G_{l_{1}}^{(1)}$$G_{1}^{(2)}$$G_{2}^{(2)}$$G_{l_{2}-1}^{(2)}$$G_{l_{2}}^{(2)}$……$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$
Figure 2: DAG augmented with a hidden root $G_{1}^{(0)}$, indicated by the
dashed arrows. The original root nodes are
$G_{1}^{(1)},\dots,G_{l_{1}}^{(1)}$.

We define GDP recursively from layer 0, the root node,

$G_{1}^{(0)}\mid\alpha_{1}^{(0)},G_{0}\sim
DP\left(\alpha_{1}^{(0)},G_{0}\right),$ (7)

where $G_{0}$ is a fixed base probability measure. Then the distribution of
the random probability measure of node $j$ in layer $k$ of DAG $D$ conditional
on the concentration parameters and random probability measures of its parent
nodes is given by,

$G_{j}^{(k)}\mid\alpha_{j}^{(k)},\\{G_{l}^{(k-1)}:l\in an^{(k,1)}(j)\\}\sim
DP\left(\alpha_{j}^{(k)},\sum_{l\in
an^{(k,1)}(j)}\pi_{jl}^{(k)}G_{l}^{(k-1)}\right),$ (8)

for $j=1,2,\dots,l_{k}$. In other words, node $j$ in layer $k$ of the DAG is
distributed according to a DP with its own concentration parameter
$\alpha_{j}^{(k)}$ and its base distribution being a weighted average of the
random probability measures of its parents in layer $k-1$ of the DAG,
$\\{G_{l}^{(k-1)}:l\in an^{(k,1)}(j)\\}$, where the weights are given by
$\\{\pi_{jl}^{(k)}:l\in an^{(k,1)}(j)\\}$, which have a unit sum $\sum_{l\in
an^{(k,1)}(j)}\pi_{jl}^{(k)}=1$. Moreover, from the Markov properties of DAG
$D$, $G_{j_{1}}^{(k)}$ and $G_{j_{2}}^{(k)}$ are conditionally independent
given their parents, $\\{G_{l}^{(k-1)}:l\in an^{(k,1)}(j_{1})\\}$ and/or
$\\{G_{l}^{(k-1)}:l\in an^{(k,1)}(j_{2})\\}$, and $G_{j}^{(k)}$ is
conditionally independent of all other random probability measures given its
Markov blanket.

We remark that HDP is a special case of the proposed GDP with a specific DAG,
fork-DAG (Figure 1(a)). Using the notations introduced, a fork-DAG is a DAG
with a unique root node and only one layer of $l_{1}$ child nodes. With this
specific DAG, the GDP is given by

$\displaystyle G_{1}^{(0)}\mid\alpha_{1}^{(0)},G_{0}\sim
DP\left(\alpha_{1}^{(0)},G_{0}\right),$ $\displaystyle
G_{j}^{(1)}\mid\alpha_{j}^{(1)},G_{1}^{(0)}\sim
DP\left(\alpha_{j}^{(1)},G_{1}^{(0)}\right),\hskip
5.69046ptj=1,2,\dots,l_{1},$

which is clearly an HDP.

###  3.2 GDP mixture model

To cluster observations that are organized into possibly non-exchangeable
groups, we use the proposed GDP in Section 3.1 as a mixing distribution of a
mixture model. Letting $j$ index the groups and $i$ index the observations
within each group, we assume that the observations
$x_{j1},x_{j2},\dots,x_{jn_{j}}$ are exchangeable within each group $j$ but
the groups may not be exchangeable. We assume that each observation within a
group is drawn independently from the mixture model (5) and $G_{j}$’s follow
the GDP (7) and (8).

###  3.3 Hyperpriors

We assign a Dirichlet prior on the weights $\\{\pi_{jl}^{(k)}\,:\,l\in
an^{(k,1)}(j)\\}$ in (8),

$\\{\pi_{jl}^{(k)}\,:\,l\in an^{(k,1)}(j)\\}\sim
Dir\left(\\{\alpha_{l}^{(k-1)}:l\in an^{(k,1)}(j)\\}\right),$ (9)

where the parameters $\\{\alpha_{l}^{(k-1)}:l\in an^{(k,1)}(j)\\}$ correspond
to the concentration parameters of the parents (generation-$1$ ancestors) of
node $j$. Since the concentration parameter of a DP relates to its precision
(inverse-variance), assuming a Dirichlet prior for the mixture weights of any
node with Dirichlet parameters proportional to the precisions of the parent
nodes is a natural choice. This gives more “weightage” to a parent node with a
higher precision as opposed to a parent node with a lower precision.

The other distributional consideration that significantly simplifies the
distribution of the random measure of any particular node is by considering a
gamma-DAG distribution on the concentration parameters $\alpha_{j}^{(k)}$’s,
which, like the distribution of $G_{j}^{(k)}$’s, also respects the same Markov
property of DAG $D$. Specifically, we assume that

$\displaystyle\alpha_{1}^{(0)}\mid\alpha_{0}\sim Gamma(\alpha_{0},1),$ (10)
$\displaystyle\alpha_{j}^{(k)}\mid\\{\alpha_{l}^{(k-1)}:l\in
an^{(k,1)}(j)\\}\sim Gamma\left(\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)},1\right),$ $\displaystyle\hskip
5.69046ptj=1,2,\dots,l_{k}.$

In other words, the concentration parameter of the root node follows a gamma
distribution with a fixed shape $\alpha_{0}$ and a unit rate. The
concentration parameter at any level of the DAG follows a conditionally gamma
distribution with the shape parameter equal to the sum of the shape parameters
of its parents. Such a choice of Gamma hyperprior on the concentration
parameters of bottom level DPs of HDP have been considered in Williamson et
al., 2013. We extend such a construction for the more general framework of our
proposed GDP. In the next section, we will see how our choice of hyperpriors
and hyperparameters leads to several compact representations of the proposed
GDP, which requires two lemmas. The first lemma is Lemma 3.1 from Sethuraman,
1994, which we state here.

######  Lemma 1 (Sethuraman, 1994).

Let $\bm{\alpha}_{1}=(\alpha_{11},\alpha_{12},\dots,\alpha_{1k})$ and
$\bm{\alpha}_{2}=(\alpha_{21},\alpha_{22},\dots,\alpha_{2k})$ be
$k$-dimensional vectors with $\alpha_{ij}>0\hskip 5.69046pt\forall\hskip
5.69046ptj=1,2,\dots,k,\hskip 5.69046pti=1,2$. Let $\bm{X}_{1}$ and
$\bm{X}_{2}$ be independent $k$-dimensional random vectors distributed as
Dirichlet distribution with parameters $\bm{\alpha}_{1}$ and
$\bm{\alpha}_{2}$, respectively. Let
$\alpha_{1\cdot}=\sum_{j=1}^{k}\alpha_{1j}$ and
$\alpha_{2\cdot}=\sum_{j=1}^{k}\alpha_{2j}$. Let $\pi$ be independent of
$\bm{X}_{1}$ and $\bm{X}_{2}$ and have a beta distribution
$Beta\left(\alpha_{1\cdot},\alpha_{2\cdot}\right)$. Then the distribution of
$\pi\bm{X}_{1}+\left(1-\pi\right)\bm{X}_{2}$ is the Dirichlet distribution
with parameter $\bm{\alpha}_{1}+\bm{\alpha}_{2}$.

The proof is provided in Section A of the Supplementary Materials for
completeness. The next lemma is an immediate extension of Lemma 1 for more
than two independent Dirichlet distributed random vectors. As the Dirichlet
distribution is a multivariate analog of the beta distribution, by considering
a Dirichlet distribution on the weights, we arrive at a similar result. This
lemma is a finite-dimensional version of Theorem 1 of Williamson et al., 2013,
which essentially states that a finite Dirichlet mixture of DPs is, in turn, a
DP with its concentration parameter being the sum of the concentration
parameters of the component DPs, and the base measure being a weighted mixture
of the corresponding mixing base measures.

######  Lemma 2.

Let $\bm{\alpha}_{1},\bm{\alpha}_{2},\dots,\bm{\alpha}_{L}$ be $k$-dimensional
vectors where $\bm{\alpha}_{i}=(\alpha_{i1},\dots,\alpha_{ik})$ with
$\alpha_{ij}>0\hskip 5.69046pt\forall\hskip 5.69046ptj=1,2,\dots,k$,
$i=1,2,\dots,L$. Let $\bm{X}_{1},\bm{X}_{2},\dots,\bm{X}_{L}$ be independent
$k$-dimensional random vectors distributed as Dirichlet distribution with
parameters $\bm{\alpha}_{1},\bm{\alpha}_{2},\dots,\bm{\alpha}_{L}$,
respectively. Let $\alpha_{i\cdot}=\sum_{j=1}^{k}\alpha_{ij},\ \
i=1,2,\dots,L$. Let $\bm{\pi}=\left(\pi_{1},\pi_{2},\dots,\pi_{L}\right)$ be
independent of $\bm{X}_{1},\bm{X}_{2},\dots,\bm{X}_{L}$ and have a Dirichlet
distribution
$Dir\left(\alpha_{1\cdot},\alpha_{2\cdot},\dots,\alpha_{L\cdot}\right)$. Then
the distribution of $\sum_{i=1}^{L}\pi_{i}\bm{X}_{i}$ is the Dirichlet
distribution with parameter $\sum_{i=1}^{L}\bm{\alpha}_{i}$.

The proof is provided in Section A of the Supplementary Materials.

##  4 Representations of the Graphical Dirichlet Process

In this section, we characterize the proposed GDP through (i) the hypergraph
representation, (ii) the stick-breaking representation, (iii) the restaurant-
type process representation, and (iv) the limit of finite mixture
representation.

###  4.1 The hypergraph representation

The GDP, along with the hyperpriors on the concentration parameters and
mixture weights, can be represented hierarchically as,

$\displaystyle\begin{aligned} &\alpha_{1}^{(0)}\mid\alpha_{0}\sim
Gamma\left(\alpha_{0},1\right),\\\ &G_{1}^{(0)}\mid\alpha_{1}^{(0)},G_{0}\sim
DP\left(\alpha_{1}^{(0)},G_{0}\right),\\\
&\alpha_{j}^{(k)}\mid\\{\alpha_{l}^{(k-1)}:l\in an^{(k,1)}(j)\\}\sim
Gamma\left(\sum_{l\in an^{(k,1)}(j)}\alpha_{l}^{(k-1)},1\right),\\\
&\\{\pi_{jl}^{(k)}\,:\,l\in an^{(k,1)}(j)\\}\mid\\{\alpha_{l}^{(k-1)}:l\in
an^{(k,1)}(j)\\}\sim Dir\left(\\{\alpha_{l}^{(k-1)}:l\in
an^{(k,1)}(j)\\}\right),\\\
&G_{j}^{(k)}\mid\alpha_{j}^{(k)},\\{G_{l}^{(k-1)}:l\in an^{(k,1)}(j)\\}\sim
DP\left(\alpha_{j}^{(k)},\sum_{l\in
an^{(k,1)}(j)}\pi_{jl}^{(k)}G_{l}^{(k-1)}\right),\end{aligned}$ (11)

for $j=1,2,\dots,l_{k}$ and $k=1,\dots,K$.

The hyperparameters of the GDP consist of the base probability measure $G_{0}$
and the concentration parameter $\alpha_{0}$. The probability measure
$G_{1}^{(0)}$ of the root node varies around the base measure $G_{0}$ with the
amount of variability governed by $\alpha_{1}^{(0)}$, which in turn is
governed by the hyperparameter $\alpha_{0}$. We now present a novel hypergraph
representation of GDP, which simplifies the graph-based distribution. The
representation follows from the gamma-DAG distribution on the concentration
parameters and standard properties of Dirichlet distribution.

######  Theorem 1 (Hypergraph Representation).

Consider a DAG $D$ that has $K$ layers and $l_{k}$ distinct nodes in layer $k$
for $k=1,\dots,K$. Under model (LABEL:gdp_equations), the distribution of the
random measure $G_{j}^{(k)}$ of node $j$ in layer $k$ of DAG $D$ can be
equivalently represented as,

$\displaystyle G_{j}^{(k)}\mid\alpha_{j}^{(k)},H_{j}^{(k,k)}\sim
DP\left(\alpha_{j}^{(k)},H_{j}^{(k,k)}\right),$ $\displaystyle
H_{j}^{(k,k)}\mid\\{\alpha_{l}^{(k-1)}:l\in
an^{(k,1)}(j)\\},H_{j}^{(k,k-1)}\sim DP\left(\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)},H_{j}^{(k,k-1)}\right),$ $\displaystyle
H_{j}^{(k,k-1)}\mid\\{\alpha_{l}^{(k-2)}:l\in
an^{(k,2)}(j)\\},H_{j}^{(k,k-2)}\sim DP\left(\sum_{l\in
an^{(k,2)}(j)}\alpha_{l}^{(k-2)},H_{j}^{(k,k-2)}\right),$
$\displaystyle\vdots$ $\displaystyle H_{j}^{(k,2)}\mid\\{\alpha_{l}^{(1)}:l\in
an^{(k,k-1)}(j)\\},G_{1}^{(0)}\sim DP\left(\sum_{l\in
an^{(k,k-1)}(j)}\alpha_{l}^{(1)},G_{1}^{(0)}\right).$

The proof is provided in the Appendix A. In words, Theorem 1 essentially
states the following. The distribution of $G_{j}^{(k)}$ is a DP with a hidden
base measure $H_{j}^{(k,k)}$ and the concentration parameter
$\alpha_{j}^{(k)}$. The hidden base measure $H_{j}^{(k,k)}$, in turn, is again
a DP with base measure $H_{j}^{(k,k-1)}$ and concentration parameter being the
sum of the concentration parameters of the generation-$1$ ancestors of
$G_{j}^{(k)}$. Recursively, the hidden base measure $H_{j}^{(k,k-1)}$ is a DP
with base measure $H_{j}^{(k,k-2)}$ and the concentration parameter being the
sum of the concentration parameters of the generation-$2$ ancestors. This
distributional pattern continues in a hierarchical fashion. Through $k-1$
hidden base measures, any node in layer $k$ can be seen to depend on the root
node $G_{1}^{(0)}$ through its ancestral relationships. We call the
representation of GDP in Theorem 1 as the hypergraph representation because
one can view $H_{j}^{(k,k-a)}$ for $a=0,\dots,k-2$ as a hypernode that
contains all the sufficient information from generation-$(a+1)$ ancestors of
$G_{j}^{(k)}$. We provide in Figure 3 an illustrative example of the
hypergraph representation showing how the hypernodes contain all the ancestral
information. From Figure 3(a), we can see that the distribution of $G_{6}$
depends on the distribution of its parents, $G_{2}$ and $G_{3}$. We refer to
$H_{2}$, consisting of $\\{G_{2},G_{3}\\}$, as a hypernode. Hypernode $H_{2}$
contains all the information about the parents of $G_{6}$. Loosely speaking,
the information of the root node $G_{1}$ (e.g., its atoms) is passed to
$G_{6}$ through $H_{2}$. Similarly, $H_{3}$, being the hypernode of
$\\{G_{3},G_{4}\\}$, contains all the information about $G_{7}$ from its
parent nodes allowing the flow of information from the root node (see Figure
3(b)). For node $G_{8}$, we have two levels of hypernodes – $H_{4}$ denotes
the first layer and consists of the parents of $G_{8}$, and $H^{*}$ denotes
the second layer and consists of generation-$2$ ancestors of $G_{8}$. Thus,
hypernodes $H_{4}$ and $H^{*}$ carry all the information from the root node
$G_{1}$ to $G_{8}$ as illustrated in Figure 3(c).

$G_{0}$$G_{1}$$G_{3}$$G_{2}$$H_{2}$$G_{4}$$G_{5}$$G_{6}$$G_{7}$$G_{8}$ (a)

$G_{0}$$G_{1}$$G_{3}$$G_{2}$$G_{4}$$H_{3}$$G_{5}$$G_{6}$$G_{7}$$G_{8}$ (b)

$G_{0}$$G_{1}$$G_{3}$$G_{2}$$G_{4}$$H^{*}$$G_{5}$$G_{6}$$G_{7}$$H_{4}$$G_{8}$
(c)

Figure 3: Illustration of hypernodes (represented by dashed ovals) of the DAG
for our motivational problem. (a) Hypernode $H_{2}$ consists of the
generation-1 ancestors (i.e., $G_{2}$ and $G_{3}$) of node $G_{6}$. (b)
Hypernode $H_{3}$ consists of the generation-1 ancestors (i.e., $G_{3}$ and
$G_{4}$) of node $G_{7}$. (c) Hypernode $H_{4}$ consists of the generation-1
ancestors (i.e., $G_{5}$, $G_{6}$, and $G_{7}$) of node $G_{8}$. Hypernode
$H^{*}$ consists of the generation-2 ancestors (i.e., $G_{2}$, $G_{3}$, and
$G_{4}$) of node $G_{8}$.

We will exploit this representation to derive the stick-breaking
representation and the limit of finite mixture representation of the proposed
GDP in the next subsections.

###  4.2 The stick-breaking representation

Given that the random measure $G_{1}^{(0)}$ of the root node is distributed as
a DP, it can be expressed using a stick-breaking representation,

$G_{1}^{(0)}=\sum_{l=1}^{\infty}\beta_{1l}^{(0)}\delta_{\phi_{l}},$ (12)

where $\phi_{l}\overset{iid}{\sim}G_{0}$ and
$\bm{\beta}_{1}^{(0)}=\left(\beta_{1l}^{(0)}\right)_{l=1}^{\infty}\sim\text{GEM}\left(\alpha_{1}^{(0)}\right)$
are mutually independent. We interpret $\bm{\beta}_{1}^{(0)}$ as a probability
measure on the positive integers. Since $G_{1}^{(0)}$ has support at the atoms
$\bm{\phi}=(\phi_{l})_{l=1}^{\infty}$, each $G_{j}^{(k)}$ necessarily has
support at these atoms as well and hence can be expressed as,

$G_{j}^{(k)}=\sum_{l=1}^{\infty}\beta_{jl}^{(k)}\delta_{\phi_{l}}.$ (13)

As with Theorem 1, the stick-breaking weights depend hierarchically on a set
of hidden weights. Letting
$\bm{\beta}_{j}^{(k)}=\left(\beta_{jl}^{(k)}\right)_{l=1}^{\infty}$ be the
stick-breaking weights for node $j$ in layer $k$ of DAG $D$ and letting
$\bm{\nu}_{j}^{(k,m)}=\left(\nu_{jl}^{(k,m)}\right)_{l=1}^{\infty},m=2,\dots,k$
be their hidden weights, we have the following corollary.

######  Corollary 1.1 (Stick-Breaking Representation).

Consider a DAG $D$ that has $K$ layers and $l_{k}$ distinct nodes in layer $k$
for $k=1,\dots,K$. The stick-breaking weights $\bm{\beta}_{j}^{(k)}$ of node
$j$ at layer $k$ of DAG $D$ can be represented as

$\displaystyle\bm{\beta}_{j}^{(k)}\mid\alpha_{j}^{(k)},\bm{\nu}_{j}^{(k,k)}\sim
DP\left(\alpha_{j}^{(k)},\bm{\nu}_{j}^{(k,k)}\right),$
$\displaystyle\bm{\nu}_{j}^{(k,k)}\mid\\{\alpha_{l}^{(k-1)}:l\in
an^{(k,1)}(j)\\},\bm{\nu}_{j}^{(k,k-1)}\sim DP\left(\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)},\bm{\nu}_{j}^{(k,k-1)}\right),$
$\displaystyle\bm{\nu}_{j}^{(k,k-1)}\mid\\{\alpha_{l}^{(k-2)}:l\in
an^{(k,2)}(j)\\},\bm{\nu}_{j}^{(k,k-2)}\sim DP\left(\sum_{l\in
an^{(k,2)}(j)}\alpha_{l}^{(k-2)},\bm{\nu}_{j}^{(k,k-2)}\right),$
$\displaystyle\vdots$
$\displaystyle\bm{\nu}_{j}^{(k,2)}\mid\\{\alpha_{l}^{(1)}:l\in
an^{(k,k-1)}(j)\\},\bm{\beta}_{1}^{(0)}\sim DP\left(\sum_{l\in
an^{(k,k-1)}(j)}\alpha_{l}^{(1)},\bm{\beta}_{1}^{(0)}\right).$

The proof of this corollary directly follows from the hypergraph
representation of Theorem 1 and is hence omitted. We call this representation
_the stick-breaking representation_ where $\bm{\nu}_{j}^{(k,k)}$ is
interpreted as a hidden probability measure on the set of positive integers
corresponding to the first hidden layer. Each hidden layer of stick-breaking
weights depend hierarchically on its previous hidden layer, denoted by
$\bm{\nu}_{j}^{(k,k-1)}$, $\bm{\nu}_{j}^{(k,k-2)}$, and so on, and finally on
the weights $\bm{\beta}_{1}^{(0)}$ of the root node.

###  4.3 The family-owned restaurant process representation

DP and HDP have the well-known Chinese restaurant process and franchise
representations. Here, we provide a culinary analog for the proposed GDP. We
refer to this process as the _family-owned restaurant process_ as it is
customary to use familial relationships to describe the relationships between
nodes in a DAG. The metaphor is as follows. An original restaurant is opened
by the ancestor of a family (the root node), which serves some dishes from a
global menu containing an infinite number of dishes. The descendants of the
ancestor open their own respective restaurants, which serve some of the dishes
already being served in the restaurants owned by their parents and possibly
some new dishes from the global menu. At each table of the original
restaurant, one dish is ordered from the menu by the first customer occupying
the table, and the dish is shared by all the other customers who sit at that
table. Any subsequent customer may either join an occupied table and share the
dish being served at that table or open a new table with a new dish from the
menu. In restaurants other than the original restaurant, however, the first
customer might choose to select a dish being served at one of the tables of
its parent restaurant or order a new dish from the menu. Since the hypergraph
representation of GDP involves hypernodes with hidden probability measures, we
introduce a notation for the number of tables serving a dish in any restaurant
and demarcate them with the notation for the number of tables serving the dish
in the hypernodes, which we refer to as hyper-restaurants.

As before, assume that there are $K$ generations in the family and there are
$l_{k}$ different restaurants in generation $k$. The restaurants correspond to
the nodes of DAG $D$. The customers coming in restaurant $j$ of generation $k$
correspond to parameters $\theta_{ji}^{(k)}$. Let
$\phi_{1},\phi_{2},\dots,\phi_{L}$ denote i.i.d. random variables distributed
according to the base distribution $G_{0}$, which are dishes from the global
menu. To maintain a count of customers and tables, we introduce two notations.
We use the notation $n_{jt}^{(k)}$ to denote the number of customers at table
$t$ in the restaurant $j$ of generation $k$ and the notation $m_{jl}^{(k)}$ to
denote the number of tables in the restaurant $j$ of generation $k$ that serve
dish $l$. Marginal counts are represented by dots at the appropriate indices.
For example, $m_{j\cdot}^{(k)}$ denotes the count of all the tables
(regardless of what dishes being served) in the restaurant $j$ of generation
$k$. We introduce the notation $\psi_{jt}^{(k,k)}$ to denote the dish served
at table $t$ in restaurant $j$ of generation $k$, chosen from the
corresponding layer-$1$ hyper-restaurant ($H_{j}^{(k,k)})$.

We integrate out random measures
$\left\\{G_{j}^{(k)},H_{j}^{(k,k)},H_{j}^{(k,k-1)},\dots,G_{1}^{(0)}\right\\}$
sequentially. First, we find the conditional distribution of
$\theta_{ji}^{(k)}$ given
$\theta_{j1}^{(k)},\theta_{j2}^{(k)},\dots,\theta_{j,i-1}^{(k)},\alpha_{j}^{(k)}$,
and $H_{j}^{(k,k)}$ with $G_{j}^{(k)}$ integrated out,

$\theta_{ji}^{(k)}\mid\theta_{j1}^{(k)},\theta_{j2}^{(k)},\dots,\theta_{j,i-1}^{(k)},\alpha_{j}^{(k)},H_{j}^{(k,k)}\sim\sum_{t=1}^{m_{j\cdot}^{(k)}}\frac{n_{jt}^{(k)}}{i-1+\alpha_{j}^{(k)}}\delta_{\psi_{jt}^{(k,k)}}+\frac{\alpha_{j}^{(k)}}{i-1+\alpha_{j}^{(k)}}H_{j}^{(k,k)},$
(14)

We let $\psi_{jt}^{(k,k-1)}$ to denote the dish served at table $t$ in the
layer-$1$ hyper-restaurant corresponding to restaurant $j$ of generation $k$,
chosen from the dishes served in the layer-$2$ hyper-restaurants
($H_{j}^{(k,k-1)})$. Integrating out the hidden measure from the current layer
$H_{j}^{(k,k)}$, the conditional distribution of $\psi_{jt}^{(k,k)}$ given
$\psi_{j1}^{(k,k-1)},\psi_{j2}^{(k,k-1)},\dots,\psi_{j1}^{(k,k)},\dots,\psi_{j,t-1}^{(k,k)},\\{\alpha_{l}^{(k-1)}:l\in
an^{(k,1)}(j)\\}$, and the hidden measure from the previous layer,
$H_{j}^{(k,k-1)}$ is given by,

$\psi_{jt}^{(k,k)}\mid\psi_{j1}^{(k,k-1)},\psi_{j2}^{(k,k-1)},\dots,\psi_{j1}^{(k,k)},\dots,\psi_{j,t-1}^{(k,k)},\\{\alpha_{l}^{(k-1)}:l\in
an^{(k,1)}(j)\\},H_{j}^{(k,k-1)}\\\
\sim{\sum}_{l=1}^{\scalebox{0.7}{$M_{j}^{{(k,1)}}$}}\frac{m_{jl}^{(k,1)}}{m_{j\cdot}^{(k,1)}+\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)}}{\delta}_{\psi_{jl}^{(k,k-1)}}+\frac{\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)}}{m_{j\cdot}^{(k,1)}+\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)}}H_{j}^{(k,k-1)},$ (15)

where the notation $m_{jl}^{(k,1)}$ denotes the number of tables in layer-1
hyper-restaurant, corresponding to restaurant $j$ of generation $k$ serving
the dish $l$. We denote by $M_{j}^{{(k,1)}}$ the number of dishes served in
the layer-$1$ hyper-restaurants and by $m_{j\cdot}^{(k,1)}$ the total number
of tables in the layer-$1$ hyper-restaurant, corresponding to the restaurant
$j$ of generation $k$. Similarly, integrating out the measure
$H_{j}^{(k,k-1)}$ and introducing the next layer of variables
$\psi_{jt}^{(k,k-2)}$, the conditional distribution of $\psi_{jt}^{(k,k-1)}$
given
$\psi_{j1}^{(k,k-2)},\psi_{j2}^{(k,k-2)},\dots,\psi_{j1}^{(k,k-1)},\dots,\psi_{j,t-1}^{(k,k-1)},\\{\alpha_{l}^{(k-2)}:l\in
an^{(k,2)}(j)\\}$, and the hidden measure from the previous layer
$H_{j}^{(k,k-2)}$ is given by,

$\psi_{jt}^{(k,k-1)}\mid\psi_{j1}^{(k,k-2)},\psi_{j2}^{(k,k-2)},\dots,\psi_{j1}^{(k,k-1)},\dots,\psi_{j,t-1}^{(k,k-1)},\\{\alpha_{l}^{(k-2)}:l\in
an^{(k,2)}(j)\\},H_{j}^{(k,k-2)}\\\
\sim{\sum}_{l=1}^{\scalebox{0.7}{$M_{j}^{{(k,2)}}$}}\frac{m_{jl}^{(k,2)}}{m_{j\cdot}^{(k,2)}+\sum_{l\in
an^{(k,2)}(j)}\alpha_{l}^{(k-2)}}{\delta}_{\psi_{jl}^{(k,k-2)}}+\frac{\sum_{l\in
an^{(k,2)}(j)}\alpha_{l}^{(k-2)}}{m_{j\cdot}^{(k,2)}+\sum_{l\in
an^{(k,2)}(j)}\alpha_{l}^{(k-2)}}H_{j}^{(k,k-2)}.$ (16)

As in the stick-breaking representation, we can recursively integrate out
hidden measures and eventually arrive at the conditional distribution of
$\psi_{jt}^{(k,2)}$ given
$\psi_{j1}^{(0)},\psi_{j2}^{(0)},\dots,\psi_{j1}^{(k,2)},\dots,\psi_{j,t-1}^{(k,2)},\\{\alpha_{l}^{(1)}:l\in
an^{(k,k-1)}(j)\\}$, and the probability measure of the root node
$G_{1}^{(0)}$,

$\psi_{jt}^{(k,2)}\mid\psi_{j1}^{(0)},\psi_{j2}^{(0)},\dots,\psi_{j1}^{(k,2)},\dots,\psi_{j,t-1}^{(k,2)},\\{\alpha_{l}^{(1)}:l\in
an^{(k,k-1)}(j)\\},G_{1}^{(0)}\\\
\sim{\sum}_{l=1}^{M_{j}^{(k,k-1)}}\frac{m_{jl}^{(k,k-1)}}{m_{j\cdot}^{(k,k-1)}+\sum_{l\in
an^{(k,k-1)}(j)}\alpha_{l}^{(1)}}{\delta}_{\psi_{jl}^{(0)}}+\frac{\sum_{l\in
an^{(k,k-1)}(j)}\alpha_{l}^{(1)}}{m_{j\cdot}^{(1)}+\sum_{l\in
an^{(k,k-1)}(j)}\alpha_{l}^{(1)}}G_{1}^{(0)},$ (17)

and the conditional distribution of $\psi_{jt}^{(0)}$ given
$\psi_{j1}^{(0)},\dots,\psi_{j,t-1}^{(0)},\alpha_{1}^{(0)}$, and the base
measure $G_{0}$,

$\psi_{jt}^{(0)}\mid\psi_{j1}^{(0)},\dots,\psi_{j,t-1}^{(0)},\alpha_{1}^{(0)},G_{0}\sim{\sum}_{l=1}^{L}\frac{m_{l}^{(0)}}{m_{\cdot}^{(0)}+\alpha_{1}^{(0)}}{\delta}_{\phi_{l}}+\frac{\alpha_{1}^{(0)}}{m_{\cdot}^{(0)}+\alpha_{1}^{(0)}}G_{0},$
(18)

where $m_{l}^{(0)}$ denotes the number of tables in the original restaurant
serving dish $l$ and $m_{\cdot}^{(0)}$ denotes the total number of tables in
the original restaurant. Note that (18) corresponds to the case where the root
node is hidden (the same as in HDP). When the root node is not hidden, a
similar formula can be derived, which is omitted for simplicity.

###  4.4 The infinite limit of finite mixture model

The GDP mixture model can be derived as the infinite limit of a finite mixture
model. Let us denote the observations and the mixture component indicator from
node $j$ in layer $k$ of DAG $D$ by $x_{ji}^{(k)}$ and $z_{ji}^{(k)}$,
respectively. Suppose $\bm{\beta}_{1}^{(0)}$ is the vector of mixing weights
for the root node. Denoting by $\bm{\beta}_{j}^{(k)}$ the mixing weights of
node $j$ in layer $k$ and by $\nu_{j}^{(k,m)}$ the corresponding mixing
weights for the hidden layer $m$, with $m=2,\dots,k$, we consider a finite
mixture version of the proposed GDP,

$\displaystyle\bm{\beta}_{1}^{(0)}\mid\alpha_{1}^{(0)}$ $\displaystyle\sim
Dir\left(\alpha_{1}^{(0)}/L,\dots,\alpha_{1}^{(0)}/L\right),$
$\displaystyle\bm{\nu}_{j}^{(k,2)}\mid\\{\alpha_{l}^{(1)}:l\in
an^{(k,k-1)}(j)\\},\bm{\beta}_{1}^{(0)}$ $\displaystyle\sim
Dir\left(\sum_{l\in
an^{(k,k-1)}(j)}\alpha_{l}^{(1)}\left(\beta_{11}^{(0)},\dots,\beta_{1L}^{(0)}\right)\right),$
$\displaystyle\vdots$
$\displaystyle\bm{\nu}_{j}^{(k,k)}\mid\\{\alpha_{l}^{(k,k-1)}:l\in
an^{(k,1)}(j)\\},\bm{\nu}_{j}^{(k,k-1)}$ $\displaystyle\sim
Dir\left(\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)}\left(\nu_{j1}^{(k,k-1)},\dots,\nu_{jL}^{(k,k-1)}\right)\right),$
$\displaystyle\bm{\beta}_{j}^{(k)}\mid\alpha_{j}^{(k)},\bm{\nu}_{j}^{(k,k)}$
$\displaystyle\sim
Dir\left(\alpha_{j}^{(k)}\left(\nu_{j1}^{(k,k)},\dots,\nu_{jL}^{(k,k)}\right)\right),$
$\displaystyle\phi_{l}\mid G_{0}$ $\displaystyle\sim G_{0},$ $\displaystyle
z_{ji}^{(k)}\mid\bm{\beta}_{j}^{(k)}$ $\displaystyle\sim\bm{\beta}_{j}^{(k)},$
$\displaystyle x_{ji}^{(k)}\mid z_{ji}^{(k)},\left(\phi_{l}\right)_{l=1}^{L}$
$\displaystyle\sim F\left(\phi_{z_{ji}^{(k)}}\right).$ (19)

The distribution of this finite mixture model approaches the GDP mixture model
as $L\rightarrow\infty$. Refer to Section B of the Supplementary Materials for
the proof. Based on this finite mixture model approximation with a large
enough truncation level $L$, we develop an efficient posterior inference
procedure of our model using a Metropolis-within-blocked-Gibbs sampler with a
specialized proposal (Director et al., 2017); see Section C of the
Supplementary Materials for details.

##  5 Simulations

Our simulations are designed to mimic the motivating application where we have
8 experimental groups, whose relationships are represented by the DAG in
Figure 4.

$1$$3$$2$$4$$5$$6$$7$$8$ Figure 4: The DAG of experimental groups.

We generated data within each of the $8$ groups from a four-component mixture
of bivariate Gaussian distributions with different covariance matrices for
each group. We drew the DP concentration parameters $\alpha_{j}$’s for the
different groups from their prior distribution (LABEL:eq:alpha_prior)
respecting the DAG in Figure 4 with $\alpha_{0}=5$. The weights of the finite
mixture model corresponding to the different groups were drawn using (4.4) and
the same DAG. The true cluster indicators of each group were sampled from a
multinomial distribution with probabilities equal to the mixture weights.
Using these true cluster indices for each group, samples were drawn from the
Gaussian distribution with the cluster-specific mean and group-specific
covariance matrix, given in Tables S2 and S3, respectively, in Section D of
the Supplementary Materials. Refer to the same section in the Supplementary
Materials for more details on our simulation strategy. In our Gibbs sampler,
the truncation level of the finite mixture model was set to $L=10$, and the
base measure for GDP, $G_{0}$, was specified as the normal-inverse-Wishart
distribution, $\mathcal{NIW}(\bm{0},0.01,\mathbb{I}_{2},2)$. Upon the
completion of the Gibbs sampler, the clusters were estimated by using the
least squares criterion (Dahl, 2006), and they were compared with the true
cluster labels for evaluation. We considered various sample sizes in each
group, which are summarized in Table 1. In all cases, we ran 15,000 iterations
of our Gibbs sampler and discarded the first 5,000 samples as burn-in.

Sample sizes | Groups
---|---
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
small | 40 | 30 | 30 | 35 | 25 | 30 | 25 | 30
moderate | 80 | 70 | 70 | 75 | 83 | 88 | 92 | 88
large | 150 | 160 | 180 | 170 | 155 | 175 | 185 | 145
unbalanced | 350 | 30 | 40 | 45 | 25 | 25 | 35 | 35
Table 1: The sample sizes for the different groups that were used to simulate
the data.

(a) Small sample size in each group

(b) Unbalanced sample sizes between groups

Figure 5: Clustering performance of GDP for different sample sizes. The colors
indicate the estimated clusters by GDP. Adjusted Rand index is reported at the
top of each panel.

The clustering results of GDP for small and unbalanced sample sizes are
visualized in Figure 5. The remaining clustering plots are shown in
Supplementary D. Across different sample sizes, the proposed GDP was able to
identify the clusters within each group with very good accuracy and was able
to link clusters across non-exchangeable groups.

We also looked at the clustering performance of GDP under a more difficult
scenario. The simulation details and clustering results are shown in
Supplementary D. Since HDP is a special case of the proposed GDP, we compared
the two methods for this difficult scenario. We also compared the clustering
performance of GDP with k-means, a widely used non-Bayesian clustering
technique. The number of clusters in k-means was taken to be the truncation
level of our GDP. All simulations were replicated $50$ times.

GDP significantly outperformed both HDP and k-means. For example, the boxplots
of Adjusted Rand indices (Hubert and Arabie, 1985; higher is better) for the
different methods are shown in Figure 6. It is evident that the Adjusted Rand
indices of GDP were almost uniformly higher than those of HDP because HDP was
not able to handle non-exchangeable groups. Similarly, the higher Adjusted
Rand indices of GDP indicated its superior clustering performance over the
k-means algorithm. Moreover, k-means algorithm does not allow sharing of
relevant clusters across the groups.

Figure 6: The boxplots of the Adjusted Rand indices for GDP, HDP, and k-means
for all sample sizes.

##  6 Real Data Analysis

With the advancement of next-generation sequencing techniques in recent years,
it is now possible to molecularly characterize individual cells, which may
provide valuable insights into complex biological systems, ranging from cancer
genomics to diverse microbial communities (Hwang et al., 2018). Colorectal
cancer is the third most common type of cancer after breast and lung cancers.
It is known that the mutation of tumor-suppressor gene Apc is an initial step
in most colorectal tumors (Morin et al., 1997). In addition, numerous studies
have been conducted to understand the effect of high-fat vs low-fat diet on
gene expressions (Jump and Clarke, 1999; Bouchard-Mercier et al., 2013; Fan et
al., 2020). We are motivated by a study that aimed to investigate how diet,
genotype, and treatment with a new cancer prevention drug (AdipoRon) against
placebo interacted to influence the expression of genes in intestinal crypt
and tumor cells. The experiments started from a baseline group where the mice
were genetically wild-type, fed with a normal diet, and treated with placebo.
Then to understand the main effects of genotype, diet, and cancer treatment on
stem cell gene expression, the experimenters introduced three new groups of
mice, each differing from the baseline group by exactly one factor (Apc knock-
out, high-fat diet, or new cancer treatment AdipoRon). To determine the two-
way interaction effects, three additional groups of mice were studied, each of
which differed from the baseline group by two factors (e.g., mice with Apc
knock-out, high-fat diet, and placebo). Lastly, for a three-way interaction,
the experimenters introduced the eighth group of mice with Apc knock-out, a
high-fat diet, and the new treatment AdipoRon. By design, these 8 experimental
groups are non-exchangeable and their relationships can be delineated by the
DAG in Figure 4. The goal of this analysis is to identify potential intestinal
molecular subtypes within each experimental group while allowing information
to be shared across these non-exchangeable groups with the proposed GDP model.
For illustration, we randomly sampled 100 cells from each of the eight groups.
The scRNA-seq data were pre-processed following standard procedure as outlined
by Hao et al. (2021) using the R package Seurat. The data was log-normalized
and scaled such that the mean expression across cells was $0$ and the variance
across cells was $1$. As a common practice in single-cell data analyses, the
uniform manifold approximation and projection (UMAP) (McInnes et al., 2018)
was used to reduce the data to two dimensions. We considered the same
truncation level, $L=10$, and the same base probability measure, $G_{0}$, as
in the simulations. We ran four parallel chains of the Gibbs sampler for
$40,000$ iterations. To monitor the convergence of the sampler, we drew the
traceplots of the log-likelihood for each of the four chains, after discarding
the initial $25,000$ samples and thinning the samples by a factor of $15$,
which indicated no lack of convergence of our sampler. We pooled the Monte
Carlo samples across different chains for posterior inference. We compared the
clustering performance with that obtained from HDP on the same data.

(a) Clustering plot for different groups by GDP.

(b) Clustering plot for different groups by HDP.

Figure 7: Clustering of the group-specific single-cell data whose dimensions
are reduced to 2 by UMAP by (a) GDP and (b) HDP.

The estimated clusters from GDP and HDP are shown in Figures 7(a) and 7(b),
respectively. As shown in Table S1 in Section C of the Supplementary
Materials, group 1 is the wild-type group receiving the placebo and a normal
diet. Each of group 2, 3, and 4 are obtained from group 1 by changing the
three factors one at a time, and hence shares some similar clusters with group
1. Group 4 is similar to the baseline group 1 but with the Apc gene knocked
out. The corresponding clustering plot (Figure 7(a)) of GDP indicates that the
Apc knock-out group seems to exhibit more heterogeneity of cells (suggesting
possibly new cellular subtypes) as compared to the wild-type group. Group 5 is
the Apc knock-out group receiving a high-fat diet and the placebo. The
clustering plot shows some resemblance with its parent groups (groups 2 and 4)
but with the absence of some parental clusters. Groups 6 and 7 show similar
clustering patterns, indicating possibly similar impact of changing the
corresponding factors from their parent groups. Groups 7 and 8 correspond to
the Apc knock-out group receiving the new treatment and fed with a normal and
high-fat diet, respectively. It can be seen that the high-fat diet group
appears to have greater molecular heterogeneity than the normal diet group.
The Figure 7(b), on the other hand, clearly shows that HDP fails to capture
meaningful clusters across the non-exchangeable groups, i.e., some points that
seemingly belong to the same cluster are assigned different labels across
groups. To quantify the difference between GDP and HDP, we computed several
internal clustering validation measures; see Liu et al., 2010 for a review of
several such measures. Table 2 compares the Calinski-Harabasz, Davies–Bouldin,
and Silhouette Index between GDP and HDP. Clearly, all of them indicate the
superior clustering performance of GDP over HDP.

| Calinski-Harabasz Index (CHI) | Davies–Bouldin Index (DBI) | Silhouette Index (SI)
---|---|---|---
GDP | 493.567 | 1.131 | 0.251
HDP | 340.981 | 1.599 | 0.095
Table 2: Different measures of internal clustering for GDP and HDP. Higher
values of CHI and SI indicates better clustering. Lower values of DBI indicate
well separated clusters.

##  7 Discussion

We have introduced the GDP as a graph-based stochastic process for modeling
dependent random measures that are linked by a DAG. We have also introduced
the corresponding infinite mixture model and presented how the GDP mixture
model can be used for clustering grouped data with non-exchangeable groups. We
provided different representations of the GDP including a novel hypergraph
representation of the original process. The posterior inference was relatively
straightforward. We illustrated our method using both simulations and an
application to a real grouped scRNA-seq dataset.

There are a few possible future directions for this work. First, it may be
possible to replace the DAG in our GDP with an undirected or chain graph. The
challenge is to define the joint distribution over a set of random measures
given the graph where the convenient DAG factorization no longer applies.
Second, it may also be possible to learn the DAG structure instead of assuming
it is known, which may require independent realizations of the GDP. Third, it
will be interesting to extend the nested DP and other Bayesian nonparametric
priors to grouped data with non-exchangeable groups.

##  Appendix A Proof of the Hypergraph Representation

We prove Theorem 1 (the hypergraph representation of the proposed GDP) of the
main manuscript in the case of our motivational problem where we have 8
groups. Note that the proof for any general DAG follows in a similar fashion
by repeated application of the two lemmas in Section 3.3 of the main
manuscript and properties of gamma and Dirichlet distributions, which,
however, requires more involved bookkeeping of the corresponding random
distributions and hence is omitted. Our proof also illustrates how the random
distribution of any particular node of the DAG is related to the root node
through a number of hidden random measures, which shows the clustering
property of our model. In our motivating example, each group corresponds to a
combination of treatment, diet, and genotype, as summarized in Table S1 in
Section C of the Supplementary Materials. The underlying DAG for the problem
is given in Figure 4 of the main manuscript where group 1 is the root node,
groups 2-4 are the layer-1 nodes, groups 5-7 are the layer-2 nodes, and group
8 is the layer-3 node. For ease of notation, instead of using $G_{1}^{(0)}$
and $\alpha_{1}^{(0)}$ to denote the random measure and the concentration
parameter of the root node, we use simply $G_{1}$ and $\alpha_{1}$ instead;
similarly for all the other nodes. Using these simplified notations, Figures
8(a) and 8(b) show the relationships among the group-specific random measures
and concentration parameters according to Figure 4 of the main manuscript.

$G_{0}$$G_{1}$$G_{3}$$G_{2}$$G_{4}$$G_{5}$$G_{6}$$G_{7}$$G_{8}$ (a)

$\alpha_{0}$$\alpha_{1}$$\alpha_{3}$$\alpha_{2}$$\alpha_{4}$$\alpha_{5}$$\alpha_{6}$$\alpha_{7}$$\alpha_{8}$
(b)

Figure 8: The DAG of the (a) random measures $G_{j}$’s and (b) concentration
parameters $\alpha_{j}$’s.

The proposed GDP mixture model for this problem is given hierarchically as,

$\displaystyle\begin{aligned} \alpha_{1}\mid\alpha_{0}&\sim
Gamma(\alpha_{0},1),&G_{1}\mid G_{0},\alpha_{1}&\sim DP(\alpha_{1},G_{0}),\\\
\alpha_{2}\mid\alpha_{1}&\sim Gamma(\alpha_{1},1),&G_{2}\mid
G_{1},\alpha_{2}&\sim DP(\alpha_{2},G_{1}),\\\ \alpha_{3}\mid\alpha_{1}&\sim
Gamma(\alpha_{1},1),&G_{3}\mid G_{1},\alpha_{3}&\sim DP(\alpha_{3},G_{1}),\\\
\alpha_{4}\mid\alpha_{1}&\sim Gamma(\alpha_{1},1),&G_{4}\mid
G_{1},\alpha_{4}&\sim DP(\alpha_{4},G_{1}),\\\
\alpha_{5}\mid\alpha_{2},\alpha_{4}&\sim
Gamma(\alpha_{2}+\alpha_{4},1),&G_{5}\mid G_{2},G_{4},\alpha_{5},\pi_{1}&\sim
DP(\alpha_{5},\pi_{1}G_{2}+(1-\pi_{1})G_{4}),\\\
&&\pi_{1}\mid\alpha_{2},\alpha_{4}&\sim Beta(\alpha_{2},\alpha_{4}),\\\
\alpha_{6}\mid\alpha_{2},\alpha_{3}&\sim
Gamma(\alpha_{2}+\alpha_{3},1),&G_{6}\mid G_{2},G_{3},\alpha_{6},\pi_{2}&\sim
DP(\alpha_{6},\pi_{2}G_{2}+(1-\pi_{2})G_{3}),\\\
&&\pi_{2}\mid\alpha_{2},\alpha_{2}&\sim Beta(\alpha_{2},\alpha_{3}),\\\
\alpha_{7}\mid\alpha_{3},\alpha_{4}&\sim
Gamma(\alpha_{3}+\alpha_{4},1),&G_{7}\mid G_{3},G_{4},\alpha_{7},\pi_{3}&\sim
DP(\alpha_{7},\pi_{3}G_{3}+(1-\pi_{3})G_{4}),\\\
&&\pi_{3}\mid\alpha_{3},\alpha_{4}&\sim Beta(\alpha_{3},\alpha_{4}),\\\
\alpha_{8}\mid\alpha_{5},\alpha_{6},\alpha_{7}&\sim
Gamma(\alpha_{5}+\alpha_{6}+\alpha_{7},1),&G_{8}\mid
G_{5},G_{6},G_{7},\alpha_{8},\bm{\gamma}&\sim
DP(\alpha_{8},\gamma_{1}G_{5}+\gamma_{2}G_{6}+\gamma_{3}G_{7}),\\\
&&\bm{\gamma}=(\gamma_{1},\gamma_{2},\gamma_{3})\mid\alpha_{5},\alpha_{6},\alpha_{7}&\sim
Dir(\alpha_{5},\alpha_{6},\alpha_{7}),\\\ \theta_{ji}\mid
G_{j}&\overset{ind}{\sim}G_{j},\\\
x_{ji}\mid\theta_{ji},G_{j},G_{0}&\overset{ind}{\sim}F(\theta_{ji}),&&i=1,\dots,n_{j},\hskip
5.69046ptj=1,\dots,8.\end{aligned}$ (20)

Now, from Theorem 1, we have the following hypergraph representation, which we
are going to prove,

$\displaystyle\begin{aligned} \alpha_{1}\mid\alpha_{0}&\sim
Gamma(\alpha_{0},1),&G_{1}\mid G_{0},\alpha_{1}&\sim DP(\alpha_{1},G_{0}),\\\
\alpha_{2}\mid\alpha_{1}&\sim Gamma(\alpha_{1},1),&G_{2}\mid
G_{1},\alpha_{2}&\sim DP(\alpha_{2},G_{1}),\\\ \alpha_{3}\mid\alpha_{1}&\sim
Gamma(\alpha_{1},1),&G_{3}\mid G_{1},\alpha_{3}&\sim DP(\alpha_{3},G_{1}),\\\
\alpha_{4}\mid\alpha_{1}&\sim Gamma(\alpha_{1},1),&G_{4}\mid
G_{1},\alpha_{4}&\sim DP(\alpha_{4},G_{1}),\\\
\alpha_{5}\mid\alpha_{2},\alpha_{4}&\sim
Gamma(\alpha_{2}+\alpha_{4},1),&G_{5}\mid H_{1},\alpha_{5}&\sim
DP(\alpha_{5},H_{1}),\\\ &&H_{1}\mid\alpha_{2},\alpha_{4}&\sim
DP(\alpha_{2}+\alpha_{4},G_{1}),\\\ \alpha_{6}\mid\alpha_{2},\alpha_{3}&\sim
Gamma(\alpha_{2}+\alpha_{3},1),&G_{6}\mid H_{2},\alpha_{6}&\sim
DP(\alpha_{6},H_{2}),\\\ &&H_{2}\mid\alpha_{2},\alpha_{3}&\sim
DP(\alpha_{2}+\alpha_{3},G_{1}),\\\ \alpha_{7}\mid\alpha_{3},\alpha_{4}&\sim
Gamma(\alpha_{3}+\alpha_{4},1),&G_{7}\mid H_{3},\alpha_{7}&\sim
DP(\alpha_{7},H_{3}),\\\ &&H_{3}\mid\alpha_{3},\alpha_{4}&\sim
DP(\alpha_{3}+\alpha_{4},G_{1}),\\\
\alpha_{8}\mid\alpha_{5},\alpha_{6},\alpha_{7}&\sim
Gamma(\alpha_{5}+\alpha_{6}+\alpha_{7},1),&G_{8}\mid H_{4},\alpha_{8}&\sim
DP(\alpha_{8},H_{4}),\\\
&&H_{4}\mid\alpha_{5},\alpha_{6},\alpha_{7},H^{*}&\sim
DP(\alpha_{5}+\alpha_{6}+\alpha_{7},H^{*}),\\\
&&H^{*}\mid\alpha_{2},\alpha_{3},\alpha_{4},G_{1}&\sim
DP(2(\alpha_{2}+\alpha_{3}+\alpha_{4}),G_{1}),\\\ \theta_{ji}\mid
G_{j}&\overset{ind}{\sim}G_{j},\\\
x_{ji}\mid\theta_{ji},G_{j},G_{0}&\overset{ind}{\sim}F(\theta_{ji}),&i=1,\dots,n_{j},\hskip
5.69046ptj=1,\dots,8.\end{aligned}$ (21)

###### Proof.

Note that the random measures $G_{2},G_{3}$, and $G_{4}$ are the layer-1
nodes. Their relationships to the root node $G_{1}$ are the same as those in
an HDP. We shall consider the relationships of the random measures of the
layer-2 and layer-3 nodes (i.e., $G_{5},G_{6},G_{7}$, and $G_{8}$) to the root
node. Let $H_{1}=\pi_{1}G_{2}+(1-\pi_{1})G_{4}$ where $G_{2}\sim
DP(\alpha_{2},G_{1})$ and $G_{4}\sim DP(\alpha_{4},G_{1})$ independently. Let
$A_{1},A_{2},\dots,A_{r}$ be a finite measurable partition of the sample space
$\Theta$. Then by the definition of DP, we have

$\displaystyle\left(G_{2}(A_{1}),G_{2}(A_{2}),\dots,G_{2}(A_{r})\right)$
$\displaystyle\sim
Dir\left(\alpha_{2}G_{1}(A_{1}),\alpha_{2}G_{1}(A_{2}),\dots,\alpha_{2}G_{1}(A_{r})\right),$
$\displaystyle\left(G_{4}(A_{1}),G_{4}(A_{2}),\dots,G_{4}(A_{r})\right)$
$\displaystyle\sim
Dir\left(\alpha_{4}G_{1}(A_{1}),\alpha_{4}G_{1}(A_{2}),\dots,\alpha_{4}G_{1}(A_{r})\right),$

which are conditionally independent given $\alpha_{2},\alpha_{4}$ and $G_{1}$.
As $\pi_{1}\sim Beta(\alpha_{2},\alpha_{4})$ independently of $G_{2}$ and
$G_{4}$, using Lemma 1, we have that, given $\alpha_{2},\alpha_{4}$ and
$G_{1}$,

$\displaystyle\pi_{1}\left(G_{2}(A_{1}),\dots,G_{2}(A_{r})\right)+(1-\pi)\left(G_{4}(A_{1}),\dots,G_{4}(A_{r})\right)$
$\displaystyle\sim
Dir((\alpha_{2}+\alpha_{4})(G_{1}(A_{1}),\dots,G_{1}(A_{r})))$
$\displaystyle\Rightarrow(H_{1}(A_{1}),\dots,H_{1}(A_{r}))\mid\alpha_{2},\alpha_{4},G_{1}$
$\displaystyle\sim
Dir((\alpha_{2}+\alpha_{4})(G_{1}(A_{1}),\dots,G_{1}(A_{r})))$
$\displaystyle\Rightarrow H_{1}\mid\alpha_{2},\alpha_{4},G_{1}$
$\displaystyle\sim DP(\alpha_{2}+\alpha_{4},G_{1})$ (22)

Thus, we have

$\displaystyle G_{5}\mid H_{1},\alpha_{5}$ $\displaystyle\sim
DP(\alpha_{5},H_{1}),$ $\displaystyle H_{1}\mid\alpha_{2},\alpha_{4},G_{1}$
$\displaystyle\sim DP(\alpha_{2}+\alpha_{4},G_{1}).$ (23)

Similarly, the other layer-2 measures $G_{6}$ and $G_{7}$ have the following
representations:

$\displaystyle G_{6}\mid H_{2},\alpha_{6}$ $\displaystyle\sim
DP(\alpha_{6},H_{2}),$ $\displaystyle H_{2}\mid\alpha_{2},\alpha_{3},G_{1}$
$\displaystyle\sim DP(\alpha_{2}+\alpha_{3},G_{1}),$ (24)

and,

$\displaystyle G_{7}\mid H_{3},\alpha_{7}$ $\displaystyle\sim
DP(\alpha_{7},H_{3}),$ $\displaystyle H_{3}\mid\alpha_{3},\alpha_{4},G_{1}$
$\displaystyle\sim DP(\alpha_{3}+\alpha_{4},G_{1}),$ (25)

where $H_{2}=\pi_{2}G_{2}+(1-\pi_{2})G_{3}$ and
$H_{3}=\pi_{3}G_{3}+(1-\pi_{3})G_{4}$.

Let $H_{4}=\gamma_{1}G_{5}+\gamma_{2}G_{6}+\gamma_{3}G_{7}$ and
$\bm{\gamma}=(\gamma_{1},\gamma_{2},\gamma_{3})\sim
Dir(\alpha_{5},\alpha_{6},\alpha_{7})$. Since $G_{5}$, $G_{6}$, and $G_{7}$
are conditionally independent given $G_{2}$, $G_{3}$, and $G_{4}$, they are
also independent given $H_{1},H_{2}$, and $H_{3}$. Therefore, we have,

$\displaystyle G_{5}\mid\alpha_{5},H_{1}$ $\displaystyle\sim
DP(\alpha_{5},H_{1}),$ $\displaystyle G_{6}\mid\alpha_{6},H_{2}$
$\displaystyle\sim DP(\alpha_{6},H_{2}),$ $\displaystyle
G_{7}\mid\alpha_{7},H_{3}$ $\displaystyle\sim DP(\alpha_{7},H_{3}).$

For any finite measurable partition $A_{1},A_{2},\dots,A_{r}$ of $\Theta$,
from Lemma 2, we have

$\displaystyle\left(H_{4}(A_{1}),\dots,H_{4}(A_{r})\right)\mid\alpha_{5},\alpha_{6},\alpha_{7},H_{1},H_{2},H_{3}$
$\displaystyle=\hskip 8.5359pt$
$\displaystyle\gamma_{1}\left(G_{5}(A_{1}),\dots,G_{5}(A_{r})\right)+\gamma_{2}\left(G_{6}(A_{1}),\dots,G_{6}(A_{r})\right)+\gamma_{3}\left(G_{7}(A_{1}),\dots,G_{7}(A_{r})\right)$
$\displaystyle\phantom{\left(G_{7}(A_{1}),\dots,\right)}\sim
Dir\left(\left(\alpha_{5}H_{1}+\alpha_{6}H_{2}+\alpha_{7}H_{3}\right)(A_{1}),\dots,\left(\alpha_{5}H_{1}+\alpha_{6}H_{2}+\alpha_{7}H_{3}\right)(A_{r})\right)$
$\displaystyle\phantom{\left(G_{7}(A_{1}),\dots,\right)}\equiv
Dir\left(\alpha^{*}\left(\left(\frac{\alpha_{5}}{\alpha^{*}}H_{1}+\frac{\alpha_{6}}{\alpha^{*}}H_{2}+\frac{\alpha_{7}}{\alpha^{*}}H_{3}\right)(A_{1}),\dots,\left(\frac{\alpha_{5}}{\alpha^{*}}H_{1}+\frac{\alpha_{6}}{\alpha^{*}}H_{2}+\frac{\alpha_{7}}{\alpha^{*}}H_{3}\right)(A_{r})\right)\right)$
$\displaystyle\phantom{\left(G_{7}(A_{1}),\dots,\right)}\equiv
Dir\left(\alpha^{*}\left(H^{*}(A_{1}),\dots,H^{*}(A_{r})\right)\right)$
$\displaystyle\phantom{\left(G_{7}(A_{1}),\dots,\right)}\Rightarrow
H_{4}\mid\alpha^{*},H^{*}\sim DP(\alpha^{*},H^{*}),$ (26)

where $\alpha^{*}=\alpha_{5}+\alpha_{6}+\alpha_{7}$ and
$H^{*}=\frac{\alpha_{5}}{\alpha^{*}}H_{1}+\frac{\alpha_{6}}{\alpha^{*}}H_{2}+\frac{\alpha_{7}}{\alpha^{*}}H_{3}$.
Note that $\alpha_{5}$, $\alpha_{6}$, and $\alpha_{7}$ are independent gamma
random variables conditionally on $\alpha_{2},\alpha_{3},\alpha_{4}$ with
shape parameters $\alpha_{2}+\alpha_{4},\alpha_{2}+\alpha_{3}$, and
$\alpha_{3}+\alpha_{4}$, respectively. Thus,

$\left(\frac{\alpha_{5}}{\alpha^{*}},\frac{\alpha_{6}}{\alpha^{*}},\frac{\alpha_{7}}{\alpha^{*}}\right)\mid\alpha_{2},\alpha_{3},\alpha_{4}\sim
Dir(\alpha_{2}+\alpha_{4},\alpha_{2}+\alpha_{3},\alpha_{3}+\alpha_{4}).$ (27)

Thus, given $G_{1},G_{2},G_{3},G_{4},\alpha_{2},\alpha_{3},\alpha_{4}$, and
from Appendices A to A, and Eq. 27, using Lemma 2 and using a similar
measurable finite partition of $\Theta$ argument, we have,

$H^{*}\mid\alpha_{2},\alpha_{3},\alpha_{4},G_{1}\sim
DP(2(\alpha_{2}+\alpha_{3}+\alpha_{4}),G_{1}),$ (28)

which completes the proof. ∎

Supplementary Materials for “Graphical Dirichlet Process for Clustering Non-
Exchangeable Grouped Data”

##  A Proof of Lemma 1 and Lemma 2

######  Lemma 1 (Sethuraman, 1994).

Let $\bm{\alpha}_{1}=(\alpha_{11},\alpha_{12},\dots,\alpha_{1k})$ and
$\bm{\alpha}_{2}=(\alpha_{21},\alpha_{22},\dots,\alpha_{2k})$ be
$k$-dimensional vectors with $\alpha_{ij}>0\hskip 5.69046pt\forall\hskip
5.69046ptj=1,2,\dots,k,\hskip 5.69046pti=1,2$. Let $\bm{X}_{1}$ and
$\bm{X}_{2}$ be independent $k$-dimensional random vectors distributed as
Dirichlet distribution with parameters $\bm{\alpha}_{1}$ and
$\bm{\alpha}_{2}$, respectively. Let
$\alpha_{1\cdot}=\sum_{j=1}^{k}\alpha_{1j}$ and
$\alpha_{2\cdot}=\sum_{j=1}^{k}\alpha_{2j}$. Let $\pi$ be independent of
$\bm{X}_{1}$ and $\bm{X}_{2}$ and have a beta distribution
$Beta\left(\alpha_{1\cdot},\alpha_{2\cdot}\right)$. Then the distribution of
$\pi\bm{X}_{1}+\left(1-\pi\right)\bm{X}_{2}$ is the Dirichlet distribution
with parameter $\bm{\alpha}_{1}+\bm{\alpha}_{2}$.

###### Proof.

Let $T_{i}\overset{ind}{\sim}Gamma\left(\alpha_{1i},\lambda\right),\hskip
5.69046pti=1,2,\dots,k$ and
$S_{i}\overset{ind}{\sim}Gamma\left(\alpha_{2i},\lambda\right),\hskip
5.69046pti=1,2,\dots,k$ independently of $T_{i}$, where $\lambda>0$. Let
$T=\sum_{i=1}^{k}T_{i}$ and $S=\sum_{i=1}^{k}S_{i}$. We know from the
reproductive property of independent gamma distributions that $T\sim
Gamma\left(\sum_{i=1}^{k}\alpha_{1i},\lambda\right)\equiv
Gamma\left(\alpha_{1\cdot},\lambda\right)$ and $S\sim
Gamma\left(\sum_{i=1}^{k}\alpha_{2i},\lambda\right)\equiv
Gamma\left(\alpha_{2\cdot},\lambda\right)$ independently of $T$. Define

$\displaystyle\bm{X}_{1}:=\left(\frac{T_{1}}{T},\frac{T_{2}}{T},\dots,\frac{T_{k}}{T}\right),\hskip
11.38092pt\bm{X}_{2}:=\left(\frac{S_{1}}{S},\frac{S_{2}}{S},\dots,\frac{S_{k}}{S}\right),\hskip
11.38092pt\text{and}\hskip 11.38092pt\pi:=\frac{T}{T+S}.$

It is easy to see that $\bm{X}_{1}\sim Dir(\alpha_{11},\dots,\alpha_{1k})$ is
independent of $\bm{X}_{2}\sim Dir(\alpha_{21},\dots,\alpha_{2k})$, and that
$\pi\sim Beta\left(\alpha_{1\cdot},\alpha_{2\cdot}\right)$. We now need to
show that $\pi$ as defined above is indeed independent of $\bm{X}_{1}$ and
$\bm{X}_{2}$ as required by the lemma. For any fixed
$\alpha_{11},\dots,\alpha_{1k}$, we have that $\sum_{i=1}^{k}T_{i}$ is a
complete and sufficient statistic for $\lambda$. Because $\bm{X}_{1}\sim
Dir(\alpha_{11},\dots,\alpha_{1k})$ is ancillary for $\lambda$, by the Basu’s
theorem (Basu, 1955), we have that $\bm{X}_{1}$ is independent of
$\sum_{i=1}^{k}T_{i}=T$. Furthermore, due to the independence of $S_{i}$ and
$T_{i}$, $i=1,\dots,k$, $\bm{X}_{1}$ is independent of $S$, and, therefore,
$\bm{X}_{1}$ is independent of $\pi=\frac{T}{T+S}$. Similarly, $\bm{X}_{2}$ is
also independent of $\pi$. Then,

$\displaystyle\pi\bm{X}_{1}+\left(1-\pi\right)\bm{X}_{2}$
$\displaystyle=\frac{T}{T+S}\left(\frac{T_{1}}{T},\frac{T_{2}}{T},\dots,\frac{T_{k}}{T}\right)+\frac{S}{T+S}\left(\frac{S_{1}}{S},\frac{S_{2}}{S},\dots,\frac{S_{k}}{S}\right)$
$\displaystyle=\left(\frac{T_{1}+S_{1}}{T+S},\frac{T_{2}+S_{2}}{T+S},\dots,\frac{T_{k}+S_{k}}{T+S}\right)\sim
Dir\left(\bm{\alpha}_{1}+\bm{\alpha}_{2}\right),$

because
$T_{i}+S_{i}\overset{ind}{\sim}Gamma(\alpha_{1i}+\alpha_{2i},\lambda)\hskip
5.69046pti=1,2,\dots,k$ and $T+S\sim
Gamma(\alpha_{1\cdot}+\alpha_{2\cdot},\lambda)$. ∎

######  Lemma 2.

Let $\bm{\alpha}_{1},\bm{\alpha}_{2},\dots,\bm{\alpha}_{L}$ be $k$-dimensional
vectors where $\bm{\alpha}_{i}=(\alpha_{i1},\dots,\alpha_{ik})$ with
$\alpha_{ij}>0\hskip 5.69046pt\forall\hskip 5.69046ptj=1,2,\dots,k$,
$i=1,2,\dots,L$. Let $\bm{X}_{1},\bm{X}_{2},\dots,\bm{X}_{L}$ be independent
$k$-dimensional random vectors distributed as Dirichlet distribution with
parameters $\bm{\alpha}_{1},\bm{\alpha}_{2},\dots,\bm{\alpha}_{L}$,
respectively. Let $\alpha_{i\cdot}=\sum_{j=1}^{k}\alpha_{ij},\ \
i=1,2,\dots,L$. Let $\bm{\pi}=\left(\pi_{1},\pi_{2},\dots,\pi_{L}\right)$ be
independent of $\bm{X}_{1},\bm{X}_{2},\dots,\bm{X}_{L}$ and have a Dirichlet
distribution
$Dir\left(\alpha_{1\cdot},\alpha_{2\cdot},\dots,\alpha_{L\cdot}\right)$. Then
the distribution of $\sum_{i=1}^{L}\pi_{i}\bm{X}_{i}$ is the Dirichlet
distribution with parameter $\sum_{i=1}^{L}\bm{\alpha}_{i}$.

###### Proof.

The proof is similar to that of Lemma 1. By noting that

$\bm{\pi}\sim
Dir\left(\alpha_{1\cdot},\alpha_{2\cdot},\dots,\alpha_{L\cdot}\right)\overset{d}{=}\left(\frac{\gamma_{1}}{\gamma},\frac{\gamma_{2}}{\gamma},\dots,\frac{\gamma_{L}}{\gamma}\right),$

where $\gamma_{i}\overset{ind}{\sim}Gamma(\alpha_{i\cdot},\lambda),\hskip
5.69046pti=1,2,\dots,L$ and $\gamma=\sum_{i=1}^{L}\gamma_{i}\sim
Gamma\left(\sum_{i=1}^{L}\alpha_{i\cdot},\lambda\right)$. The remaining proof
follows from standard properties of Dirichlet distributions and mimics the
proof of Lemma 1. ∎

##  B Proof of the Infinite Limit of Finite Mixture Model

The GDP mixture model can be derived as the infinite limit of a finite mixture
model. Let us denote the observations and the mixture component indicator from
node $j$ in layer $k$ of DAG $D$ by $x_{ji}^{(k)}$ and $z_{ji}^{(k)}$,
respectively. Let $\bm{\beta}_{1}^{(0)}$ be the vector of mixing weights for
the root node. Denoting by $\bm{\beta}_{j}^{(k)}$ the mixing weights of node
$j$ in layer $k$ and by $\nu_{j}^{(k,m)}$ the corresponding mixing weights for
the hidden layer $m$, with $m=2,\dots,k$, we have

$\displaystyle\begin{aligned} \bm{\beta}_{1}^{(0)}\mid\alpha_{1}^{(0)}&\sim
Dir\left(\alpha_{1}^{(0)}/L,\dots,\alpha_{1}^{(0)}/L\right),\\\
\bm{\nu}_{j}^{(k,2)}\mid\\{\alpha_{l}^{(1)}:l\in
an^{(k,k-1)}(j)\\},\bm{\beta}_{1}^{(0)}&\sim Dir\left(\sum_{l\in
an^{(k,k-1)}(j)}\alpha_{l}^{(1)}\left(\beta_{11}^{(0)},\dots,\beta_{1L}^{(0)}\right)\right),\\\
&\vdots\\\ \bm{\nu}_{j}^{(k,k)}\mid\\{\alpha_{l}^{(k,k-1)}:l\in
an^{(k,1)}(j)\\},\bm{\nu}_{j}^{(k,k-1)}&\sim Dir\left(\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)}\left(\nu_{j1}^{(k,k-1)},\dots,\nu_{jL}^{(k,k-1)}\right)\right),\\\
\bm{\beta}_{j}^{(k)}\mid\alpha_{j}^{(k)},\bm{\nu}_{j}^{(k,k)}&\sim
Dir\left(\alpha_{j}^{(k)}\left(\nu_{j1}^{(k,k)},\dots,\nu_{jL}^{(k,k)}\right)\right),\\\
\phi_{l}\mid G_{0}&\sim G_{0},\\\
z_{ji}^{(k)}\mid\bm{\beta}_{j}^{(k)}&\sim\bm{\beta}_{j}^{(k)},\\\
x_{ji}^{(k)}\mid z_{ji}^{(k)},\left(\phi_{l}\right)_{l=1}^{L}&\sim
F\left(\phi_{z_{ji}^{(k)}}\right).\end{aligned}$ (1)

###### Proof.

Consider the random probability measure

$G_{1}^{(0)\,L}=\sum_{l=1}^{L}\beta_{1l}^{(0)}\delta_{\phi_{l}}.$

Ishwaran and Zarepour, 2002 shows that for every measurable function $g$,
integrable with respect to $G_{0}$, we have, given $\alpha_{1}^{(0)}$, as
$L\rightarrow\infty$

$\int g(\theta)dG_{1}^{(0),L}(\theta)\overset{\mathcal{D}}{\rightarrow}\int
g(\theta)dG_{1}^{(0)}(\theta).$

Further, consider

$\displaystyle
G_{j}^{(k)\,L}=\sum_{l=1}^{L}\beta_{jl}^{(k)}\delta_{\phi_{l}},$
$\displaystyle
H_{j}^{(k,m)\,L}=\sum_{l=1}^{L}\nu_{jl}^{(k,m)}\delta_{\phi_{l}},\hskip
11.38092ptm=2,\dots,k.$

Let $(A_{1},\dots,A_{r})$ be a measurable partition of the sample space
$\Theta$. Let $K_{t}=\\{l=1,\dots,L:\phi_{l}\in A_{t}\\},\ t=1,\dots,r$, where
$r\leq L$. Assuming that $G_{0}$ is non-atomic, the $\phi_{l}$’s are distinct
with probability one, implying that any partition of $\\{1,\dots,L\\}$
corresponds to some partition of $\Theta$. Thus, as
$\bm{\beta}_{j}^{(k)}\mid\alpha_{j}^{(k)},\bm{\nu}_{j}^{(k,k)}\sim
Dir\left(\alpha_{j}^{(k)}\left(\nu_{j1}^{(k,k)},\dots,\nu_{jL}^{(k,k)}\right)\right)$,
from the properties of Dirichlet distribution, we have,

$\left(G_{j}^{(k)\,L}(A_{1}),\dots,G_{j}^{(k)\,L}(A_{r})\right)=\left(\sum_{l\in
K_{1}}\beta_{jl}^{(k)},\dots,\sum_{l\in K_{r}}\beta_{jl}^{(k)}\right)\sim
Dir\left(\alpha_{j}^{(k)}\sum_{l\in
K_{1}}\nu_{jl}^{(k,k)},\dots,\alpha_{j}^{(k)}\sum_{l\in
K_{r}}\nu_{jl}^{(k,k)}\right).$

Thus,

$\displaystyle G_{j}^{(k)\,L}\mid\alpha_{j}^{(k)},H_{j}^{(k,k)\,L}\sim
DP\left(\alpha_{j}^{(k)},H_{j}^{(k,k)\,L}\right).$

Similarly,

$\displaystyle H_{j}^{(k,k)\,L}\mid\\{\alpha_{l}^{(k-1)}:l\in
an^{(k,1)}(j)\\},H_{j}^{(k,k-1)\,L}\sim DP\left(\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)},H_{j}^{(k,k-1)\,L}\right),$ $\displaystyle
H_{j}^{(k,k-1)\,L}\mid\\{\alpha_{l}^{(k-2)}:l\in
an^{(k,2)}(j)\\},H_{j}^{(k,k-2)\,L}\sim DP\left(\sum_{l\in
an^{(k,2)}(j)}\alpha_{l}^{(k-2)},H_{j}^{(k,k-2)\,L}\right),$
$\displaystyle\vdots$ $\displaystyle
H_{j}^{(k,2)\,L}\mid\\{\alpha_{l}^{(1)}:l\in
an^{(k,k-1)}(j)\\},G_{1}^{(0)\,L}\sim DP\left(\sum_{l\in
an^{(k,k-1)}(j)}\alpha_{l}^{(1)},G_{1}^{(0)\,L}\right).$

By letting $L\rightarrow\infty$, the marginal distribution that this finite
mixture model induces on the observations,
$\bm{x}_{j}^{(k)}=(x_{j1}^{(k)},x_{j2}^{(k)},\dots)$, approaches the proposed
GDP mixture model. ∎

##  C Finite Mixture Model Approximation and Posterior Inference

The posterior inference of the proposed GDP mixture model is carried out using
a blocked Gibbs sampler. For concreteness, we will present the finite mixture
model approximation of the GDP for our motivating example and posterior
inference based on this approximation. In our motivating application, we have
8 experimental groups. Each group corresponds to a combination of treatment,
diet, and genotype; see Table S1 where we use binary indicators to denote the
genotype, the two levels of diet, and the two treatment regimes. The design of
the experiments naturally introduces dependencies among the experimental
groups, which are represented by the DAG in Figure S1, where group 1 is the
root node, groups 2-4 are the layer-1 nodes, groups 5-7 are the layer-2 nodes,
and group 8 is the layer-3 node. For ease of notation, instead of using
$G_{1}^{(0)}$ and $\alpha_{1}^{(0)}$ to denote the random measure and the
concentration parameter of the root node, we use simply $G_{1}$ and
$\alpha_{1}$ instead; similarly for all the other nodes.

$1$$3$$2$$4$$5$$6$$7$$8$ Figure S1: The DAG of experimental groups. |  Group Diet Treatment Genotype $1$ 0 0 0 $2$ 1 0 0 $3$ 0 1 0 $4$ 0 0 1 $5$ 1 0 1 $6$ 1 1 0 $7$ 0 1 1 $8$ 1 1 1 Table S1: Each experimental group corresponds to a combination of diet, treatment, and genotype. Diet = 1 corresponds to high-fat diet and 0 corresponds to normal diet, Treatment = 1 corresponds to AdipoRon and 0 corresponds to no therapy, Genotype = 1 corresponds to Apc knock-out and 0 corresponds to wild type.
---|---

Recall that from the main text, the finite truncation of the infinite mixture
model representation is given by,

$\displaystyle\begin{aligned} \bm{\beta}_{1}^{(0)}\mid\alpha_{1}^{(0)}&\sim
Dir\left(\alpha_{1}^{(0)}/L,\dots,\alpha_{1}^{(0)}/L\right),\\\
\bm{\nu}_{j}^{(k,2)}\mid\\{\alpha_{l}^{(1)}:l\in
an^{(k,k-1)}(j)\\},\bm{\beta}_{1}^{(0)}&\sim Dir\left(\sum_{l\in
an^{(k,k-1)}(j)}\alpha_{l}^{(1)}\left(\beta_{11}^{(0)},\dots,\beta_{1L}^{(0)}\right)\right),\\\
&\vdots\\\ \bm{\nu}_{j}^{(k,k)}\mid\\{\alpha_{l}^{(k,k-1)}:l\in
an^{(k,1)}(j)\\},\bm{\nu}_{j}^{(k-1)}&\sim Dir\left(\sum_{l\in
an^{(k,1)}(j)}\alpha_{l}^{(k-1)}\left(\nu_{j1}^{(k,k-1)},\dots,\nu_{jL}^{(k,k-1)}\right)\right),\\\
\bm{\beta}_{j}^{(k)}\mid\alpha_{j}^{(k)},\bm{\nu}_{j}^{(k,k)}&\sim
Dir\left(\alpha_{j}^{(k)}\left(\nu_{j1}^{(k,k)},\dots,\nu_{jL}^{(k,k)}\right)\right),\\\
\phi_{l}\mid G_{0}&\sim G_{0},\\\
z_{ji}^{(k)}\mid\bm{\beta}_{j}^{(k)}&\sim\bm{\beta}_{j}^{(k)},\\\
x_{ji}^{(k)}\mid z_{ji}^{(k)},\left(\phi_{l}\right)_{l=1}^{L}&\sim
F\left(\phi_{z_{ji}^{(k)}}\right).\end{aligned}$ (2)

Using the simplified notations for the group-specific random measures and
concentration parameter, from the finite truncation of the infinite mixture
model representation from Eq. (2), we have, for this motivating problem,

$\displaystyle\alpha_{1}\mid\alpha_{0}$ $\displaystyle\sim
Gamma(\alpha_{0},1),$ $\displaystyle\bm{\beta}_{1}\mid\alpha_{1}$
$\displaystyle\sim Dir(\alpha_{1}/L,\dots,\alpha_{1}/L),$
$\displaystyle\alpha_{2}\mid\alpha_{1}$ $\displaystyle\sim
Gamma(\alpha_{1},1),$
$\displaystyle\bm{\beta}_{2}\mid\alpha_{2},\bm{\beta}_{1}$ $\displaystyle\sim
Dir(\alpha_{2}\beta_{11},\dots,\alpha_{2}\beta_{1L}),$
$\displaystyle\alpha_{3}\mid\alpha_{1}$ $\displaystyle\sim
Gamma(\alpha_{1},1),$
$\displaystyle\bm{\beta}_{3}\mid\alpha_{3},\bm{\beta}_{1}$ $\displaystyle\sim
Dir(\alpha_{3}\beta_{11},\dots,\alpha_{3}\beta_{1L}),$
$\displaystyle\alpha_{4}\mid\alpha_{1}$ $\displaystyle\sim
Gamma(\alpha_{1},1),$
$\displaystyle\bm{\beta}_{4}\mid\alpha_{4},\bm{\beta}_{1}$ $\displaystyle\sim
Dir(\alpha_{4}\beta_{11},\dots,\alpha_{4}\beta_{1L}),$
$\displaystyle\alpha_{5}\mid\alpha_{2},\alpha_{4}$ $\displaystyle\sim
Gamma(\alpha_{2}+\alpha_{4},1),$
$\displaystyle\bm{\beta}_{5}\mid\alpha_{5},\bm{\nu}_{1}$ $\displaystyle\sim
Dir(\alpha_{5}\nu_{11},\dots,\alpha_{5}\nu_{1L}),$
$\displaystyle\bm{\nu}_{1}\mid\alpha_{2},\alpha_{4},\bm{\beta}_{1}$
$\displaystyle\sim
Dir\left((\alpha_{2}+\alpha_{4})(\beta_{11},\dots,\beta_{1L})\right),$
$\displaystyle\alpha_{6}\mid\alpha_{2},\alpha_{3}$ $\displaystyle\sim
Gamma(\alpha_{2}+\alpha_{3},1),$
$\displaystyle\bm{\beta}_{6}\mid\alpha_{6},\bm{\nu}_{2}$ $\displaystyle\sim
Dir(\alpha_{6}\nu_{21},\dots,\alpha_{6}\nu_{2L}),$
$\displaystyle\bm{\nu}_{2}\mid\alpha_{2},\alpha_{3},\bm{\beta}_{1}$
$\displaystyle\sim
Dir\left((\alpha_{2}+\alpha_{3})(\beta_{11},\dots,\beta_{1L})\right),$
$\displaystyle\alpha_{7}\mid\alpha_{3},\alpha_{4}$ $\displaystyle\sim
Gamma(\alpha_{3}+\alpha_{4},1),$
$\displaystyle\bm{\beta}_{7}\mid\alpha_{7},\bm{\nu}_{3}$ $\displaystyle\sim
Dir(\alpha_{7}\nu_{31},\dots,\alpha_{7}\nu_{3L}),$
$\displaystyle\bm{\nu}_{3}\mid\alpha_{3},\alpha_{4},\bm{\beta}_{1}$
$\displaystyle\sim
Dir\left((\alpha_{3}+\alpha_{4})(\beta_{11},\dots,\beta_{1L})\right),$
$\displaystyle\alpha_{8}\mid\alpha_{5},\alpha_{6},\alpha_{7}$
$\displaystyle\sim Gamma(\alpha_{5}+\alpha_{6}+\alpha_{7},1),$
$\displaystyle\bm{\beta}_{8}\mid\bm{\nu}_{4},\alpha_{8},$ $\displaystyle\sim
Dir(\alpha_{8}\nu_{41},\dots,\alpha_{8}\nu_{4L}),$
$\displaystyle\bm{\nu}_{4}\mid\alpha_{5},\alpha_{6},\alpha_{7},\bm{\eta}$
$\displaystyle\sim
Dir\left((\alpha_{5}+\alpha_{6}+\alpha_{7})(\eta_{1},\dots,\eta_{L})\right),$
$\displaystyle\bm{\eta}\mid\alpha_{2},\alpha_{3},\alpha_{4},\bm{\beta}_{1}$
$\displaystyle\sim
Dir\left(2(\alpha_{2}+\alpha_{3}+\alpha_{4})(\beta_{11},\dots,\beta_{1L})\right),$
$\displaystyle z_{ji}\mid\bm{\beta}_{j}$
$\displaystyle\overset{ind}{\sim}Cat(1:L,\bm{\beta}_{j}),$ $\displaystyle
x_{ji}\mid z_{ji},(\phi_{l})_{l=1}^{L}$
$\displaystyle\overset{ind}{\sim}F(\phi_{z_{ji}}),\hskip 5.69046pt$
$\displaystyle i=1,\dots,n_{j},\hskip 5.69046ptj=1,\dots,8.$ (3)

With the above distributional structure, Gibbs sampling is straightforward. We
use $\pi(.)$ and $\pi(.\mid-)$ to denote the prior distribution and the
conditional distribution, respectively, of the parameter specified in the
argument. The full conditional distribution for the atoms is given by,

$\pi(\\{\phi_{l}\\}_{l=1}^{L}\mid-)\propto\prod_{l=1}^{L}\left[\left\\{\prod_{j=1}^{8}\prod_{i=1}^{n_{j}}F(x_{ji}\mid\phi_{l})^{\mathds{1}(z_{ji}=l)}\right\\}\pi(\phi_{l})\right].$
(4)

The full conditional distributions for the latent cluster labels are given by,

$P(z_{ji}=l\mid-)\propto\beta_{jl}F(x_{ji}\mid\phi_{l}),\hskip
8.5359ptl=1,\dots,L,\ \ i=1,\dots,n_{j}\ \ j=1,\dots,8.$ (5)

The full conditional distribution for the stick-breaking weights is given by,

$\pi(\bm{\beta}_{1}\mid-)\propto\frac{\prod_{l=1}^{L}\beta_{1l}^{m_{1l}+\frac{\alpha_{1}}{L}}\left\\{\beta_{2l}^{\alpha_{2}}\beta_{3l}^{\alpha_{3}}\beta_{4l}^{\alpha_{4}}\nu_{1l}^{\alpha_{2}+\alpha_{4}}\nu_{2l}^{\alpha_{2}+\alpha_{3}}\nu_{3l}^{\alpha_{3}+\alpha_{4}}\eta_{l}^{2(\alpha_{2}+\alpha_{3}+\alpha_{4})}\right\\}^{\beta_{1l}}}{\prod_{l=1}^{L}\left\\{\Gamma\left((\alpha_{2}+\alpha_{4})\beta_{1l}\right)\Gamma\left((\alpha_{2}+\alpha_{3})\beta_{1l}\right)\Gamma\left((\alpha_{3}+\alpha_{4})\beta_{1l}\right)\Gamma\left(2(\alpha_{2}+\alpha_{3}+\alpha_{4})\beta_{1l}\right)\right\\}}\\\
\times\frac{1}{\prod_{l=1}^{L}\left\\{\Gamma\left(\alpha_{2}\beta_{1l}\right)\Gamma\left(\alpha_{3}\beta_{1l}\right)\Gamma\left(\alpha_{4}\beta_{1l}\right)\right\\}},$
(6)

where $m_{1l}=\sum_{i=1}^{n_{1}}\mathds{1}(z_{1i}=l)$, $l=1,\dots,L$. The full
conditionals for $\bm{\beta}_{j},\ \ j=2,\dots,8$, are in closed form,

$\pi(\bm{\beta}_{j}\mid-)\sim Dir(\bm{m}_{j}+\alpha_{j}\bm{\beta}_{1}),\ \
\text{where $\bm{m}_{j}=(m_{j1},\dots,m_{jL})$ and
$m_{jl}=\sum_{i=1}^{n_{j}}\mathds{1}(z_{ji}=l)$, $l=1,\dots,L$}.$ (7)

By letting $\bm{B}(\bm{a})$ to denote the multivariate beta function, i.e.,
for a $L$-dimensional vector $\bm{a}=(a_{1},\dots,a_{L})$ with $a_{i}>0$, we
have,

$\bm{B}(\bm{a})=\frac{\prod_{l=1}^{L}\Gamma(a_{l})}{\Gamma(\sum_{l=1}^{L}a_{l})},$

where $\Gamma(\cdot)$ is the gamma function. Then the full-conditional
distribution of the hidden weights are given by,

$\pi(\bm{\nu}_{1}\mid-)\propto\frac{1}{\bm{B}(\alpha_{5}\bm{\nu}_{1})}\prod_{l=1}^{L}\left\\{\beta_{5l}^{\alpha_{5}\nu_{1l}}\nu_{1l}^{(\alpha_{2}+\alpha_{4})\beta_{1l}-1}\right\\},$
(8)
$\pi(\bm{\nu}_{2}\mid-)\propto\frac{1}{\bm{B}(\alpha_{6}\bm{\nu}_{2})}\prod_{l=1}^{L}\left\\{\beta_{6l}^{\alpha_{6}\nu_{2l}}\nu_{2l}^{(\alpha_{2}+\alpha_{3})\beta_{1l}-1}\right\\},$
(9)
$\pi(\bm{\nu}_{3}\mid-)\propto\frac{1}{\bm{B}(\alpha_{7}\bm{\nu}_{3})}\prod_{l=1}^{L}\left\\{\beta_{7l}^{\alpha_{7}\nu_{3l}}\nu_{3l}^{(\alpha_{3}+\alpha_{4})\beta_{1l}-1}\right\\},$
(10)
$\pi(\bm{\nu}_{4}\mid-)\propto\frac{1}{\bm{B}(\alpha_{8}\bm{\nu}_{4})}\prod_{l=1}^{L}\left\\{\beta_{4l}^{\alpha_{8}\nu_{4l}}\nu_{4l}^{(\alpha_{5}+\alpha_{6}+\alpha_{7})\eta_{l}-1}\right\\},$
(11)
$\pi(\bm{\eta}\mid-)\propto\frac{1}{\bm{B}((\alpha_{5}+\alpha_{6}+\alpha_{7})\bm{\eta})}\prod_{l=1}^{L}\left\\{\eta_{l}^{2(\alpha_{2}+\alpha_{3}+\alpha_{4})\beta_{1l}-1}\nu_{4l}^{(\alpha_{5}+\alpha_{6}+\alpha_{7})\eta_{l}}\right\\}.$
(12)

The full conditionals for the concentration parameters are given by,

$\pi(\alpha_{1}\mid-)\propto\frac{e^{-\alpha_{1}}\alpha_{1}^{\alpha_{0}-1}\alpha_{2}^{\alpha_{1}}\alpha_{3}^{\alpha_{1}}\alpha_{4}^{\alpha_{1}}}{\\{\Gamma(\alpha)\\}^{3}\bm{B}((\alpha_{1}/L,\dots,\alpha_{1}/L)}\prod_{l=1}^{L}\beta_{1l}^{\frac{\alpha_{1}}{L}}$
(13)
$\pi(\alpha_{2}\mid-)\propto\frac{e^{-\alpha_{2}}\alpha_{2}^{\alpha_{1}-1}\alpha_{5}^{\alpha_{2}}\alpha_{6}^{\alpha_{2}}\left[\prod_{l=1}^{L}\left\\{\beta_{2l}^{\beta_{1l}}\nu_{1l}^{\beta_{1l}}\nu_{2l}^{\beta_{1l}}\eta_{l}^{2\beta_{1l}}\right\\}^{\alpha_{2}}\right]\Gamma(\alpha_{2})\Gamma(2(\alpha_{2}+\alpha_{3}+\alpha_{4}))}{\prod_{l=1}^{L}\left\\{\Gamma(\alpha_{2}\beta_{1l})\Gamma((\alpha_{2}+\alpha_{4})\beta_{1l})\Gamma((\alpha_{2}+\alpha_{3})\beta_{1l})\Gamma(2(\alpha_{2}+\alpha_{3}+\alpha_{4})\beta_{1l})\right\\}},$
(14)
$\pi(\alpha_{3}\mid-)\propto\frac{e^{-\alpha_{3}}\alpha_{3}^{\alpha_{1}-1}\alpha_{6}^{\alpha_{3}}\alpha_{7}^{\alpha_{3}}\left[\prod_{l=1}^{L}\left\\{\beta_{3l}^{\beta_{1l}}\nu_{2l}^{\beta_{1l}}\nu_{3l}^{\beta_{1l}}\eta_{l}^{2\beta_{1l}}\right\\}^{\alpha_{3}}\right]\Gamma(\alpha_{3})\Gamma(2(\alpha_{2}+\alpha_{3}+\alpha_{4}))}{\prod_{l=1}^{L}\left\\{\Gamma(\alpha_{3}\beta_{1l})\Gamma((\alpha_{2}+\alpha_{3})\beta_{1l})\Gamma((\alpha_{3}+\alpha_{4})\beta_{1l})\Gamma(2(\alpha_{2}+\alpha_{3}+\alpha_{4})\beta_{1l})\right\\}},$
(15)
$\pi(\alpha_{4}\mid-)\propto\frac{e^{-\alpha_{4}}\alpha_{4}^{\alpha_{1}-1}\alpha_{5}^{\alpha_{4}}\alpha_{7}^{\alpha_{4}}\left[\prod_{l=1}^{L}\left\\{\beta_{4l}^{\beta_{1l}}\nu_{1l}^{\beta_{1l}}\nu_{3l}^{\beta_{1l}}\eta_{l}^{2\beta_{1l}}\right\\}^{\alpha_{4}}\right]\Gamma(\alpha_{4})\Gamma(2(\alpha_{2}+\alpha_{3}+\alpha_{4}))}{\prod_{l=1}^{L}\left\\{\Gamma(\alpha_{4}\beta_{1l})\Gamma((\alpha_{2}+\alpha_{4})\beta_{1l})\Gamma((\alpha_{3}+\alpha_{4})\beta_{1l})\Gamma(2(\alpha_{2}+\alpha_{3}+\alpha_{4})\beta_{1l})\right\\}},$
(16)
$\pi(\alpha_{5}\mid-)\propto\frac{e^{-\alpha_{5}}\alpha_{5}^{\alpha_{2}+\alpha_{4}-1}\alpha_{8}^{\alpha_{5}}\left[\prod_{l=1}^{L}\left\\{\beta_{5l}^{\nu_{1l}}\nu_{4l}^{\eta_{l}}\right\\}^{\alpha_{5}}\right]\Gamma(\alpha_{5})}{\prod_{l=1}^{L}\left\\{\Gamma(\alpha_{5}\nu_{1l})\Gamma((\alpha_{5}+\alpha_{6}+\alpha_{7})\nu_{l})\right\\}},$
(17)
$\pi(\alpha_{6}\mid-)\propto\frac{e^{-\alpha_{6}}\alpha_{6}^{\alpha_{2}+\alpha_{3}-1}\alpha_{8}^{\alpha_{6}}\left[\prod_{l=1}^{L}\left\\{\beta_{6l}^{\nu_{2l}}\nu_{4l}^{\eta_{l}}\right\\}^{\alpha_{6}}\right]\Gamma(\alpha_{6})}{\prod_{l=1}^{L}\left\\{\Gamma(\alpha_{6}\nu_{2l})\Gamma((\alpha_{5}+\alpha_{6}+\alpha_{7})\nu_{l})\right\\}},$
(18)
$\pi(\alpha_{7}\mid-)\propto\frac{e^{-\alpha_{7}}\alpha_{7}^{\alpha_{3}+\alpha_{4}-1}\alpha_{8}^{\alpha_{7}}\left[\prod_{l=1}^{L}\left\\{\beta_{7l}^{\nu_{3l}}\nu_{4l}^{\eta_{l}}\right\\}^{\alpha_{7}}\right]\Gamma(\alpha_{7})}{\prod_{l=1}^{L}\left\\{\Gamma(\alpha_{7}\nu_{3l})\Gamma((\alpha_{5}+\alpha_{6}+\alpha_{7})\nu_{l})\right\\}},$
(19)
$\pi(\alpha_{8}\mid-)\propto\frac{e^{-\alpha_{8}}\alpha_{8}^{\alpha_{5}+\alpha_{6}+\alpha_{7}-1}\left[\prod_{l=1}^{L}\beta_{8l}^{\alpha_{8}\nu_{4l}}\right]\Gamma(\alpha_{8})}{\prod_{l=1}^{L}\Gamma(\alpha_{8}\nu_{4l})}.$
(20)

Note that the full conditionals of $\alpha_{j},\ j=1,\dots,8$,
$\bm{\beta}_{1},\bm{\nu}_{1},\bm{\nu}_{2},\bm{\nu}_{3},\bm{\nu}_{4}$, and
$\bm{\eta}$ are not standard distributions that have direct samplers. We adopt
a Metropolis-within-Gibbs strategy to sample from their corresponding full
conditional distributions. Since $\alpha_{j}$’s are real-valued, sampling
using a Metropolis step is straightforward. However, the main bottleneck in
sampling are the weights
$\bm{\beta}_{1},\bm{\nu}_{1},\bm{\nu}_{2},\bm{\nu}_{3},\bm{\nu}_{4}$ and
$\bm{\eta}$, which have a complex structure on the simplex. To mitigate this
problem, we use the SALTSampler (Director et al., 2017) for which the
implementation is publicly available as an R package.

##  D Simulation details

Our simulations are designed to mimick the motivating application where we
have 8 experimental groups. See Table S1 for our experimental design
represented in terms of binary indicators denoting the levels of diet,
treatment, and genotype. The corresponding DAG is given in Figure S1.

For our simulation study, we generated data within each of the $8$ groups from
a four-component mixture of bivariate Gaussian distributions with different
covariance matrices for each group. Taking $\alpha_{0}=5$, we drew the
concentration parameters for the different groups $\alpha_{j}$’s, the mixture
model weights, $\bm{\beta}_{j}$’s, $\bm{\nu}_{j}$’s, and $\bm{\eta}_{j}$, and
the true cluster indicators $z_{ji}$’s for each of the different groups using
(C). Given the cluster indicators, the data were generated from the Gaussian
distribution with the true cluster-specific means $\phi_{l}$’s given in Table
S2 and the group-specific covariance matrices given in Table S3. Note that
within each group, the same covariance matrix was used for all clusters.

Cluster | Mean
---|---
1 | (-2, -5)
2 | (0, 0)
3 | (-3, 3)
4 | (3, -3)
Table S2: True cluster-specific means. Group | Covariance
---|---
1 | $\begin{bmatrix}0.8&0.3\\\ 0.3&0.8\end{bmatrix}$
2 | $\begin{bmatrix}0.85&0.25\\\ 0.25&0.85\end{bmatrix}$
3 | $\begin{bmatrix}1&0.1\\\ 0.1&1\end{bmatrix}$
4 | $\begin{bmatrix}0.8&-0.1\\\ -0.1&0.8\end{bmatrix}$
5 | $\begin{bmatrix}0.8&-0.2\\\ -0.2&0.9\end{bmatrix}$
6 | $\begin{bmatrix}0.8&0\\\ 0&0.8\end{bmatrix}$
7 | $\begin{bmatrix}0.75&0.25\\\ 0.25&0.75\end{bmatrix}$
8 | $\begin{bmatrix}1.1&0.1\\\ 0.1&1.1\end{bmatrix}$
Table S3: True covariance matrices for different groups.

In our Gibbs sampler, the truncation level of the finite mixture model was set
to $L=10$, and the base measure for GDP, $G_{0}$, was specified as the normal-
inverse-Wishart distribution, $\mathcal{NIW}(\bm{0},0.01,\mathbb{I}_{2},2)$.
Upon the completion of the Gibbs sampler, the clusters were estimated by using
the least squares criterion (Dahl, 2006), and they were compared with the true
cluster labels for evaluation. We considered a variety of sample sizes as well
as a case with very imbalanced design, which are summarized in Table S4. In
all cases, we ran 15,000 iterations of our Gibbs sampler and discarded the
first 5,000 samples as burn-in.

Group | Sample sizes
---|---
| small | moderate | large | unbalanced
$1$ | 40 | 80 | 150 | 350
$2$ | 30 | 70 | 160 | 30
$3$ | 30 | 70 | 180 | 40
$4$ | 35 | 75 | 170 | 45
$5$ | 25 | 83 | 155 | 25
$6$ | 30 | 88 | 175 | 25
$7$ | 25 | 92 | 185 | 35
$8$ | 30 | 88 | 145 | 35
Table S4: The sample sizes for the different groups that were used to simulate
the data.

We presented the results of clustering for small sample sizes and unbalanced
sample sizes in the main document. Figure S2 shows the results of clustering
for moderate and large sample sizes in each group.

(a) Moderate sample size in each group

(b) Large sample size in each group

Figure S2: Clustering performance of GDP for additional sample sizes. The
colors indicate the estimated clusters by GDP. Adjusted Rand index is reported
at the top of each panel.

We further considered the case, wherein the simulation scenario was difficult
with possibly overlapping clusters. We generated data within each of the $8$
groups from a ten-component mixture of bivariate Gaussian distributions with
different covariance matrices for each group. The choice of mixture model
weights for the first fours groups are summarized in Table S5.

Group |  Mixture weights $\bm{\beta}_{j}$
---|---
1 | $(0.100,0.100,0.100,0.100,0.100,0.100,0.100,0.100,0.100,0.100)^{\top}$
2 | $(0.167,0.167,0.167,0.167,0.167,0.056,0.056,0.056,0.000,0.000)^{\top}$
3 | $(0.095,0.095,0.095,0.000,0.000,0.143,0.143,0.143,0.143,0.143)^{\top}$
4 | $(0.030,0.030,0.030,0.182,0.182,0.182,0.182,0.182,0.000,0.000)^{\top}$
Table S5: True group-specific mixture model weights.

The mixture weights for all other groups were taken to be the mean of the
mixture weights of their parent, e.g., the mixture weight for group 5 was the
mean of the mixture weights of groups 2 and 4. The true cluster indicators
$z_{ji}$’s for each of the different groups were drawn using (C) and the true
mixture weights. Given the cluster indicators, the data were generated from
the Gaussian distribution with the true cluster-specific means $\phi_{l}$’s
given in Table S6 and the group-specific covariance matrices given in Table
S3.

Cluster | Mean
---|---
1 | (-2.5, 0)
2 | (0, 0)
3 | (2.5, 0)
4 | (2.5, -2.5)
5 | (-3, -3)
6 | (2, 2)
7 | (-2, 5)
8 | (5, 8)
9 | (-5, -8)
10 | (8, -8)
Table S6: True cluster-specific means.

In our Gibbs sampler, the truncation level of the finite mixture model was set
to $L=20$, the hyperparameter $\alpha_{0}$ was taken to be 1, and the base
measure for GDP, $G_{0}$, was specified as the normal-inverse-Wishart
distribution, $\mathcal{NIW}(\bm{0},0.01,\mathbb{I}_{2},2)$. Upon the
completion of the Gibbs sampler, the clusters were estimated by using the
least squares criterion (Dahl, 2006), and they were compared with the true
cluster labels for evaluation. We again considered a variety of sample sizes
as summarized in Table S4. In all cases, we ran 25,000 iterations of our Gibbs
sampler and after discarding the first 15,000 samples as burn-in, considered
thinning of the samples by a factor 5. The clustering results are shown in
Figure S4. Clearly, GDP was able to identify the overlapping clusters within
each group and link them across groups for all simulation scenarios with
reasonable accuracy as measured by Adjusted Rand indices for each group (shown
in the plots).

(a) Small sample size in each group

(b) Moderate sample size in each group

(a) Large sample size in each group

(b) Unbalanced sample size in each group

Figure S4: Clustering performance of GDP for various sample sizes and
difficult simulation scenario. The colors indicate the estimated clusters by
GDP. Adjusted Rand index is reported at the top of each panel.

In the main manuscript, we reported the boxplot of Adjusted Rand indices for
50 replicates. Further investigation regarding the choice of $\alpha_{0}$
revealed no significant impact in clustering performance. Figure S5 shows that
boxplot of Adjusted Rand indices for 50 replicates, comparing GDP, HDP, and
k-means with $\alpha_{0}$ taken to be 6. In all situations, GDP uniformly out-
performed the other two methods.

Figure S5: The boxplots of the Adjusted Rand indices for GDP, HDP, and k-means
for all sample sizes. In all simulations $\alpha_{0}$ was taken to be 6.

##  E Real Data Analysis plots

We present the traceplots (Figure S7) of the log-likelihood for each of the
four parallel chains of our sampler, corresponding to the real data analysis
of the main document, after discarding the initial $25,000$ samples and
thinning the samples by a factor of $15$. The traceplots indicate the presence
of local modes, necessitating the need to concatenate posterior samples across
these chains for more efficient and reliable inference.

(a) GDP (a) HDP Figure S7: Traceplots of log-likelihood for each group for (a)
GDP and (b) HDP.

## References

  * Antoniak (1974) Antoniak, C. E.  Mixtures of dirichlet processes with applications to bayesian nonparametric problems.  _The Annals of Statistics_ , 2(6):1152–1174, 1974\.  ISSN 00905364.  URL http://www.jstor.org/stable/2958336.
  * Barrios et al. (2013) Barrios, E., Lijoi, A., Nieto-Barajas, L. E., and Prünster, I.  Modeling with normalized random measure mixture models.  _Statistical Science_ , 28(3):313–334, 2013.
  * Basu (1955) Basu, D.  On statistics independent of a complete sufficient statistic.  _Sankhyā: The Indian Journal of Statistics (1933-1960)_ , 15(4):377–380, 1955.  ISSN 00364452.  URL http://www.jstor.org/stable/25048259.
  * Bouchard-Mercier et al. (2013) Bouchard-Mercier, A., Paradis, A.-M., Rudkowska, I., Lemieux, S., Couture, P., and Vohl, M.-C.  Associations between dietary patterns and gene expression profiles of healthy men and women: a cross-sectional study.  _Nutrition Journal_ , 12(1):24, 2013.  doi: 10.1186/1475-2891-12-24.  URL https://doi.org/10.1186/1475-2891-12-24.
  * Camerlenghi et al. (2019) Camerlenghi, F., Dunson, D. B., Lijoi, A., Prünster, I., and Rodríguez, A.  Latent Nested Nonparametric Priors (with Discussion).  _Bayesian Analysis_ , 14(4):1303 – 1356, 2019\.  doi: 10.1214/19-BA1169.  URL https://doi.org/10.1214/19-BA1169.
  * Cifarelli and Regazzini (1978) Cifarelli, D. and Regazzini, E.  Problemi statistici non parametrici in condizioni di scambiabilita parziale e impiego di medie associative.  Technical report, Tech. rep., Quaderni Istituto Matematica Finanziaria dell’Universita di Torino, 1978.
  * Dahl (2006) Dahl, D. B.  Model-based clustering for expression data via a dirichlet process mixture model.  _Bayesian Inference for Gene Expression and Proteomics_ , 2006.
  * Dahl et al. (2017) Dahl, D. B., Day, R., and Tsai, J. W.  Random partition distribution indexed by pairwise information.  _Journal of the American Statistical Association_ , 112(518):721–732, 2017.  doi: 10.1080/01621459.2016.1165103.  URL https://doi.org/10.1080/01621459.2016.1165103.  PMID: 29276318.
  * Das et al. (2023) Das, S., Niu, Y., Ni, Y., Mallick, B. K., and Pati, D.  Blocked gibbs sampler for hierarchical dirichlet processes, 2023.
  * De Blasi et al. (2013) De Blasi, P., Favaro, S., Lijoi, A., Mena, R. H., Prünster, I., and Ruggiero, M.  Are gibbs-type priors the most natural generalization of the dirichlet process?  _IEEE transactions on pattern analysis and machine intelligence_ , 37(2):212–229, 2013.
  * de Finetti (1938) de Finetti, B.  Sur la condition d’equivalence partielle.  _Actual. Sci. Ind._ , 739:5–18, 1938.
  * Dey et al. (2022) Dey, D., Datta, A., and Banerjee, S.  Graphical gaussian process models for highly multivariate spatial data.  _Biometrika_ , 109(4):993–1014, 2022.
  * Director et al. (2017) Director, H. M., Gattiker, J., Lawrence, E., and Wiel, S. V.  Efficient sampling on the simplex with a self-adjusting logit transform proposal.  _Journal of Statistical Computation and Simulation_ , 87(18):3521–3536, 2017.  doi: 10.1080/00949655.2017.1376063.  URL https://doi.org/10.1080/00949655.2017.1376063.
  * Dunson and Herring (2005) Dunson, D. and Herring, A.  Bayesian latent variable models for mixed discrete outcomes.  _Biostatistics (Oxford, England)_ , 6:11–25, 02 2005.  doi: 10.1093/biostatistics/kxh025.
  * Escobar and West (1995) Escobar, M. D. and West, M.  Bayesian density estimation and inference using mixtures.  _Journal of the American Statistical Association_ , 90(430):577–588, 1995.  doi: 10.1080/01621459.1995.10476550.  URL https://www.tandfonline.com/doi/abs/10.1080/01621459.1995.10476550.
  * Fan et al. (2020) Fan, X., Yao, H., Liu, X., Shi, Q., Lv, L., Li, P., Wang, R., Tang, T., and Qi, K.  High-fat diet alters the expression of reference genes in male mice.  _Frontiers in Nutrition_ , 7, 2020.  ISSN 2296-861X.  doi: 10.3389/fnut.2020.589771.  URL https://www.frontiersin.org/articles/10.3389/fnut.2020.589771.
  * Ferguson (1973) Ferguson, T. S.  A bayesian analysis of some nonparametric problems.  _The Annals of Statistics_ , 1(2):209–230, 1973\.  doi: 10.1214/aos/1176342360.  URL https://doi.org/10.1214/aos/1176342360.
  * Gelfand et al. (2005) Gelfand, A. E., Kottas, A., and MacEachern, S. N.  Bayesian nonparametric spatial modeling with dirichlet process mixing.  _Journal of the American Statistical Association_ , 100(471):1021–1035, 2005.  ISSN 01621459.  URL http://www.jstor.org/stable/27590632.
  * Griffin and Steel (2006) Griffin, J. E. and Steel, M. F. J.  Order-based dependent dirichlet processes.  _Journal of the American Statistical Association_ , 101(473):179–194, 2006.  ISSN 01621459.  URL http://www.jstor.org/stable/30047448.
  * Gu and Dunson (2023) Gu, Y. and Dunson, D. B.  Bayesian pyramids: identifiable multilayer discrete latent structure models for discrete data.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 85(2):399–426, 2023.
  * Hao et al. (2021) Hao, Y., Hao, S., Andersen-Nissen, E., III, W. M. M., Zheng, S., Butler, A., Lee, M. J., Wilk, A. J., Darby, C., Zagar, M., Hoffman, P., Stoeckius, M., Papalexi, E., Mimitou, E. P., Jain, J., Srivastava, A., Stuart, T., Fleming, L. B., Yeung, B., Rogers, A. J., McElrath, J. M., Blish, C. A., Gottardo, R., Smibert, P., and Satija, R.  Integrated analysis of multimodal single-cell data.  _Cell_ , 2021.  doi: 10.1016/j.cell.2021.04.048.  URL https://doi.org/10.1016/j.cell.2021.04.048.
  * Hjort et al. (2010) Hjort, N. L., Holmes, C., Müller, P., and Walker, S. G.  _Bayesian nonparametrics_ , volume 28.  Cambridge University Press, 2010.
  * Hubert and Arabie (1985) Hubert, L. and Arabie, P.  Comparing partitions.  _Journal of Classification_ , 2(1):193–218, 1985\.  doi: 10.1007/BF01908075.  URL https://doi.org/10.1007/BF01908075.
  * Hwang et al. (2018) Hwang, B., Lee, J. H., and Bang, D.  Single-cell rna sequencing technologies and bioinformaticspipelines.  _Experimental & Molecular Medicine_, 50(8):1–14, 2018.  doi: 10.1038/s12276-018-0071-8.  URL https://doi.org/10.1038/s12276-018-0071-8.
  * Iorio et al. (2004) Iorio, M. D., Müller, P., Rosner, G. L., and MacEachern, S. N.  An anova model for dependent random measures.  _Journal of the American Statistical Association_ , 99(465):205–215, 2004.  doi: 10.1198/016214504000000205.  URL https://doi.org/10.1198/016214504000000205.
  * Iorio et al. (2006) Iorio, M. D., Johnson, W. O., Müller, P., and Rosner, G. L.  A ddp model for survival regression.  2006\.
  * Ishwaran and James (2001) Ishwaran, H. and James, L. F.  Gibbs sampling methods for stick-breaking priors.  _Journal of the American Statistical Association_ , 96(453):161–173, 2001.  doi: 10.1198/016214501750332758.  URL https://doi.org/10.1198/016214501750332758.
  * Ishwaran and Zarepour (2002) Ishwaran, H. and Zarepour, M.  Exact and approximate sum representations for the dirichlet process.  _The Canadian Journal of Statistics / La Revue Canadienne de Statistique_ , 30(2):269–283, 2002.  ISSN 03195724.  URL http://www.jstor.org/stable/3315951.
  * Jump and Clarke (1999) Jump, D. B. and Clarke, S. D.  Regulation of gene expression by dietary fat.  _Annual Review of Nutrition_ , 19(1):63–90, 1999\.  doi: 10.1146/annurev.nutr.19.1.63.  URL https://doi.org/10.1146/annurev.nutr.19.1.63.  PMID: 10448517.
  * Kleinman and Ibrahim (1998) Kleinman, K. P. and Ibrahim, J. G.  A semiparametric bayesian approach to the random effects model.  _Biometrics_ , 54(3):921–938, 1998.  ISSN 0006341X, 15410420.  URL http://www.jstor.org/stable/2533846.
  * Lijoi et al. (2022) Lijoi, A., Prünster, I., and Rebaudo, G.  Flexible clustering via hidden hierarchical dirichlet priors.  _Scandinavian Journal of Statistics_ , feb 2022.  doi: 10.1111/sjos.12578.  URL https://doi.org/10.1111%2Fsjos.12578.
  * Liu et al. (2010) Liu, Y., Li, Z., Xiong, H., Gao, X., and Wu, J.  Understanding of internal clustering validation measures.  In _2010 IEEE International Conference on Data Mining_ , pages 911–916, 2010.  doi: 10.1109/ICDM.2010.35.
  * MacEachern (1999) MacEachern, S. N.  Dependent nonparametric processes.  In _ASA Proceedings of the Section on Bayesian Statistical Science_ , Alexandria, VA, 1999. American Statistical Association.
  * MacEachern (2000) MacEachern, S. N.  Dependent dirichlet processes.  Technical report, Department of Statistics, The Ohio State University, 2000.
  * MacEachern and Müller (1998) MacEachern, S. N. and Müller, P.  Estimating mixture of dirichlet process models.  _Journal of Computational and Graphical Statistics_ , 7(2):223–238, 1998.  ISSN 10618600.  URL http://www.jstor.org/stable/1390815.
  * Mallick and Walker (1997) Mallick, B. K. and Walker, S. G.  Combining information from several experiments with nonparameter priors.  _Biometrika_ , 84(3):697–706, 1997.  ISSN 00063444.  URL http://www.jstor.org/stable/2337589.
  * McInnes et al. (2018) McInnes, L., Healy, J., and Melville, J.  Umap: Uniform manifold approximation and projection for dimension reduction, 2018.  URL https://arxiv.org/abs/1802.03426.
  * Morin et al. (1997) Morin, P. J., Sparks, A. B., Korinek, V., Barker, N., Clevers, H., Vogelstein, B., and Kinzler, K. W.  Activation of β-catenin-tcf signaling in colon cancer by mutations in β-catenin or apc.  _Science_ , 275(5307):1787–1790, 1997.  doi: 10.1126/science.275.5307.1787.  URL https://www.science.org/doi/abs/10.1126/science.275.5307.1787.
  * Müller et al. (2015) Müller, P., Quintana, F. A., Jara, A., and Hanson, T.  _Bayesian nonparametric data analysis_ , volume 1.  Springer, 2015.
  * Müller et al. (2004) Müller, P., Quintana, F., and Rosner, G.  A method for combining inference across related nonparametric bayesian models.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 66(3):735–749, 2004.  ISSN 13697412, 14679868.  URL http://www.jstor.org/stable/3647503.
  * Nieto-Barajas and Contreras-Cristán (2014) Nieto-Barajas, L. E. and Contreras-Cristán, A.  A bayesian nonparametric approach for time series clustering.  _Bayesian Analysis_ , 9(1):147–170, 2014.
  * Pearl (2009) Pearl, J.  _Causality_.  Cambridge university press, 2009.
  * Pitman (2002) Pitman, J.  Poisson–dirichlet and gem invariant distributions for split-and-merge transformations of an interval partition.  _Comb. Probab. Comput._ , 11(5):501–514, sep 2002.  ISSN 0963-5483.  doi: 10.1017/S0963548302005163.  URL https://doi.org/10.1017/S0963548302005163.
  * Quintana et al. (2020) Quintana, F. A., Mueller, P., Jara, A., and MacEachern, S. N.  The dependent dirichlet process and related models, 2020.  URL https://arxiv.org/abs/2007.06129.
  * Ren et al. (2008) Ren, L., Dunson, D. B., and Carin, L.  The dynamic hierarchical dirichlet process.  In _Proceedings of the 25th international conference on machine learning_ , pages 824–831, 2008.
  * Rodríguez et al. (2008) Rodríguez, A., Dunson, D. B., and Gelfand, A. E.  The nested dirichlet process.  _Journal of the American Statistical Association_ , 103(483):1131–1154, 2008.  doi: 10.1198/016214508000000553.  URL https://doi.org/10.1198/016214508000000553.
  * Sethuraman (1994) Sethuraman, J.  A constructive definition of dirichlet priors.  _Statistica Sinica_ , 4(2):639–650, 1994.  ISSN 10170405, 19968507.  URL http://www.jstor.org/stable/24305538.
  * Srebro and Roweis (2005) Srebro, N. and Roweis, S.  Time-varying topic models using dependent dirichlet processes.  _UTML, TR# 2005_ , 3, 2005.
  * Teh et al. (2006) Teh, Y. W., Jordan, M. I., Beal, M. J., and Blei, D. M.  Hierarchical dirichlet processes.  _Journal of the American Statistical Association_ , 101(476):1566–1581, 2006.  doi: 10.1198/016214506000000302.  URL https://doi.org/10.1198/016214506000000302.
  * Williamson et al. (2013) Williamson, S., Dubey, A., and Xing, E.  Parallel markov chain monte carlo for nonparametric mixture models.  In _International Conference on Machine Learning_ , pages 98–106. PMLR, 2013.
  * Zhang et al. (2010) Zhang, J., Song, Y., Zhang, C., and Liu, S.  Evolutionary hierarchical dirichlet processes for multiple correlated time-varying corpora.  In _Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 1079–1088, 2010.