Weighted stochastic block model.

We propose a weighted stochastic block model (WSBM) which extends the stochastic block model to the important case in which edges are weighted. We address the parameter estimation of the WSBM by use of maximum likelihood and variational approaches, and establish the consistency of these estimators. The problem of choosing the number of classes in a WSBM is addressed. The proposed model is applied to simulated data and an illustrative data set.

Introduction

Networks are used in many scientific disciplines to represent interactions among objects of interest. For example, in the social sciences, a network typically represents social ties between actors. In biological sciences, a network can represent interactions between proteins.

The stochastic block model (SBM) (Holland et al. 1983; Snijders and Nowicki 1997) is a popular generative model which partitions vertices into latent classes. Conditional on the latent class allocations, the connection probability between two vertices depends only on the latent classes to which the two vertices belong. Many extensions of the class SBM have been proposed which include the degree correlated SBM (Karrer and Newman 2011; Peng and Carvalho 2016), mixed membership SBM (Airoldi et al. 2008) and overlapping SBM (Latouche et al. 2011).

The SBM and many of its variants are usually restricted to Bernoulli networks. However, many binary networks are produced after applying a threshold to a weighted relationship (Ghasemian et al. 2016) which results in the loss of potentially valuable information. Although most of the literature has focused on binary networks, there is a growing interest in weighted graphs (Barrat et al. 2004; Newman 2004; Peixoto 2018).

In particular, a number of clustering methods have been proposed for weighted graphs including algorithm based and model based methods. Algorithm based methods for clustering of weighted graphs can be further divided into two classes: algorithms which do not explicitly optimize any criteria (Pons and Latapy 2005; von Luxburg 2007) and those directly optimize a criterion (Clauset et al. 2004; Stouffer and Bascompte 2011). Model based methods (Mariadassou et al. 2010; Aicher et al. 2013, 2015; Ludkin 2020) attempt to take into account the random variability in the data. A recent review of graph clustering methods is given by (Leger et al. 2014).

Mariadassou et al. (2010) presents a Poisson mixture random graph model for integer valued networks and proposes a variational inference approach for parameter estimation. The model can account for covariates via a regression model. In Zanghi et al. (2010), a mixture modelling framework is considered for random graphs with discrete or continuous edges. In particular, the edge distribution is assumed to follow an exponential family distribution. Aicher et al. (2013) proposed a general class of weighted stochastic block model for dense graphs where edge weights are assumed to be generated according to an exponential family distribution. In particular, their construction produces complete graphs, in which every pair of vertices is connected by some real-valued weight. Since most real-world networks are sparse, the constructed model cannot be applied directly. To address this shortcoming, Aicher et al. (2015) extends the work of Aicher et al. (2013) and models the edge existence using a Bernoulli distribution and the edge weights using an exponential family distribution. The contributions of edge-existence distribution and edge-weight distribution in the likelihood function are then combined via a simple tuning parameter. However, their construction does not result in a generative model and it is not obvious how to simulate network observations from the proposed model. More recently, Ludkin (2020) presents a generalization of the SBM which allows artbitrary edge weight distributions and proposes a reversible jump Markov chain Monte Carlo sampler for estimating the parameters and the number of blocks. However, the use of continuous probability distribution to model the edge weights implies that the resulting graph is complete whereby every edge is present. This assumption is unrealistic for many applications whereby a certain proportion of the real-valued edges is 0. Haj et al. (2020) presents a binomial SBM for weighted graphs and proposes a variational expectation maximization algorithm for parameter estimation.

In this paper, we propose a weighted Stochastic Block model (WSBM) with gamma weights which aims to capture the information of weights directly using a generative model. Both maximum likelihood estimation and variational methods are considered for parameter estimation where consistency results are derived. We also address the problem of choosing the number of classes using the Integrated Completed Likelihood (ICL) criteria (Biernacki et al. 2000). The proposed models and inference methodology are applied to an illustrative data set.

Model specification

In this section, we present the weighted stochastic block model in detail and introduce the main notations and assumptions.

We let denote the set of directed weighted random graphs where is the set of countable vertices, is the set of edge-existence adjacency matrix, and is the set of weighted adjacency matrix. Given a random adjacency matrix , if an edge exists from vertex i to vertex j and otherwise. The associated weighted random adjacency matrix is given by: for , if , , and otherwise. Let be a probability measure on .

Generative model

