An empirical Bayes approach to stochastic blockmodels and graphons: shrinkage estimation and model selection

Results on SBM: homogeneous block connectivity

We designed a constrained SBM that generates affiliation networks, i.e., two vertices within the same community connect with probability λ, and from different communities with probability ε < λ. We also added a parameter ρ ∈ (0,1] to control the sparsity of the graph. The corresponding true connectivity matrix is

$${\mathrm{\Theta }}^{\ast }=\rho {\left(\begin{array}{cccc}\lambda & \epsilon & \cdots & \epsilon \\ \epsilon & \lambda & \cdots & ⋮\\ ⋮& & \ddots & \epsilon \\ \epsilon & \cdots & \epsilon & \lambda \end{array}\right)}_{K^{\ast }\times K^{\ast }},$$where K* is the number of communities.

To generate dense graphs (model 1), we set λ = 0.9, ε = 0.1, and ρ = 1. We generated graphs with n = 200 vertices and the number of communities K* ∈ {10, 11, …, 18}. For each choice of K*, we generated 100 networks independently. For each network, all the nodes were randomly divided into K* clusters with equal probability 1/K*, and then connected according to the connectivity matrix Θ* and their cluster labels. Note that the simulated node clusters had very different sizes, ranging between 7 and 35, due to the high variance in block size.

We also used λ = 0.9, ε = 0.1 and ρ = 0.2 to generate sparse graphs (model 2), while keeping K* = 10 but changing the network size n ∈ {200, 250, 300,350, 400, 450}. For each network size n, we followed the same procedure as in model 1 and generated 100 networks independently. Here “sparse” refers to a lower edge density around 0.035, which is 20% of the graphs generated in model 1.

For a simulated graph, we applied the variational Bayes algorithm (Latouche, Birmele & Ambroise, 2012) with an input number of clusters K = 1, …, 20, from which we obtained K communities and a Bayesian estimate ${\hat{\mathrm{\Theta }}}_{VBEM}(K)$ of the connecting probabilities among the K × K blocks. Given the estimated communities by the variational Bayes algorithm, we found ${\hat{\mathrm{\Theta }}}_{MLE}(K)$ as in (4) and our empirical Bayes estimate ${\hat{\mathrm{\Theta }}}_{EB}(K)$ as in (11) and compared them to the VBEM estimate. As the estimates are functions of K, so are their MSEs as defined in (18). Let MSE_MLE(K) be the mean squared error of the MLE by plugging ${\hat{\mathrm{\Theta }}}_{MLE}(K)$ into (18), where each element $\hat{\mathrm{\Theta }}{}_{ij}^{\prime }$ is given by ${\hat{\mathrm{\Theta }}}_{MLE}(K)$ and the partition Z. Then we define $\tilde{K}$ as the number of clusters that minimizes the MSE of the MLE, i.e.,

(22) $$\tilde{K}=\underset{K}{\arg min}MS{E}_{MLE}(K)$$over the input range of K. For the 100 graphs generated under the same matrix Θ*, they share the same K* while each one of them defines a corresponding $\tilde{K}$. Both K* and $\tilde{K}$ were used in our comparisons on model selection criteria for the number of blocks. In particular, for a general graphon, K* may not be clearly defined and in such a case, $\tilde{K}$ serves as the reference for comparison.

For dense graphs (model 1), as shown in Fig. 2, we compared the MSEs (18) of the three estimates of Θ to the true connectivity matrix and presented the ratio of the MSE of our EB estimate to the MSEs of the MLE and VBEM estimate. For dense stochastic block models, the accuracy of MLE and that of VBEM were close, whereas EB gave better estimates for almost all K values, i.e., MSE ratios were smaller than 100%. We see a significantly smaller MSE ratio when K is close to K*, especially when K* is relatively small. For example, the MSE ratios EB/MLE and EB/VBEM were lower than 10% at K = K* when K* = 10, …, 15. When K* went bigger, such as K* = 17, 18 in the simulation, the $\tilde{K}$ for most of the graphs was less than K*, and the MSE ratios reached a minimum level at some K < K*, which was slightly above 50%.

Table 1 presents the model selection results on the simulated dense graphs from model 1, where we define E_K* and $E_{\tilde{K}}$ as the average deviation of the selected number of blocks $\hat{K}$ from K* and from $\tilde{K}$ respectively, i.e.,

