A Bayesian Nonparametric Latent Space Approach to Modeling Evolving Communities in Dynamic Networks

The evolution of communities in dynamic (time-varying) network data is a prominent topic of interest. A popular approach to understanding these dynamic networks is to embed the dyadic relations into a latent metric space. While methods for clustering with this approach exist for dynamic networks, they all assume a static community structure. This paper presents a Bayesian nonparametric model for dynamic networks that can model networks with evolving community structures. Our model extends existing latent space approaches by explicitly modeling the additions, deletions, splits, and mergers of groups with a hierarchical Dirichlet process hidden Markov model. Our proposed approach, the hierarchical Dirichlet process latent position clustering model (HDP-LPCM), incorporates transitivity, models both individual and group level aspects of the data, and avoids the computationally expensive selection of the number of groups required by most popular methods. We provide a Markov chain Monte Carlo estimation algorithm and apply our method to synthetic and real-world networks to demonstrate its performance.


Introduction
Many naturally occurring networks contain discrete changes in community structure. When high school students move across the country for college, old friendship groups often dissolve, leading the way for new friendship groups to form. After World War II, the Eastern and Western Blocs emerged and dominated the network of global alliances. However, after the fall of the Soviet Union, these blocs reshuffled into new political alliances. When exposed to external stimuli, regions of the brain activate before becoming dormant. By identifying these community-level phase changes, we can gain valuable insight into the rich processes that generate dynamic (longitudinal or time-varying) networks.
In this work, we address this problem of inferring discrete changes in a network's community structure. See Figure 1 for a concrete example. In this case, a dynamic Figure 1: Adjacency matrices from a dynamic sequence of networks that exhibit an evolving community structure. At t = 4 the groups split from two into six groups, then at t = 7 the groups merge into four communities. The rows and columns of each adjacency matrix have been sorted according to the ground truth labels. These networks were simulated according to the procedure described in Section 5.2.
network's community structure contains two groups during the first three time points, splits into six groups at the fourth time point, and then merges into four groups during the seventh time point. These additions, deletions, splits, and mergers of groups in dynamic networks is our primary object of interest. For brevity, we refer to these group-level dynamics as a network's evolving community structure. Note that this does not refer to individual actors (or nodes) moving between existing communities. While we will model such actor-level dynamics, our primary goal is to infer changes in community structure that occur on a larger macro scale.
To infer an evolving community structure in dynamic networks, we adopt the latent space approach to network modeling. Originally developed in Hoff et al. (2002), latent space models (LSMs) embed actors within a Euclidean space (the distance model) or a hypersphere (the projection model). Closeness in the latent space increases the probability that the two actors form an edge in the observed network. For this reason, one interprets the proximity of two actors in the latent space as an indication that they have similar characteristics. The LSM's popularity stems from the intuitive meaning of the embeddings and its ability to naturally incorporate desirable sociological features such as homophily, reciprocity, and transitivity. The distance model was extended to undirected dynamic networks by Sarkar and Moore (2006) and directed dynamics networks by Sewell and Chen (2015). The latent space model was first applied to the community detection problem in the case of a single static network by Handcock et al. (2007). Their proposed latent position clustering model (LPCM) uses a Gaussian mixture model to cluster the latent positions embedded according to the Euclidean distance model. The idea is that a cluster of actors in the latent space will correspond to a densely connected community in the network. This work was extended to dynamic networks by Sewell and Chen (2017). Their approach allows the actors to move between a static number of groups at each point in time. They accomplish this by endowing the actor's latent trajectories with autoregressive hidden Markov model (AR-HMM) dynamics.
The LPCM is an excellent model for community detection in dynamic networks, but it has two key shortcomings that make it inadequate for inferring an evolving community structure: (1) the number of communities must be set a priori so that uncertainty in the community structure cannot be assessed in a fully Bayesian way, and (2) the number of communities must remain static over time. The first shortcoming is important since multiple clusterings often fit the data well, which means we need a way to compare uncertainty across these partitions. The second shortcoming implies that existing LPCMs are unable to allow the number of communities to change over time so that inference over evolving communities is not possible. The current inability of LPCMs to handle time-varying groups in longitudinal networks motivates our investigation into a suitable Bayesian methodology.
In this work, we propose a Bayesian nonparametric model that addresses the shortcomings outlined in (1) and (2). As in Sewell and Chen (2017), we model the actor's trajectories through the latent space as a collection of independent and identically distributed (iid) AR-HMMs with Gaussian emissions. However, unlike their approach, we place a hierarchical Dirichlet process (HDP) prior on the Markov chain's transition distributions to allow for inference over the number of communities. Furthermore, we no longer assume that the hidden label Markov chain is homogeneous. Instead, we allow for time-inhomogeneous transition distributions. This modification provides us with the flexibility to model the additions, deletions, splits, and mergers of communities. Since our approach combines the advantages of the HDP with the LPCM for modeling dynamic networks, we refer to our model as the hierarchical Dirichlet process latent position clustering model (HDP-LPCM). To the best of our knowledge, there is no other latent space approach that accomplishes our modeling criteria.
The rest of this paper is organized as follows. Section 2 reviews other Bayesian nonparametric approaches for modeling evolving communities. In Section 3, we describe the hierarchical Dirichlet process latent position clustering model (HDP-LPCM), our proposed Bayesian nonparametric method for modeling evolving communities in dynamic networks. In Section 4, we elaborate on the details of our Metropolis-Hastings within Gibbs estimation procedure as well as our methodology for posterior summarization. We demonstrate the empirical performance of our proposed method through various simulation studies in Section 5. In Section 6, we use our model to analyze two realworld dynamic networks: a network of military alliances during the Cold War (Gibler, 2009b) and a narrative network constructed from the Game of Thrones television series (Beveridge and Chemers, 2018). Section 7 contains a discussion.

Related Methods
Although modeling evolving communities in dynamic networks is new to the latent space methodology, the statistical modeling of varying communities has a long history in the broader network literature. For brevity, we focus on Bayesian nonparametric approaches because they are most related to our work. We restrict the discussion to the three most popular varieties of statistical network models: stochastic block models, mixedmembership stochastic block models, and latent feature models.
For stochastic block models, Kemp et al. (2006) introduced the infinite relational model (IRM), which uses a Dirichlet process prior to infer the number of blocks in the traditional stochastic block model. Ishiguro et al. (2010) extended this work to longitudinal networks in their dynamic infinite relational model (dIRM). Similar to our work, the dIRM assumes the block memberships form a time-inhomogeneous Markov chain. In the mixed-membership literature, Fan et al. (2015) introduced the dynamic infinite mixedmembership stochastic block model. Their model also utilizes the hierarchical Dirichlet process framework to extend the original mixed-membership model (Airoldi et al., 2008;Fu et al., 2009;Ho et al., 2011) to a dynamic setting. In contrast to our method, they do not assume an HMM structure and only use the HDP to re-sample the mixed-membership vectors at every time step. Finally, in the latent feature modeling literature, Kim and Leskovec (2013) introduced the nonparametric multi-group membership model. They utilize a distance-dependent Indian buffet process (dd-IBP) to assign each actor a latent binary feature vector. Due to their nonparametric prior, they can infer an evolving number of binary features.
While some of these competing approaches have addressed shortcomings (1) and (2), LPCMs have many appealing advantages that make our extension worthwhile. Foremost is the LPCM's ability to capture network structures on multiple scales. Unlike the block models described above, which assume connections are purely governed by a global clustering, latent space models allow for heterogeneous connectivity patterns due to an actor's local position. Furthermore, unlike other methods, the latent space provides an interpretable visualization of the entire network. In practice, these visualizations are crucial for allowing domain experts to make qualitative assessments and critiques of the statistical methodology.