We now describe the procedure of generating a sample of random graph (V, X, Y) with n vertices from .The generative framework described above is a straightforward extension of the binary stochastic block model whereby a positive weight is generated according to a gamma distribution for each edge. In particular, (X, Z) is a realization of the binary directed SBM. The gamma distribution is chosen due to its flexibility in the sense that, depending on the value of its shape parameter, it can represent distributions of different shapes.

Let be the vector of latent block allocations for the vertices, and set with . For each vertex , draw its block label from a multinomial distribution

Let be a matrix with entries in [0, 1]. Conditional on the block allocations , the entries for of the edge-existence adjacency matrix is generated from independently a Bernoulli distribution

Let and be matrices with entries taking values in the positive reals. Conditional on the latent block allocations and edge-existence adjacency matrix X, the weighted adjacency matrix is generated independently from where denotes the gamma distribution and is the Dirac delta function.

The log-likelihood of the observations and is given bywhere the sum is over all possible latent block allocations, is the probability of latent block allocation , andis the complete data log-likelihood, where is the gamma probability density function with shape parameter a and rate parameter b,

Assumptions

We present several assumptions needed for identifiability and consistency of maximum likelihood estimates. The following four assumptions were presented in Celisse et al. (2012) and are needed in this paper.

Assumption 1

For every , there exists such that

Assumption 2

There exists such that for all

Assumption 3

There exists such that for all ,

Assumption 4

There exists and such that for all , for all ,where and is any realized block allocation under the WSBM.

(A1) requires that no two classes have the same connectivity probabilities. If this assumption is violated, the resulting model has too many classes and is non-identifiable. (A2) requires that the connectivity probability between any two classes strictly lies within a closed subset of the unit interval. Note that this assumption is slightly more restrictive compared to assumption 2 of Celisse et al. (2012) in that we do not consider the boundary cases where . The boundary cases require special treatment and are not pursued in this paper. (A3) ensures that no class is empty with high probability while (A4) is the empirical version of (A3). We note that (A4) is satisfied asymptotically under the generative framework in Sect. 2.1 since the block allocations are generated according to a multinomial distribution.

In addition to the four assumptions above, we also have the following constraints on the gamma parameters.

Assumption 5

For every , there exists such that(A5) requires that no two classes have the same weight distribution. This assumption is the exact counterpart of (A1).

The log-likelihood function (1) contains degeneracies that prevent the direct estimation of parameters . To see this, we note that the probability density function of a gamma distribution is given byBy Stirling’s formula, we haveSetting ,Therefore, letting while keeping , we have . One can therefore show that the log-likelihood function is unbounded above. To avoid likelihood degeneracy, we compactify the parameter space. That is, we restrict the parameter space to a compact subset which contains the true paraemters. Therefore, we have the following assumption.

Assumption 6

There exists and such that for all ,With this assumption, it is easy to see that the log-likelihood function is bounded for any sample size.

Identifiability

Sufficient conditions for identifiability of binary SBM with two classes have been first obtained by Allman et al. (2009). Celisse et al. (2012) show that the SBM parameters are identifiable up to a permutation of class labels under the conditions that has distinct coordinates and . The condition on is mild since the set of vectors violating this assumption has Lebesgue measure 0. The identifiability of weighted SBM is more challenging where the only known result (Section 4 of Allman et al. (2011)) requires all entries of to be distinct. We note that the assumptions in the previous section are not necessarily sufficient but are necessary to ensure that the identifiability of the parameters.

Asymptotic recovery of class labels

We study the posterior probability distribution of the class labels given the random adjacency matrix and weight matrix , which is denoted by . Since and are random, is also random.

Let be the true conditional distribution of which depends on the true parameters . We study the convergence rate of towards 1 with respect to .

The matrices are permutation-invariant if one permutes both its rows and columns according to some permutation . Let be the matrices defined byand define the setTwo vectors of class labels z and are equivalent if there exists such that for all i. We let [z] denote the equivalence class of z and will omit the square-brackets in the equivalence class notation as long as no confusion arises.

The following result extends Theorem 3.1 of Celisse et al. (2012) to the case of WSBM.

Theorem 1

Under assumptions (A1)–(A6), for every ,uniformly with respect to , and for some depending only on but not on . Here with . Furthermore, can be replaced by under assumptions (A1)–(A3) and (A5)–(A6).

Maximum likelihood estimation of WSBM parameters