(23) $$E_{K^{\ast }}={\displaystyle \frac{1}{M}\sum _{t=1}^M|{\hat{K}}_t-K^{\ast }|,E_{\tilde{K}}={\displaystyle \frac{1}{M}\sum _{t=1}^M|{\hat{K}}_t-{\tilde{K}}_t|,}}$$where t ∈ {1,…,M} is the index of the graphs generated under the same Θ*, ${\hat{K}}_t$ is the estimated number of clusters by a model selection criterion, and ${\tilde{K}}_t$ is the $\tilde{K}$ defined by (22) for the t-th graph. When K* was small, such as 10 ≤ K* ≤ 13, ${\mathcal{J}}_{VBEM}$ and ${\mathcal{J}}_{EB}$ gave the same results, and both accurately selected $\hat{K}=K^{\ast }$ as the optimal number of blocks. As K* increased, ${\mathcal{J}}_{EB}$ outperformed ${\mathcal{J}}_{VBEM}$, and was comparable to ${\mathcal{J}}_{CVRP}$ in terms of E_K*. In fact, for a limited graph size n = 200 here, the average number of vertices in each block will be smaller as K* increases, making it hard for small communities to be detected. Therefore, $\tilde{K}$ may better reflect the number of clusters that fit well the observed network. Considering this, we see ${\mathcal{J}}_{EB}$ had both smaller E_K* and $E_{\tilde{K}}$ than ${\mathcal{J}}_{VBEM}$ in general, which indicates the superiority of our model selection method. ${\mathcal{J}}_{CVRP}$ showed relatively stable performance in terms of E_K* and $E_{\tilde{K}}$, but the results were not satisfactory for small K*. In summary, from the simulation results on dense graphs (model 1), EB has demonstrated the highest estimation accuracy, especially when the clustering algorithm finds the true number of communities, and the EB model selection criterion generally selects the best model.

Detecting the true number of blocks for a sparse graph (model 2) is harder because of fewer edge connections in a block. Thus, we fixed K* = 10 and varied the network size n from 200 to 450. In terms of estimation accuracy, Fig. 3 shows that our EB estimate had better performance than MLE in almost all the cases (except when K = 1 under which the two estimates were identical), and the MSE ratio kept decreasing as K increased. In particular, for K = K* = 10, the MSE ratio of EB over MLE was about 95%. If the number of blocks is overestimated (say K > 15), the MSE ratio can drop to <90%. When compared to VBEM, for a small network size n and a small number of blocks K, EB estimates can be slightly less accurate (<5% increase in MSE), but as K increases and becomes close to K*, the MSE ratio decreases to the same level as that of EB over MLE. As reported in Table 2, for all the cases ${\mathcal{J}}_{EB}$ achieved the best model selection performance with the smallest E_K* and $E_{\tilde{K}}$ among the three methods. This highlights the usefulness of our model selection criterion for the more challenging sparse graph settings.

More detailed results for both models 1 and 2 in this simulation study can be found in the Supplemental Material.

Results on SBM: heterogeneous block connectivity

In this section, we show how the performance changes when heterogeneous block connectivity probabilities are used. We consider the following connectivity matrix

$${\mathrm{\Theta }}^{\ast }=\rho {\left(\begin{array}{cccc}{\lambda }_1& {\epsilon }_{12}& \dots & {\epsilon }_{1K^{\ast }}\\ {\epsilon }_{21}& {\lambda }_2& \dots & ⋮\\ ⋮& & \ddots & {\epsilon }_{(K^{\ast }-1)K^{\ast }}\\ {\epsilon }_{K^{\ast }1}& \dots & {\epsilon }_{K^{\ast }(K^{\ast }-1)}& {\lambda }_{K^{\ast }}\end{array}\right)}_{K^{\ast }\times K^{\ast }},$$where the values are sampled from uniform distributions.

Similar to model 1, to generate dense graphs (model 1s), we set ρ = 1, and drew λ_i ∼ U(0.5, 0.9) for i ∈ {1, …, K*} and ε_ij ∼ U(0.3, 0.5) for i, j ∈ {1, …, K*}, i ≠ j. We generated graphs with n = 200 vertices and the number of communities K* ∈ {10, 11, …, 18}. For each choice of K*, we generated 100 networks independently with parameters sampled from the above uniform distributions. For each network, the nodes were randomly divided into K* clusters with equal probability 1/K*, and then connected according to the connectivity matrix Θ* and their cluster labels. Overall, EB outperformed MLE and VBEM with respect to MSE. Different from model 1 where the smallest ratios of EB/MLE and EB/VBEM were observed at K = K* for most of the simulations, the MSE ratio of EB over MLE and VBEM decreases smoothly with the increase in K (Fig. 4). In terms of model selection, EB was better than VBEM when K* ≥ 14 and comparable to VBEM with smaller K, although the improvement was slightly less substantial. The detailed results are reported in Table 3.

