Dynamic degree-corrected blockmodels for social networks: a nonparametric approach

TL;DR

This paper introduces a nonparametric Bayesian method for modeling social networks with degree heterogeneity and community detection, adaptable to both static and dynamic networks, using Dirichlet processes and Gibbs sampling.

Contribution

It proposes a flexible, data-driven nonparametric approach combining degree correction and community detection in social networks, extendable to dynamic settings.

Findings

Effective community detection in real-world networks.

Automatic determination of number of communities and clusters.

Demonstrated applicability to dynamic social networks.

Abstract

A nonparametric approach to the modeling of social networks using degree-corrected stochastic blockmodels is proposed. The model for static network consists of a stochastic blockmodel using a probit regression formulation and popularity parameters are incorporated to account for degree heterogeneity. Dirichlet processes are used to detect community structure as well as induce clustering in the popularity parameters. This approach is flexible yet parsimonious as it allows the appropriate number of communities and popularity clusters to be determined automatically by the data. We further discuss some ways of extending the static model to dynamic networks. We consider a Bayesian approach and derive Gibbs samplers for posterior inference. The models are illustrated using several real-world benchmark social networks.

Figures

Equations

$$(G(A_1),\dots ,G(A_r))\sim \text{Dirichlet}(\alpha G_0(A_1),\dots ,\alpha G_0(A_r))$$

$${\theta }_i|G\stackrel{\text{iid}}{\sim }G\text{where}G|\alpha ,G_0\sim \text{DP}(\alpha ,G_0).$$

$$y_{ij}|p_{ij}\stackrel{\text{indep}}{\sim }\text{Bernoulli}(p_{ij}),$$

$${\mu }_{ij}={\theta }_i+{\theta }_j+\sum _{k=1}^K{\beta }_k^{\ast }\text{mathbbm}1\{z_i=z_j\},$$

$${\theta }_i|G$$

$$G$$

$${\beta }_i|H$$

$$H$$

$$y_{t,ij}|p_{t,ij}\stackrel{\text{indep}}{\sim }\text{Bernoulli}(p_{t,ij}).$$

$${\mu }_{t,ij}={\theta }_{it}+{\theta }_{jt}+\sum _{k=1}^K{\beta }_k^{\ast }\text{mathbbm}1\{z_i=z_j\}.$$

$${\theta }_{it}|G$$

$$G$$

$${\mu }_{t,ij}=\eta y_{t-1,ij}\text{mathbbm}1\{t>1\}+{\theta }_i+{\theta }_j+\sum _{k=1}^K{\beta }_k^{\ast }\text{mathbbm}1\{z_i=z_j\},$$

$$p(c|\alpha )p(\alpha )p(z|\nu )p(\nu )p({\theta }^{\ast })p({\beta }^{\ast })\prod _{i<j}p(y_{ij}|{\zeta }_{ij})p({\zeta }_{ij}|c_i,c_j,{\theta }^{\ast },z_i,z_j,{\beta }^{\ast }),$$

$$\text{P}(z_i=k|z_{-i},\nu )=\frac{m_{-i,k}}{n-1+\nu }\text{for k∈z−i},\text{P}(z_i\ne z_j\text{ for all }j\ne i|z_{-i},\nu )=\frac{\nu }{n-1+\nu },$$

$$\text{P}(c_i=\ell |c_{-i},\alpha )=\frac{n_{-i,\ell }}{n-1+\alpha }\text{for ℓ∈c−i},\text{P}(c_i\ne c_j\text{ for all }j\ne i|c_{-i},\alpha )=\frac{\alpha }{n-1+\alpha }.$$

$$\begin{gathered}\text{P}(z_{i}=k|\text{rest})=am_{-i,k}\exp\Big{\{}\beta_{k}^{*}\sum_{j\neq i:\;z_{j}=k}(\zeta_{ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}})-\frac{m_{-i,k}}{2}{\beta_{k}^{*}}^{2}\Big{\}}\\ \text{for $k\in z_{-i}$ and }\text{P}(z_{i}\neq z_{j}\text{ for all }j\neq i|\text{rest})=a\nu,\end{gathered}$$

$$\begin{gathered}\text{P}(c_{i}=\ell|\text{rest})=bn_{-i,\ell}\exp\Big{\{}\theta^{*}_{\ell}\sum_{j\neq i}(\zeta_{ij}-\theta^{*}_{c_{j}}-Z_{ij}^{T}\beta^{*})-\frac{n-1}{2}{\theta^{*}_{\ell}}^{2}\Big{\}}\\ \text{for $\ell\in c_{-i}$ and }\text{P}(c_{i}\neq c_{j}\text{ for all }j\neq i|\text{rest})=b\alpha\frac{\sigma_{c}}{\sigma_{\theta}}\exp\Big{\{}\frac{\mu_{c_{i}}^{2}}{2\sigma_{c}^{2}}\Big{\}},\end{gathered}$$

$$p(c|\alpha )p(\alpha )p(z|\nu )p(\nu )p({\theta }^{\ast })p({\beta }^{\ast })\prod _{t=1}^T\prod _{i<j}p(y_{t,ij}|{\zeta }_{t,ij})p({\zeta }_{t,ij}|c_{it},c_{jt},{\theta }^{\ast },z_i,z_j,{\beta }^{\ast }),$$

$$\text{P}(c_{it}=\ell |c_{-it},\alpha )=\frac{n_{-it,\ell }}{nT-1+\alpha }\text{for ℓ∈c−it},$$

$$\text{P}(c_{it}\text{ is not equal to any value in }{c}_{-it}|c_{-it},\alpha )=\frac{\alpha }{nT-1+\alpha },$$

$$\begin{gathered}\text{P}(z_{i}=k|\text{rest})=am_{-i,k}\exp\Big{\{}\beta_{k}^{*}\sum_{j\neq i:z_{j}=k}\sum_{t}(\zeta_{t,ij}-\theta^{*}_{c_{it}}-\theta^{*}_{c_{jt}})-\tfrac{Tm_{-i,k}}{2}{\beta_{k}^{*}}^{2}\Big{\}}\\ \text{for $k\in z_{-i}$ and }\text{P}(z_{i}\neq z_{j}\text{ for all }j\neq i|\text{rest})=a\nu,\end{gathered}$$

$$\begin{gathered}\text{P}(c_{it}=\ell|\text{rest})=bn_{-it,\ell}\exp\Big{\{}\theta^{*}_{\ell}\sum_{j\neq i}(\zeta_{t,ij}-\theta^{*}_{c_{jt}}-Z_{ij}^{T}\beta^{*})-\tfrac{n-1}{2}{\theta^{*}_{\ell}}^{2}\Big{\}}\\ \text{for $\ell\in c_{-i}$ and }\text{P}(c_{it}\neq\text{any value in $c_{it}$}|\text{rest})=b\alpha\frac{\sigma_{c}}{\sigma_{\theta}}\exp\Big{\{}\frac{\mu_{c_{it}}^{2}}{2\sigma_{c}^{2}}\Big{\}},\end{gathered}$$

