Non-parametric Bayesian modeling of complex networks

A simple and very powerful approach to modeling structure in a complex network is to use a mixture model, leading to a Bayesian version of the socalled stochastic blockmodel [42]. In a mixture model, the observations are assumed to be distributed according to a mixture of $K$ components belonging to some parametric family. Conditioned on knowing which mixture components generated each datum, the observations are assumed independent. In a mixture model for network data, each node belongs to a single mixture component, and since each edge is associated with two nodes, its likelihood will depend on two components. Thus, the likelihood of the network will take the following form,

where the product ranges over all node pairs, and the parameters are given by $\theta=\big{\{}\{z_{i}\}_{i=1}^{N},\phi\big{\}}$ where $z_{i}$ indicates which mixture component the $i$th node belongs to and $\phi$ denotes any further parameters. In the most simple setting, each term in the likelihood could be a Bernoulli distribution (a biased coin flip),

such that $\phi_{k,\ell}$ denotes the probability of an edge between two nodes in group $k$ and $\ell$. To finish the specification of the model, we must define prior distributions for the mixture component indicators $z$ as well as the link probabilities $\phi$. Starting with $\phi$, a natural choice would be independent $\mathrm{Beta}$ distributions for each pair of components,

where the parameters for example can be set to $a=b=1$ to yield a uniform distribution. A natural choice for $z$ would be a $K$-dimensional categorical distribution,

parameterized by $\pi=\{\pi_{k}\}_{k=1}^{K}$ where $\sum_{k=1}^{K}\pi_{k}=1$. How, then, should $\pi$ be chosen? We could for example set each of these parameters to a fixed value, e.g. $\pi_{k}=\tfrac{1}{K}$, but this would be a strong prior assumption specifying that the mixture components have the same number of members on average. A more flexible option would be to define a hierarchical prior, where $\pi$ is generated from a Dirichlet distribution,

where $\mathrm{B}(\alpha)$ is the multinomial beta function, which can be expressed using the gamma function,

Since each component a priori is equally likely, we select the concentration parameters to be equal to each other, $\alpha_{1}=\dots=\alpha_{K}=\tfrac{A}{K}$ such that the scale of the distribution is $\sum_{k=1}^{K}\alpha_{k}=A$. This results in a joint prior over $z$ and $\pi$ given by

where $n_{k}$ denotes the number of $z_{i}$’s with the value $k$.

As we are not particularly interested in the mixture component probabilities $\pi$ (they can be considered nuisance parameters) we can compute the effective prior over $z$ by marginalizing over $\pi$ which has a closed form expression due to the conjugacy between the Dirichlet and Categorical distributions (i.e., the posterior distribution of $\pi$ has the same functional form as the prior),

This resulting effective prior distribution is known as a multivariate Pólya distribution.

Furthermore, the link probabilities $\phi$ can also be considered nuisance parameters, and can also be marginalized analytically due to the conjugacy between the Beta and Bernoulli distributions,

where the product ranges over all pairs of components and $m_{k,\ell}$ and $\bar{m}_{k,\ell}$ denote the number of links and non-links between nodes in component $k$ and $\ell$ respectively.

In the previous section we specified a parametric Bayesian mixture model for complex networks. In the following we move to the non-parametric setting in which the number of mixture components is allowed to be countably infinite. First, consider what happens when the number of components is much larger than the number of nodes in the graph. In that situation, many of the components will not have any nodes assigned to them; in fact, no more than $N$ components can be non-empty, corresponding to the worst case situation where each node has a component of its own. To handle the situation with an infinite number of components, we can not explicitly represent the components but, as we will show in the following, we need only an explicit representation of the finite number of non-empty components.

As we defined the model so far we have introduced $K$ labelled mixture components. This means that if we for example have $N=5$ nodes and $K=4$ components, we assign a separate probability to, say, the configurations $\{1,2,1,4,2\}$ and $\{3,4,3,2,4\}$ even though they correspond to the same clustering of the network nodes. A better choice is to specify the probability distribution directly over the equivalence class of partitions of the network nodes. Since we have $K$ labels in total to choose from, there are $K$ possible labellings for the first component, $K-1$ for the second, etc. resulting in a total of