We also used ρ = 0.2 and the same setting for λ and ε as above to generate sparse graphs (model 2s), while keeping K* = 10 but changing the network size n ∈ {200, 250, 300,350, 400, 450}. For each network size n, we followed the same procedure as in model 1s and generated 100 networks independently. The results are similar to the homogeneous case (model 2). The detailed results are provided in the Supplemental Material.

Results on graphon model

Following the same design as in Latouche & Robin (2016), we choose a graphon function

$$W(x,y)=\rho {\lambda }^2(xy{)}^{\lambda -1}$$with two parameters $\lambda \le 1/\sqrt{\rho }$. Here, ρ controls the sparsity of the graph, as the expected number of edges is proportional to ρ, and λ controls the concentration of the degrees, so that more edges will concentrate on fewer nodes if λ is large. We chose ρ ∈ {10⁻¹, 10^−1.5, 10⁻² } and λ ∈ {2, 3, 5}, and simulated graphs of size n = 100 (model 3) and of size n = 316 (≈10^2.5) (model 4). For each network, we used SBM approximation with the number of clusters K = 1, 2, …, 10. Using (22), we also defined $\tilde{K}$ as the number of blocks that minimizes the mean squared error (19) of the MLE, i.e., $\tilde{K}$ is the number of communities that best fits the observed network. Figure 5 shows the graphon function for some values of (ρ, λ). The parameter ρ controls the scale of the function, and thus the grophon functions reach the maximum height when ρ = 10⁻¹. While not shown in the figure, for ρ = 10^−1.5 or 10⁻² the functions are scaled down and have lower values. Meanwhile, λ controls the concentration of the function, such that a graphon defined by a higher value of λ shows a highly concentrated peak as for λ = 5 in the figure.

The MSE ratios between our EB estimate and the other two competing methods, MLE and VBEM, are shown in Fig. 6 for graphs of size n = 316. The results for n = 100 are similar and relegated to the Supplemental Material. In general, our EB method achieved higher accuracy with smaller MSEs than the other two methods. For most cases, our EB estimate was more accurate than the MLE, with the MSE ratios between 60% and 100%. Compared to VBEM, our EB estimate achieved substantially smaller MSEs with ratios below 20%. For both graph sizes, the improvement of the EB method over the other two competitors was especially significant when the graph was sparse (ρ small). In such a case, fewer connections between nodes are observed in a network, and there is a high probability to have zero edge within the cluster. For blocks with lower connectivity, MLE tends to underestimate their connectivity, while shrinkage helps the situation by borrowing information from other blocks.

The model selection results are reported in Table 4. Since the true number of communities under the graphon model is not clearly defined, we used $\tilde{K}$ as the ground-truth to evaluate model selection performance. For both n = 100 and n = 316, the mean absolute deviation $E_{\tilde{K}}$ (23) of the $\hat{K}$ selected by our criterion ${\mathcal{J}}_{EB}$ was either the smallest or was very close to the smallest value among the three methods. While EB and VBEM were generally comparable, CVRP showed unstable performance as its $E_{\tilde{K}}$ could be much larger than the other two methods in some cases (such as ρ = 10⁻¹ and ρ = 10^−1.5). See Supplemental Material for more detailed results.

To expand the scope of this study, we further compared the performance of our EB method on graphons with a non-Bayesian approach. A commonly used algorithm is network histogram approximation (NHA) developed by Olhede & Wolfe (2014). The authors showed the universality of graphon approximation through regular stochastic block model and introduced an automatic bandwidth selection rule to select the best block model to represent graphon functions. The method fist divides degree-sorted vertices into equal-sized groups and selects the histogram bandwidth that maximizes the likelihood under an SBM. Given the automatically selected histogram bandwidth, the model parameters are estimated by the MLE in (4). In the comparison, we substitute the MLE estimate with our EB estimate to see if it can improve the accuracy.

In our simulation we used NHA to estimation the graphon functions in model 3 and 4 with ρ = 10⁻¹, using the suggested parameter (c = 4 in Olhede & Wolfe (2014)) to select the NHA bandwidth. NHA did not work on the sparser cases since too many nodes have a degree of zero. Table 5 shows that EB indeed improved the graphon estimation by NHA as well. The MSE for each set of parameters shown in the table is the average results from 100 networks. For λ = 2 in which case the graphon function has lower variability, EB outperformed MLE substantially. For larger λ’s, the two methods had comparable accuracy.