$$p(y,\zeta,z,\beta^{*},\nu,c,\theta^{*},\alpha)=p(c|\alpha)p(\alpha)p(z|\nu)p(\nu)p(\theta^{*})p(\beta^{*})p(\eta)\\ \times\prod_{i<j}\Big{\{}\Big{[}\prod_{t\geq 1}p(y_{t,ij}|\zeta_{t,ij})\Big{]}p(\zeta_{1,ij}|c_{i},c_{j},\theta^{*},z_{i},z_{j},\beta^{*})\Big{[}\prod_{t\geq 2}p(\zeta_{t,ij}|c_{i},c_{j},\theta^{*},z_{i},z_{j},\beta^{*},\eta,y_{t-1,ij})\Big{]}\Big{\}},$$

$$\begin{gathered}\text{P}(z_{i}=k|\text{rest})=am_{-i,k}\exp\Big{\{}\beta_{k}^{*}\sum_{t}\sum_{j\neq i:z_{j}=k}(\tilde{\zeta}_{t,ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}}),-\frac{Tm_{-i,k}}{2}{\beta_{k}^{*}}^{2}\Big{\}}\\ \text{for $k\in z_{-i}$ and }\text{P}(z_{i}\neq z_{j}\text{ for all }j\neq i|\text{rest})=a\nu,\end{gathered}$$

$$\begin{gathered}\text{P}(c_{i}=\ell|\text{rest})=bn_{-i,\ell}\exp\Big{\{}\theta^{*}_{\ell}\sum_{t}\sum_{j\neq i}(\tilde{\zeta}_{t,ij}-\theta^{*}_{c_{j}}-Z_{ij}^{T}\beta^{*})-\tfrac{T(n-1)}{2}{\theta^{*}_{\ell}}^{2}\Big{\}}.\\ \text{for $\ell\in c_{-i}$ and }\text{P}(c_{i}\neq c_{j}\text{ for all }j\neq i|\text{rest})=b\alpha\frac{\sigma_{c}}{\sigma_{\theta}}\exp\Big{\{}\frac{\mu_{c_{i}}^{2}}{2\sigma_{c}^{2}}\Big{\}},\end{gathered}$$

$$\sum _{i<j}|\text{mathbbm}{1}_{\{c_i^{\ast }=c_j^{\ast }\}}-S_{ij}|,$$

$$p({\zeta }_{ij}|\text{rest})$$

$$\displaystyle\propto\mathbbm{1}\{\zeta_{ij}>0\}^{y_{ij}}\mathbbm{1}\{\zeta_{ij}\leq 0\}^{1-y_{ij}}\exp\Big{\{}-\frac{1}{2}\big{[}\zeta_{ij}^{2}-2\zeta_{ij}(\theta^{*}_{c_{i}}+\theta^{*}_{c_{j}}+Z_{ij}^{T}\beta^{*})\big{]}\Big{\}}.$$

$$\therefore\text{P}(z_{s}=k|\text{rest})=a^{\prime}m_{-s,k}\exp\Big{\{}-\frac{1}{2}\sum_{(i,j)\in\mathcal{P}_{s}}(\zeta_{ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}}-\beta_{k}^{*}\mathbbm{1}\{z_{i}=z_{j}=k\})^{2}\Big{\}}\;\text{for $k\in z_{-s}$},$$

$$\text{P}(z_s\ne z_j\text{ for all }j\ne s|\text{rest})$$

$$\displaystyle=a^{\prime}\nu\int\exp\Big{\{}-\frac{1}{2}\sum_{(i,j)\in\mathcal{P}_{s}}(\zeta_{ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}}-\beta_{k}^{*}\mathbbm{1}\{z_{i}=z_{j}=k\})^{2}\Big{\}}\frac{1}{\sqrt{2\pi}\sigma_{\beta}}\exp\Big{\{}-\frac{{\beta_{k}^{*}}^{2}}{2\sigma_{\beta}^{2}}\}\;d\beta_{k}^{*}$$

$$\displaystyle=a^{\prime}\nu\exp\Big{\{}-\frac{1}{2}\sum_{(i,j)\in\mathcal{P}_{s}}(\zeta_{ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}})^{2}\Big{\}},$$

$$\text{P}(c_i|\text{rest})$$

Full text

1 Introduction

Social networks play a central role in the dissemination of information (Westerman et al., 2014), formation of alliances (Gulati, 1998), transmission of disease (Cauchemez et al., 2011) and many other areas. It is thus important to understand the underlying structure of a social network and the behavioral patterns in the interactions. A common characteristic of social networks is that they often exhibit community structure, where certain groups of nodes (representing the social actors) are more densely connected within each group than across groups. The community structure may be present due to various factors such as similar interests, social stature or physical locations. Studying the nodal attributes associated with the communities can provide a greater understanding of the network topology, behavior patterns and network dependent processes such as epidemic spreading. However, identifying the community structure in a network can be challenging as the number of communities is typically unknown and the communities can vary in size and rate of interaction. Moreover, results can be distorted if the broad degree distributions often observed in real networks are not taken into account (Karrer and Newman, 2011). In this article, we propose a non-parametric approach to community detection in social networks by using independent Dirichlet processes (DP Ferguson, 1973) to capture the blockstructure in the social network and induce clustering in the activity level of nodes.

The partitioning of nodes into structurally equivalent groups, such that nodes in the same group relate with other nodes in the exact same way, was first discussed by Lorrain and White (1971), followed by White et al. (1976). Building upon the work of Holland and Leinhardt (1981) and Fienberg and Wasserman (1981), (Holland et al., 1983) generalized this deterministic concept and formulated stochastic blockmodels to allow for data variability. In a stochastic blockmodel, the nodes of a network are partitioned into groups and the distribution of ties between the nodes depends only on the group membership of the nodes and the probabilities of interactions between different groups. The stochastic blockmodel is generative and a wide variety of network structures, such as community, hierarchical or core-periphery, can be produced through different choices of the probability matrix. In a priori blockmodeling, exogenous actor attribute data are used to partition the nodes, while the discovery of blockstructures from relational data is referred to as a posteriori blockmodeling (Wasserman and Anderson, 1987). Snijders and Nowicki (1997) studied a posteriori blockmodeling for undirected networks when there are only two groups and derived procedures for finding the blockstructure using both maximum likelihood and Bayesian estimation. Nowicki and Snijders (2001) extend their approach to directed valued networks where the number of classes is fixed and address the nonidentifiability problem of the class labels. Handcock et al. (2007) consider a different clustering approach based on latent space models (Hoff et al., 2002), which posits that the probability of a tie is dependent on the positions of the actors in some unobserved space and decreases with distance.

The stochastic blockmodel has been extended in many ways. To overcome the restriction that each actor can only belong to one group, Airoldi et al. (2008) develop mixed membership stochastic blockmodels (MMSB), where each node is associated with a membership vector describing the probability of the node belonging to each of the groups. Each node can also assume different group membership when interacting with different nodes. Latouche et al. (2011) considers overlapping stochastic blockmodels, where each node can belong simultaneously to multiple groups with independent probabilities. The infinite relational model introduced by Kemp et al. (2006) allows the number of groups to be determined automatically by the data by drawing the membership vector from a Chinese restaurant process (CRP Pitman, 2006). A brief review on the CRP is given in Section 2.

