A Class of Temporal Hierarchical Exponential Random Graph Models for Longitudinal Network Data

TL;DR

This paper introduces a new class of statistical models called THERGM for analyzing longitudinal networks with community structure, capturing evolving clusters and network features over time.

Contribution

It proposes a novel two-stage estimation strategy combining latent space clustering with temporal ERG models for dynamic network analysis.

Findings

Effective in identifying evolving communities

Good fit to simulated data

Accurate link prediction results

Abstract

As a representation of relational data over time series, longitudinal networks provide opportunities to study link formation processes. However, networks at scale often exhibits community structure (i.e. clustering), which may confound local structural effects if it is not considered appropriately in statistical analysis. To infer the (possibly) evolving clusters and other network structures (e.g. degree distribution and/or transitivity) within each community, simultaneously, we propose a class of statistical models named Temporal Hierarchical Exponential Random Graph Models (THERGM). Our generative model imposes a Markovian transition matrix for nodes to change their membership, and assumes they join new community in a preferential attachment way. For those remaining in the same cluster, they follow a specific temporal ERG model (TERGM). While a direct MCMC based Bayesian estimation is computational infeasible, we propose a two-stage strategy. At the first stage, a specific dynamic latent space model will be used as the working model for clustering. At the second stage, estimated memberships are taken as given to fit a TERG model in each cluster. We evaluate our methods on simulated data in terms of the mis-clustering rate, as well as the goodness of fit and link prediction accuracy.

Equations

$${\mathbb{P}}_{\theta }\{Y=y\}=\text{exp}\left({\theta }^{\prime }S(y)-\psi (\theta )\right)$$

$$\text{logit}\left(P_{\theta }\{Y_{i,j}=1|Y_{i,j}^c\}\right)={\theta }^{{}^{\prime }}{c}_{i,j}$$

$$\mathbb{P}(Y=y|M=m)=\prod _{k=1}^K\mathbb{P}(Y_k=y_k|M=m)\prod _{l=1}^{k-1}\mathbb{P}(Y_{kl}=y_{kl}|M=m)$$

$${\mathbb{P}}_{{\theta }_k}(Y_k=y_k|M=m)=exp\{{\theta }_k^{\prime }{S}_k(y_k)-{\psi }_k({\theta }_k)\}$$

$$\mathbb{P}(Y_{kl}=y_{kl}|M=m)=\prod _{i\in A_k,j\in A_l}\mathbb{P}(Y_{ij}=y_{ij}|M_i=m_i,M_j=m_j)$$

$$\mathbb{P}(Y^1,Y^2,\dots ,Y^t|Y^0)=\mathbb{P}(Y^t|Y^{t-1})\mathbb{P}(Y^{t-1}|Y^{t-2})\cdots \mathbb{P}(Y^1|Y^0).$$

$$\mathbb{P}(Y^t|Y^{t-1},\mathbf{\theta })=\exp \left\{{\mathbf{\theta }}^{\prime }S(Y^t,Y^{t-1})-\psi (\mathbf{\theta },Y^{t-1})\right\}$$