The HDP Latent Position Clustering Model
We consider binary relational data between n individuals recorded over T time periods. These relations are collected in a sequence of n × n binary adjacency matrix Y t , t = 1, . . . , T , where the entries y ijt indicate the presence (y ijt = 1) or absence (y ijt = 0) of an edge between individuals i and j at time t. For clarity, we only consider undirected random graphs without self-loops, so that Y t is a symmetric matrix with zeros on the diagonal. The extension of our model to directed networks or networks with weighted edges is straightforward and discussed in Section 7.
To each of the n individuals, we associate a latent (unobserved) position that may vary through time in a p-dimensional Euclidean latent space. We represent the latent position of individual i at time t with the vector-valued random variable X i t ∈ R p . In addition, we collect a snapshot of all individual latent positions at time t in the n × p matrix X t = (X 1 t , . . . , X n t ) . Like the traditional LPCM, we assign each actor to a latent group at each time point. Note that their assignment may change over time. However, unlike the LPCM, we assume that the number of groups changes over time to accommodate the network's evolving community structure. We use G to denote the total number of groups in the network over all observational periods. We refer to a group that contains at least one actor at time t as an active group at time t. We use G t ⊆ {1, . . . , G} to indicate the set of all active groups at time t. Note that under the assumption of an evolving community structure, G t is a random set and may grow or shrink over time. We encode the latent group membership of individual i at time t with the integer-valued random variable Z i t ∈ G t . The collection of group assignments for all individuals at time t is summarized by the n-dimensional vector Z t = (Z 1 t , . . . , Z n t ) . Common to the latent space literature, we assume that the latent community labels only play a role in the distribution of the latent positions. Formally, we assume that the formation of edges in the dynamic networks Y 1:T is conditionally independent of the actor labels Z 1:T given the latent positions X 1:T , i.e., Y 1:T ⊥ Z 1:T | X 1:T . Note that throughout the remainder of this work, we will use the notation A 1:K to refer to the sequence (A 1 , A 2 , . . . , A K ). This allows us to decompose the joint probability distribution as follows: P (Y 1:T , X 1:T , Z 1:T ) = P (Y 1:T | X 1:T ) P (X 1:T , Z 1:T ). (3.1) This independence assumption says that the probability of a tie is solely determined by the underlying latent positions of the actors. This decomposition is consistent with the latent space clustering idea described earlier. Specifically, the notion that an underlying cluster of actors in the latent space results in observed communities. We believe such a generative model is natural for modeling communities in networks. For example, friend groups in social networks often form due to the similar interests or characteristics of their members.

The Euclidean Distance Model
As in Hoff et al. (2002), we posit that the probability of an edge forming between actors increases as the Euclidean distance between actors decreases. Let d ijt = X i t − X j t 2 denote the Euclidean distance between actors i and j at time t. A conventional link between the conditional probability of forming an edge and d ijt is the logistic regression model: , (3.2) where η ijt is a linear predictor that depends on the distances d ijt and the vector ψ holds any additional parameters. A sequence of conditional independence assumptions results in the full network likelihood. First, we assume that the longitudinal networks are conditionally independent given the latent positions, i.e., Y t ⊥ Y s | X 1:T for all s, t ∈ {1, . . . , T }. Second, we posit that edges between actors form independently conditioned on their latent positions at each time point. Under these assumptions, the likelihood of the adjacency matrices factors as a product over the networks at each time point and the set of dyads D: In this work, we focus on undirected networks without self-loops for which 3.4) and the linear predictor takes the form (3.5) so that ψ = {β 0 }. Various versions of this likelihood have been proposed, such as in Sarkar and Moore (2006), Krivitsky et al. (2009), andChen (2015). The intercept parameter β 0 is sometimes referred to as the abundance. This is because higher values of β 0 result in a higher probability of forming edges.