Karrer and Newman (2011) note that the stochastic blockmodel often yield poor fits to real-world networks whose degree distributions are much broader than that generated by the stochastic blockmodel. To account for heterogeneity in the degree of nodes, they propose degree-corrected stochastic blockmodels, which modify the probability of a tie between node $i$ in group $g_{i}$ and node $j$ in group $g_{j}$ from $\omega_{g_{i}g_{j}}$ to $\theta_{i}\theta_{j}\omega_{g_{i}g_{j}}$, where $\omega_{rs}$ denotes the probability of a tie between group $r$ and $s$ while $\theta_{i}$ measures the activity level or “popularity” of node $i$. Estimates of the parameters are derived using maximum-likelihood and they demonstrate that the degree-corrected blockmodel leads to improved community detection. Gopalan et al. (2013) consider a related “assortative MMSB with node popularities” model, that considers a logit link and extends the MMSB to incorporate node popularities. A stochastic variational inference algorithm (Hoffman et al., 2013) is developed for posterior inference. Peng and Carvalho (2014) consider degree-corrected stochastic blockmodels using a Bayesian approach and a logistic regression formulation. Posterior inference is obtained via data augmentation with latent Pólya-Gamma variables and a canonically mapped centroid estimator that addresses label non-identifiability.

In this article, we focus on degree-corrected stochastic blockmodels for community detection in undirected social networks using a non-parametric Bayesian approach. The static model is formulated using probit regression and Dirichlet processes are used to capture the communities in the network and induce clustering among the popularity parameters. This approach is highly flexible yet parsimonious as it does not require the number of communities to be fixed in advance and instead allows the appropriate number of communities and popularity clusters to be determined automatically by the data. Our model integrates the approach of Kemp et al. (2006) who uses the CRP to detect community structure and that of Ghosh et al. (2010), who use the DP to induce clustering among the “productivity” and “attractiveness” parameters of a variation of the $p_{1}$ model (Holland and Leinhardt, 1981) and a social relations model (Gill and Swartz, 2007). While Ghosh et al. (2010) implements Bayesian analysis using WinBUGS software (Lunn et al., 2000) by considering a truncated DP (Ishwaran and Zarepour, 2000), we derive a Gibbs sampler for posterior inference. We also discuss several ways in which the static model can be extended to dynamic networks and illustrate the applicability of the proposed models using benchmark social networks.

This article is organized as follows. We review the Dirichlet process briefly in Section 2 and describe the proposed models in Section 3. We first present a model for static networks and then discuss extensions of this model to dynamic networks. In Section 4, we describe how posterior inference for the proposed model can be obtained using Gibbs samplers. In Section 5, we use the proposed models to analyze three real-world social networks. We conclude with a discussion of future research directions in Section 6.

2 The Dirichlet process

The Dirichlet process (DP, Ferguson, 1973) is widely used in Bayesian nonparametric models, particularly in DP mixture models as a prior over distributions. Let $(\Theta,\mathcal{B})$ be a measurable space with $G_{0}$ a probability measure on the space and $\alpha$ a positive real number. A random probability measure $G$ is distributed as a DP with base distribution $G_{0}$ and concentration parameter $\alpha$, written $G\sim\text{DP}(\alpha,G_{0})$, if

[TABLE]

for every finite measurable partition $A_{1},\dots,A_{r},$ of $\Theta$. The base distribution $G_{0}$ is the mean of the DP and $\alpha$ describes the concentration of mass around $G_{0}$. The larger $\alpha$ is, the more the DP will concentrate mass around $G_{0}$. Suppose the random variables $\{\theta_{i}|i=1,\dots,n\}$ are assigned the DP prior $G$. This is denoted as

[TABLE]

Note that $G$ is discrete with probability one even when $G_{0}$ is continuous. The stick-breaking construction of Sethuraman (1994) shows that $G=\sum_{l=1}^{\infty}\pi_{l}\delta_{\theta^{*}_{l}}$, where $\delta_{\theta}$ is a point mass concentrated at $\theta$ and $\theta^{*}_{l}\overset{\text{iid}}{\sim}G_{0}$. The random weights are defined by $\pi_{l}=V_{l}\prod_{j=1}^{l-1}(1-V_{j})$ where $V_{l}\overset{\text{iid}}{\sim}\text{Beta}(1,\alpha)$. The construction of the weights $\{\pi_{l}\}$ can be interpreted as starting with a stick of unit length and recursively breaking off a proportion $V_{l}$ of the remaining stick length. The stick-breaking construction shows clearly that $G$ is a discrete probability distribution. This implies that $\{\theta_{i}\}$ generated from $G$ will have non-negligible probability of having the same value. Thus, the DP will induce clustering among $\{\theta_{i}\}$ such that within each cluster, the $\theta_{i}$’s will assume the same value.

Another metaphor on the DP is provided by the Chinese restaurant process (CRP, Aldous, 1985), which likens $\theta_{i+1}$ to the $(i+1)^{\text{th}}$ customer entering a Chinese restaurant with infinitely many tables. The customer either selects an occupied table with probability proportional to the number of customers sitting there or a new table with probability proportional to $\alpha$. The CRP illuminates the “rich gets richer” phenomenon of the DP where larger clusters grow faster.

3 Non-parametric models for social networks

Let $N=\{1,\dots,n\}$ be the set of $n$ actors of interest and $y=[y_{ij}]$ be a $n\times n$ adjacency matrix where $y_{ij}$ is an indicator of a link from actor $i$ to actor $j$. In this article, we focus on undirected networks without self-links and multiple links. Hence $y$ is symmetric and the diagonal elements of $y$ are zeros. When the network of interest is observed at multiple (discrete) time points, $T$, we let $y_{t}=[y_{t,ij}]$ be the $n\times n$ adjacency matrix representing the state of the network at a time $t$ for $t=1,\dots,T$.

3.1 Static model

First we introduce a model for a static network $y$ that aims to detect community structure while incorporating heterogeneity through node-specific popularity parameters. For $1\leq i<j\leq n$, we assume that

[TABLE]

and introduce latent variable $\zeta_{ij}|\mu_{ij}\overset{\text{indep}}{\sim}N(\mu_{ij},1)$ where

[TABLE]