We briefly summarize a few key observations from the simulation studies. It is seen that EB estimates had smaller MSEs than the other two methods in most of the cases above. For the dense SBM (model 1), the accuracy of EB estimate was much higher. The relative low variance in connectivity across different blocks led to higher degree of shrinkage and information sharing among the EB estimates. For the sparse SBM (model 2), heterogenous SBM (model 1s and 2s) and graphon models (model 3 and 4), EB showed moderate improvements over the two competing methods in general. When the graph is sparse, EB can be much more accurate than VBEM, as shown in Fig. 6. As for model selection, EB generally selected the number of clusters $\hat{K}$ that was closer to K* and $\tilde{K}$ in all the models above, which demonstrates the usefulness of our hierarchical model for deriving likelihood-based model selection criterion.

Alternative clustering and running time comparison

Our results and numerical comparisons were conducted to demonstrate the uniform accuracy improvement: By varying the input number of clusters so some cluster results could be very inaccurate, our EB estimates reached smaller MSEs for almost all the clustering results. To further demonstrate this point, we also applied our EB estimates after spectral clustering. As shown in Fig. 7, our method improved the parameter estimation accuracy as well: Under the same simulation setting as in Figs. 2 and 3, the EB/MLE MSE ratio shows a similar pattern to the results of the previous simulation in SBMs.

The computation of our EB method is only the maximization of the likelihood (8, 9). The objective is the sum of two separate functions. Thus, we just need to maximize two bi-variate functions, regardless of the problem size (n, K). In general, the computation time is negligible compared to the graph clustering step. Table 6 reports the average running times (in seconds) of spectral clustering (T_C) and our EB estimation (T_E) by BFGS for various network size n and number of communities K, on a single 2.6 GHz Intel i7 core. It is seen from the table that for large problems (n = 10,000, K = 500), the running time of EB estimation step was less than 1% of the runtime of spectral clustering.

Political blogs

First we consider the French political blogosphere network from Latouche, Birmelé & Ambroise (2011). The network is made of 196 vertices connected by 2,864 edges. It was built from a single day snapshot of political blogs automatically extracted on October 14th, 2006 and manually classified by the “Observatoire Presidentiel” project (Zanghi, Ambroise & Miele, 2008). In this network, nodes correspond to hostnames and there is an edge between two nodes if there is a known hyperlink from one hostname to the other. The four main political parties that are present in the data set are the UMP (french republican), liberal party (supporters of economic-liberalism), UDF (moderate party), and PS (french democrat). However, in the dataset annotated by Latouche, Birmelé & Ambroise (2011) there are K* = 11 different node labels in total, since they considered analysts as well as subgroups of the parties. The test statistic by the method in Yan et al. (2012) yielded a p-value of 0.08 according to the bootstrap distribution suggested by the authors, which indicates that the regular SBM is a fair representation of this network compared to DCSBM.

Given the known community memberships, we constructed a connectivity matrix Θ* = (θ*_ab)_K*×K* with entries

(24) $${\theta }_{ab}^{\ast }=X_{ab}^B/n_{ab},a,b\in \{1,\dots ,K^{\ast }\},$$where $X_{ab}^B$ is the number of edges observed in block (a, b), $n_{ab}=|Z_{true}^{-1}(a)|\cdot |Z_{true}^{-1}(b)|$ for a ≠ b and $n_{aa}=|Z_{true}^{-1}(a)|\cdot (|Z_{true}^{-1}(a)|-1)/2$, and K* is the true number of communities. Then the MSE (18) between an estimate $\hat{\mathrm{\Theta }}(K)$ and $\mathrm{\Theta }$* (24) were used as an accuracy metric to compare estimated connectivity matrices, where K is the input number of clusters.

We also used test data likelihood as another comparison metric. We randomly sampled 70% of the nodes, denoted by V_o, as observed training data, and estimated a connectivity matrix $\hat{\mathrm{\Theta }}=({\hat{\theta }}_{ij}{)}_{K^{\ast }\times K^{\ast }}$ from their edge connections and true memberships. Denote by V_t the test data nodes not used in the estimation. Recall that X_ij is the (i, j)th element in the adjacency matrix of the network. Then test data likelihood ${\mathcal{L}}_{test}$ was calculated according to (2) given the $\hat{\mathrm{\Theta }}$ estimated by a method,