labellings corresponding to the same partitioning, where $\bar{K}$ is the number of non-empty components. Thus, defining a parameter $\bar{z}$ that holds the partitioning of the network nodes, we have

Since $\bar{z}$ represents partitions rather than labels it can be finitely represented. We can now simply let the number of components go to infinity by computing the limit of the prior distribution for $\bar{z}$,

The details involved in computing this limit can be found in [15] and [40]. The limiting stochastic process is known as a Chinese restaurant process (CRP) [2] (for an introduction to the CRP, see [14]). Compactly, we may write

In summary, the generative process for the infinite relational model can be expressed as follows:

The network nodes are partitioned according to a Chinese restaurant process; a probability of linking between each pair of node clusters is simulated from a Beta distribution; and each link in the network is generated according to a Bernoulli distribution depending on which clusters the pair of nodes belong to.

Identically, in the notation of exponential random graph models the likelihood in Eq. (29) can be expressed as

where the sufficient statistics are the counts of links between each pair of clusters, $s(X,z)=\{m_{k,\ell}\}$, and the natural parameter is the log odds of links between each pair of clusters, $\theta(\phi)=\left\{\log\tfrac{\phi_{k,\ell}}{1-\phi_{k,\ell}}\right\}$.

In Gibbs sampling the variables are iteratively sampled from their conditional distribution, and repeating this process the samples will eventually approximate the posterior distribution. We iteratively sample the partition assignments, $\bar{z}_{n}$, from their conditional distribution,

where $\bar{z}^{\setminus n}$ denotes all partition assignments except $\bar{z}_{n}$. An expression for this conditional distribution can be found by considering which terms in the likelihood and prior will change when node $n$ is assigned to a different partition. For the prior in Eq. (24) we have

where $n^{-}_{k}$ is the number of nodes associated with component $k$ not counting node $n$. Adding node $n$ to an existing component increases the argument of the corresponding Gamma function by one, effectively multiplying the prior by $n^{-}_{k}$, whereas adding the node to a new cluster increases $\bar{K}$ by one, effectively multiplying the prior by $A$. For the likelihood in Eq. (22), adding node $n$ to partition $k$ effectively multiplies the likelihood by

where $m^{\setminus n}_{k,\ell}$ and $\bar{m}^{\setminus n}_{k,\ell}$ denote the number of links and non-links between nodes in component $k$ and $\ell$, not counting any links from node $n$, and $r_{n,\ell}$ is the number of links from node $n$ to any nodes in component $\ell$.

In order to perform Gibbs sampling we can now simply consider each node in turn; for each partition (including a new, empty partition) compute the product of Eq. (35) and (36); normalize to yield a categorical distribution over partitions; and sample a new $\bar{z}_{n}$ according to this distribution. The final algorithm is summarized in Fig. 5. The result after running the Gibbs sampler for $2T$ iterations is a set of samples of $\bar{z}$, where usually the first half is discarded for burn in. This yields a final ensemble $\{\bar{z}^{(t)}:t\in 1,\dots,T\}$ approximately sampled from the posterior.

In the algorithm outlined in Figure 5 it can be observed that there are two loops: One over the $T$ simulated samples and one over the $N$ nodes in the network. In each run of the inner loop, a node is assigned to a cluster by the Gibbs sampler. In the following we consider the number of clusters $K$ a constant (although of course it will vary depending on the network data), and examine how the computational complexity of the algorithm depends on the number of nodes and edges in the network.

In the code in Figure 5 the variables M0, M1 and m, which hold the counts of nonlinks, links, and nodes, are re-computed in each iteration. In a more sensible implementation, these quantities would be precomputed and efficiently updated during the Gibbs sampling.

Evaluating the probability of assigning a node to each cluster then requires the computation of the vector r which holds the count of links from node $n$ to each of the clusters. The time complexity of this computation is on the order of the node degree. Looping over the nodes gives a total time complexity of $\mathrm{O}(L)$ where $L$ is the number of edges in the graph. To calculate the probabilities of assigning the nodes to the clusters for all $N$ Gibbs samples requires $2K^{2}N$ evaluations of the (logarithm of the) Beta function so the time complexity of this computation is $\mathrm{O}(N)$. As a result, since in general $L>N$, the total computational complexity of the Gibbs sampler for the IRM model is $\mathrm{O}(L)$. Figure 6 demonstrates that this linear scaling is observed in practice when analyzing networks of varying numbers of nodes and edges.