and $y_{ij}|\zeta_{ij}=1$ if $\zeta_{ij}>0$ and 0 if $\zeta_{ij}\leq 0$. Here we are considering a probit link where $\Phi^{-1}(p_{ij})=\mu_{ij}$ and $\Phi(\cdot)$ denotes the cumulative distribution function of the standard normal. The parameter $\theta_{i}$ represents the popularity or activity level of actor $i$, $K\leq n$ denotes the total number of groups or communities in the network and $z_{i}\in\{1,\dots,K\}$ represents the group membership of actor $i$. The coefficient $\beta_{k}^{*}$ measures the rate of interaction in the $k$th community. Members within a community are assumed to interact with each other at a common rate. A high $\beta_{k}^{*}$ indicates a tight or close-knit community where members interact at a high rate while a low $\beta_{k}^{*}$ indicates a group with little interaction. The third term on the right-hand side of (1) resembles a stochastic blockmodel where non-diagonal entries of the probability matrix are set to a common value (not necessarily zero). In (1), the probability of interaction, $p_{ij}$, between actors $i$ and $j$ depends on their individual popularities as well as the interaction rate of their community if they belong to the same community. An interaction between actors from different communities is driven only by their popularities. Thus the presence of a link can be explained by homophily in terms of community membership or popularities and the popularity parameters $\{\theta_{i}\}$ and community assignments $\{z_{i}\}$ are competing to explain the observed network.

For model parsimony, a DP is used to induce clustering among the popularity parameters $\{\theta_{i}\}$. We assume

[TABLE]

where the base distribution $G_{0}$ is $N(0,\sigma_{\theta}^{2})$ and $\alpha\sim\text{Gamma}(a_{\alpha},b_{\alpha})$. Let $\theta^{*}=[\theta_{1}^{*},\dots,\theta_{L}^{*}]^{T}$ denote the set of unique values among $\{\theta_{1},\dots,\theta_{n}\}$ and let $c_{i}$ indicate the latent class associated with $\theta_{i}$ so that $\theta_{i}=\theta_{c_{i}}^{*}$.

To detect the communities in the network, we consider another DP, $H$, which is independent of $G$. We introduce a $\beta_{i}$ for each actor $i$ where $\beta_{i}=\beta_{z_{i}}^{*}$ and assume

[TABLE]

where $H_{0}$ is $N(0,\sigma_{\beta}^{2})$ and $\nu\sim\text{Gamma}(a_{\nu},b_{\nu})$. Let $\beta^{*}=[\beta_{1}^{*},\dots,\beta_{K}^{*}]^{T}$ be the set of unique values among $\{\beta_{1},\dots,\beta_{n}\}$.

In this non-parametric approach, the number of clusters $L$ among $\{\theta_{i}\}$ and the number of communities $K$ are not fixed in advance. Instead, they are random and to be inferred from the data. The prior distribution of $L$ depends on the concentration parameter $\alpha$, with a larger $\alpha$ implying a larger $L$ a priori. To avoid overfitting and for greater interpretability, $\alpha$ will typically be small relative to $n$. Uncertainty about $L$ can be expressed by placing a prior on $\alpha$ and we consider a Gamma prior here. The relation between $K$ and $\nu$ is similar.

Next, we propose some ways of extending the static model to model dynamic networks. Suppose we observe networks $y_{t}=[y_{t,ij}]$ for $t=1,\dots,T$. For the dynamic models below, we assume that for $t=1,\dots,T$, $1\leq i<j\leq n$,

[TABLE]

As before we consider the probit link function and introduce the latent variable $\zeta_{t,ij}|\mu_{t,ij}\sim N(\mu_{t,ij},1)$ for 1$\leq i<j\leq n$, $t=1,\dots,T$ such that $y_{t,ij}|\zeta_{t,ij}=1$ if $\zeta_{t,ij}>0$ and 0 if $\zeta_{t,ij}\leq 0$. Thus, $p_{t,ij}=\Phi(\mu_{t,ij})$.

3.2 Dynamic model 1

Dynamic model I assumes that the community memberships remain unchanged over time but the popularities of the actors can vary with time. This assumption is appropriate for data where the communities arise due to factors that do not or are unlikely to vary drastically over time, for instance, gender, race, physical locations and job positions. In such cases, the changes in ties may be attributed to variations in the activity levels of individual nodes. For $1\leq i<j\leq n$ and $t=1,\dots,T$, let

[TABLE]

In resemblance of the static model, we assume that the $\{\theta_{it}\}$ are independent and induce clustering among them using a DP,

[TABLE]

where $G_{0}$ is $N(0,\sigma_{\theta}^{2})$ and $\alpha\sim\text{Gamma}(a_{\alpha},b_{\alpha})$. For this model, let $\theta^{*}=[\theta_{1}^{*},\dots,\theta_{L}^{*}]^{T}$ denote the set of unique values among $\{\theta_{11},\dots,\theta_{nT}\}$ and $c_{it}$ indicate the latent class associated with $\theta_{it}$ so that $\theta_{it}=\theta_{c_{it}}^{*}$ for $i=1,\dots,n$, $t=1,\dots,T$. The $\{\beta_{k}^{*}\}$ and $\{z_{i}\}$ are modeled using a DP as described in (3).

3.3 Dynamic model II

Dynamic model II extends the static model by allowing the tie between nodes $i$ and $j$ at time $t$ to depend on the existence of the tie at the previous time point. It assumes that the popularities and community memberships of the actors remain unchanged over time. For $1\leq i<j\leq n$ and $t=1,\dots,T$, let

[TABLE]

where $\eta\sim N(0,\sigma_{\eta}^{2})$. The coefficient $\eta$ can be interpreted as a measure of the persistence of ties in the network once they are formed. A positive $\eta$ implies that a tie between two actors is more likely to be present at time $t$ if a tie was present at time $t-1$ than if it were not, conditional on their popularities and community memberships. On the other hand, a negative $\eta$ would imply that a tie is more likely to be present at time $t$ if the tie was absent at the previous time point than if it were present. The parameters $\{\theta_{i}\}$ are modeled as in (2) and $\{z_{i}\}$ and $\{\beta_{k}^{*}\}$ are modeled as in (3). The popularities and communities inferred from this model smooths out the noise in the data and provide an overview of the behavior of actors over time.

4 Posterior inference

We use Gibbs samplers to derive posterior inference for the proposed models. To obtain the updates in the Gibbs sampler, we derive the posterior distribution of each variable conditional on the rest. Detailed derivations are given in the Appendix. Sampling from the DP is performed using the methods described in Neal (2000) while the concentration parameters $\alpha$ and $\nu$ are sampled using the method described in Escobar and West (1995).

For ease in representation, we introduce the following notations. Let $Z_{ij}$ be a binary vector of length $K$ where the $k$th element is 1 if $z_{i}=z_{j}=k$, 0 otherwise, and $Z=[Z_{12},Z_{13},\dots,Z_{(n-1),n}]^{T}$ be a $n(n-1)/2\times K$ matrix. We also define $\zeta_{ij}=\zeta_{ji}$ and $\zeta_{t,ij}=\zeta_{t,ji}$ for $1\leq j<i\leq n$ and $t=1,\dots,T$. Let $\mathcal{S}_{m}=\{(i,j)|i<j,c_{i}=c_{j}=m\}$, $\mathcal{S}_{t,m}=\{(i,j)|i<j,c_{it}=c_{jt}=m\}$, $\mathcal{P}_{m}=\{(i,j)|j\neq i,c_{i}=m,c_{j}\neq m\}$, and $\mathcal{P}_{t,m}=\{(i,j)|j\neq i,c_{it}=m,c_{jt}\neq m\}$. We use $\text{TN}(x|\mu,\sigma,a,b)$ denote the truncated normal distribution with density $\frac{1}{\sigma}\phi(\frac{x-\mu}{\sigma})/(\Phi(\frac{b-\mu}{\sigma})-\Phi(\frac{a-\mu}{\sigma}))$, where $\phi(\cdot)$ denotes the density of the standard normal. In the algorithms presented below, we use $K$ and $L$ to represent the current number of communities and popularity clusters respectively at each iteration and $\beta^{*}=[\beta^{*}_{1},\dots,\beta_{K}^{*}]$ and $\theta^{*}=[\theta^{*}_{1},\dots,\theta_{L}^{*}]$ to represent the states currently associated with the clusters.