(25) $$\begin{array}{cc}{\mathcal{L}}_{test}& =\prod _{i\in V_{}o,j\in V_{}t}{\hat{\theta }}_{z_i{z}_j}^{X_{ij}}{(1-{\hat{\theta }}_{z_i{z}_j})}^{1-X_{ij}}\times \prod _{k<j\in V_{}t}{\hat{\theta }}_{z_j{z}_k}^{X_{jk}}{(1-{\hat{\theta }}_{z_j{z}_k})}^{1-X_{jk}},\end{array}$$where z_i, z_j, z_k are the known ground truth labels (ground truth) of the nodes. Note that X_ij ∈ {0, 1} is the edge connection between a vertex i in the training data and a vertex j in the test data, while X_jk is the edge connection between two vertices j and k in the test data. We repeated this procedure 100 times independently to find the distribution of test data likelihood ${\mathcal{L}}_{test}$ across random sample splitting of the n nodes into V_o and V_t.

We applied VBEM to detect communities with an input number of clusters K = 1, 2, …, 20. The MSE ratios of EB over the other two competing methods were calculated and plotted against K in Fig. 8A. It is clear that EB achieved smaller MSE than the other two methods for all values of K. When K was close to or greater than K* = 11, EB provided more accurate estimates than both MLE and VBEM with smaller MSEs. Figure 8B shows the box-plot of test data log-likelihood values across 100 random sample splitting. From the box-plots, we see that the test data likelihood of EB was significantly higher than the other two estimates. These comparisons confirm that EB estimates were more accurate than the other two competing methods in terms of both metrics. In terms of model selection, Fig. 8C plots the three model selection criteria, ${\mathcal{J}}_{CVRP}$, ${\mathcal{J}}_{VBEM}$, ${\mathcal{J}}_{EB}$, over the input range of K. All three model selection criteria have been standardized to [0, 1] with a higher value indicating a better model, such that the best model is selected by the maximizer of each criterion. Accordingly, CVRP, VBEM and EB estimated $\hat{K}=1$, 12 and 10, respectively, while the true K* = 11. The $\hat{K}$ by VBEM and EB were both reasonably close to the ground-truth, while CVRP did not work well in this case. Figures 8D and 8E show the distributions of the shrinkage values η_ab at K = K*. We see that the diagonal blocks had higher shrinkage. Around 70% of the η_ab’s were around 0, which means that most blocks had a similar estimate to the MLE, while a few blocks with large η_ab borrowed information from shrinkage and increased the estimation accuracy.

Email network

The Email-Eu-core network (Eucore) is a directed network generated using email data from a large European institute, consisting of incoming and outgoing communications between members of the institute from 42 departments. Leskovec & Krevl (2014) organized the data and labeled which department each individual node belongs to, i.e., the “ground-truth” community memberships. The network has n = 1,005 nodes and 25,571 directed edges, which we converted to undirected ones by removing their orientations. Although the test of Yan et al. (2012) suggested rejection of the hypothesis that a regular SBM is the true underlying model, our results on this network still show the improvement brought by EB assuming a regular SBM. We leave the generalization of our EB estimate to DCSBM as a future direction.

We applied VBEM to detect communities with an input number of clusters K = 1, 2, …, 50. The MSE ratios of EB over the other two competing methods were calculated and plotted against K in Fig. 9A. Similarly, EB achieved smaller MSE than the other two methods for all values of K. The MSE ratios ranged from 60% to 90%. When the input number of communities K was close to or greater than K* = 42, the improvement of EB over the competing methods became more substantial. Figure 9B shows higher test data likelihood of EB than the other two estimates. Figure 9C shows the values of three model selection criteria for k ∈ {1, …, 50}. The three methods, CVRP, VBEM and EB, gave estimates $\hat{K}=1$, 39 and 39, respectively. The $\hat{K}$ by VBEM and EB were both reasonably close to the ground-truth of K* = 42, while the performance of CVRP was much worse on this dataset. Moreover, ${\mathcal{J}}_{VBEM}$ is relatively flat around the estimated $\hat{K}$, while the curve of ${\mathcal{J}}_{VBEM}$ shows a higher sensitivity. From the distributions of the shrinkage values η_ab in (d) and (e), we see η ≥ 0.3 for a good number of diagonal and off-diagonal blocks, which led to substantial shrinkage and better performance than the MLE.