For comparison, Monte Carlo maximum likelihood inference in exponential random graph model based on endogenous network statistics requires the simulation of random networks from the ERGM distribution, which is inherently an $\mathrm{O}(N^{2})$ operation. We should note though, that in practice we would not expect the time complexity of the IRM to scale linearly in the number of edges, since the number of clusters most likely would increase with the size of the network and since the number of required iterations of the Gibbs sampler might also go up.

Say we are interested in some function $f(\bar{z})$ that depends on the model. We can now compute the posterior distribution of this quantity,

approximated by an empirical average over the posterior samples. For example, the approach can be used to compute the posterior distribution over the number of components in the mixture model or other quantities of interest.

Missing data is easily handled in the Bayesian framework, simply by leaving out the terms in the likelihood corresponding to unobserved links. If we observe only a part of the network and are interested in predicting the presence or absence of an unobserved link between two nodes, we can simply compute the posterior predictive distribution of the missing link,

Here $X$ denotes the observed part of the network, and $\bar{z}^{(t)}$ is simulated from the posterior distribution where only the observed part of the network is conditioned on. Inserting $p(x_{i,j}|\bar{z},X)=\int p(x_{i,j}|\theta,\bar{z})p(\theta|\bar{z},X)d\theta$ yields

where

Predicting missing links can be used to compare different models: A number of links can be excluded when fitting the models which can then be compared by assessing their ability to predict the held-out links. Since the links in a network are highly correlated and because many networks exhibit a highly imbalanced distribution of links and nonlinks, care must be taken in choosing a hold out test set in an appropriate way. If the test set is chosen to balance the number of links and nonlinks, its distribution will not correspond to the full network which makes the absolute link prediction results difficult to interpret. Thus, although indicative of a model’s predictive performance, this approach is perhaps best suited for the relative comparison of different models. If, on the other hand, several examples of full networks are available, a whole network can used as test data making the absolute link prediction results directly interpretable.

Finally, we might be interested in examining how well our model describes the data to assess if the model is appropriate for the data at hand or if a more suitable model should be constructed. A principled approach to achieving this is posterior predictive checking. First, an ensemble of replicated networks from the posterior predictive distribution is generated from,

which as before can be approximated using samples of $\bar{z}$ simulated from the posterior using Eq. (41). Now, the idea is to compare characteristics of the observed network, such as the degree distribution, clustering coefficient, and characteristic path length, with the posterior predictive distribution of these properties, approximated by the empirical distribution over the ensemble of replicated networks. If the model fits well, the observed characteristic of the network should be quite likely under the posterior predictive distribution, whereas a large discrepancy indicates model mismatch. Posterior predictive checking is useful for model critisism, i.e., for exploring lack of fit as opposed to testing whether the model is correct. Discovering network characteristics for which the model does not fit the data well can inspire to the development of more sophisticated models; however, even a simple model which does not fit the data in all respects can be useful.

In a directed network, the links have an associated direction, so that they point from one node to another. A directed network can be represented by an asymmetric adjacency matrix, and the directionality of links between groups can be modelled through the parameter $\phi$ by the existence of asymmetric interactions between the groups such that $\phi_{k,\ell}\neq\phi_{\ell,k}$. This double the number of link probability parameters $\phi$. The rest of the model is unaffected, except for the likelihood which must now be evaluated not for each pair of nodes but for each ordered pair of nodes. This extension of the infinite relational model assigns different probabilities to links in each direction between each pair of clusters, but has only a single parameter for the link probability within each cluster—thus, directionality is not modelled within clusters.

A bipartite network is defined as a set of links between two disjoint sets of nodes, possibly with different cardinality. The adjacency matrix for a bipartite network can thus be non-square. We can then use two independent Chinese restaurant processes to model the clustering of the two sets of nodes,