For the static model, the joint distribution $p(y,\zeta,z,\beta^{*},\nu,c,\theta^{*},\alpha)$ is given by

[TABLE]

where $\zeta=\{\zeta_{11},\dots,\zeta_{n-1,n}\}$, $c=\{c_{1},\dots,c_{n}\}$ and $z=\{z_{1},\dots,z_{n}\}$.. Note that $p(c|\alpha)$ and $p(z|\nu)$ are defined as (Neal, 2000):

[TABLE]

where $z_{-i}=z\backslash z_{i}$, $c_{-i}=c\backslash c_{i}$, $m_{-i,k}=\sum_{z_{j}\in z_{-i}}\mathbbm{1}\{z_{j}=k\}$ and $n_{-i,\ell}=\sum_{c_{j}\in c_{-i}}\mathbbm{1}\{c_{j}=\ell\}$. The Gibbs sampler for the static model is outlined in Algorithm 1. In step 2, suppose that the number of distinct values in $z_{-i}$ is $K^{\prime}$. In the update, $z_{i}$ can either assume one of these $K^{\prime}$ distinct values or a new value not assumed by any $z_{j}\in z_{-i}$. (5) describe these probabilities and $a$ is a constant that ensures these $K^{\prime}+1$ probabilities sum to one. The same idea applies to (6)–(10),

For dynamic model I, the joint distribution $p(y,\zeta,z,\beta^{*},\nu,c,\theta^{*},\alpha)$ is given by

[TABLE]

where $\zeta=\{\zeta_{1,11},\dots,\zeta_{T,n-1,n}\}$, $c=\{c_{11},\dots,c_{nT}\}$ and $z=\{z_{1},\dots,z_{n}\}$. Note that $p(c|\alpha)$ is defined as

[TABLE]

where $c_{-it}=c\backslash c_{it}$ and $n_{-it,\ell}$ be the number of indicators in $c_{-it}$ that are equal to $\ell$. The definition of $p(z|\nu)$ remains as in (4). The Gibbs sampler for dynamic model I is outlined in Algorithm 2.

For dynamic model II, the joint distribution is given by

[TABLE]

where $\zeta=\{\zeta_{1,11},\dots,\zeta_{T,n-1,n}\}$, $c=\{c_{1},\dots,c_{n}\}$ and $z=\{z_{1},\dots,z_{n}\}$. Note that $p(c|\alpha)$ and $p(z|\nu)$ are as defined in (4). The Gibbs sampler for dynamic model II is outlined in Algorithm 3.

We code Algorithms 1–3 in Julia and all experiments are run on an Intel Core i5 CPU @ 3.30GHz, 8.0GB RAM. We note that it is also possible to use software such as OpenBUGS to obtain posterior inference for the proposed models by considering a truncated DP approach Ishwaran and Zarepour (2000). However, we observe that the runtime in OpenBUGS is significantly longer than Julia especially as the number of nodes increases. For the examples, we initialize multiple MCMC chains from random starting points and make use of diagnostic plots to check for convergence.

4.1 Cluster Analysis

Given the sample of clusterings from the MCMC output, we can assess clustering by computing the posterior similarity matrix $S$, which is a $n\times n$ symmetric matrix whose $(i,j)$ entry contains the posterior probability that actors $i$ and $j$ belong to the same cluster. This probability is estimated by the proportion of times actors $i$ and $j$ cluster together and it is not affected by the problem of “label-switching” (labels associated with clusters may change during MCMC runs, see, e.g. Stephens, 2000) or the number of clusters varying across iterations.

We can also compute a single (hard) clustering estimate, for instance, by using the maximum a posteriori (MAP) approach or methods based on the posterior similarity matrix or Rand index (see discussion in Fritsch and Ickstadt, 2009). Here we consider the Binder’s loss function (Binder, 1978), which is defined as the total number of disagreements between the estimated and true clustering among all pairs of actors. The R package mcclust provides a function, minbinder, that can be used to find the clustering $c^{*}=[c_{1}^{*},\dots,c_{n}^{*}]$ that minimizes the posterior expectation of this loss. The posterior expected loss can be written as

[TABLE]

where the sum is taken over all possible pairs of actors and $S_{ij}$ is the $(i,j)$ entry of the posterior similarity matrix.

5 Applications