$$Y_{i,j}^{(t)}=\{\begin{array}{ll}1& \text{if there is an edge between vertices }i\text{ and }j\\ 0& \text{otherwise}\end{array}$$

$$\mathbb{P}(\{Y{\}}_1^t,\{M{\}}_1^t|Y^{(0)},M^{(0)})=$$

$$\mathbb{P}(Y^{(1)}|Y^{(0)},M^{(1)},M^{(0)})\mathbb{P}(\{M{\}}_1^t|Y^{(0)})$$

$$\mathbb{P}(Y^{(t)}|Y^{(t-1)},M^{(t)},M^{(t-1)})=$$

$$\mathbb{P}(Y_{k,new}^{(t,t-1)}|M^{(t)},M^{(t-1)})$$

$$\prod _{l=1}^{k-1}\mathbb{P}(Y_{kl}^{(t)}|M^{(t)})$$

$$\mathbb{P}(Y_k^{(t,t-1)}|Y_k^{(t-1,t)},{\theta }_k)=\exp \left\{{{\mathbf{\theta }}_k}^{\prime }S(Y_k^{(t,t-1)},Y_k^{(t-1,t)})-\psi ({\mathbf{\theta }}_k,Y^{(t-1,t)})\right\}$$

$$\mathbb{P}(\{M{\}}_0^t)=\mathbb{P}(M^{(t)}|M^{(t-1)})\cdots \mathbb{P}(M^{(1)}|M^{(0)})\mathbb{P}(M^{(0)}).$$

$$P(M^t|M^{t-1})=\prod _i\left(\eta {\delta }_{m_i^t,m_i^{t-1}}+(1-\eta )q_{m_i^t}\right),$$

$$P(Y=y|Z,X,\beta )=\prod _{i\ne j}P(Y_{i,j}=y_{i,j}|x_{i,j},z_i,z_j,\beta )$$

$$\text{logit}\{P(Y_{i,j}=1|x_{i,j},z_i,z_j,\beta )\}={\beta }_0^T{x}_{i,j}-{\beta }_1|z_i-z_j|$$

$$z_i\stackrel{\text{iid}}{\sim }\sum _{k=1}^K{\lambda }_k{\text{MVN}}_d({\mu }_k,{\sigma }_k^2{\text{I}}_d)$$

$$P(\{Y^t,Z_t,M_t{\}}_1^T)=\prod _{t=1}^T\mathbb{P}(Y^t|Z^t)\times \mathbb{P}(\{Z_t,M_t{\}}_1^T)$$

$$M_i^t|M_i^1,\dots ,M_i^{t-1}=M_i^t|M_i^{t-1}$$

$$Z_i^t|Z_i^1,\dots ,Z_i^{t-1},M_i^1,\dots ,M_i^t=Z_i^t|Z_i^{t-1},M_i^t$$

Full text

22footnotetext:

1 Introduction

Relational data over a series of time points are becoming more available as the information technology advances, for example, [salathe2010high] study flu transmission by using sensors to record human contact in a high school. One of the most attractive topics in network research is community detection or clustering (see [fortunato2010community] for a review). Speaking of scientific literatures, communities could be corresponding to research fields. Those fields may or may not be pre-specified, take the “network research community” as example, it has seen mixed contributions and gaining momentums from scholars in Physics, Computer Science, Sociology, Statistics and other disciplines [newman2011structure, wasserman1994social, jackson2008social]. Multiple snapshots, or even at the finest grain as a stream of relational events, provide more information than a static network. It provides us opportunities to capture the trajectories of evolving network structures, with applications to trend analysis [leskovec2007patterns] and link prediction [sarkar2014].

Besides the community structure, there could be various types of network effects going on locally. For example, the common sense that “my friends’ friends are more likely to be my friends” and “people with same interests are more likely to be friends” are typically referred to as transitivity and homophily, respectively. With a long history of being studied by Sociologists [granovetter1973strength, mcpherson2001birds], such kind of network structures can be flexibly specified and therefore estimated/tested by a powerful statistical tool called Exponential Random Graph Models (ERGM) [frank1986markov, wasserman1996logit]. We should have reasonable suspect that transitivity or homophily acts equally in different communities and it would be natural to construct a statistical model that consider local structural effects as well as clustering, which represents a more macro level structure. However, current works for large-scale network structure analysis have mostly been limited to base on the Stochastic Block Models (SBM) [holland1983stochastic, snijders1997block], which assumes no dependence among ties within the same block (cluster). SBM is adopted for its simplicity in the hope that it approximates well from the perspective of an efficient clustering estimator. Nice statistical properties have been established for computationally effective algorithmic methods such as spectral clustering, under both the static [rohe2011spectral, qin2013regularized, jin2015score] and dynamic [han2014consistent] settings. While efforts are being made to extend this often over-simplified family of models to a boarder class [wan2015class], ERGM is generally considered having degeneracy problem [handcock2003assessing] and too expensive to compute the MCMC approximated MLE [snijders2002markov]. In fact, significant improvement on model specification and estimation technique has been made [snijders2006new, hunter2006curved, hummel2012improving] for the past decade. More importantly, the longitudinal extensions of ERGM rarely suffer from the model degeneracy problem [hanneke2010temporal]. It can be intuitively understood as that structural effects became much more clear when the changes are explicitly modeled, rather than being very likely blurred in a cumulative snapshot.

In this paper, we propose a class of Temporal Hierarchical Exponential Random Graph Models (THERGM) to fully utilize the information contained in longitudinal network data, which may greatly improve the overall goodness of fit and link prediction accuracy. It seems that we pay the price of a non-scalable inference procedure in order to gain the modeling power of ERGM, but fortunately, [leskovec2009size] suggests that roughly 100 nodes is probably a universal "natural size" of a meaningful cluster, including social, information, Internet and biological networks. While this is the right scale that current ERGM implementation [handcock2008statnet] could be routinely carried out, it is still computationally infeasible to directly combine with the community structure in a straightforward Bayesian approach. To tackle this problem, we propose a two-stage strategy that uses a Dynamic Latent Space Model [hoff2002latent, sewell2014dynamic, sewell2016community] to do the clustering at the first stage. That working model should be corresponded to the particular ERG Model you choose as specifically as possible. For example, a distance measure or some covariates should be included if the (hypothesized) true model contains transitivity or homophily effect, respectively. A Latent Space Model that accounts for all dependencies among relational ties plus an appropriate model of community evolution could be expected to recover the clustering change trajectory asymptotically. Then at the second stage, we take the clustering estimates at each time point as given, and further fit a temporal ERG model for actors remaining in each community. In this way, we harness the describing power of more complex models (THERGM versus dynamic SBM) in a timely manner by having the sizes after decomposition under control. For researchers who only care about the inference of (temporal) ERGM, our approach could provide correctness, rather than feasibility, as it could be easily confounded by the hierarchical structure if there are unknown clusters.

The rest of paper is organized as following: In Section 2, we briefly state the form of Exponential Random Graph Models. In Section 3 we propose our new Hierarchical ERGM for temporal clustered network data. In Section 4, we describe the two-phase strategy for inference. In Section 5, various simulation settings for our generative models are used to evaluate the performance of the proposed two-phase strategy. We summarize contributions and discuss limitations as well as future work in Section 6.

2 Related Works

2.1 Exponential Random Graph Models

Exponential Random Graph Models is a family of statistical models taking the following form of the probability functions

[TABLE]

with arbitrary statistics $S(y)$ for researchers to specify structures of scientific interests. The interpretation of parameters $\theta$ is typically based on the log odds ratio of forming a tie, conditional on the rest of the graph since:

[TABLE]

where $Y^{c}_{i,j}=\{Y_{u,v}|\text{ for all }u<v,(u,v)\neq(i,j)\}$ represents all other ties except $Y_{i,j}$, $c_{i,j}=S\left(y^{(ij1)}\right)-S\left(y^{(ij0)}\right)$ is the change statistic with $y^{(ij0)}$ and $y^{(ij1)}$ denoting the adjacency matrices with the $(i,j)$th element equal to $1$ and [math] while all others are the same as $y$.

2.2 Hierarchical ERGM

[schweinberger2015local] proposed a hierarchical extension of ERGM by introducing the local dependence that breaks down the dependence of random graph $Y$ into subgraphs. Assume there is a partition of the vertices V into $K\geq 2$ non-empty finite subsets $V_{1},\dots,V_{K}$, such that the within- and between- neighbourhood subgraphs $Y_{k,l}$ given membership $M$ satisfy

[TABLE]

The within-neighborhood probability measures $\mathbb{P}_{k,k}$ take specific ERGM forms as

[TABLE]

whereas the between-neighborhood probability measures $\mathbb{P}_{k,l}$ induce independence between subgraphs, and the between-neighborhood ties are assumed to be independent

[TABLE]

2.3 Temporal ERGM

[hanneke2010temporal] proposed a temporal extension of ERGM by making a Markov assumption on the network from one time step to the next. Specifically, $Y^{t}$ is independent of $Y^{1},\ldots,Y^{t-2}$ given $Y^{t-1}$. Taking the first observation $Y^{0}$ as given, the joint distribution can be factorized:

[TABLE]

Given the Markov assumption, ERGM is generalized for evolving networks assuming $Y^{t}|Y^{t-1}$ admits an ERGM representation. The conditional PDF takes the following form:

[TABLE]

The normalizing constant $\psi$ now also depends on $Y^{t-1}$, it is still intractable and the same MCMC approximated MLE techniques apply.

[krivitsky2014jrssb] further extended the discrete Temporal ERGM to a separate parameterization of incidence from duration, motived by the reality that social processes and factors of ties forming are very likely not the same as those of ties dissolving.

2.4 Other Temporal Models

Another popular class of models for dynamic networks is Stochastic Actor-Oriented Models (SAOM) [snijders2001statistical, snijders2010introduction]. [leifeld2015theoretical] explains how SAOM has an ERGM as its limiting distribution under a specific but rarely applied rate function. TERGM outperforms SAOM substentially when the DSP of SAOM was not met so TERGM is more robust in the sense of model mis-specification. [kim2013nonparametric] models the birth and death of individual groups via a distance dependent Indian Buffet Process and capture the evolution of node group memberships via a Factorial Hidden Markov model.

3 Temporal Hierarchical ERGM

We begin with notations. Consider a temporal series of networks $\{{G}^{(1)},\ldots,{G}^{(T)}\}$, where ${G}^{(t)}\equiv\{V^{(t)},E^{(t)}\}$ represents the network observed at time $t$, consisting a set of vertices $V^{(t)}=\{1,\dots,n^{(t)}\}$ and edges $E^{(t)}=\{(i,j)|i,j\in V^{(t)}\}$. We assume a fixed number of $K$ communities(or clusters, neighborhoods, blocks), and each vertex $i$ has a membership(or color) $m_{i}^{(t)}=k$ for $k\in{1,\dots,K}$, which may vary at a different time $t^{\prime}$. The distribution of communities is also allowed to change, as it reflects the overall trend of the system, for example, the rising of “big data” in Statistics. In this paper, $G^{(t)}$ is restricted to a binary adjacency matrix $Y^{(t)}=\{Y_{i,j}^{(t)}\}_{1\leq i\neq j\leq n}$ where

[TABLE]

Most of the literature assume that $V=V^{(t)}$ for all $t$ is a fixed set over time, while any dynamically changing network in real world would have nodes entering or leaving the system. It is because of that network structural analysis are concerning about the effects of existing ties on the formation (or dissolution) of future relations, which are implicitly among the same set of nodes. If the nodes change too volatile, it is the reason(s) of entering (leaving) should be modeled. For example, during the first 3 months from the launch of facebook.com, what really matters is how to attract people to open accounts, rather than introducing friends’ friends. Here we further assume this situation applies to the community level. That is to say, for subsets $V^{(t-1)}_{k}$ and $V^{(t)}_{k}$ of vertices belonging to community $k$ at any two consecutive time points, the remaining part $V^{(t-1)}_{k}\cap V^{(t)}_{k}$ is at least a significant portion of the union $V^{(t-1)}_{k}\cup V^{(t)}_{k}$. Together with a finite and pre-fixed number of clusters $K$, we are targeting a gradually evolving system. For streaming graph data, see for example [aggarwal2010stream].

3.1 Our model

Now we propose a class of Temporal Hierarchical Exponential Random Graph Models (THERGM), by assuming 1. $Y^{t}$ can be decomposed by the current clustering $M^{t}$ at any time $t$; 2. Markov property of the temporal dependency on $Y^{t-1}$ and also implicitly $M^{t-1}$:

[TABLE]

Denote $V^{t}_{k,remain}\equiv V^{t-1}_{k}\cap V^{t}_{k}$ as the set of remaining nodes in group $k$ during the time $t-1$ and $t$, so that ${G}^{t-1,t}_{k}\equiv\{V^{t}_{k,remain},E^{t-1,t}_{k}\}$ and ${G}^{t,t-1}_{k}\equiv\{V^{t}_{k,remain},E^{t,t-1}_{k}\}$ are the subgraphs spanned by the common nodes at time $t-1$ and $t$, respectively. $Y^{t,t-1}_{k}$ and $Y^{t-1,t}_{k}$ are the corresponding adjacency matrices:

[TABLE]

Note the first item in the RHS of equation 3.1 accounts for the ties associated with remaining nodes in group $k$, which takes the temporal ERGM form:

[TABLE]

Note that it is implicitly depending on $M^{(t)}$ and $M^{(t-1)}$. The second item accounts for the ties associated with nodes newly joined group $k$ from time $t-1$ to time $t$, which can be assumed to follow a preferential attachment process [barabasi1999emergence]. The third item accounts for the between-neighborhood ties, which are assumed to be independent as in equation 5.

To complete the THERGM specification, a sensible layer of community evolution is needed. Similar to longitudinal network observations, Markovian assumptions are made to simplify the joint distribution

[TABLE]

and there are at least three ways to do it:

• •

[yang2011detecting, sewell2016community] use a stationary transition matrix $\{\beta_{hk}\}$, to denote the probability of a node in cluster $h$ at time $t-1$ to become a member of cluster $k$ at next step $t$.

• •

[xing2010aoas, ho2011evolving] impose a logistic normal prior with a State Space model to model the trajectory.

• •

[ghasemian2015detectability] propose that with probability $\eta$, each node keeps its label from one time step to the next, otherwise it chooses a new label $m_{i}^{t}$ from the prior $q_{r}$. The transition probability for community memberships is

[TABLE]

We adopt the first way of modeling, and implicitly assume that the community evolution does not depend on the realized networks. See latent feature propagation model of [heaukulani2013dynamic], to capture the phenomenon that observed social interactions from the past can also affect future unobserved structure. The proposed dynamic process can be illustrated in Figure 1:

4 Model inference

A straightforward full Bayesian approach to estimate THERGM parameters is not practical. At the clustering level, the sample space of membership $M$ is $K^{n}$ where $K$ is the number of blocks and $n$ is the number of nodes. At the within-block TERGM part, the sample size of edge variable $Y_{k}$ is $2^{n_{k}\choose 2}$ where $n_{k}$ is the number of nodes in $k$th block (each pair of nodes can have a link present or absent in an undirected binary network). So both parts need MCMC or other sampling methods to do the approximation, directly combining them makes the problem intractable. In this section, we propose a two-stage strategy to tackle this problem.

4.1 Clustering with a working model

Now, we temporarily jump out of the ERGM paradigm and take a look at another class of statistical models that mainly embracing the classical latent factor idea. Instead of explicitly modeling dependence, the Latent Structure Models (LSM) postulate latent nodal variables $Z$ and conditional independence of $Y_{i,j}$ given those variables $Z_{i}$ and $Z_{j}$.

[hoff2002latent] introduced the concept of unobserved "social space" within which each node has a position so that a tie is independent of all others given the unobserved positions of the pair it connects to:

[TABLE]

where $X$ are observed covariates, and $\beta$ and $Z$ are parameters and positions to be estimated. [handcock2007cluster] took a subclass Distance Models where the probability of a tie is modeled as a function of some measure of distance between the latent space positions of two nodes:

[TABLE]

with restriction of $\sqrt{\frac{1}{n}\sum_{i}{|z_{i}|}^{2}}=1$ for the identification purpose. Then they imposed a finite mixture of multivariate Gaussian distribution for $z_{i}$ to represent clustering:

[TABLE]

where non-negative $\lambda_{k}$ is the probability that an individual belongs to the $k$th group, with $\sum\limits_{k=1}^{K}\lambda_{k}=1$.

[sewell2016community] proposed a dynamic extension of LSM, assuming at each time point, each node belongs to one of a fixed number $K$ of clusters, and the membership may change over time. Denote the latent position of node $i$ at time $t$ as $Z_{i}^{t}$ and the corresponding membership as $M_{i}^{t}$. Let $Z^{t}=({Z_{1}^{t}}^{\prime},\ldots,{Z_{n}^{t}}^{\prime})$ and $M_{t}=(M_{1}^{t},\ldots,M_{n}^{t})$, with the additional assumption that given the latent positions $Z^{t}$, $Y^{t}$ and $Y^{s}$, $s\neq t$, are conditionally independent, the decomposition follows:

[TABLE]

Two assumptions are made on the latent positions and the cluster assignments. First, the cluster assignments are assumed to follow a Markov process, i.e.,

[TABLE]

Second, given the current cluster assignment and all previous cluster assignments and latent positions, we assume the current latent positions depend only on the previous latent position and the current cluster assignments, i.e.,

[TABLE]

4.2 Identification and label switching

The cluster label is identified up to a permutation. We adopted the hamming error definition to define a correct membership, which is the minimal clustering errors over a permutation $\pi$ of the set $\{1,2,\ldots,K\}$. See section $2.6$ of [jin2015score] for details. Here we want to make explicitly an assumption that the label switching problem between two consecutive time steps can be solved by the permutation of the hamming error. This assumption could be easily satisfied when the temporal changes are not very dramatic.

5 Experiments

5.1 Simulation setting

We use Dynamic Stochastic Block Model (DSBM) [matias2017dynamic] as a comparison. The dissolving rate is set to $10$ time steps. The initial density is set to be $0.1$ and kept roughly stable if without the transitivity effect. Other parameters are … All the simulations are carried out using the R package Simulator [bien2016simulator].

One example of the simulated longitudinal networks, with some summary statistics are shown in Figure 2:

5.2 Transition matrix estimate

In a lot of applications, the overall trend of the communities change, for example, the rising of a certain group. The following river plot is an illustration of how nodes “flow” from one “river” into another. The height of each bar, in other words, the width of each river, represents the relative size of that community, illustrated in Figure 3.

5.3 Evolving membership estimate

We investigate the effects of sample size (per cluster) and the transitivity strength under four scenarios, which are the combination of two factors: the time transition rate (slow v.s. quick) and within cluster density (easy when the difference of within- and between- cluster density is large). The results in Figure 4 indicate that when it is slow and easy, all methods perform well, while under quicker or harder situations, the dynamic latent space model perform better as it is a better working model.

If we dig into the longitudinal mis-clustering rates, we can see that DSBM has a non-negligible probability that error rate could not be controlled, while DLSM can as shown in Figure 5.

5.4 Cluster specific ERGM estimate

The TERGM estimates conditioned on the membership estimates in the first step is carried out by [krivitsky2014jrssb]. Results illustrated in Figure 6.

5.5 Goodness of fit

To answer the natural question that if my working model is already good enough, why bother fitting within-cluster TERGMs, we show Figure 7.

5.6 Link prediction

As in a lot of machine learning applications, the out-of-sample prediction accuracy is the most trustworthy metric. So we conducted a one time step out-of-sample prediction and calculated the Area Under the ROC curves (AUC) in Figure 8.

6 Discussion and future work

As the experiments show, the clustering performance is not perfect, and the mis-clustered nodes have huge impact on the second stage TERGM inference. We are currently working on an augmented Bayesian method to jointly infer the clustering and local structures. The tool we use is the newly developed probabilistic programming language Edward [tran2016edward, tran2017deep], which integrates the Google Tensorflow [abadi2016tensorflow] framework with its capability to utilize GPU.