For the binary SBM, consistency of parameter estimation have been shown for profile likelihood maximization (Bickel and Chen 2009), spectral clustering method (Rohe et al. 2011), method of moments approach (Bickel et al. 2011), method based on empirical degrees (Channarond et al. 2012), and others (Choi et al. 2012). Consistency of both maximum likelihood estimation and variational approximation method are established in Celisse et al. (2012) and Bickel et al. (2013) where asymptotic normality is also established in Bickel et al. (2013). Abbe (2018) reviews recent development in the stochastic block model and community detections.

Ambroise and Matias (2012) proposes a general class of sparse and weighted SBM where the edge distribution may exhibit any parametric form and studies the consistency and convergence rates of various estimators considered in their paper. However, their model requires the edge existence parameter to be constant across the graph, that is, for all q, l, or can be modelled as where is the indicator function and . Furthermore, they also assume that conditional on the block assignments, the edge weight is modelled using a parametric distribution with a single parameter . They further impose the restriction thatThese assumptions are more restrictive than those imposed in this paper. Jog and Loh (2015) studies the problem of characterizing the boundary between success and failure of MLE when edge weights are drawn from discrete distributions. More recently, Brault et al. (2020) studies the consistency and asymptotic normality of the MLE and variational estimators for the latent block model which is a generalization of the SBM. However, the model considered in Brault et al. (2020) is restricted to the dense setting and requires the observations in the data matrix to be modelled by univariate exponential family distributions.

This section addresses the consistency of the MLE of WSBM. In particular, we extend the results obtained in the pioneering paper of Celisse et al. (2012) to the case of weighted graphs. Our proof closely follows the proof of consistency of the MLE in Celisse et al. (2012). The MLE consistency proof of and require different treatments since there are edges but only n vertices. The following result established the MLE consistency of .

Theorem 2

Assume that assumptions (A1), (A2), (A3), (A5), (A6) hold. Let us define the MLE of byThen for any metric on ,

Under additional assumption on the rate of convergence of the estimators of , consistency of can be established.

Theorem 3

Let denote the MLE of and assume that , , and , thenfor any metric d in .

Variational estimators

Direct maximization of the log-likelihood function is intractable except for very small graphs since it involves a sum over terms. In practice, approximate algorithms such as Markov Chain Monte Carlo (MCMC) and variational inference algorithms are often used for parameter inference. For the SBM, both MCMC and variational inference approaches have been proposed (Snijders and Nowicki 1997; Daudin et al. 2008). Variational inference algorithms have also been developed for mixed membership SBM (Airoldi et al. 2008), overlapping SBM (Latouche et al. 2011), and the weighted SBM proposed in Aicher et al. (2013). This section develops a variational inference algorithm for the WSBM which can be considered a natural extension of the algorithm proposed in Daudin et al. (2008) for the SBM.

The variational method consists in approximating by a product of n multinomial distributions. Let denote a set of product multinomial distributionswhereFor any , the variational log-likelihood is defined byHere denotes the Kullback-Leibler divergence between two probability distributions, which is nonnegative. Therefore, provides a lower bound on the log-likelihood function. We have thatwhere denotes the probability mass function of a Bernoulli distribution with parameter p, and recall that denotes the density function of a gamma distribution .

The variational algorithm works by iteratively maximizing the lower bound with respect to the approximating distribution , and estimating the model parameters. Maximization of with respect to consists of solvingwhere can be replaced by plug-in estimates. This has a closed-form solution given byConditional on , the variational estimators of are found by solving,Closed-form updates for and exist and are given byOn the other hand, updates for and do not have a closed form since the maximum likelihood estimators of the two parameters of a gamma distribution do not have closed forms. However, using the fact that a gamma distribution is a special case of a generalized gamma distribution, Ye and Chen (2017) derived simple closed-form estimators for the two parameters of gamma distribution. The estimators were shown to be strongly consistent and asymptotically normal. For , let us define the quantitiesthen the updates for are given byWe obtain the variational estimators by computing (3), (4), (5), (6), (7) until convergence.

We now address the consistency of the variational estimators derived above. The following two propositions are the counterpart of Theorem 2 and 3 for variational estimators. We omit the proof since they follow similar arguments as the proof of Corollary 4.3 and Theorem 4.4 of Celisse et al. (2012).

Proposition 1

Assume that assumptions (A1), (A2), (A3), (A5) and (A6), and let be the variational estimators defined above. Then for any distance on ,

Proposition 2