We investigate the performance of the static and dynamic models on three well-known social network datasets and compare the fitted models with results obtained previously by other authors. The first is a karate club network studied by Zachary (1977) from 1970 to 1972, the second is a dolphins social network (Lusseau et al., 2003) and the third is a dataset collected by Kapferer (1972) at an African clothing factory in Zambia. These datasets are available at the UCI Network Data Repository (https://networkdata.ics.uci.edu/).

5.1 Karate club network

This dataset contains 78 undirected friendship links among 24 members, which are constructed based on interactions outside of club activities. Due to disputes over the price of karate lessons, the club was divided informally into two factions, led by the karate instructor “Mr Hi” (actor 1) and the president “John A.” (actor 34) respectively (these names are pseudonyms). During the study, the club eventually split into two separate clubs when Mr Hi was fired for trying to raise lesson fees unilaterally and his supporters left to join the new club formed by Mr Hi. All members joined clubs following their own factions except actor 9, who crossed factions to join Mr Hi’s club because he was only three weeks away from a test for black belt at the time of the split and he could not bear to give up his rank.

We fit the static model to this dataset using Algorithm 1. Three chains were run in parallel, each consisting of 40,000 iterations with the first 30,000 discarded as burn-in. The total runtime is 172 seconds. A thinning factor of 5 was applied and the remaining 6000 samples were used for posterior inference. We set $a_{\nu}=b_{\nu}=a_{\alpha}=b_{\alpha}=5$ and $\sigma_{\theta}^{2}=\sigma_{\beta}^{2}=1$. Figure 1 shows the posterior distributions of the number of communities ($K$), the number of popularity clusters ($L$), and the DP concentration parameters $\alpha$ and $\nu$.

The mode of $K$ is 3 and that of $L$ is 4. The fitted model is quite parsimonious with a relatively small number of clusters for both popularity and community. Figure 2 shows the posterior similarity matrices for the clusterings according to community (left) and popularity (right). Each element in the matrix represents the proportion of times that the actors concerned belong to the same cluster.

Figure 3 plots the posterior mean of $\theta_{i}$ against the degree for each actor. While the factional leaders, Mr Hi (actor 1) and John A. (actor 34), and a few other actors {2, 3, 33} have high popularity, the rest of the members have much lower activity levels generally.

Using the similarity matrices, we compute hard clustering estimates using Binder’s loss function. There are three communities, one of which contains a single node {3} and three popularity clusters. Figure 4 shows plots of the karate club network where nodes of the same color belong to the same cluster and singletons are not colored. We run Algorithm 1 again, fixing $z$ and $c$ to obtain estimates of $\beta^{*}$ and $\theta^{*}$ for these clusterings. The conditional posterior mean and standard deviation (in brackets) of these parameters for each cluster are shown in the legend of Figure 4.

There are three clusters for the popularity parameters $\{\theta_{i}\}$, the first contains {Mr Hi, John A., 3}, the second contains {2, 33} and the third contains all remaining members. For the communities, we note that the $\beta_{k}^{*}$ for groups 1 and 3 are strongly positive, indicating a high interaction rate within each group. The posterior mean and standard deviation of $\beta_{k}^{*}$ for the singletons necessarily equal that of the prior distribution. Note that group 3 corresponds exactly to the faction led by John A. as concluded in Zachary (1977) while group 1 together with the singleton {3} correspond to the faction led by Mr Hi. From the posterior probability matrix, actor 3 has a posterior probability of about 0.4 of being clustered together with members in group 1 (Mr Hi’s faction) and a probability of about 0.05 of being clustered together with members in group 3 (John A.’s faction). It is thus reasonable to combine actor 3 with group 1. Hence, our proposed static model is able to identify members in the factions accurately.

Incidentally, if we drop $\{\theta_{i}\}$ from the static model and consider just the blockmodel, we obtain five clusters, four of which are singletons: {1}, {3}, {33}, {34} and the fifth cluster contains all other members. This result is similar to the phenomenon discussed in Karrer and Newman (2011), who note that the non-degree-corrected blockmodel with $K=2$ splits the network into high-degree and low-degree nodes instead of by factions, while the degree-corrected version splits it according to factions albeit with one misclassification. In addition, Bickel and Chen (2009) observe that the non-degree-corrected blockmodel with $K=4$ splits the network according to factions correctly after the merging of sub-communities. These observations highlight the importance of accounting for degree variation in blockmodels as well as the difficulties in determining an appropriate number of clusters. Our static model tries to address these issues using a non-parametric approach via the automatic clustering structures induced by the DP. We observed that the clusters identified by the static model can be sensitive to the DP concentration parameters in some cases. For example, if we adopt a more conservative prior, say by setting $a_{\nu}=b_{\nu}=a_{\alpha}=b_{\alpha}=10$, then we obtain three communities, the first corresponds to the faction led by John A., the second contains {5, 6, 7, 11, 17} and the third contains all remaining members. Here, the second cluster emerges as one with a higher interaction rate than the third. However, merging the second and third clusters still yields Mr Hi’s faction.

While the clustering estimates return hard partitions of the network which are easy to interpret, the posterior similarity matrices reveal finer details regarding the degree of affiliation of actors towards the clusters that they are assigned to in the hard split. For the posterior similarity matrix for popularities,there are two main blocks but the partitioning among actors {1, 34, 3, 33, 2} is not so straightforward. For the posterior similarity matrix for communities, actor 10 is assigned to the cluster led by John A., but he has a somewhat lower posterior probability ($\sim$ 0.5) of being together with the other members in this cluster than the rest, and also has some posterior probability ($\sim$ 0.2) of being in the same cluster as members in Mr Hi’s faction.

5.2 Dolphins social network

Lusseau et al. (2003) constructed an undirected social network describing the associations among a community of 62 bottlenose dolphins living off Doubtful Sound, New Zealand after observing them for seven years from 1994–2001. This dataset has been widely studied in community detection, see for instance, Lusseau and Newman (2004) and Cao et al. (2015). In this network, the nodes represent dolphins and the ties represent higher than expected frequency of being sighted together. Of the 62 dolphins, 33 are males, 25 are females and the gender of the remaining 4 are unknown.

We apply Algorithm 1 to this network, using 15,000 iterations with a burn-in of 5000 iterations and a thinning factor of 5 in each chain. Three chains were run in parallel and the total runtime is 250 seconds. We set $a_{\nu}=b_{\nu}=a_{\alpha}=b_{\alpha}=10$ and $\sigma_{\theta}^{2}=\sigma_{\beta}^{2}=1$. The marginal posterior distributions of $K$, $L$, $\nu$ and $\alpha$ are shown in Figure 5. The posterior of $K$ is concentrated on larger values as compared to $L$ and $K$ has a mode of 7 while the mode of $L$ is 2.

The posterior similarity matrices in Figure 6 show the community and popularity clustering structure in this network. Around five communities can be seen in the matrix on the left while the right matrix shows faint outlines of two clusters.

Next we use Binder’s loss to obtain clustering estimates for the community structure and popularity clusterings based on the MCMC samples. This yields 16 communities and a single popularity cluster. Of the 16 communities, 9 are singletons so there are essentially only 7 communities. We run Algorithm 1 again, fixing $c$ and $z$ to obtain estimates of $\beta^{*}$ and $\theta^{*}$ for these clusterings. The estimate of $\theta^{*}$ is $-0.92\pm 0.03$. Figure 7 shows the observed dolphins social network where the nodes are labeled with the names of the dolphins, and males, females and dolphins of unknown gender are represented using squares, circles and triangles respectively. Nodes of the same color belong to the same community while the singletons are not colored. From Figure 7, most of the singletons can be regarded as peripherals (e.g. Zig, TR82, Quasi, MN23); they have few links and lie at the margins of the network. While some of them can clearly be pushed into certain clusters, others such as SN89 lie at the edge of different groups. The estimate of $\beta_{k}^{*}$ is indicative of the rate of interaction for each group $k$ and this is shown in the legend along with the standard deviation in brackets. Groups 1–3 and 6–7 represent close-knit communities while groups 4–5 have low within-group interaction rates.

Previously, Lusseau and Newman (2004) studied the community structure of this dolphins network by using a clustering algorithm which is based on removing links with high “betweeness” measures to extract the groupings (Girvan and Newman, 2002). They also investigated the role that gender and age homophily played in the formation of communities. They concluded that there are 2 main communities and 4 sub-communities; the first sub-community matches group 1 exactly, the second matches group 2 together with the singletons (Zig, TR82, Quasi, MN23), the third matches groups 7 and 4 combined and the fourth matches groups 3 and 6 combined plus the singletons TR120, TR88, TSN83, Zipfel and SN89. We note that the posterior similarity matrix does suggests some of these combinations. Thus, the communities detected by Algorithm 1 agree largely with the results of Lusseau and Newman (2004) and also that of Cao et al. (2015). In addition, Figure 7 also provides some evidence of assortative mixing by sex. For example, group 6 consists almost entirely of females while groups 2 and 7 are composed of mainly males.

5.3 Kapferer’s tailor shop network

Kapferer (1972) collected data on the interactions among 39 workers in a tailor shop in Zambia, Southern Africa, from June 1965 to February 1966, and he examined how these social networks relate to major events taking place in the factory. The workers’ duties can be classified into eight categories: head tailor (worker number 19), cutter (16), line 1 tailor (1–3, 5–7, 9, 11–14, 21, 24), button machiner (25–26), line 3 tailor (8, 15, 20, 22–23, 27–28), ironer (29, 33, 39), cotton boy (30–32, 34–38) and line 2 tailor (4, 10, 17–18). These positions require different levels of skills and some like the head tailor, cutter, line 1 tailors and button machiners were perceived as having more prestige. Here we focus on the symmetric “sociational” networks (based on convivial interactions) recorded at two time points, the first was before an aborted strike and the second was after a successful strike for higher wages. The network at the second time point (223 edges) is much denser than the first (158 edges) as the workers strive to be more united (thereby expanding their social relations) in their efforts to change the wage system. This dataset has been widely studied, for instance, by Mitchell (1989) and Nowicki and Snijders (2001) using block structures and Thiemichen et al. (2016) using Bayesian exponential random graph models.

5.3.1 Dynamic model I

First, we fit dynamic model I to the data using Algorithm 2. Dynamic model I assumes that the communities remain constant over time and that the emergence or dissolution of ties are due to changes in the activity level of individual actors. The hyperparameters are set as $a_{\nu}=b_{\nu}=a_{\alpha}=b_{\alpha}=10$ and $\sigma_{\theta}^{2}=\sigma_{\beta}^{2}=1$. We use three parallel chains, each with 15,000 iterations and the first 5000 iterations are discarded as burn-in. The total runtime is 139 seconds. A thinning factor of 5 was applied and posterior inferences are based on the remaining 6000 iterations. The posterior distributions of $K$, $L$, $\nu$ and $\alpha$ are shown in Figure 8. The mode of $K$ is 6 and the mode of $L$ is 4.

Next we compute the posterior similarity matrices and use Binder’s function to obtain hard clustering estimates. This yields nine communities and three popularity clusters. Of the nine communities, five are singletons so there are essentially only four communities. We run Algorithm 2 again, fixing $z$ and $c$ to obtain estimates of $\beta^{*}$ and $\theta^{*}$ for these clusterings. The results are shown in Figure 9, and the mean and standard deviation (in brackets) of $\beta^{*}$ and $\theta^{*}$ are reported for each group.

The first row shows the four communities which are constant across the two time points (the singletons (19–21, 25–26) are not colored). The shapes of the nodes represent the positions of the workers as explained in the legend. The plots indicate a high degree of job homophily in the communities even though these social networks are constructed based on casual interactions (“general conversation, the sharing of gossip and the enjoyment of a drink together”, Kapferer, 1972). In particular, groups 1 and 2 consists of workers with jobs perceived to be of higher prestige: cutter, line 1 and line 2 tailors, group 3 consists of line 3 tailors and group 9 consists of all the ironers and cotton boys. The estimates of $\beta_{k}^{*}$ are strongly positive, indicating a high interaction rate within each group.

There are three popularity clusters with increasing means, $-1.41$ (group 1), $-0.46$ (group 2) and 0.57 (group 3). Thus, we can consider the three clusters as representing “low”, “average” and “high” popularity. Actors 19 (head tailor) and 16 (cutter) are the only two actors with high popularity at $t=1$, and they maintained high popularity at $t=2$. This is not surprising since they are regarded by Kapferer (1972) to be in “supervisory” positions and play critical roles in the operation of the factory. From the barplot (left) in Figure 10, the number of workers with low popularity decreased from $t=1$ to $t=2$ while the number with average or high popularity increased. This reflects the efforts of the workers in expanding social ties after the first unsuccessful strike.

Examining the results more closely using the barplot (right) in Figure 10, the proportion of workers with low and average popularity actually remained unchanged over the two time points for the ironers, cotton boys and line 2 tailors (positions with lower prestige). The changes in popularity arise mainly from line 1 tailors, button machiners and line 3 tailors. In particular, two line 1 tailors, {21, 24} and a button machiner {25} moved from average to high popularity. These observations are consistent with the analysis of Kapferer (1972), who noted that line 1 tailors made a strong attempt to expand their links after the first unsuccessful strike as they stand to benefit the most from the change in wage system. Mitchell (1989) also noted that the button machiner, Meshak (actor 25) played a crucial role in the unfolding events at the factory and for the latter part was regarded as a supervisor by the factory owner.

5.3.2 Dynamic model II

Next, we fit dynamic model II to the data using Algorithm 3. In this model, the parameter $\eta$ provides an indication of the persistence of ties formed. The probability that a tie is formed at any time point depends on whether a tie exists at the previous time point as well as the community membership of the nodes and their popularities. Using 15,000 iterations, discarding the first 5000 as burn-in and applying a thinning factor of 5 for each of three independent chains, the total runtime is 106 seconds. Figure 11 shows the posterior distributions of $K$, $\nu$, $L$, $\alpha$ and $\eta$ based on 6000 MCMC samples. The modes of $K$ and $L$ are both 6. The posterior mean of $\eta$ is 0.58 and its posterior mass is concentrated on positive values. This indicates that a tie is likely to persist at the second time point given that it existed at the first time point. Figure 12 shows the posterior similarity matrices. We note that the block structures are not clear-cut.

Figure 13 shows the hard clusterings computed using Binder’s function and the estimates of $\beta^{*}$ and $\theta^{*}$ for these clusterings. The communities detected are largely similar to that of Dynamic model I except for changes to the assignment of individuals {14, 16, 19, 21}. A new “community” consisting of $\{19,21\}$ appeared. However, the $\beta_{k}^{*}$ estimate for this group is 0.54 with a large standard deviation of 0.8. Thus, this is not truly a “community” in the sense that there is a high interaction rate between the actors.

The number of popularity clusters increased from three in dynamic model I to six in model II. In model II, the popularity of an actor summarizes his activity level across all time points. Figure 14 shows how the mean of $\theta_{i}$ varies with the degree of an actor at each time point. The head tailor and the cutter have significantly higher popularity than the other workers, followed by actor 24 (Ibrahim) and actors in popularity group 2. We note that group 2 includes several individuals who play significant roles in the factory’s social relationships (Kapferer, 1972). These include Lyashi (11), who tried to win followers in support of his view of the factory structure, Hastings (13), who took on many supervisory duties of the cutter at time 2, Meshak (25), who was regarded as a leader by the factory owner, and Mubanga (34), an influential figure among unskilled workers.

6 Conclusion and future work

We present a non-parametric Bayesian approach for detecting communities in social networks, using degree-corrected stochastic blockmodels. In the proposed static model, the number of communities and popularity clusters does not have to be fixed in advance and is inferred from the data automatically through the use of the DP. For the karate club network and the dolphins social network, we find that the static model returns sensible results although there is some sensitivity to the DP concentration parameters. The inferred popularity clusters also summarizes the popularities of the actors and helps in the identification of key players in the network. We discuss two extensions of the static model to dynamic networks. Dynamic model I enables the study of the change in activity level of actors over the entire duration while dynamic model II provides a measure of the persistence of links formed in the network. While the Gibbs samplers are feasible for small networks, they do not scale well to large networks and more efficient methods of estimation, such as variational approximation methods, can be developed.

7 Acknowledgments

Linda Tan is supported by the National University of Singapore Overseas Postdoctoral Fellowship.

Appendix A Updates of static model

• •

For $i<j$,

[TABLE]

• •

For $s=1,\dots,n$, $p(z_{s}|\text{rest})\propto p(z_{s}|z_{-s},\nu)\prod_{(i,j)\in\mathcal{P}_{s}}p(\zeta_{ij}|\theta^{*}_{c_{i}},\theta_{c_{j}}^{*},\beta_{z_{s}}^{*})$.

[TABLE]

[TABLE]

where $a^{\prime}$ is a normalizing constant that ensures the probabilities sum to one. Hence we can simplify the expressions to that in (5).

• •

\!\begin{aligned} p(\beta^{*}|\text{rest})&\propto\exp\Big{\{}-\frac{1}{2}\sum_{i<j}(\zeta_{ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}}-Z_{ij}^{T}\beta^{*})^{2}\Big{\}}\exp\Big{(}-\frac{{\beta^{*}}^{T}\beta^{*}}{2\sigma_{\beta}^{2}}\Big{)}\\ &\propto\exp\Big{\{}-\frac{1}{2}\Big{(}{\beta^{*}}^{T}\Big{(}Z^{T}Z+\frac{1}{\sigma_{\beta}^{2}}\Big{)}\beta^{*}-2{\beta^{*}}^{T}\sum_{i<j}Z_{ij}(\zeta_{ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}})\Big{\}}.\\ \end{aligned}

