Dynamic stochastic blockmodels: Statistical models for time-evolving networks

Dynamic sto c hastic blo c kmo dels: Statistical mo dels for time-ev olving net w orks Kevin S. Xu ? and Alfred O. Hero I II Departmen t of Electrical Engineering and Computer Science, Univ ersity of Michigan, Ann Arbor, MI, USA {xukevin, hero}@umich.edu Abstract. Significan t efforts hav e gone into the dev elopment of statis- tical mo dels for analyzing data in the form of netw orks, such as so cial net works. Most existing work has focused on mo deling static net works, whic h represen t either a single time snapshot or an aggregate view o ver time. There has b een recent interest in statistical mo deling of dynamic networks , whic h are observ ed at m ultiple points in time and offer a ric her represen tation of many complex phenomena. In this paper, we prop ose a state-space mo del for dynamic netw orks that extends the well-kno wn sto chastic blockmo del for static net works to the dynamic setting. W e then prop ose a procedure to fit the mo del using a modification of the extended Kalman filter augmen ted with a lo cal search. W e apply the pro cedure to analyze a dynamic so cial netw ork of email communication. Keyw ords: dynamic net work, sto c hastic blo c kmo del, state-space model 1 In tro duction Man y complex ph ysical, biological, and so cial phenomena are naturally repre- sen ted by net works. T remendous efforts ha ve b een dedicated to analyzing net- w ork data, which has led to the developmen t of many formal statistical mo dels for netw orks. Most research has fo cused on static netw orks, which either repre- sen t a single time snapshot of the phenomenon being inv estigated or an aggregate view o ver time. As such, statistical models for static net w orks ha ve a long history in statistics and so ciology among other fields [2]. How ever, most complex phe- nomena, including so cial b eha vior, are time-v arying, whic h has led researchers to consider dynamic, time-ev olving netw orks. In this pap er, we consider dynamic net works represented by a sequence of snapshots of the netw ork at discrete time steps. W e characterize suc h net works using a set of unobserved time-varying states from whic h the observ ed snapshots are derived. W e propose a state-space mo del for dynamic netw orks that com bines t wo t yp es of statistical mo dels: a static mo del for the individual snapshots and a temp oral mo del for the evolution of the states. The netw ork snapshots are mo d- eled using the stochastic blockmodel [5], a simple parametric mo del commonly ? Curren t affiliation: 3M Corp orate Researc h Lab oratory , St. Paul, MN, USA 2 Kevin S. Xu and Alfred O. Hero I II used in the analysis of static so cial net works. The state evolution is mo deled b y a stochastic dynamic system. Using a Cen tral Limit Theorem appro ximation, w e dev elop a near-optimal pro cedure for fitting the prop osed mo del in the on-line setting where only past and present netw ork snapshots are av ailable. The infer- ence pro cedure in volv es a mo dification of the extended Kalman filter, which is used for state tracking in many applications [3], augmented with a lo cal search strategy . W e apply the prop osed pro cedure to analyze a dynamic so cial netw ork of email comm unication and predict future email activity . 2 Related w ork Sev eral statistical mo dels for dynamic netw orks ha ve previously been prop osed b y extending a static mo del to the dynamic setting in a similar fashion to our prop osed mo del [2]. Two such mo dels include temp oral extensions of the expo- nen tial random graph mo del [1] and laten t space model [13]. More closely related to the state-space mo del we prop ose are several temp oral extensions of sto c has- tic blo c kmo dels (SBMs). SBMs divide no des in the netw ork into m ultiple classes and generate edges indep enden tly with probabilities θ ab dep enden t on the class mem b erships a, b of the no des [5]. Y ang et al. [15] propos e a dynamic SBM in- v olving a transition matrix that sp ecifies the probabilit y that a no de in class i at time t switches to class j at time t + 1 for all i, j, t and fit the mo del using Gibbs sampling and simu lated annealing. Ho et al. [4] prop ose a temp oral extension of a mixed-membership version of the SBM using linear state-space mo dels for the class membership vectors of no de clusters. One ma jor difference b etw een [4, 15] and this pap er is that w e treat the edge probabilities θ ab as time-varying states , while [4, 15] treat them as time-in v ariant parameters. In addition, our mo del allows for a simpler inference pro cedure using a Central Limit Theorem appro ximation. W e demonstrate the imp ortance of the time-v arying states for analysis of a dynamic so cial netw ork in Section 5. 3 Static sto c hastic blo c kmo dels W e first in tro duce notation and summarize the static sto c hastic blockmodel (SSBM), whic h w e use as the static model for the individual net work snap- shots. W e represent a dynamic netw ork by a time-indexed sequence of graphs, with W t = [ w t ij ] denoting the adjacency matrix of the graph observed at time step t . w t ij = 1 if there is an edge from no de i to no de j at time t , and w t ij = 0 otherwise. W e assume that the graphs are directed, i.e. w t ij 6 = w t j i in general, and that there are no self-edges, i.e. w t ii = 0. W ( s ) denotes the set of all snapshots up to time s , { W s , W s − 1 , . . . , W 1 } . The notation i ∈ a indicates that no de i is a member of class a . | a | denotes the n umber of no des in class a . The classes of all no des at time t is given b y a vector c t with c t i = a if i ∈ a at time t . W e denote the submatrix of W t corresp onding to the relations betw een no des in class a and class b by W t [ a ][ b ] . W e denote the v ectorized equiv alent of a matrix X , Dynamic sto c hastic blo c kmo dels 3 i.e. the v ector obtained by simply stac king columns of X on top of one another, b y x . Doubly-indexed subscripts such as x ij denote entries of matrix X , while singly-indexed subscripts suc h as x i denote en tries of the vectorized equiv alent x . Consider a snapshot at an arbitrary time step t . An SSBM is parameterized b y a k × k matrix Θ t = [ θ t ab ], where θ t ab denotes the probabilit y of forming an edge b et w een a no de in class a and a no de in class b , and k denotes the n umber of classes. The SSBM decomp oses the adjacency matrix in to k 2 blo c ks, where each block is asso ciated with relations b et ween no des in tw o classes a and b . Each blo c k corresp onds to a submatrix W t [ a ][ b ] of the adjacency matrix W t . Th us, given the class membership vector c t , eac h entry of W t is an indep enden t realization of a Bernoulli random v ariable with a blo ck-dependent parameter; that is, w t ij ∼ Bernoulli  θ t c t i c t j  . SBMs are used in t wo settings: 1. The a priori blo c kmo deling setting, where class memberships are known or assumed, and the ob jectiv e is to estimate the matrix of e dge pr ob abilities Θ t . 2. The a p osteriori blo c kmo deling setting, where the ob jectiv e is to simultane- ously estimate Θ t and the class mem b ership v ector c t . Since eac h entry of W t is indep enden t, the lik eliho o d for the SBM is given by f  W t ; Φ t  = Y i 6 = j  θ t c i c j  w t ij  1 − θ t c i c j  1 − w t ij = exp ( k X a =1 k X b =1  m t ab log  θ t ab  +  n t ab − m t ab  log  1 − θ t ab  ) , (1) where m t ab = P i ∈ a P j ∈ b w t ij denotes the n umber of observe d edges in block ( a, b ), and n t ab = ( | a || b | a 6 = b | a | ( | a | − 1) a = b (2) denotes the num b er of p ossible edges in blo c k ( a, b ) [6]. The parameters are given b y Φ t = Θ t in the a priori setting, and Φ t = { Θ t , c t } in the a p osteriori setting. In the a priori setting, a sufficien t statistic for estimating Θ t is the matrix Y t of blo ck densities (ratio of observ ed edges to p ossible edges within a blo ck) with en tries y t ab = m t ab /n t ab . Y t also happens to b e the maxim um-likelihoo d estimate of Θ t , whic h can b e sho wn [6] by setting the deriv ative of the logarithm of (1) to 0. Estimation in the a p osteriori setting is more in volv ed, and many metho ds ha ve b een prop osed, including Gibbs sampling [8], labe l-switc hing [6, 16], and sp ectral clustering [12]. The label-switching metho ds use a heuristic for solving the combinatorial optimization problem of maximizing the likelihoo d (1) ov er the set of p ossible class memberships, which is to o large to p erform an exhaustive searc h. 4 Kevin S. Xu and Alfred O. Hero I II 4 Dynamic sto c hastic blo c kmo dels W e prop ose a state-space mo del for dynamic netw orks that consists of a tem- p oral extension of the static sto c hastic blo ckmodel. First we presen t the mo del and inference procedure for a priori blo c kmo deling, and then we discuss the ad- ditional steps necessary for a p osteriori blo c kmo deling. The inference pro cedure is on-line, i.e. the state estimate at time t is formed using only observ ations from time t and earlier. 4.1 A priori blo c kmo dels In the a priori SSBM setting, Y t is a sufficien t statistic for estimating Θ t as discussed in Section 3. Thus in the a priori dynamic SBM setting, we can equiv- alen tly treat Y t as the observ ation rather than W t . The entries of W t [ a ][ b ] are indep enden t and iden tically distributed (iid) Bernoulli ( θ t ab ); th us by the Cen tral Limit Theorem, the sample mean y t ab is appro ximately Gaussian with mean θ t ab and v ariance ( σ t ab ) 2 = θ t ab (1 − θ t ab ) /n t ab , where n t ab w as defined in (2). W e assume that y t ab is indeed Gaussian for all ( a, b ) and p osit the linear observ ation mo del Y t = Θ t + Z t , where Z t is a zero-mean iid Gaussian noise matrix with v ariance ( σ t ab ) 2 for the ( a, b )th entry . In the dynamic setting where past snapshots are av ailable, the observ ations w ould b e given by the set Y ( t ) . The set Θ ( t ) can then b e viewed as states of a dynamic system that is generating the noisy observ ation sequence. W e complete the mo del b y sp ecifying a mo del for the state evolution o ver time. Since θ t ab is a probability and m ust b e b ounded b et ween 0 and 1, we instead w ork with the matrix Ψ t = [ ψ t ab ] where ψ t ab = log( θ t ab ) − log (1 − θ t ab ), the logit of θ t ab . A simple mo del for the state evolution is the random walk ψ t = ψ t − 1 + v t , where ψ t is the v ector representation of the matrix Ψ t , and v t is a random v ector of zero-mean Gaussian entries, commonly referred to as pro cess noise, with co v ariance matrix Γ t . The entries of the pro cess noise vector are not necessarily indep enden t or identically distributed (unlike the entries of Z t ) to allow for states to evolv e in a correlated manner. The observ ation mo del can then b e written in terms of ψ t as 1 y t = h  ψ t  + z t , (3) where the function h : R k 2 → R k 2 is defined by h i ( x ) = 1 / (1 + e − x i ), i.e. the logistic function applied to each entry of x . W e denote the cov ariance matrix of z t b y Σ t , which is a diagonal matrix 2 with entri es given b y ( σ t ab ) 2 . A graphical represen tation of the proposed model for the dynamic netw ork is shown in Fig. 1. 1 Note that w e hav e conv erted the block densities Y t and observ ation noise Z t to their resp ectiv e vector representations y t and z t . 2 The indices ( a, b ) for ( σ t ab ) 2 are conv erted into a single index i corresponding to the v ector representation z t . Dynamic sto c hastic blo c kmo dels 5 Logistic SBM . . . Fig. 1. Graphical representation of the prop osed mo del. The rectangular b o xes denote observ ed quantities, and the ov als denote unobserved quantities. The logistic SBM refers to applying the logistic function to each en try of Ψ t to obtain Θ t then generating W t using Θ t and c t . T o p erform inference on this mo del, we assume the initial state is Gaussian distributed, i.e. ψ 0 ∼ N  µ 0 , Γ 0  , and that { ψ 0 , v 1 , . . . , v t , z 1 , . . . , z t } are mu- tually indep enden t. If (3) w as linear in ψ t , then the optimal estimate of ψ t in terms of minim um mean-squared error would b e giv en by the Kalman filter [3]. Due to the non-linearity , we apply the extended Kalman filter (EKF), which linearizes the dynamics ab out the predicted state and provides an near-optimal estimate of ψ t . The predicted state under the random walk mo del is simply ˆ ψ t | t − 1 = ˆ ψ t − 1 | t − 1 with cov ariance R t | t − 1 = R t − 1 | t − 1 + Γ t . Let J t denote the Ja- cobian of h ev a luated at the predicted state ˆ ψ t | t − 1 . The EKF up date equations are as follo ws [3]: Near-optimal Kalman gain: K t = R t | t − 1  J t  T h J t R t | t − 1  J t  T + Σ t i − 1 P osterior state estimate: ˆ ψ t | t = ˆ ψ t | t − 1 + K t h y t − h  ˆ ψ t | t − 1 i P osterior estimate cov ariance: R t | t =  I − K t J t  R t | t − 1 The p osterior state estimate ˆ ψ t | t pro vides a near-optimal fit to the mo del at time t given the observed sequence W ( t ) . How to c ho ose the hyperparameters  µ 0 , Γ 0 , Σ t , Γ t  in an optimal manner is b eyond the scop e of this pap er and is discussed in [14, c hap. 5]. 4.2 A p osteriori blo c kmo dels In many applications, the class memberships c t are not known a priori and must b e estimated along with Ψ t . This can b e done using lab el-switching metho ds [6, 16], but rather than maximizing the likelihoo d, we maximize the p osterior state density given the entire sequence of observ ations W ( t ) up to time t to accoun t for the prior information. This is done by alternating b etw een lab el- switc hing and applying the EKF. 6 Kevin S. Xu and Alfred O. Hero I II The p osterior state density is given by f  ψ t | W ( t )  ∝ f  W t | ψ t , W ( t − 1)  f  ψ t | W ( t − 1)  . (4) By the conditional indep endence of current and past observ ations given the cur- ren t state, W ( t − 1) drops out of the first term in (4). It can thus b e obtained simply b y substituting h ( ψ t ) for θ t in (1). The second term in (4) is equiv alent to f  ψ t | y ( t − 1)  b ecause the class mem b erships at all previous time steps hav e already b een estimated. By applying the Kalman filter to the linearized tem- p oral mo del [3], f  ψ t | y ( t − 1)  ∼ N  ˆ ψ t | t − 1 , R t | t − 1  . Thus the logarithm of the p osterior density is given by log f  ψ t | W ( t )  = c − 1 2  ψ t − ˆ ψ t | t − 1  T  R t | t − 1  − 1  ψ t − ˆ ψ t | t − 1  + k X a =1 k X b =1  m t ab log  h  ψ t ab  +  n t ab − m t ab  log  1 − h  ψ t ab  , (5) where c is a constan t term indep enden t of ψ t that can b e ignored 3 . W e use the log-p osterior (5) as the ob jective function for lab el-switc hing. W e find that a simple lo cal search (hill climbing) algorithm [11] initialized using the estimated class memberships at the previous time step suffices, b ecause only a small fraction of no des change classes betw een time steps in most applications. A t the initial time step, w e employ the sp ectral clustering algorithm of Sussman et al. [12] for the SSBM as the initialization. 5 Application to Enron email net w ork W e demonstrate the prop osed pro cedure on a dynamic so cial netw ork con- structed from the Enron corpus [9, 10], which consists of ab out 0 . 5 million email messages b et ween 184 Enron employ ees from 1998 to 2002. W e place directed edges betw een emplo yees i and j at time t if i sends at least one email to j during w eek t . Each time step corresponds to a 1-week in terv al. W e make no distinction b et w een emails sen t “to”, “cc”, or “b cc”. In addition to the email data, the roles of most of the employ ees within the company (e.g. CEO, president, manager, etc.) are a v ailable, whic h we use as classes for a priori blockmodeling. Employ ees with unkno wn roles are placed in an “others” class. 5.1 State tracking W e b egin by examining the temp oral v ariation of the states, which we refer to as state tr acking . Recall that the states Ψ t corresp ond to the logit of the edge probabilities Θ t . W e first apply the a priori EKF to obtain the state estimates ˆ ψ t | t and their v ariances (the diagonal of R t | t ). Applying the logistic function, we can then obtain the estimated edge probabilities ˆ Θ t | t with confidence in terv als. 3 A t the initial time step, ˆ ψ 1 | 0 = µ 0 and R 1 | 0 = Γ 0 + Γ 1 . Dynamic sto c hastic blo c kmo dels 7 Recipient class Sender class 1 2 3 4 5 6 7 1 2 3 4 5 6 7 0 0.1 0.2 0.3 0.4 0.5 0.6 (a) W eek 59: a normal week Recipient class Sender class 1 2 3 4 5 6 7 1 2 3 4 5 6 7 0 0.1 0.2 0.3 0.4 0.5 0.6 (b) W eek 89: CEO Skilling resigns Fig. 2. Estimated edge probabilit y matrices for tw o selected w eeks. En try ( i, j ) denotes the estimated probability of an edge from class i to class j . Classes are as follo ws: (1) directors, (2) CEOs, (3) presidents, (4) vice-presidents, (5) managers, (6) traders, and (7) others. Notice the increase in the probability of edges from CEOs during the week of Skilling’s resignation. 0 20 40 60 80 100 120 0 0.2 0.4 0.6 Week Edge probability (a) Presidents to presidents 0 20 40 60 80 100 120 0 0.05 0.1 0.15 Week Edge probability (b) Others to others Fig. 3. A priori EKF estimated edge probabilities ˆ θ t | t ab (solid lines) with 95% confidence in terv als (shaded region) for selected a, b b y week. An increase in edge probabilities b et w een Enron presidents (a) o ccurs prior to a similar increase b et ween those in other roles (b) suggesting insider knowledge. 8 Kevin S. Xu and Alfred O. Hero I II Examining the temporal v ariation of ˆ Θ t | t rev eals some in teresting trends. F or example, a large increase in the probabilities of edges from CEOs is found at w eek 89. This is the w eek in whic h CEO Jeffrey Skilling resigned and is confirmed to b e the cause of the increased probabilities b y examining the con tent of the emails. Fig. 2 shows a comparison of the matrix ˆ Θ t | t during a normal week and during the w eek Skilling resigned. Another in teresting trend is highlighted in Fig. 3, where the temp oral v ari- ation of tw o selected edge probabilities ov er the en tire data trace with 95% confidence interv als is shown. Edge probabilities b et ween Enron presiden ts show a steady increase as Enron’s financial situation worsens, hinting at more fre- quen t and widespread insider discussions, while emails b et ween others (not of one of the six known roles) b egin to increase only after Enron falls under federal in vestigation. A key observ ation from this analysis is the imp ortance of mo deling the edge probabilities as time-v arying states, as opp osed to time-inv arian t parameters as in [4, 15]. Indeed the temp or al variation of the edge probabilities is what rev eals the internal dynamics of this time-ev olving social net work. F urthermore, the temp oral mo del pro vides estimates with less uncertaint y than the static SBM, with 95% confidence in terv als that are 24% narrow er on av erage. 5.2 Dynamic link prediction Next w e turn to the task of dynamic link prediction, which differs from static link prediction [7] b ecause the link predictor m ust simultaneously predict the new edges that will b e formed at time t + 1, as well as the current edges (as of time t ) that will disapp ear at time t + 1, from the observ ations W ( t ) . The latter task is not addressed b y most static link prediction metho ds in the literature. Since the SBM assumes sto c hastic equiv alence b etw een no des in the same class, the EKF alone is only a goo d predictor of the block densities Y t , not the edges themselv es. How ever, the EKF can b e combined with a predictor that operates on individual edges to form a link predictor. A simple individual- lev el predictor is the exp onen tially-weigh ted moving a verage (EWMA) given by ˆ W t +1 = λ ˆ W t + (1 − λ ) W t . Using a conv ex com bination of the EKF and EWMA predictors, w e obtain a b etter link predictor that incorp orates b oth blo c k-level c haracteristics (through the EKF) and individual-lev el c haracteristics (through the EWMA). This can b e seen from the receiver op erating c haracteristic (ROC) curv es in Fig. 4. The a p osteriori EKF slightly outp erforms the a priori EKF b ecause the a p osteriori EKF finds a b etter fit to the dynamic SBM via a b etter assignmen t of no des to classes than the a priori (assumed) assignment. 6 Conclusion This pap er prop oses a statistical mo del for dynamic net works that utilizes a set of unobserved time-v arying states to c haracterize the dynamics of the net work. The prop osed mo del extends the well-kno wn sto c hastic blo c kmo del for static Dynamic sto c hastic blo c kmo dels 9 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 False positive rate True positive rate A priori EKF + EWMA A posteriori EKF + EWMA EWMA Fig. 4. Comparison of ROC curv es for link prediction on Enron data. T rue p ositiv e rate denotes the fraction of actual edges that are correctly predicted, and false p ositive rate denotes the fraction of non-edges that are predicted to b e edges. The conv ex com- bination of either EKF with the EWMA outp erforms the EWMA alone by accounting for blo c k-level characteristics. net works to the dynamic setting can b e used for either a priori or a p osteriori blo c kmo deling. The main contribution of the pap er is a near-optimal on-line in- ference pro cedure for the prop osed mo del using a mo dification of the extended Kalman filter, augmen ted with a local search. W e applied the proposed inference pro cedure to the Enron email netw ork and disco vered some interesting trends when w e examined the estimated states. One such trend w as a steady increase in emails b et ween Enron presidents as Enron’s financial situation worsened, while emails b et w een other employ ees remained at their baseline levels un til Enron fell under federal inv estigation. In addition, the prop osed pro cedure show ed promis- ing results for predicting future email activity . W e b eliev e the prop osed mo del and inference pro cedure can b e applied to rev eal the internal dynamics of many other dynamic net works. Ac knowledgmen ts. This work w as partially supp orted by the Army Research Office grant W911NF-12-1-0443. Kevin Xu was partially supp orted by an aw ard from the Natural Sciences and Engineering Researc h Council of Canada. References [1] Ahmed, A., Xing, E.P .: Recov ering time-v arying netw orks of dep endencies in so cial and biological studies. Pro c. Nat. Acad. Sci. 106(29), 11878–11883 (2009) 10 Kevin S. Xu and Alfred O. Hero I II [2] Golden b erg, A., Zheng, A.X., Fien b erg, S.E., Airoldi, E.M.: A survey of sta- tistical net work mo dels. F ound. T rends Mach. Learn. 2(2), 129–233 (2010) [3] Ha ykin, S.: Kalman filtering and neural netw orks. Wiley-Interscience (2001) [4] Ho, Q., Song, L., Xing, E.P .: Evolving cluster mixed-mem b ership blo ck- mo del for time-v arying net works. In: Pro ceedings of the 14th Int. Conf. Artif. In tell. Statist. (2011) [5] Holland, P .W., Laskey , K.B., Leinhardt, S.: Sto chastic blo c kmo dels: First steps. So c. Netw. 5(2), 109–137 (1983) [6] Karrer, B., Newman, M.E.J.: Sto c hastic blockmodels and communit y struc- ture in net works. Phys. Rev. E 83, 016107 (2011) [7] Lib en-No w ell, D., Klein b erg, J.: The link-prediction problem for so cial net- w orks. J. Am. So c. Inf. Sci. 58(7), 1019–1031 (2007) [8] No wicki, K., Snijders, T.A.B.: Estimation and prediction for sto c hastic blo c kstructures. J. Am. Stat. Asso c. 96(455), 1077–1087 (2001) [9] Prieb e, C.E., Conroy , J.M., Marchette, D.J., P ark, Y.: Scan statistics on Enron graphs. Comput. Math. Organ. Theory 11(3), 229–247 (2005) [10] Prieb e, C.E., Conroy , J.M., Marchette, D.J., P ark, Y.: Scan statistics on Enron graphs (2009), http://cis.jhu.edu/ ~ parky/Enron/enron.html [11] Russell, S.J., Norvig, P .: Artificial in telligence: A mo dern approach. Pren tice Hall, 2nd edn. (2003) [12] Sussman, D.L., T ang, M., Fishkind, D.E., Prieb e, C.E.: A consis- ten t adjacency spectral em b edding for sto c hastic blo ckmodel graphs. arXiv:1108.2228v3 [stat.ML] (2012) [13] W estveld, A.H., Hoff, P .D.: A mixed effects model for longitudinal relational and net w ork data, with applications to in ternational trade and conflict. Ann. Appl. Stat. 5(2A), 843–872 (2011) [14] Xu, K.S.: Computational metho ds for learning and inference on dynamic net works. Ph.D. thesis, Universit y of Mic higan (2012) [15] Y ang, T., Chi, Y., Zhu, S., Gong, Y., Jin, R.: Detecting communities and their ev olutions in dynamic so cial netw orks—a Bay esian approac h. Mach. Learn. 82(2), 157–189 (2011) [16] Zhao, Y., Levina, E., Zh u, J.: Consistency of communit y detection in net- w orks under degree-corrected stochastic blo c k models. The Annals of Statis- tics (in press) (2012)