Assume that the variational estimators converge at rate 1/n to , respectively, and assumptions (A1), (A2), (A3), (A5), (A6) hold. We havewhere d denotes any distance between vectors in .

Note a stronger assumption on the convergence rate 1/n of is assumed for Proposition 2 compared to in Theorem 3. The same assumption is also used in Theorem 4.4. of Celisse et al. (2012).

Choosing the number of classes

In real world applications, the number of classes is typically unknown and needs to be estimated from the data. For the SBM, a number of methods have been developed to determine the number of classes, including log-likelihood ratio statistic (Wang and Bickel 2017), composite likelihood (Saldaña et al. 2017), exact integrated complete data likelihood (Côme and Latouche 2015) and Bayesian framework (Yan 2016) based methods. Model selection for variants of SBM have also been investigated (Latouche et al. 2014).

We apply the integrated classification likelihood (ICL) criterion developed by Biernacki et al. (2000) to choose the number of classes for the WSBM. The ICL is an approximation of the complete-date integrated likelihood. The ICL criterion for SBM have been derived by Daudin et al. (2008) under the assumptions that the prior distribution of factorizes and a non-informative Dirichlet prior on . Here we follow the approach of Daudin et al. (2008) to derive an approximate ICL for the WSBM.

Let denote the model with Q blocks, the ICL criterion is an approximation of the complete-data integrated likelihood:where is the prior distribution of the parameters. Assuming a non-informative Jeffreys prior on , a Stirling approximation to the gamma function, and finally a BIC approximation to the conditional log-likelihood function, the approximate ICL can be derived. For a model with Q blocks, and assuming a non-informative Jeffreys prior on , the approximate ICL criterion is:where is the estimate of . The derivation of above follows exactly the same lines as the proof of Proposition 8 of Daudin et al. (2008).

Simulation

We validate the theoretical results developed in previous sections by conducting simulation studies. In particular, we investigate how fast parameter estimates of WSBM converge to their true values and the accuracy of posterior block allocations. Additionally, we investigate the performance of ICL in choosing the number of blocks.

Experiment 1 (two-class model)

For each fixed number of nodes n, 50 realizations of WSBM are generated based on the fixed parameter setting given in (8 and 9). The variational inference algorithm derived in Sect. 5 is then applied to estimate the model parameters and class allocations.

We can see from Table 1 that the estimated model parameters converge to their true values as the number of nodes increases while the posterior class assignment is accurate across any number of nodes. Table 2 shows the ICL criterion tends to select the correct number of classes, especially when the number of nodes is large.

Table 1

Convergence analysis of posterior class allocations and parameter estimates under the two-class model

n	∑i=1nI_z_i=zi∗/n	||θ-θ^∗||_2	||π-π^∗||_F	||α-α^∗||_F	||β-β^∗||_F
25	1	0.031	0.057	0.855	0.434
50	1	0.029	0.035	0.654	0.322
100	1	0.027	0.012	0.256	0.070
200	1	0.025	0.006	0.116	0.046
500	1	0.021	0.002	0.053	0.024

Table 2

Frequency of choosing Q blocks by ICL under different number of nodes n under the two-class model

n | Q	1	2	3	4	5
25	0	44	3	3	0
50	0	48	2	0	0
100	0	50	0	0	0
200	0	50	0	0	0
500	0	50	0	0	0

Experiment 2 (three-class model)

50 network realizations are obtained under the three-class model with parameter values given in (10 and 11) for a range of values n.We can see from Table 3 that the estimated parameters converge to their true values quickly as the number of nodes increases. The ICL criterion tends to overestimate the number of classes when the number of nodes is small, but consistently selects the correct model when the number of nodes is large (Table 4).

Table 3

Convergence analysis of posterior class allocations and parameter estimates under the three-class model

n	∑i=1nI_z_i=zi∗/n	||θ-θ^∗||_2	||π-π^∗||_F	||α-α^∗||_F	||β-β^∗||_F
25	0.961	0.116	0.178	7.86	136.72
50	1	0.039	0.039	1.256	7.866
100	1	0.033	0.026	0.481	1.426
200	1	0.014	0.024	0.228	1.011
500	1	0.003	0.006	0.197	0.222

Table 4

Frequency of choosing Q blocks by ICL under different number of nodes n under the three-block model

n|Q	1	2	3	4	5
25	0	3	37	8	2
50	0	0	43	7	0
100	0	0	50	0	0
200	0	0	49	1	0
500	0	0	50	0	0