Note that $Z^{T}Z=\sum_{i<j}Z_{ij}Z_{ij}^{T}$ is a $K\times K$ diagonal matrix where the $k$th diagonal element counts the number of pairs of $(z_{i},z_{j})$ that assume a common value $k$.

• •

For $i=1,\dots,n$,

[TABLE]

[TABLE]

• •

For $m=1,\dots,L$,

[TABLE]

Appendix B Updates of dynamic model I

• •

For $t=1,\dots,T$, $i<j$,

[TABLE]

• •

For $s=1,\dots,n$, $p(z_{s}|\text{rest})\propto p(z_{s}|z_{-s},\nu)\prod_{t}\prod_{(i,j)\in\mathcal{P}_{s}}p(\zeta_{t,ij}|c_{it},c_{jt},\theta^{*},\beta_{k}^{*})$.

[TABLE]

for $k\in z_{-s}$ and

[TABLE]

where $a^{\prime}$ is a normalizing constant to ensure probabilities sum to one. Hence we can simplify the expressions to (7).

• •

\begin{aligned} p(\beta^{*}|\text{rest})&\propto\exp\Big{\{}-\frac{1}{2}\sum_{t}\sum_{i<j}(\zeta_{t,ij}-\theta^{*}_{c_{it}}-\theta^{*}_{c_{jt}}-Z_{ij}^{T}\beta^{*})^{2}\Big{\}}\exp\Big{(}-\frac{{\beta^{*}}^{T}\beta^{*}}{2\sigma_{\beta}^{2}}\Big{)}\\ &\propto\exp\Big{\{}-\frac{1}{2}\Big{(}{\beta^{*}}^{T}TZ^{T}Z\beta^{*}-2{\beta^{*}}^{T}\sum_{i<j}Z_{ij}\sum_{t}(\zeta_{t,ij}-\theta^{*}_{c_{it}}-\theta^{*}_{c_{jt}})+\frac{{\beta^{*}}^{T}\beta^{*}}{\sigma_{\beta}^{2}}\Big{)}\Big{\}}.\end{aligned}