An AR-HMM for Latent Space Dynamics
Recall that our goal is to endow the dynamics of the latent space with a probabilistic structure that allows for inference over an evolving collection of communities. To accomplish this goal, we must place an adequately flexible joint distribution over the latent space, P (X 1:T , Z 1:T ). To begin, we adopt the assumptions taken by Sewell and Chen (2017). Specifically, we assume that an actor's community assignments form a Markov chain, i.e., where D = denotes equality in distribution. Similarly, we assume an actor's latent position follows a Markov process with an additional dependence on the current community assignment, i.e., These two assumptions allow us to factor the marginal density of an individual actor's trajectory. In particular, we conclude that an actor's trajectory follows an independent autoregressive hidden Markov model (AR-HMM): where p(Z i 1 ) is the initial distribution over actor i's initial community assignment, p(Z i t | Z i t−1 ) is the transition distribution between actor i's label assignment at time t − 1 and time t, and p(X i t | X i t−1 , Z i t ) is the emission distribution for actor i's latent position at time t, which we allow to depend on the previously emitted value X i t−1 . To properly model evolving communities, we make an important departure from previous dynamic latent space models. In particular, we expand the probabilistic model to include time-inhomogeneous Markov chains where the transition matrix p(Z i t | Z i t−1 ) can vary over time. This is in contrast with previous methods, which assume that the probability for an actor to transition from community i to community j is the same at all times. We argue that time-inhomogeneous transitions are essential characteristics of evolving communities. For example, the addition of a group requires the transition matrices to add a non-zero probability to transition into that new group. Furthermore, once the new group is added to the network, an actor's probability to transition into that group may approach a different steady-state than when it initially appeared. Although models that utilize homogeneous Markov chains may infer changing groups through the inclusion of empty clusters, their inferences will certainly be biased due to smoothing over sudden changes in group structure.
Tractable inference over these time-inhomogeneous latent trajectories is made possible by tying together the transition distributions and emission distributions of every actor. We collect the transition probabilities into a single row-stochastic transition matrix Π t for each time step t = 2, . . . , T . Each entry is defined as (Π t ) jk = π t jk = p(Z i t = k | Z i t−1 = j), which is the same for all n actors. In what follows, it will be useful to isolate the j-th row of Π t in the vector π t j . Furthermore, we refer to {π t j } as the set of rows of the transition matrix at time t. For simplicity, we assume that all actors share a common initial state Z i 0 = 0. This allows us to view the initial distribution of the Markov chain as a special transition distribution denoted by π 1 0 = (π 1 01 , . . . , is equal for all actors. Conditioned on the initial and transition probability matrices as well as the parameters of the emission distributions θ, the data generating process of all the actor trajectories is characterized by a collection of iid AR-HMMs. In this work, we model the joint density over all actor trajectories as where I p is the p × p identity matrix, and N (X | µ, Σ) is the normal density with mean vector µ and covariance Σ evaluated at the point X. Note that in this case θ = {{µ g , σ 2 g } G g=1 , λ}. Like the clustering model in Handcock et al. (2007), the communities are modeled as a multivariate normal distribution in the latent space with mean location µ g , and spherical covariance σ 2 g I p . As in the longitudinal clustering approach of , the mean position X i t is equal to λµ g + (1 − λ)X i t−1 where λ ∈ (0, 1). This is a blend between the actor's previous position and the current assigned group location. Consequently, λ is called the blending coefficient.
To complete the model, we must specify a prior on the initial distribution π 1 0 , the transition probabilities Π 2:T , and the parameters of the emission distribution θ = {{µ g , σ 2 g } G g=1 , λ}. A parametric approach is taken in Sewell and Chen (2017). Their model restricts itself to time-homogeneous Markov chains and assumes a fixed number of communities throughout time. As previously explained, their model is inadequate for inference over an evolving community structure. They are neither able to quantify uncertainty in the number of groups nor account for time-inhomogeneous transitions. For this reason, we turn to the flexibility afforded by a Bayesian nonparametric approach.

Background: Dirichlet Processes and the HDP-HMM
The Dirichlet process (DP) is a distribution over discrete probability measures: where β ∼ GEM(γ) denotes the sticking-breaking process (Sethuraman, 1994): From now on we will use the notation G 0 ∼ DP(γ, H) to indicate draws from a DP with concentration parameter γ and base measure H. Since draws from a DP are discrete with probability one, the DP cannot be used as a general nonparametric prior over continuous densities. To extend the DP to continuous density estimation, one often uses G 0 as a mixing measure over some parametric class of distributions f θ . This construction is known as the DP mixture model. The sampling process begins by drawing indicator variables z i iid ∼ β, and generating observations as y i ∼ f θz i . Note that in this work, we adopt the convention that when β represents a distribution over the natural numbers, i.e., K j=1 β j = 1, then we use z ∼ β to mean z ∼ K j=1 β j δ(z, j), where δ(i, j) is the Kronecker delta.
The hierarchical Dirichlet process (HDP) is a distribution over a collection of discrete probability measures that share a common set of atoms: The group-specific distributions of the HDP, π j , are often used as priors over the rows of an infinite dimensional transition matrix, i.e., (Π) jk = π jk for j ∈ N. This formulation allows for the construction of the hierarchical Dirichlet process hidden Markov model (HDP-HMM) (Teh et al., 2006), which is a natural Bayesian nonparametric extension of the Bayesian HMM (Robert et al., 2000). The sampling mechanism of the HDP-HMM proceeds as follows: one samples the hidden states sequentially as z t ∼ π z t−1 , and the observations are linked to the global parameters via y t ∼ f θz t . Note that the sharing of atoms induced by the HDP prior allows the Markov chain to utilize a single global set of parameters, θ k , at all time points. Despite its popularity the original HDP-HMM struggles to model Markov chains with long state durations. To remedy this issue Fox et al. (2011b) introduced the sticky hierarchical Dirichlet process (sticky HDP): where δ j is a vector of zeros except for a single one at the j-th index. Analogous to the HDP, the sticky HDP can be used as a prior over the transition matrices of an HMM. In this case, the extra stickiness parameter κ > 0 biases the process towards self-transitions. As a result, Fox et al. (2011b) found that the corresponding sticky HDP-HMM better models the longer state durations found in real world applications. For this reason, we chose to use the sticky variant of the HDP as a prior over the transition matrices in Equation (3.9).

An HDP Prior for Time-Inhomogeneous Markov Chains
We now present our extension to the latent position clustering model that can infer an evolving community structure in dynamic networks. To accomplish this objective, we place the following Bayesian nonparametric prior on the transition matrices in Equation (3.9): (3.14) In this work, we take H to be a multivariate normal distribution and F to be an inverse gamma distribution, i.e., µ g For the initial distribution, we set κ = 0 and allocate it a separate concentration parameter α 0 to distinguish it from the other transition distributions. Note that, unlike the traditional HDP-HMM, our model re-samples the rows of the transition matrix, π t g , at each time step. This difference allows us to infer time-inhomogeneous Markov chains. As we previously argued, this time-inhomogeneity is crucial for modeling evolving communities. In short, our model extends the dynamic Euclidean distance model proposed in Sewell and Chen (2017) by adding a nonparametric HDP prior over time-inhomogeneous transition matrices. For this reason, we refer to our model as the hierarchical Dirichlet process latent position clustering model (HDP-LPCM). A graphical representation of the HDP-LPCM is depicted in Figure 2.
Although our prior borrows heavily from the sticky HDP-HMM formulation, there are three key differences between the state-space model utilized in the HDP-LPCM and the traditional sticky HDP-HMM. First, the emitted values X i t have an autoregressive dependence. This modification was originally explored outside the network setting by Fox et al. (2011a) to model switching linear dynamical systems with Gaussian observations. However, in our case, the observations are a sequence of binary adjacency matrices. The second difference is that the HDP-HMM contains a single state sequence Z 1:T , while we infer multiple state sequences Z i 1:T for i = 1, . . . , n from a single observed sequence of networks. The final departure from the HDP-HMM is the re-sampling of the rows of the transition matrix {π t g } at every time step. This difference introduces time-inhomogeneity into the Markov chain, which is crucial for modeling the additions, deletions, splits, and mergers of groups found in evolving communities.
In addition to properly modeling an evolving community structure, the other important accomplishment of our model is its ability to infer the number of communities from the data in a fully Bayesian way. A typical parametric approach (specifying |G t | at each time step) would require comparing a combinatorial amount of models. This task is Figure 2: The graphical model for the HDP-LPCM. The dependence on α and κ is only displayed for the last transition matrix, although it is present in all transition matrices from t = 2, . . . , T . This change is for clarity only. The HDP-LPCM consists of a collection of actor specific iid time-inhomogeneous AR-HMMs whose shared state space is inferred by an HDP. The latent positions are linked to the observed network through the Euclidean distance model of Hoff et al. (2002).
computationally infeasible for even a small number of groups and time steps. Through a nonparametric prior, we naturally incorporate model selection into our inference procedure, which avoids the computationally expensive model selection step found in many LPCMs (Handcock et al., 2007;Sewell and Chen, 2017).

The Weak-Limit Approximation
For inference, we utilize the weak-limit approximation to the HDP proposed in Fox et al. (2011b). The approximation replaces the infinite dimensional DPs with finite L-dimensional Dirichlet distributions as follows: (3.15) Practically, the weak limit approximation transforms an infinite transition matrix into a finite L × L matrix. The parameter L gives us control over the approximation with the accuracy increasing as L goes to infinity (Ishwaran and Zarepour, 2000). In practice, one often sets L larger than some a priori upper limit on the number of communities.

Estimation
We take a fully Bayesian approach to estimation. In what follows, we describe a Markov chain Monte Carlo (MCMC) method to sample from the posterior of the HDP-LPCM described in Section 3. We implement a Metropolis-Hastings within Gibbs MCMC scheme with the goal of identifying an evolving community structure consistent with the inferred posterior distribution.

Blocked Metropolis-Hastings within Gibbs Sampler
As outlined in Section 3, the joint distribution over all the variables (Equation (3.1)) factors as Equation (3.15) .
(4.1) To complete the model we assign the following priors: where N (0,1) (µ, σ 2 ) indicates the normal distribution with mean µ and variance σ 2 truncated to the range (0, 1). We can realize these samples by following a Gibbs sampling algorithm in which we iteratively sample from the appropriate conditional distributions of X 1:T , Z 1:T , π 1 0 , Π 2:T , β, µ 1:L , σ 2 1:L , λ, and β 0 . Due to our choice of priors, most conditional distributions are conjugate so that the Gibbs updates are derived in standard fashion. See Section S.1 of the supplementary materials for the details on these conditional distributions. Metropolis-Hastings (MH) steps are necessary for the latent positions and the intercept. In both cases, we use a normal random walk proposal. In all experiments, we tune the proposal step sizes using an initial tuning phase so that the proposed moves have a 25% -40% acceptance rate. The only samplers that need special care are the block sampler for the actor labels Z 1:T , the sampler for the global prior β, as well as the additional hyperparameter sampling schemes. We describe these samplers in the next sections. The full Metropolis-Hastings within Gibbs sampler is outlined in Algorithm 1.
In the following sections, we use dot notation to indicate summations over an index, e.g., for a random variable w ab , w ·b = a w ab , w a· = b w ab and w ·· = a b w ab . In addition, we use n kjt to denote the number of actors that transitioned from group k to group j at time t and n kt to denote the number of actors in group k at time t.
2. Update latent positions X 1:T via MH with a normal random walk proposal.
3. Update β 0 via MH with a normal random walk proposal.

Sampling Z 1:T
Although previous work on dynamic latent space clustering models (Sewell and Chen, 2017) update Z 1:T coordinate wise, we chose to utilize a blocked Gibbs sampler. In particular, we use a variant of the forward-backward algorithm for AR-HMMs (Rabiner, 1989) to jointly sample an actor's group assignments Z i 1:T . We made this change to increase the mixing rate of the sampler. Previous work has shown that sampling Z i 1:T as a block as opposed to coordinate wise can lead to drastically faster mixing rates in HMMs. This is especially the case when the elements of Z i 1:T are highly related in the posterior distribution, which is likely due to our stickiness assumption (Scott, 2002). The full sampling algorithm is detailed in Algorithm 2. A derivation of the sampler is provided in Section S.2 of the supplementary materials. A similar sampler is found in Fox et al. (2011a); however, they do not consider multiple state sequences with time-inhomogeneous transition probabilities as we do in this work.

Sampling β
We follow the sampling strategy outlined in Teh et al. (2006) and Fox et al. (2011b). These samplers use the Chinese restaurant franchise (CRF) metaphor for the hierarchical Dirichlet process to derive closed-form conditional distributions. In the CRF metaphor, each group j ∈ {1, . . . , J} in Equation (3.12) is associated with a restaurant. The entire collection of restaurants (or franchise) shares a common menu of dishes θ 1:K , where K is the total number of draws from the base measure H. The metaphor describes the probabilistic mechanism used by the HDP to partition a collection of individuals (or customers) pre-assigned to J restaurants into K clusters that are identified by their common dish assignments.
For completeness and to properly describe the auxiliary variable samplers in the next section, we present a brief description of the CRF metaphor. As customers arrive at a restaurant, the host seats them at an occupied table with probability proportional to the number of people already seated there and assigns them to a new table with a probability proportional to α. When the first customer arrives at a table, a single dish from the global menu is served. The probability to receive an existing dish is proportional to the number of times that dish is served across all restaurants, and a new dish is prepared with probability proportional to γ. For the sticky HDP, the dish assignment process is modified to account for the stickiness parameter κ. In this case, each restaurant has a specialty dish. Once the waiter receives a table's order, the waiter has a probability ρ = κ/(α + κ) to override that order and serve the restaurant's specialty dish instead.
In the HDP-LPCM, we identify each row of the transition matrices, {π t g }, with a restaurant in the franchise. In other words, each pair (g, t) indexes the restaurants in the CRF metaphor in a way that a bijective map exists between restaurant j → (g, t). Furthermore, the customers assigned to restaurant j correspond to all actors assigned to group g at time t−1. With these mappings in mind, we can directly apply the re-sampling algorithm for β described in Fox et al. (2011b) with a few straightforward modifications to account for the multiple hidden state sequences.

Sampling m,m, and w
The sampler for β and the hyperparameter samplers presented in the next section require the introduction of three auxiliary variables m,m, and w. These variables keep track of certain statistics of the CRF sampling process that when made available preserve conjugacy in the model. The components of these variables are as follows: m jks is the number of tables at restaurant j that ordered dish k at time s,m jks is the number of tables at restaurant j that were served dish k at time s, and w jts is a binary variable indicating whether table t at restaurant j at time s had their dish choice overridden by the waiter. Note that m jks andm jks may differ due to the waiter overriding a table's order. The algorithm derived to sampler these auxiliary variables is similar to the one found in Fox et al. (2011b) with modifications that account for the increased number of restaurants. We present our auxiliary variable sampler in Algorithm 3. Given n jks , the number transitions from state j to state k at time s, and ρ = κ/(α + κ), sample the auxiliary variables m,m, w as follows: 1. For each k ∈ {1, . . . , L}, set m 0k1 =m 0k1 = 0, and n = 0. For each customer in restaurant 0 eating dish k at time 1, that is for i = 1, . . . , n 0k1 , sample Increment n by one, and if x = 1 increment m 0k1 andm 0k1 by one.
Increment n by one, and if x = 1 increment m jks by one.
Set the number of tables in restaurant j ordering dish k at time s to: m jks = m jks , j = k, m jjs − w j·s , j = k.
Algorithm 3: Sampler for auxiliary variables needed for the update of the global transition distribution β. Note that we do not need to incorporate overridden dishes for m 0k1 because we draw π 1 0 from a DP(α 0 , β), which has no stickiness parameter κ.

Sampling Hyperparameters
The HDP-LPCM has many hyperparameters that must be set by the practitioner. If one has no a priori knowledge of appropriate hyperparameter settings, it can be advantageous to incorporate hyperparameter exploration into the sampling procedure. In this section, we detail additional hyperparameter samplers as well as our initial choices of hyperparameter values used in the applications of the following sections. The full hyperparameter sampler is presented in Algorithm S.2 in Section S.3 of the supplementary materials. In our experience, the hyperparameters that influence inference the most are τ 2 , b, γ, α 0 , α, and κ. The prior variance of the group means τ 2 controls the scale of the latent space. The prior scale of the group variances b controls the size of the clusters in the network. The concentration parameters γ, α 0 , and α influence the a priori number of clusters in the network with higher values indicating more groups. Finally, larger values of κ result in a priori higher probabilities of self-transitions. With this intuition about the role of these parameters, we proceed to describe our choice of hyperparameter priors and initial values.
Determining the scale of the latent space is a difficult task a priori due to the unobserved nature of the space. For this reason, we place the following priors on τ 2 and b: where Γ -1 (a, b) indicates an Inverse Gamma distribution with shape parameter a and scale parameter b, and Γ(c, d) represents a Gamma distribution with shape and scale parameters c and d respectively. Under these priors, the updates for the variance of the cluster means τ 2 and the scale parameter for the group variances b are conjugate. We use a heuristic similar to the latentnet package (Krivitsky and Handcock, 2008) to set the remaining hyperparameter values. In particular, we expect the size of the latent space to grow with the number of actors n. For this reason, we chose a τ and b τ such that τ 2 has mean E [τ 2 ] = n 2/p /50 and standard deviation 4E [τ 2 ], where n is the number of nodes and p is the dimension of the latent space. In addition, we chose b, c, and d such that the mode of σ 2 is τ 2 with a standard deviation of 4E [b]. Furthermore, the initial values for τ 2 and b are set to their prior mean values and a = 2. The idea behind this choice is to set a prior that can generate group shapes ranging from the scale of the entire latent space to a much smaller scale.
To update the hyperparameters of the HDP prior γ, α 0 , α, and κ, we utilize a straightforward modification of the samplers in Teh et al. (2006) and Fox et al. (2011b). Both samplers employ or extend the sampler developed in Escobar and West (1995), which we describe in Algorithm S.1 in Section S.3 of the supplementary materials. For all simulations and data applications, we place the following disperse priors on the concentration parameters: (4.4) Note that under these priors, E [ρ] = 0.8, which matches our large prior belief in selftransitions. The initial values are γ = 1, α = α 0 = 1, and κ = 4. In our experience, the sampler is insensitive to these initial values.

Initialization Scheme
Due to the high dimensionality of our model, we can greatly reduce the number of iterations required to reach convergence by choosing good initial values for the model parameters. To initialize the latent positions X (1) 1:T and the intercept β (1) 0 , we used the MAP estimate from a short MCMC chain (1,000 iterations) of the model with no clustering of Sewell and Chen (2015). Note that we modified their sampler to use the undirected likelihood with the linear predictor given in Equation (3.5). Due to the block sampling of the node labels Z 1:T and the fact that they are sampled before the other clustering parameters, the sampler is not affected by the initial values of the node labels. We initialized the group locations and shapes (µ 1:L , σ 2 1:L ) (1) by running a k-means algorithm for longitudinal data (Genolini and Falissard, 2010) on the initial latent positions with k = L and using the empirical mean and variance estimates of the groups as initial values. We set the initial blending coefficient λ (1) = 0.9 to encourage early clustering. We sampled the remaining parameters β (1) , (π 1 0 ) (1) , and Π (1) 2:T from their respective priors. We set the prior parameters on the intercept to µ β 0 = β (1) 0 and σ 2 β 0 = 2. Finally, we set the prior parameters on the blending coefficient to µ λ = 0.9 and σ 2 λ = 0.01.

Posterior Summarization
The selection of a single point estimate that adequately summarizes the full posterior is a challenging and open problem in Bayesian nonparametric statistics. A simple solution uses the posterior mode. However, in Bayesian nonparametric models, the MCMC chain often only visits a single partition once, which makes frequency estimates unreliable. Another approach is to choose the partition corresponding to the MAP (maximum a posteriori) estimate. A downside of the MAP approach is that it does not marginalize over the uncertainty in the partitions. As a result, MAP clustering estimates tend to over-fit the data. An ideal summarization methodology should take into account the clustering uncertainty suggested by the posterior samples.
In this work, we take a decision theoretic approach to posterior summarization. Specifically, we select the partition that minimizes the posterior expectation of an appropriately chosen loss function over dynamic clusterings. In the case of static clustering, two popular choices of loss functions are Binder's loss (Binder, 1978) and the variation of information (Meilȃ, 2007). Their use for posterior summarization was advocated by Lau and Green (2012) and Wade and Ghahramani (2018) respectively. Note that minimizing Binder's loss is equivalent to maximizing the Rand index, another popular measure of clustering performance. Both loss functions have the desirable property that they are metrics over the space of clusterings. Furthermore, the optimization of the expected losses only depends on the posterior co-occurrence probabilities, so avoids complications due to label-switching.
To extend this approach to dynamic clusterings, we minimize the posterior expected time-averaged variation of information (VI): arg min which averages the static VI's at each time-step. Similar to Wade and Ghahramani (2018), this is accomplished by minimizing the following lower-bound: arg min (4.6) Notice that this expression only depends on the posterior co-occurrence probabilities, which can be estimated from the posterior samples. We perform this optimization by searching over all posterior samples from the MCMC chain (after an appropriate burnin). Possible ties are broken by selecting the sample with the highest likelihood (Equation (3.3)). Note that this minimization does not require the partitions to be visited by the Markov chain; however, by restricting the search to the sampled partitions, we have access to the associated latent space for later visualization of the network.

Non-Identifiability of the Latent Positions
Since latent position models depend on the distance between actors, it is clear that they are invariant to translations, rotations, and reflections of the latent space. Any posterior inference that utilizes these positions must correct for such a non-identifiability. We take the approach of Procrustes matching commonly employed in the literature. This involves post-processing the samples by rotating and translating them to match a reference layout. For a detailed description of this procedure, see the original work by Hoff et al. (2002). We use the sample chosen by the procedure in Section 4.7 as our reference layout.

Simulation Study
We designed a simulation study to assess the HDP-LPCM's performance on synthetic networks with known community dynamics that mimic those found in real-world networks. We considered two scenarios: a sequence of networks with (1) a single static number of groups active over all time steps, and (2) a number of groups that change over time because the groups have either merged or split. We labeled these scenarios the time-homogeneous and time-inhomogeneous simulations, respectively. We expect both scenarios to occur in real-world dynamic networks, so it is important that the HDP-LPCM adequately handles both regimes.
We evaluated the performance of our method in four different ways. To measure the goodness-of-fit of the HDP-LPCM on the training data, we calculated the in-sample area under the receiver operating characteristic curve (AUC) based on the model's edge predictions. To evaluate the performance of the node level clustering, we used two popular metrics: the adjusted Rand index (ARI) and the variation of information (VI) described in Section 4.7. For the ARI, a value closer to one indicates a good clustering. For the VI, a value closer to zero indicates a good clustering. To measure the HDP-LPCM's ability to quantify uncertainty in the number of clusters, we recorded the posterior probabilities of |G t | at each time point. To assess the ability of our posterior summarization methodology (Section 4.7) to select the correct number of communities, we recorded the number of active communities selected by the model at each time step and compared to the ground truth value.
All models were estimated using a Markov chain consisting of 50,000 samples. We used the initial 5,000 samples to tune the step sizes of the Metropolis-Hastings samplers. We discarded the following 10,000 samples as burn-in leaving 35,000 samples for inference. Furthermore, we fixed the truncation level of the weak-limit approximation (Equation (3.15)) to L = 10 in all experiments. Finally, we repeated each simulation 50 times. For each repetition, we simulated a network starting with a different random seed.

Time-Homogeneous Simulation
This simulation contains a single set of groups and transition matrices common to all time points. We used a simulation scheme similar to the one found in Sewell and Chen (2017) but modified for undirected networks. The networks were generated according to the sampling mechanisms described in Sections 3.1 and 3.2. We fixed the total number of groups to G = 6 at each time point. We chose the blending coefficient λ = 0.8, the dimension of the latent space p = 2, and the intercept β 0 = 1.0. We set the group locations to (µ 1 , . . . , µ 6 ) = −1.5 1.5 −3 3 0 0 0 0 0 0 −2 2 . (5.1) We drew the group shapes σ g from a Γ -1 (6, 0.5 × 10 −1 ) distribution, and the initial distribution π 1 0 from a six-dimensional Dirichlet(10, . . . , 10). The six rows of the transition matrix, π t g for g = 1, . . . , 6 and t = 2, . . . , 6, were chosen to be proportional to We set the constant in the above equation equal to 20, which yields self-transition probabilities ranging from 0.83 to 0.89. The latent actor positions, X 1:T , and group assignments, Z 1:T , were drawn from Equation (3.9). Finally, the adjacency matrices were generated according to Equation (3.2) with the linear predictor η ijt given by Equation (3.5). We simulated 50 dynamic networks with T = 6 time points and n = 120 actors each. The results of these simulations are displayed in the first row of Table 1. The in-sample AUC is 0.842, which indicates that the model fits the observed data well. Furthermore,  the VI is 0.0324 and the adjusted Rand index is 0.990. Since the values are close to zero and one respectively, we can conclude that the clustering performance is very good. Next, we evaluate the model's ability to quantify the uncertainty in the number of clusters. Boxplots of the posterior probabilities of |G t | over the 50 simulations are displayed in Figure 3. In most cases, |G t | = 6 is correctly identified with the highest posterior probability. Finally, we assess the ability of our posterior summarization method to select a partition with the correct number of communities. The number of groups selected by this method at each time step is displayed in Figure 4. Our methodology is able to correctly select the six group partitioning in all fifty simulations.

Time-Inhomogeneous Simulation
We designed this simulation to test our model's ability to detect changes in group structure by allowing the number of groups to vary over time. We generated 50 networks with T = 9 time points and n = 120 actors each. There are G = 6 groups in the networks overall. The simulations begin with two groups, these two groups split into six groups at t = 4, and then the six groups merge into four groups at t = 7. Note that the adjacency matrices generated from this procedure are displayed in the introductory example ( Figure  1).
Once again, the networks were simulated according to the sampling mechanisms described in Sections 3.1 and 3.2. This simulation differs from the time-homogeneous simulation in how we specified the transition matrices π t g . We drew the initial distribution π 1 0 from a two-dimensional Dirichlet(10, 10). We chose the rows of the transition matrix, π t g for g ∈ G t−1 and t = 2, . . . , 9, proportional to 1 {1∈Gt} where G t indicates the set of active groups at time t. At t = 1, 2, 3 there are two groups so that G 1 = G 2 = G 3 = {1, 2}. From t = 4, 5, 6 these two groups split into six groups so that G 4 = G 5 = G 6 = {1, 2, 3, 4, 5, 6}. Finally, these six groups merge into four groups at t = 7, 8, 9 so that G 7 = G 8 = G 9 = {3, 4, 5, 6}. For all time points, we set the constant in Equation (5.3) to 20 except for t = 4. At t = 4, we set the constant equal to 1 so that the nodes were more evenly distributed among the six newly created groups. The results of the simulation are displayed in the second row of Table 1. The median in-sample AUC is 0.850, which indicates that the model fits the observed data well. Furthermore, the median VI is 0.0292 and the median ARI is 0.990. Since the values are close to zero and one, respectively, we conclude that the clustering performance is very good. We note that the standard deviations of these metrics are much higher in this simulation. To make sense of this discrepancy, we display boxplots of the metrics in Figure 5. We see that the VI and the ARI have a few large outliers for some simulations; however, the majority of simulations are clustered near the mean values, which indicates a good clustering.
Next, we evaluate the ability of the HDP-LPCM to quantify the changing macrolevel community structure. The posterior probabilities for |G t | are displayed in Figure  6. Recall that the ground truth is |G t | = 2 for t = 1 − 3, |G t | = 6 for t = 4 − 6, and |G t | = 4 for t = 7 − 9. For most time points the correct number of groups has the highest posterior probability. The exception is t = 7, which is when the six groups merge into four groups. As is common in Dirichlet process mixture models, the model tends to create a few small clusters, which can inflate the group count. This is the case here, where the extra clusters at t = 7 tend to be composed of less than five nodes. Once again we assess our posterior summarization method's ability to select a clustering with the correct number of communities. The number of groups selected by this method at each time step is displayed in Figure 7. For the majority of simulations, our method selects the correct number of clusters. The method struggles the most at t = 7, typically because of a few small clusters. Overall, the clustering is correct and we feel confident the HDP-LPCM is able to adequately infer evolving communities. Figure 6: Boxplots of the posterior probabilities of |G t | over the 50 time-inhomogeneous simulations. The data was generated with |G t | = 2 for t = 1 − 3, |G t | = 6 for t = 4 − 6, and |G t | = 4 for t = 7 − 9. In most cases, the HDP-LPCM is able to assign high posterior probability to the correct number of communities. Figure 7: The number of groups predicted at each time step by the VI estimator for the time-inhomogeneous simulation. The data was generated with |G t | = 2 for t = 1 − 3, |G t | = 6 for t = 4 − 6, and |G t | = 4 for t = 7 − 9. Our posterior summarization method is able to select the correct number of communities in the majority of the simulations.

Real Data Application
In this section, we demonstrate the utility of our proposed HDP-LPCM on a variety of real-world dynamic networks with an evolving community structure. We include two applications: inferring changing international military alliances during the Cold War and detecting dynamic plotlines in the television series Game of Thrones. Also, in Section S.4 of the supplementary materials, we show that the HDP-LPCM corroborates many previous findings concerning Sampson's monastery network (Sampson, 1969), a standard pedagogical example in the network literature.

International Military Alliances
We begin by using the HDP-LPCM to examine changes in international military alliances during the first three decades of the Cold War (1950 -1979). We use the Formal Alliances (v4.1) dataset curated as part of the Correlates of War Project (Gibler, 2009b). The raw dataset records all formal alliances -mutual defense pacts, non-aggression treaties, and ententes -among nations between 1816 and 2012. The goal of our analysis is to uncover the competing political blocs that defined the Cold War period in history and to determine any points in time where that alliance structure changed.
For this analysis, we use the yearly dyadic dataset, which records an undirected edge between two nations if there is a formal alliance between them during that year. To simplify the analysis, we discretized the dynamic networks into five year chunks from 1950 -1979. We removed the nodes with a degree less than two from each network to focus on the larger political blocs found in the dataset. In addition, we required a nation to have at least one alliance during 1950 -1979. We binarized the relations so that a connection between nations i and j at time t means they had at least one active alliance during those five years. This preprocessing resulted in T = 6 undirected binary networks that each contains n = 107 actors.
We fit our proposed model with a truncation level L = 25 to the international military alliances networks using 50,000 iterations for tuning, 50,000 iterations for burn-in, which left a remaining 400,000 iterations for inference. The trace plots, autocorrelation functions (ACFs), and marginal densities are displayed in Figure S.5 in Section S.5 of the supplementary materials. Visual inspection of the trace plots indicates the algorithm has converged.
The alluvial diagram (Figure 8) of the partition selected with the procedure described in Section 4.7 shows that the HDP-LPCM estimates six communities overall with five communities active during all six time points. Note that the model uses group 2 to collect isolated nodes and small intermittent alliances. This interpretation is supported by the fact thatσ 2 = 132, which is 28 times larger than the second largest group shapeσ 1 = 4.67, and it encompasses most of the latent space. Thus we consider groups 1, 3, 4, 5, 6, and 7 as the major blocs inferred by the model. Note that these blocs strongly dictate the network's dynamics with an inferred blending coefficient λ = 0.994. The estimated latent spaces for the years 1950 -1954 and 1960 -1964 are depicted in Figure 9 and Figure 10. The latent spaces for the remaining years are included in the supplementary materials. The static communities (groups 1, 3, 4, 5, and 6) coincide with long term regional alliances during the Cold War. Group 1 corresponds to the Western Bloc, consisting primarily of nations that are a part of the North Atlantic Treaty Organization (NATO) and the Western European Union (WEU). The competing Eastern Bloc is represented by group 5, which consisting of the Soviet Union and its satellite states such as East Germany, Czechoslovakia, and Poland. Group 3 consists of the Latin American coalition of the Organization of American States (OAS), founded in 1948 to oppose socialism. Group 4 consists of member nations of the Southeast Asia Treaty Organization (SEATO), which is the Asian equivalent of NATO. Group 6 is the Arab League (at the time the League of Arab States) formed in 1945 to protect the interests of Arab countries. As depicted in the alluvial diagram, an interesting finding of the HDP-LPCM is that there is little exchange of nations between these groups over this time period. This finding bolsters the claim that "once the 'cold war' confrontation became apparent ... many nations cast their lot with either the American or Soviet bloc" (Small and Singer, 1969).
The evolving community structure (birth of group 7) is a result of the emergence of Africa as a world power. In particular, seventeen African nations gained their independence in 1960 alone. As a result of these newly independent nations, a large number of regional alliances formed in the early 1960s (Gibler, 2009a). This is reflected in the HDP-LPCM by the emergence of group 7 at t = 3, which encompasses the newly independent African nations. Fewer African alliances formed over the next decade. Our model reflects this fact by only including the Union of Central African States, which consists of former French colonies, in group 7 at t = 4 and 5. Finally, the Economic Community of West Figure 9: Latent space of international military alliances for the years 1950 -1954. The group means µ g are denoted by a +, and the group shapes σ g are displayed as twostandard deviation ellipses. The names of each nation and node number are annotated. The undirected edges are also displayed. For clarity, all unconnected nodes and group 2's two-standard deviation ellipses are removed. Figure 10: Latent space of international military alliances for the years 1960 -1964. The group means µ g are denoted by a +, and the group shapes σ g are displayed as twostandard deviation ellipses. The names of each nation and node number are annotated. The undirected edges are also displayed. For clarity, all unconnected nodes and group 2's two-standard deviation ellipses are removed.
African States (ECOWAS) formed in 1975, which resulted in the re-introduction of many western African states into group 7 at t = 6. The dynamic nature of the HDP-LPCM is essential in revealing the emergence of Africa as a world power and demonstrates the importance of incorporating community evolution in latent space network modeling.

Character Interactions in Game of Thrones
In this section, we study the networks of character interactions in the television series Game of Thrones. The goal of our analysis is to use community detection to pinpoint coherent dynamic plotlines within the series. This is an interesting case study because the series' narrative is known for its many dynamic stories and characters who freely move between them. As such, we expect various groups of characters to form and die out across the series' lifetime and actors to freely move between these groups.
We utilize the networks compiled by Beveridge and Chemers (2018), who parsed fangenerated scripts found on the user curation site Genius. The original dataset consists of weighted character-character interactions split up between the eight seasons of the television series. The weight of an edge equals the number of interactions during a given season. Since this is primarily a pedagogical example, we restrict the dataset to the first four seasons of the show. To remove minor characters, we only keep interactions that occur greater than or equal to 10 times each season. The final result is T = 4 binary undirected networks with a total of n = 165 actors each.
We fit the HDP-LPCM to this dataset with a truncation level of L = 25 using 50,000 iterations for tuning, 50,000 iterations for burn-in, which left a remaining 400,000 iterations for inference. The trace plots of the unnormalized log-posterior, the intercept β 0 , and the blending coefficient λ are displayed in Figure S.10 in Section S.5 of the supplementary materials. Visual inspection of the trace plots indicates the model has converged.
The alluvial diagram ( Figure 11) and associated latent space ( Figure 12) reveal a dynamic group structure. The model infers six groups overall; however, only three groups are active during the first season of the show. During the second season, group 2 splits off into groups 3, 4, 6, and 7. In addition, group 5 is created out of characters from group 6. Note that the model uses group 6 to collect inactive characters, so we exclude it from further analysis. After the change in season two, the network's group structure remains constant for seasons three and four. The latent space of seasons three and four are included in the supplementary materials. Furthermore, since the blending coefficient λ = 0.974, we conclude that the groups drive the evolution of the network. Overall, the ability of the HDP-LPCM to infer an evolving group structure is crucial for properly understanding the network's dynamics.
The inferred groups and their dynamics are consistent with the storylines in the Game of Thrones series. During season one, there are two active groups. Group 2 consists of all characters on Westeros, while group 1 centers around Deanery Targaryen's story on Essos. Starting at season two, group 2 contains characters revolving around Arya Stark's story arc, House Lannister, and other characters at King's Landing. Group 3 contains Bran Stark's group of companions (Hodor, Rickon, Osha, Meera, and Jojen) as they travel north of the wall. The plotlines related to the remainder of House Stark are contained in group 7. This group is composed of both Robb Stark's contingent as well as those related to Theon Greyjoy, who currently rules Winterfell the ancestral home of the Starks. Group 4 pertains to the story north of the wall. It most notably contains Jon Snow and Sam Tarly. Lastly, group 5 revolves around the story of Stannis Baratheon and his quest to regain the Iron throne, which was introduced in season two. In conclusion, our model recovers narratively coherent groups within the Game of Thrones character interaction networks and identifies their temporal dynamics within the storyline.

Discussion
In this article, we proposed the hierarchical Dirichlet process latent position clustering model (HDP-LPCM) for dynamic networks. This is the first, to our knowledge, latent position model that can detect evolving community structures. In addition, we demonstrated that that HDP-LPCM still performs well when there is a static group structure. To accomplish our modeling goals, we used Bayesian nonparametric priors to provide simultaneous inference over the number of communities at each time point in the network and the dynamics of each actor's latent position. Furthermore, our MCMC inference procedure has computational advantages over existing approaches. In particular, we avoided the BIC approximations of Handcock et al. (2007) and the computationally expensive process of estimating a large number of models.
Season 1 Season 2 Figure 12: Latent space of the Game of Thrones character interaction network for season 1 and season 2. The group means µ g are denoted by a +, and the group shapes σ g are displayed as two-standard deviation ellipses. The names of each character and node number are annotated. The undirected edges are also displayed. For clarity, all unconnected nodes and group 6's two-standard deviation ellipses are removed.
In this work, we focused solely on community detection for binary undirected networks; however, our methodology can easily be extended to other network types. For example, both directed and weighted networks can be accommodated through the minor modifications of the likelihood presented in . Due to the conditional dependence structure of the HDP-LPCM, one only needs to modify the Metropolis-Hastings steps (i.e., steps 2 and 3) in Algorithm 1 to use the new network likelihood.
Common to most latent distance models, our proposed MCMC estimation method is time-intensive for large networks due to the quadratic scaling of the procedure. To alleviate this issue, one could utilize the case-control likelihood of Raftery et al. (2012). This method approximates the full likelihood by sub-sampling the unconnected edges, which allows for inference that scales linearly with the number of nodes in the network.
To further decrease the model's runtime, one could create a variational Bayes or stochastic variational Bayes (Hoffman et al., 2013) algorithm. However, the lack of conjugacy in the distance model means such an algorithm is non-trivial to implement. For this reason, we leave these extensions to future work.
Although our model can infer the number of communities from the data, it does not provide inference for the dimension of the latent space p. In the numerical studies, we fixed p = 2, which allowed for intuitive visualizations as well as being adequate for the networks we considered in this work. Furthermore, we believe such visualizations allow for human-in-the-loop model scrutiny, which is a primary strength of the latent position model over other modeling approaches such as latent feature models (Miller et al., 2009). Of course, selecting an appropriate dimension of the latent space is crucial if one is focused on tasks such as link prediction. A straightforward approach is to use a model selection criteria such as BIC or the Watanabe-Akaike information criterion (WAIC) (Watanabe, 2010). Another possibility is to put a prior over the dimensionality such as in Durante and Dunson (2014). Regardless, the selection of p in latent space models remains an interesting and open question for future research.
Despite these possible extensions, we believe the HDP-LPCM has a broad range of applicability. As demonstrated on the real-world networks in this paper, applications include the detection of changing alliances structures as well as inference regarding narrative plotlines. We believe our methodology can be built upon to provide new tools for understanding the sequential evolution of communities in dynamic networks.

Supplementary Materials for "A Bayesian Nonparametric Latent Space Approach to Modeling Evolving Communities in Dynamic Networks"
Joshua Daniel Loyal and Yuguo Chen

S.1 Conditional Distributions Used in the Gibbs Sampler
In this section, we present the conditional distributions used in the Metropolis-Hastings within Gibbs sampler for the HDP-LPCM. In what follows, let I kt = {i : Z i t = k} be the set of nodes in group k at time t. Furthermore, let n kt = |I kt | be the number of actors in group k at time t, and let n k· = T t=1 n kt . We denote conditioning on everything but the parameters of interest as | ·. The full conditionals are Note that we set the following priors on the HDP parameters: Given a previous set of hyperparameters ( , α ( −1) , and κ ( −1) , update the current parameters as follows: 1. Initialize current parameters to the values of the ( − 1)th iteration.
7. Update ρ: where we exclude the considered dishes at s = 1 because the initial restaurant does not include a stickiness parameter.

S.4 Sampson's Monastery
In this section, we apply the HDP-LPCM to dynamic networks of social relations between 18 monks in an isolated New England monastery (Sampson, 1969 friend during time t. The result is T = 3 undirected networks with n = 18 actors. Since this dataset is commonly analyzed using latent position and latent cluster models, we will compare our results to findings of previous analyses to critique our model's performance.
In addition, we extend past analyses by providing a measure of stability for the network's group structure. In particular, we aim to answer the question of whether the three factions were present over all observational periods or split-off during a particular time point.
We fit the HDP-LPCM using 15,000 iterations for tuning, 20,000 iterations for burnin, which left a final chain of length 165,000 iterations for inference. The trace plots of the unnormalized log-posterior, the intercept β 0 and the blending coefficient λ are displayed in Figure S.1. We also provide ACF plots and a kernel density estimate of the marginal posterior distributions. Visual inspection of the trace plots indicates that the algorithm has converged.
Our analysis corroborates many of the previous findings about the group structure of Sampson's monastery. The posterior mean of the blending coefficient λ is 0.869, which confirms Sampson's assertion that a strong group effect drives the dynamics of the monk's social network. Furthermore, previous static and dynamic analysis widely agree that a 3 or 4 component model adequately describes the network (Handcock et al., 2007;Krivitsky et al., 2009;Ryan et al., 2017). The posterior distribution over the number of groups at each time step G t is displayed in Figures S.2(a) -S.2(c). These distributions confirm that 3 groups is the most probable number of groups at each point in time.
A distinctive advantage of our approach compared to other latent space models is the ability to quantify the stability of the factions within the monastery. To do this, we follow the decision theoretic procedure described in Section 4.7 to select a representative  three group structure during all time points. Furthermore, we observe that the monk's group memberships do not change over time. As further confirmation of our method, we note that the HDP-LPCM assigns all monks to the same communities as Sampson's ground truth labels. Overall, our method is not only consistent with previous latent space analysis of Sampson's monastery network, but also expands the literature through its indication that the group structure is static over the observed time period. . The group means µ g are denoted by a +, and the group shapes σ g are displayed as two-standard deviation ellipses. The names of each monk and node number are annotated. The undirected edges are also displayed.

S.5 Additional Figures
In this section, we include the remaining figures for the two datasets analyzed in the real data applications of the main text (Section 6). This includes the trace plots for the MCMC algorithm as well as additional latent space visualizations. The trace plots of the unnormalized log-posterior, the intercept β 0 and the blending coefficient λ for the HDP-LPCM fit to the international military alliances network are displayed in Figure S.5. The latent spaces for 1955-1959, 1965-1969, 1970-1974, and 1975-1979 are displayed in Figures S.6, S.7, S.8, S.9 respectively.  The group means µ g are denoted by a +, and the group shapes σ g are displayed as twostandard deviation ellipses. The names of each nation and node number are annotated. The undirected edges are also displayed. For clarity, all unconnected nodes and group 2's two-standard deviation ellipses are removed. The group means µ g are denoted by a +, and the group shapes σ g are displayed as twostandard deviation ellipses. The names of each nation and node number are annotated. The undirected edges are also displayed. For clarity, all unconnected nodes and group 2's two-standard deviation ellipses are removed.   The group means µ g are denoted by a +, and the group shapes σ g are displayed as two-standard deviation ellipses. The names of each character and node number are annotated. The undirected edges are also displayed. For clarity, all unconnected nodes and group 6's two-standard deviation ellipses are removed. Figure S.12: Latent space of the Game of Thrones character interaction network for season 4. The group means µ g are denoted by a +, and the group shapes σ g are displayed as two-standard deviation ellipses. The names of each character and node number are annotated. The undirected edges are also displayed. For clarity, all unconnected nodes and group 6's two-standard deviation ellipses are removed.