Computational complexity

The computational complexity of the variational estimators derived in Sect. 5 scales as , which has the same complexity as the variational algorithm developed by Daudin et al. (2008). Therefore, the algorithm may be prohibitively expensive for networks with more than 1000 nodes. The estimated computational time for the two-class model in Sect. 7.1 and for the three-class model in Sect. 7.2 for various values of n are shown in Fig. 1. The estimated computing time at each n is the average running time of the variational algorithm over 20 replications.

Fig. 1

Estimated computing time for the two-class and three-class models

Application: Washington bike data set

We apply the WSBM to analyse the Washington bike sharing scheme data set.1 Information with respect to start stations and end stations of trips as well as length (travel time) of trips are available in the data set. We select a time window of one week staring from January 10th, 2016 and construct the adjacency matrix X and weight matrix Y as follows:The resulting network consists of 370 nodes with an average out-degree of 36.14. The average total length of trips between any pair of stations is 42.15 minutes. We apply the ICL criterion to select the number of classes for the WSBM. For each number of classes Q, the variational inference algorithm is fitted to the network 20 times with 20 random initializations and the highest value of ICL is recorded, and the six-class model is chosen (Table 5). Each bike station is plotted on the map in Fig. 2 where its colour represents the estimated class assignment. We observe that bike stations in class 6 (colored in brown) tend to be concentrated in the central area of Washington wheraas stations in class 3 (colored in red) tend to be located further from the center. Figure 2 shows some spatial effect in the class assignment of bike stations whereby stations that are close in distance tend to be in the same cluster, with the exception of class 3 (colored in red). One potential extension of the model is to take into account the spatial locations of the bike stations as covariates.

Table 5

Model selection for the Washington Bike dataset using ICL criterion

Q	ICL
1	− 107157.24
2	− 93505.61
3	− 91688.33
4	− 90813.17
5	− 90300.90
6	− 89061.38
7	− 89326.35

Fig. 2

Bike stations. Class 1: blue. Class 2: green. Class 3: red. Class 4: cyan. Class 5: black. Class 6: brown (color figure online)

if there is trip starting from station i and finishing at station j.

is the total length of trips (in minutes) from station i to station j.

The estimated class proportions shown in 12 indicate that class 3 has the largest number of stations whereas class 4 has the smallest number of stations. The estimated shows that within class connectivity is generally higher compared to between class connectivity. We further observe that the connection probabilities between bike stations in class 1, 4, and 6 are substantially higher. Interestingly, we observe a near symmetry in the matrix indicating that the probability of having a trip from a station in class k to another station in class l is similar with the probability of having a trip from class l to class k.

The estimated densities of travel time between each pair of classes are shown in Fig. 3. We observe that the majority of the estimated densities have mode and mean close to 0, particularly for the estimated densities in the diagonal of Fig. 3. This implies that the total travel times between stations in the same class are quite short. In comparison, the total travel time between stations in different classes tend to be longer. This is reasonable as the distance between bike stations in different classes tend to be longer which in turn requires longer travel time.

Fig. 3

Estimated densities of total travel time between stations

Discussion

This paper proposes a weighted stochastic block model (WSBM) for networks. The proposed model is an extension of the stochastic block model. A variational inference strategy is developed for parameter estimation. Asymptotic properties of maximum likelihood estimators and variational estimators are derived, and the problem of choosing the number of classes is addressed by using an ICL criteria. Simulation studies are conducted to evaluate the performance of variational estimators and the use of ICL to determine the number of classes. The proposed model and inference methods are an illustrative data set.

It is straightforward to extend the WSBM to allow node covariates. Let be the covariates for each node , and let be the covariates for each pair of nodes, . The edge probability between a pair of nodes i, j with node i in block q and node j in block l can be modelled aswhere and . Conditional on block assignments and existence of an edge between a pair of nodes i, j, can be modelled as a Gamma random variable with mean and variance wherewhere and .

Many possible future extensions are possible. First, it is desirable to investigate further theoretical properties of maximum likelihood and variational estimators of WSBM parameters such as asymptotic normality of the estimators. Furthermore, some of the assumptions imposed in this work in order to ensure consistency of the estimators maybe relaxed. Moreover, the number of blocks is assumed to be fixed in the asymptotic analysis of the estimators. It would be interesting to allow the number of blocks to grows as the number of nodes grows.