and change the likelihood to (cf. Eq. (29)

This latter parameterization is also useful for the modeling of directed networks when the groupings of the nodes may be different for the rows and columns of the adjacency matrix.

In a weighted network, each edge has a (scalar) weight associated with it. Depending on the type of weights, the Bernoulli likelihood can be changed to some other suitable distribution: For example, if the weights are positive integers [36], a Poisson distribution could be employed,

where $\lambda$ is the rate parameter for the edge weights, playing the role of $\phi$ in the Bernoulli model, c.f. Eq. (29). As a prior over $\lambda$, the typical choice is a Gamma distribution, replacing the Beta priors for $\phi$. If the weights are real numbers [18] an observational model based on a Normal distribution might be appropriate,

Here, we have two sets of parameters, $\mu$ and $\sigma^{2}$, denoting the means and variances of the edge weights between nodes in groups $i$ and $j$. Again, appropriate priors for $\mu$ and $\sigma^{2}$ should be selected.

Sometimes the data consists of multiple observations of networks on the same set of nodes (see [26, 34, 36, 3]). The only required change to the model is that the likelihood should be evaluated as the product of the likelihoods for each observed network. It can then either be assumed that the clustering structure as well as the link probabilities are equal across the multiple networks, that the clustering structure is shared but the link probabilities only shared according to an additional clustering of the multiple graphs, or that the clustering is shared but each network has an individual set of link probabilities,$\phi$. When the link probabilities are analytically marginalized this leads to three different expressions for the marginal likelihood.

To demonstrate the non-parametric Bayesian modeling framework in practise, we analyzed three real networks: {LaTeXdescription}

Zachary’s Karate club is an undirected unweighted network of friendships between 34 members of a karate club at a US university in the 1970s [59]. A total of 74 undirected links between the members of the Karate club are observed. In the analysis the standard IRM model was used.

The only complete connectome currently recorded of an organism is the directed integer weighted network of the 8,799 connections between the 302 neurons of the Caenorhabditis Elegans. The network has been compiled in [57]. In the analysis the weighted IRM model with a Poisson likelihood and Gamma priors was used.

The drugs and side effects network is a bipartite network on marketed medicines and their recorded adverse drug reactions extracted from public documents and package inserts. The network currently consists of 996 drugs and 4,199 side effects with 100,049 unweighted links between drugs and side effects [29]. In the analysis the bipartite IRM model was used. These three networks in turn represent three important complex network application domains within social science, neuroscience and bio-informatics.

The parameters of the models were inferred by Markov chain Monte Carlo sampling such that 250 iterations were used as burn in for the sampler and 250 iterations for drawing samples from the posterior. To improve mixing the data was analyzed based on 5 randomly initialized runs. In addition to a Gibbs sampler (as described in Fig. 5) the socalled split-merge sampler described in [26, 24] was also employed. The hyper-parameters for the Beta distribution we set to $a=b=1$. The posterior distribution of the number of components was computed. For assessing model fit by prediction of missing links, excluded 10% of links and an equivalent number of non-links in the analysis of the Zachary’s Karate Club data and 5% of links and an equivalent number of non-links in the two larger Connectome and Drug-side effects networks. For posterior predictive checking, we stored every 25th posterior sample and generated 20 replicated networks for each sample for each of the five random initializations. From the ensemble of these networks the distribution over the network characteristics; degree mean, degree standard deviation, characteristic path length and clustering coefficient were calculated and compared to the true values of these quantities for the actual network.

The results of the modeling is given in Fig. 7. The figure illustrates the network as well as a permutation of the networks adjacency matrix. The nodes are color coded according to the partition given by the sample with highest posterior likelihood across the five random initializations. From the permuted adjacency matrices it can be seen that the nodes of the networks have been grouped into clusters that share similar patterns of interactions, defining regions of network homogeneities. These blocks are color coded according to the expected value of the corresponding group interactions using a logarithmic gray scale. On the right, the posterior distribution of the number of components is shown as well as the models performance in predicting held-out links. The link prediction performance is quantified by the area under curve (AUC) of the receiver operator characteristic (ROC) [34]. In addition, the results of the posterior predictive checking of the models ability to account for the mean and standard deviation of the degree distribution, characteristic path length and clustering coefficient are given. Since the degree is explicitly modelled in the IRM model, this posterior predictive check serves only as a sanity check: The IRM should by definition get this right except for a small bias due to the prior.

From these results it can be seen that the IRM model accounts well for all the considered characteristics in the Zachary Karate Club network but that it poorly accounts for the degree standard deviation, the clustering coefficient, and the characteristic path length of the connectome of C. Elegans. As expected, the average degree falls within the lower tails of the simulated distributions for all the estimated models, but they underestimate the standard deviation of the node degree of both the connectome and drugs and side effect networks. This highlights a deficiency of the IRM model, namely that it does not explicitly model the degree distribution.

While the IRM is adept in identifying blocks of homogeneous network regions with the $\phi$ (nuisance) parameter specifying the density of each of these blocks, it does not explicitly account for microscale properties such as triangles and node degree. Hence, the clustering coefficient, characteristic path length and the standard deviation of the degree distribution is not well accounted for by the model as is evident in the posterior predictive checks. Despite these limitations, the infinite relational model does well account for mesoscale structure in the networks as quantified by its ability to predict links: For all the three networks the infinite relational model is able to predict links significantly better than random guessing.

Apart from being able to predict links the IRM model has made the structure of the networks substantially more comprehensible by reducing the complex network of pairwise interactions to a much smaller number of groups (defined by $\bar{z}$) and their interactions (defined by $\phi$). The IRM can therefore be considered an efficient framework for compressing a large complex network to a smaller network constituting consistent patterns of interactions between groups of nodes which can substantially facilitate in the understanding of mesoscale network patterns. For example, the analysis of the Zachary’s Karate Club network, the infinite relational model reveals six groups of club members including two large groups and two singletons (actually the posterior has support for five to eight groups, so these other configurations should also be considered in the interpretation of the results). It is known from the literature [59] that the karate club later split into two fractions, corresponding to the two large groups, led by the president and the instructor which are the two singletons.

Next, we compare the IRM model to several other methods on a set of social networks derived from a study of intra-organizational relations: {LaTeXdescription}

This set of undirected networks [9] consists of two types of relations defined on the same set of nodes corresponding to employees in a consultancy company. Links in the first network signifies employees who iteract whereas links in the second network signifies that either of the employees thinks that the other has expertise in an area important to her. The networks were generated by thresholding and symmetrizing the original directed weighted networks [9]. The two networks are highly correllated since employees would be expected to interact frequently with colleagues with important expertise.

We used the first of the two networks for training and examined the model fit by assessing the posterior predictive distribution of the node degree distribution. We fit an IRM model as well as two other non-parametric Bayesian network models, the infinite multiple membership relational model (IMRM) and the Bayesian community detection model (BCD) which are discussed further in the sequel. These models were fit using MCMC with $10,000$ rounds of Gibbs sampling where the first half of the samples were discarded for burn-in. Furthermore, we fit an exponential random graph model (ERGM) using the network statistics sociality and gwdegree [35] as well as a latent position and cluster model (ERGMM) [28] using a latent space of dimension four and six latent clusters (varying these parameters gave similar results). To compare how well the models fit the data we plotted the posterior predictive distribution of the degree distribution (see Figure 8). The results show that the two most flexible models, the ERGM and the IMRM fit the data very well in terms of reproducing the degree distribution. The fit of the IRM and BCD models which are both simple latent cluster models is less good: Both models appear to underestimate the number of nodes with a high degree, i.e., employees interacting with more than 15 colleagues. The ERGMM model on the other hand appears to overestimate the number of nodes with degrees around 15–20.

Next, we compared the models’ predictive performance by evaluating their ability to predict links in the second network (see Figure 8). Here, all models except the ERGM performed on par, suggesting that the inclusion of latent variables in the model is beneficial for this task.

Although, the IRM model is very flexible in terms of the structure it is able to account for, specialized non-parametric latent class methods have been proposed that specifically aim at extracting specific types of network structures. These models can be characterized by the restrictions which they impose on the between-class interactions $\phi$.

In [21] the $\phi$ matrix is constrained to only include two parameters, a within-group link probability $\rho_{w}$ and a between-group link probability $\rho_{b}$ such that

In [36] the within-group link probabilities are individual for each group but between-group probabilities are shared for all combinations of groups,

Both of the models mentioned above are inspired by the notion of communities defined as

“the organization of vertices in clusters, with many edges joining vertices of the same cluster and comparatively few edges joining vertices of different clusters.” [11]

This definition is used explicitly in [36] forming the Bayesian Community Detection (BCD) method. The BCD is based on the following non-parametric generative model that strictly enforces community structure by constraining the diagonal elements of the $\phi$ matrix to be larger than the off-diagonal elements. The generative model for BCD is given by

According to the model the probability of a link between communities $k$ and $\ell$ is strictly smaller than $w_{k,\ell}$ defined as the minimum over the two communities of some number $\gamma_{k}$ times the within-community probability $\phi_{k}$. This is enforced by generating the between-group probabilities according to an incomplete Beta distribution (BetaInc). The parameters $\gamma_{k}$ define a relative gap between link probabilities within and between communities, such that $\gamma_{k}=1$ says that there should be fewer links (on average) between than within, and $\gamma_{k}=0$ says that no links can be generated from nodes in community $k$ to other communities. The gap parameter $\gamma_{k}$ can in turn be learned from data and used to define the extend to which networks are community structured.

In [23] the IRM model was extended to handle irrelevant data entries by letting these entries constitute a separate noise cluster forming the subset infinite relational model (SIRM). The generative model for SIRM can be written as

For each node, the binary variable $r_{i}=0$ indicates that the node belongs to the noise cluster. For all pairs of nodes $(i,j)$ not in the noise cluster the model is identical to the IRM model; however, links between pairs of nodes of which at least one is in the noise cluster are generated with a shared probability $\rho$.

All the above extensions can potentially improve on identification of structure in complex networks by substantially reducing the parameters space of the within and between group interaction matrix $\phi$ compared to the IRM model. The above extensions are illustrated in Figure 9.

In [34] the binary matrix factorization model [32] based on an IBP is considered for network data. The following generative model embodies the latent feature relational model (LFRM)

where $\sigma[x]$ is a sigmoid function such as the the logit or probit. This model is inspired by the IRM in its parameterization but the model admits the nodes to belong to multiple groups, i.e., for each node to possess multiple features.

In [43] the infinite latent attribute model (ILAM) is proposed in which each of the nodes have a number of associated binary feature, and within each feature the nodes belong to an individual subcluster. The model can be summarized by the generative process,

For each feature $m$, the nodes that possess that feature are clustered according to a CRP. Here, $s$ is a bias term and $c^{(m)}_{i}$ is the cluster assignment of the $i$th node in the $m$th latent feature.

Both the LFRM and ILAM have been demonstrated to perform better than the IRM on a variety of link-prediction tasks [34, 43]. An important property of these models is that they allow for the membership of nodes in one group to inhibit the probability of linking to nodes in other groups as $\phi$ may include negative (i.e. antagonistic) elements. This property may indeed be an important reason for the models’ superior link prediction performance over IRM.

In [37, 38] the infinite multiple-membership relational model (IMRM) was proposed. Here the probability of observing a link between vertex $i$ and $j$ is generated independently given the (multiple) groups that vertex $i$ and $j$ belongs to and their interactions $\phi$. The generative model for the IMRM is given by

If, for example, node $i$ possesses feature $k$ and node $j$ possess feature $\ell$ the quantity $\phi_{k,\ell}$ denotes the probability of a link being generated between node $i$ and $j$ on account of that pair of features. The expression $1-\prod_{k,\ell}(1-\phi_{k,\ell})^{z_{i,k}z_{j,\ell}}$ defines the probability of observing a links between vertex $i$ and $j$ as the total probability of one or more of the pairs of features possessed by the two nodes to independently generate the link. This construction is referred to as a “noisy or process”. Notably, the IRM model is recovered when nodes belong to one and only one group.

Contrary to the LFRM and ILAM, the IMRM scales computationally in the number of observed links in the network rather than the number of potential links in the network which admits large scale analysis (see [38] for the details). However, scalability comes at the price of not being able to model antagonistic interactions between groups as for LFRM and ILAM. The LFRM and IMRM are illustrated in Figure 10.

The IBP is useful for defining non-parametric representations of binary latent variable models and both the LFRM and ILAM can be considered non-parametric latent variable models within the exponential random graph formulation. One approach for model order selection within framework of exponential random graph models is to impose sparse priors. The IBP can here be considered a non-parametric sparse prior for latent variable modeling in general as also proposed for factor analysis in [27]. As such, the IBP works in a similar manner as a slab-and-spike type prior, where a feature is either present or not according to the IBP while its contribution if present can be drawn separately. This can be used to extend existing sparse latent variable models within the exponential random graph model framework to form non-parametric models.

For instance, a non-parametric version of a latent factor model [20] can be defined by the following generative process using the IBP as a non-parametric sparsity promoting prior.

In [7] the perhaps most simple non-parametric model for hierarchical organization is proposed. This model imposes a uniform prior over all binary trees, which in the following we refer to as $\mathrm{UBT}$. The probability of generating a link between two nodes is defined by a parameter located at the level of their nearest common ancestor in the binary tree. A model for network with $N$ nodes thus has $N-1$ such parameters associated with each of the internal nodes in the tree. The generative model for the hierarchical random graph is given by

where $t_{i,j}$ denotes the index the nearest common ancestral node of vertex $i$ and $j$. In [47] a related generative model for binary hierarchies is proposed where each edge in the tree has an associated weight that defined the propensity in which the network complies with the given split.

One way to view the hierarchical random graph models is by first considering the top level of the hierarchy. Here the set of nodes is split into two partitions, and a single parameter is assigned to model the probability of observing a link between nodes in the two partitions. Next, the process continues recursively on the two partitions until each node is in a partition for itself. This framework was generalized and extended to the Mondrian process [48] which can be seen as a distribution over a $k$-dimensional tree. Used as a prior in a non-parametric Bayesian model of a bipartite network, at the top level the Mondrian process splits either of the two sets of nodes (chosen by random) into two partitions and continues this random bisectioning of the nodes until a stopping criterion is met. Parameters are then assigned to model the probability of links between each of the resulting pairs.

In [19] the uniform prior over binary trees of [7] where replaced by a uniform prior over multifurcating trees and the leafs of the trees rather than terminating at each vertex of the graph terminate at the levels of clusters generated from a CRP based on the following generative model

Here $K_{\bar{z}}$ denotes the number of clusters in $\bar{z}$ and $\mathrm{UT}$ defines a uniform prior over multifurcating trees. A benefit of this model is that it can be used to detect the presence of hierarchical structure as it includes the IRM model in its hypothesis space defined by a split at the root of the tree directly into all K clusters (i.e. forming a flat hierarchy). The model of [7] can on the other hand be considered the special case where the CRP only generates singleton clusters while the tree structure is strictly binary. As the leafs terminates in clusters rather than singletons the complexity of the model is in general substantially reduced compared to the models of [7, 47, 48] while the CRP defines the level at which to terminate the tree.

In [52] the Gibbs fragmentation tree was used as prior over multifurcating trees terminating at the vertex level of the network according to the following generative model

The Gibbs fragmentation tree is closely related to the two parameter nested Chinese restaurant process [2] differing in explicitly accounting for the occurrence in the nested CRP of trivial non-splits. The Gibbs fragmentation tree has several attractive properties. It is i) exchangeable in that the distribution does not depend on the labelling of the leaf nodes, ii) Markovian in that a subtree of the full tree is in turn a Gibbs fragmentation tree, and iii) consistent in that marginalizing over all leafs not considered in the subtree has the same distribution as only considering the Gibbs fragmentation tree of the subtree, see also [31, 52]. Apart from these attractive properties, the Gibbs fragmentation tree gives explicit control of the prior over multifurcating trees by its two parameters $\alpha$ and $\beta$, that makes it possible to bias the model toward deep vs. flat hierarchies. The probability of a given Gibbs fragmentation tree can be calculated using a simple recursive formula, see also [31, 52].