• •

For $i=1,\dots,n$, $t=1,\dots,T$,

[TABLE]

[TABLE]

• •

For $m=1,\dots,L$,

[TABLE]

Appendix C Updates of dynamic model II

• •

For $t=1,\dots,T$, $1\leq i<j\leq n$,

[TABLE]

• •

\begin{aligned} p(z_{s}|\text{rest})&\propto p(z_{s}|z_{-s},\nu)\prod_{t}\prod_{(i,j)\in\mathcal{P}_{s}}p(\zeta_{t,ij}|\theta^{*}_{c_{i}},\theta_{c_{j}}^{*},\eta,y,\beta_{z_{s}}^{*})\\ &\propto p(z_{s}|z_{-s},\nu)\exp\Big{\{}-\frac{1}{2}\sum_{t}\sum_{(i,j)\in\mathcal{P}_{s}}(\tilde{\zeta}_{t,ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}}-\beta_{k}^{*}\mathbbm{1}\{z_{i}=z_{j}=k\})^{2}\Big{\}}\end{aligned}

For $k\in z_{-s}$,

[TABLE]

and

[TABLE]

where $a^{\prime}$ is a normalizing constant to ensure probabilities sum to one. Hence we can simplify the expressions to (9)

• •

\begin{aligned} p(\beta^{*}|\text{rest})&\propto\exp\Big{\{}-\frac{1}{2}\sum_{t}\sum_{i<j}(\tilde{\zeta}_{t,ij}-\theta^{*}_{c_{it}}-\theta^{*}_{c_{jt}}-Z_{ij}^{T}\beta^{*})^{2}\Big{\}}\exp\Big{(}-\frac{{\beta^{*}}^{T}\beta^{*}}{2\sigma_{\beta}^{2}}\Big{)}\\ &\propto\exp\Big{\{}-\frac{1}{2}\Big{(}{\beta^{*}}^{T}TZ^{T}Z\beta^{*}-2{\beta^{*}}^{T}\sum_{i<j}Z_{ij}\sum_{t}(\tilde{\zeta}_{t,ij}-\theta^{*}_{c_{it}}-\theta^{*}_{c_{jt}})+\frac{{\beta^{*}}^{T}\beta^{*}}{\sigma_{\beta}^{2}}\Big{)}\Big{\}}.\end{aligned}

• •

For $i=1,\dots,n$,

[TABLE]

[TABLE]

• •

For $m=1,\dots,L$,

[TABLE]

• •

\begin{aligned} p(\eta&|\text{rest})\propto\exp\Big{\{}-\frac{1}{2}\sum_{t=2}^{T}\sum_{i<j}(\zeta_{t,ij}-\eta y_{t-1,ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}}-\beta^{T}Z_{ij})^{2}-\frac{\eta^{2}}{2\sigma_{\eta}^{2}}\Big{\}}\\ &\propto\exp\Big{\{}-\frac{\eta^{2}}{2}\Big{(}\frac{1}{\sigma_{\eta}^{2}}+\sum_{t=2}^{T}\sum_{i<j}y_{t-1,ij}^{2}\Big{)}+\eta\sum_{t=2}^{T}\sum_{i<j}y_{t-1,ij}(\zeta_{t,ij}-\theta^{*}_{c_{i}}-\theta^{*}_{c_{j}}-\beta^{T}Z_{ij})\Big{\}}.\end{aligned}