Community Detection with Dependent Connectivity

TL;DR

This paper introduces a new community detection method that models dependent connectivities within communities, improving accuracy over traditional models by capturing edge correlations without requiring likelihood specification.

Contribution

It proposes a novel approach using Bahadur representation to incorporate within-community edge dependence, addressing limitations of existing models like SBM.

Findings

The method reduces estimation bias.

It accelerates convergence of algorithms.

It outperforms variational EM in simulations.

Abstract

In network analysis, within-community members are more likely to be connected than between-community members, which is reflected in that the edges within a community are intercorrelated. However, existing probabilistic models for community detection such as the stochastic block model (SBM) are not designed to capture the dependence among edges. In this paper, we propose a new community detection approach to incorporate within-community dependence of connectivities through the Bahadur representation. The proposed method does not require specifying the likelihood function, which could be intractable for correlated binary connectivities. In addition, the proposed method allows for heterogeneity among edges between different communities. In theory, we show that incorporating correlation information can lower estimation bias and accelerate algorithm convergence. Our simulation studies show that the proposed algorithm outperforms the popular variational EM algorithm assuming conditional independence among edges. We also demonstrate the application of the proposed method to agricultural product trading networks from different countries.

Figures

Tables

Algorithm 1

Step 1: Input an initial membership probability for each node: ${\alpha }_{iq}^{\left(0\right)},\mathrm{\hspace{0.33em}1}\le i\le N,\mathrm{\hspace{0.33em}1}\le q\le K$

through spectral clustering on individual sample networks

Step 2: At the $s$th iteration, given ${\{{\beta }_{ql}^{\left(s-1\right)},{\rho }_q^{\left(s-1\right)}\}}_{q,l=1}^K$ and ${\left\{{\alpha }_i^{\left(s-1\right)}\right\}}_{i=1}^N$ from the $\left(s-1\right)$th iteration:

(i) Maximization: block-wise update ${\beta }_{ql}^{\left(s\right)}$ and ${\rho }_q^{\left(s-1\right)},q,l=1,\cdots ,K$;

(a) Obtain ${\beta }_{ql}^{\left(s\right)}$ through GEE with current membership as working correlation

(b) Update ${\rho }_q^{\left(s\right)}$ through the method of moments estimator given current membership

(ii) Expectation: given ${\{{\beta }_{ql}^{\left(s\right)},{\rho }_q^{\left(s\right)}\}}_{q,l=1}^K$ , update ${\left\{{\alpha }_i^{\left(s\right)}\right\}}_{i=1}^N$:

${\alpha }_{iq}^{\left(s\right)}=\frac{{\alpha }_{iq}^{\left(s-1\right)}\tilde{L}\left(𝐘|{\alpha }_{-i}^{\left(s-1\right)},Z_{iq}=1\right)}{{\sum }_{k=1}^K{\alpha }_{ik}^{\left(s-1\right)}\tilde{L}\left(𝐘|{\alpha }_{-i}^{\left(s-1\right)},Z_{ik}=1\right)},i=1,\cdots ,N,q=1,\cdots ,K.$

Step 3: Iterate until .

Step 4: Obtain the membership $z_i$ of clusters by

${\left\{{\alpha }_i^{\left(s\right)}\right\}}_{i=1}^N$: $z_i={\max }_k\{{\alpha }_{i1}^{\left(s\right)},\cdots ,{\alpha }_{iK}^{\left(s\right)}\},i=1,\cdots ,N$.

Table 2. Table 1: Adjusted Rand Index between estimated membership and true membership for networks with two communities and weak marginal signal averaging on 50 replicates.

		Unbalanced community			Balanced community
		$M=20$	$M=40$	$M=60$	$M=20$	$M=40$	$M=60$
$\rho =0$	VEM	0.38	0.41	0.48	0.31	0.28	0.28
	Bahadur2nd	0.36	0.41	0.47	0.32	0.29	0.29
	Bahadur4th	0.35	0.37	0.47	0.30	0.29	0.30
$\rho =0.3$	VEM	0.34	0.34	0.34	0.28	0.28	0.28
	Bahadur2nd	0.94	0.98	0.99	0.96	0.99	1.00
	Bahadur4th	0.96	0.99	1.00	0.99	0.99	1.00
$\rho =0.6$	VEM	0.34	0.34	0.34	0.29	0.28	0.28
	Bahadur2nd	0.96	0.99	0.99	0.97	1.00	1.00
	Bahadur4th	0.99	1.00	1.00	0.99	1.00	1.00

Table 3. Table 2: Adjusted Rand Index between estimated membership and true membership for networks with two communities and strong marginal signal averaging on 50 replicates.

		Unbalanced community			Balanced community
		$M=20$	$M=40$	$M=60$	$M=20$	$M=40$	$M=60$
$\rho =0$	VEM	0.78	0.92	0.98	0.76	0.90	0.97
	Bahadur2nd	0.73	0.91	0.97	0.77	0.92	0.98
	Bahadur4th	0.69	0.86	0.95	0.72	0.92	0.98
$\rho =0.3$	VEM	0.78	0.81	0.83	0.68	0.79	0.84
	Bahadur2nd	0.99	0.99	1.00	0.98	1.00	1.00
	Bahadur4th	0.99	0.99	1.00	0.99	1.00	1.00
$\rho =0.6$	VEM	0.78	0.89	0.83	0.84	0.92	0.88
	Bahadur2nd	0.99	1.00	1.00	0.99	1.00	1.00
	Bahadur4th	0.99	1.00	1.00	0.99	1.00	1.00

Table 4. Table 3: Estimation of within-community parameter β 11 = 1 subscript 𝛽 11 1 \beta_{11}=1 for networks with two communities and weak marginal signal.

		Unbalanced community			Balanced community
		$M=20$	$M=40$	$M=60$	$M=20$	$M=40$	$M=60$
$\rho =0$	VEM	${0.56}_{0.42}$	${0.59}_{0.29}$	${0.58}_{0.20}$	${0.64}_{0.32}$	${0.57}_{0.16}$	${0.64}_{0.18}$
	Bahadur2nd	${0.57}_{0.42}$	${0.58}_{0.30}$	${0.57}_{0.21}$	${0.61}_{0.28}$	${0.57}_{0.16}$	${0.66}_{0.20}$
	Bahadur4th	${0.52}_{0.42}$	${0.55}_{0.28}$	${0.57}_{0.19}$	${0.58}_{0.27}$	${0.58}_{0.18}$	${0.65}_{0.19}$
$\rho =0.3$	VEM	${0.49}_{0.30}$	${0.50}_{0.17}$	${0.52}_{0.14}$	${0.58}_{0.24}$	${0.58}_{0.18}$	${0.59}_{0.12}$
	Bahadur2nd	${0.81}_{0.48}$	${0.84}_{0.32}$	${0.89}_{0.27}$	${0.95}_{0.24}$	${0.93}_{0.16}$	${0.92}_{0.14}$
	Bahadur4th	${0.85}_{0.47}$	${0.83}_{0.31}$	${0.89}_{0.27}$	${0.96}_{0.24}$	${0.93}_{0.16}$	${0.93}_{0.14}$
$\rho =0.6$	VEM	${0.56}_{0.22}$	${0.54}_{0.20}$	${0.52}_{0.15}$	${0.61}_{0.27}$	${0.61}_{0.16}$	${0.60}_{0.14}$
	Bahadur2nd	${1.01}_{0.42}$	${1.04}_{0.35}$	${1.00}_{0.29}$	${0.95}_{0.31}$	${1.00}_{0.19}$	${0.96}_{0.15}$
	Bahadur4th	${0.99}_{0.25}$	${1.05}_{0.15}$	${1.01}_{0.13}$	${0.97}_{0.31}$	${1.01}_{0.19}$	${0.97}_{0.16}$

Table 5. Table 4: Estimation of within-community parameter β 22 = 1.5 subscript 𝛽 22 1.5 \beta_{22}=1.5 for networks with two communities and weak marginal signal.

		Unbalanced community			Balanced community
		$M=20$	$M=40$	$M=60$	$M=20$	$M=40$	$M=60$
$\rho =0$	VEM	${1.43}_{0.43}$	${1.42}_{0.34}$	${1.45}_{0.26}$	${1.18}_{0.40}$	${0.94}_{0.16}$	${0.94}_{0.15}$
	Bahadur2nd	${1.50}_{0.39}$	${1.49}_{0.31}$	${1.45}_{0.25}$	${1.21}_{0.42}$	${0.93}_{0.21}$	${0.97}_{0.22}$
	Bahadur4th	${1.56}_{0.37}$	${1.49}_{0.30}$	${1.46}_{0.23}$	${1.19}_{0.47}$	${0.94}_{0.24}$	${0.96}_{0.22}$
$\rho =0.3$	VEM	${1.31}_{0.23}$	${1.40}_{0.11}$	${1.37}_{0.11}$	${1.05}_{0.21}$	${0.92}_{0.16}$	${0.92}_{0.16}$
	Bahadur2nd	${1.56}_{0.19}$	${1.50}_{0.10}$	${1.49}_{0.09}$	${1.48}_{0.22}$	${1.45}_{0.19}$	${1.44}_{0.14}$
	Bahadur4th	${1.55}_{0.19}$	${1.50}_{0.09}$	${1.49}_{0.09}$	${1.48}_{0.22}$	${1.45}_{0.19}$	${1.45}_{0.14}$
$\rho =0.6$	VEM	${1.46}_{0.16}$	${1.43}_{0.16}$	${1.38}_{0.13}$	${1.16}_{0.21}$	${1.09}_{0.21}$	${1.06}_{0.22}$
	Bahadur2nd	${1.73}_{0.29}$	${1.60}_{0.15}$	${1.52}_{0.12}$	${1.73}_{0.28}$	${1.60}_{0.29}$	${1.64}_{0.15}$
	Bahadur4th	${1.69}_{0.25}$	${1.60}_{0.15}$	${1.52}_{0.13}$	${1.73}_{0.26}$	${1.61}_{0.29}$	${1.64}_{0.15}$

Table 6. Table 5: Estimation of within-community parameter β 12 = 0 subscript 𝛽 12 0 \beta_{12}=0 for networks with two communities and weak marginal signal.

		Unbalanced community			Balanced community
		$M=20$	$M=40$	$M=60$	$M=20$	$M=40$	$M=60$
$\rho =0$	VEM	${0.52}_{0.35}$	${0.57}_{0.24}$	${0.47}_{0.22}$	${0.22}_{0.31}$	${0.39}_{0.14}$	${0.41}_{0.11}$
	Bahadur2nd	${0.51}_{0.32}$	${0.58}_{0.23}$	${0.48}_{0.21}$	${0.23}_{0.30}$	${0.41}_{0.16}$	${0.39}_{0.15}$
	Bahadur4th	${0.51}_{0.29}$	${0.63}_{0.22}$	${0.48}_{0.20}$	${0.25}_{0.28}$	${0.40}_{0.17}$	${0.41}_{0.13}$
$\rho =0.3$	VEM	${0.68}_{0.24}$	${0.68}_{0.13}$	${0.69}_{0.10}$	${0.42}_{0.14}$	${0.35}_{0.12}$	${0.40}_{0.10}$
	Bahadur2nd	$-{0.02}_{0.25}$	${0.00}_{0.15}$	${0.00}_{0.11}$	${0.03}_{0.20}$	$-{0.05}_{0.16}$	$-{0.02}_{0.12}$
	Bahadur4th	$-{0.02}_{0.24}$	${0.00}_{0.14}$	${0.00}_{0.11}$	${0.03}_{0.18}$	$-{0.06}_{0.16}$	${0.03}_{0.12}$
$\rho =0.6$	VEM	${0.72}_{0.17}$	${0.71}_{0.11}$	${0.70}_{0.09}$	${0.41}_{0.18}$	${0.45}_{0.11}$	${0.48}_{0.11}$
	Bahadur2nd	$-{0.05}_{0.17}$	$-{0.03}_{0.13}$	${0.02}_{0.11}$	${0.00}_{0.19}$	${0.01}_{0.12}$	${0.03}_{0.12}$
	Bahadur4th	$-{0.04}_{0.17}$	$-{0.03}_{0.13}$	$-{0.02}_{0.11}$	$-{0.02}_{0.18}$	${0.00}_{0.12}$	${0.03}_{0.11}$

Table 7. Table 6: Performance comparison given misspecified inter-community correlation with balanced community and weak marginal signal averaging on 50 replicates.

		$\sigma =0.5$			$\sigma =1.5$
		$M=20$	$M=40$	$M=60$	$M=20$	$M=40$	$M=60$
$\rho =0.3$	VEM	0.28	0.28	0.29	0.28	0.28	0.29
	Bahadur2nd	0.90	0.99	1.00	0.76	0.99	0.99
	Bahadur4th	0.96	1.00	1.00	0.87	0.98	1.00
$\rho =0.6$	VEM	0.28	0.28	0.29	0.28	0.28	0.29
	Bahadur2nd	0.94	0.99	1.00	0.87	0.99	1.00
	Bahadur4th	0.99	1.00	1.00	0.94	0.99	1.00

Table 8. Table 7: Performance comparison given misspecified inter-community correlation with unbalanced community and weak marginal signal averaging on 50 replicates.

		$\sigma =0.5$			$\sigma =1.5$
		$M=20$	$M=40$	$M=60$	$M=20$	$M=40$	$M=60$
$\rho =0.3$	VEM	0.32	0.33	0.33	0.33	0.33	0.33
	Bahadur2nd	0.89	0.98	0.99	0.89	0.95	0.97
	Bahadur4th	0.95	0.99	0.99	0.93	0.94	0.94
$\rho =0.6$	VEM	0.34	0.33	0.34	0.33	0.33	0.33
	Bahadur2nd	0.91	0.96	0.98	0.91	0.95	0.94
	Bahadur4th	0.95	0.96	0.97	0.92	0.93	0.92

Table 9. Table 8: Clustering of nations in the agricultural products trading networks given 4 communities

	VEM	Bahadur4th
Community 1	Brazil, Denmark, Finland, Ireland	Austria, Denmark, Finland, Ireland, Poland
	Lebanon, Russia, Sweden, Switzerland	Russia, Sweden, Switzerland, Turkey
	Turkey, Ukraine, Argentina, Israel	Bulgaria, Croatia, Czech, Greece, Hungary
	Mexico, Norway, Portugal, Chile	Israel, Lithuania, Norway, Portugal
	South Africa, Qatar	Romania, Slovakia, Slovenia, Ukraine
Community 2	Australia, Canada, Hong Kong, Mainland	Brazil, Hong Kong, Taiwan, Indonesia
	Taiwan, India, Indonesia, Malaysia	Lebanon, Philippines, Korea, Argentina
	Japan, Philippines, Korea, Singapore	Mexico, Chile, New Zealand
	Thailand, U.S.A, New Zealand	South Africa, Qatar
Community 3	Belgium, France, Germany, Italy	Belgium, France, Germany, Italy
	Netherlands, Spain, United Kingdom	Netherlands, Spain, United Kingdom
Community 4	Austria, Poland, Bulgaria, Croatia	Australia, Canada, Mainland, India
	Czech, Greece, Hungary, Lithuania	Japan, Malaysia, Singapore
	Romania, Slovakia, Slovenia	Thailand, U.S.A

Equations

$$z_i\sim Multinominal(1,{\alpha }_i),$$

$$f_{ql}(Y_{ij}^m):=P(Y_{ij}^m|z_i=q,z_j=l)\sim Bern({\mu }_{ql}),\text{for}i,j\in \{1,\cdots ,N\},q,l=1,\cdots ,K,$$

$$\displaystyle\mu_{ql}^{m}=exp(\beta_{ql}x^{m}_{ij})/\big{\{}1+exp(\beta_{ql}x^{m}_{ij})\big{\}},$$

$$logP(\mathbf{Y};\mathbf{Z})=\sum _{m=1}^M\sum _{q=1}^K\sum _{i=1}^N{Z}_{iq}log{\alpha }_q+\sum _{m=1}^M\sum _{q,l=1}^K\sum _{i<j}^N{Z}_{iq}{Z}_{jl}{f}_{ql}(Y_{ij}^m;{\beta }_{ql}).$$

$$\displaystyle E_{R(\bm{Z},\bm{\tau})}\big{\{}logP(\bm{Y;Z})\big{\}}=-$$

$$\sum _{m=1}^M\sum _{q,l=1}^K\sum _{i<j}^N{\tau }_{iq}{\tau }_{jl}{f}_{ql}(Y_{ij}^m).$$

$$\displaystyle logL_{ind}(\bm{Y|Z})=\displaystyle\frac{1}{M}\sum_{m=1}^{M}\sum_{q,l=1}^{K}\sum_{{i<j}}^{N}Z_{iq}Z_{jl}\Big{\{}y^{m}_{ij}log\;\mu_{ql}+(1-y^{m}_{ij})log\;(1-\mu_{ql})\Big{\}}.$$

$$\hat{\mathbf{Z}}=\underset{\mathbf{Z}}{\mathrm{argmax}}log{P}_{ind}(\mathbf{Y}|\mathbf{Z}).$$

$$\sum _{i<j;u<v}^N{Z}_{ik}{Z}_{jk}{Z}_{uk}{Z}_{vk}{\hat{y}}_{ij}^m{\hat{y}}_{uv}^m\ge 0,$$

$$P(Y_1=y_1,\cdots ,Y_T=y_T)$$

$$\displaystyle\sum_{1\leq j_{1}<j_{2}<j_{3}\leq T}\rho_{j_{1}j_{2}j_{3}}\hat{y}_{j_{1}}\hat{y}_{j_{2}}\hat{y}_{j_{3}}+\cdots+\rho_{12\cdots T}\hat{y}_{1}\hat{y}_{2}\cdot\cdot\hat{y}_{T}\big{]},$$

$${\mu }_j=E(Y_j),{\hat{y}}_j=\frac{y_j-E(y_j)}{\sqrt{E(y_j)(1-E(y_j))}},$$

$${\rho }_{j_1{j}_2}=E({\hat{y}}_{j_1}{\hat{y}}_{j_2}),{\rho }_{j_1{j}_2{j}_3}=E({\hat{y}}_{j_1}{\hat{y}}_{j_2}{\hat{y}}_{j_3}),\cdots ,{\rho }_{12\cdots T}=E({\hat{y}}_1{\hat{y}}_2\cdot \cdot {\hat{y}}_T).$$

$$\text{mbox}log\tilde{L}(\mathbf{Y}|\mathbf{Z})=$$

$$\displaystyle+\sum_{m=1}^{M}\mbox{log}\Big{\{}1+\sum_{k=1}^{K}\frac{\rho_{k}}{2}\!\max\{\!\sum_{\begin{subarray}{c}i<j;u<v\\ (i,j)\neq(u,v)\end{subarray}}^{N}Z_{ik}Z_{jk}Z_{uk}Z_{vk}\hat{y}_{ij}^{m}\hat{y}_{uv}^{m},0\}\Big{\}}\Bigg{\}},$$

$$\displaystyle\mbox{log}L_{cor}(\bm{Y|Z})=\displaystyle\frac{1}{M}\Bigg{\{}\sum_{m=1}^{M}\mbox{log}\Big{\{}1+\sum_{k=1}^{K}\frac{\rho_{k}}{2}\max\{\!\sum_{\begin{subarray}{c}i<j;u<v\\ (i,j)\neq(u,v)\end{subarray}}^{N}Z_{ik}Z_{jk}Z_{uk}Z_{vk}\hat{y}_{ij}^{m}\hat{y}_{uv}^{m},0\}\Big{\}}\Bigg{\}}.$$

$$\frac{\tilde{L}(\mathbf{Y}|\bm{Z_{-i}},Z_{iq}=1)}{\tilde{L}(\mathbf{Y}|\bm{Z_{-i}},Z_{ik}=1)}=exp\Big{\{}\mbox{log}\tilde{L}(\mathbf{Y|\bm{Z_{-i}}},Z_{iq}=1)-\mbox{log}\tilde{L}(\mathbf{Y|\bm{Z_{-i}}},Z_{ik}=1)\Big{\}}.$$

$$\frac{\tilde{L}(\mathbf{Y}|{\mathbf{Z}}_{-\mathbf{i}},Z_{iq}=1)}{\tilde{L}(\mathbf{Y}|{\mathbf{Z}}_{-\mathbf{i}},Z_{ik}=1)}=\frac{L_{ind}(\mathbf{Y}|{\mathbf{Z}}_{-\mathbf{i}},Z_{iq}=1)}{L_{ind}(\mathbf{Y}|{\mathbf{Z}}_{-\mathbf{i}},Z_{ik}=1)}\frac{L_{cor}(\mathbf{Y}|{\mathbf{Z}}_{-\mathbf{i}},Z_{iq}=1)}{L_{cor}(\mathbf{Y}|{\mathbf{Z}}_{-\mathbf{i}},Z_{ik}=1)},$$

$$\displaystyle E\big{(}\hat{Y}^{m}_{i_{1}j_{1}}\hat{Y}^{m}_{i_{2}j_{2}}\hat{Y}^{m}_{i_{3}j_{3}}\hat{Y}^{m}_{i_{4}j_{4}}\big{)}\geq E\big{(}\hat{Y}^{m}_{i_{1}j_{1}}\hat{Y}^{m}_{i_{2}j_{2}})E\big{(}\hat{Y}^{m}_{i_{3}j_{3}}\hat{Y}^{m}_{i_{4}j_{4}})=\rho_{i_{1}j_{1}i_{2}j_{2}}\rho_{i_{3}j_{3}i_{4}j_{4}}.$$

$$\frac{{\rho }_k}{2}\sum _{\begin{array}{c}i<j,u<v\\ (i,j)\ne (u,v)\end{array}}^N{Z}_{ik}{Z}_{jk}{Z}_{uk}{Z}_{vk}{\hat{y}}_{ij}^m{\hat{y}}_{uv}^m={\rho }_k\sum _{s<t}^{N_0}{\gamma }_s^m{\gamma }_t^m.$$

$$\displaystyle\sum_{\begin{subarray}{c}s_{1}<t_{1},s_{2}<t_{2}\\ (s_{1},t_{1})\neq(s_{2},t_{2})\end{subarray}}^{N_{0}}\!\!\!\!\!\!\!\!\!\!\!\frac{E(\gamma_{s_{1}}^{m}\gamma_{t_{1}}^{m}\gamma_{s_{2}}^{m}\gamma_{t_{2}}^{m})}{2}\gamma_{s_{1}}^{m}\gamma_{t_{1}}^{m}\gamma_{s_{2}}^{m}\gamma_{t_{2}}^{m}\geq\!\!\!\!\!\!\!\!\!\!\!\sum_{\begin{subarray}{c}s_{1}<s_{2},t_{1}<t_{2}\\ (s_{1},t_{1})\neq(s_{2},t_{2})\end{subarray}}^{N_{0}}\!\!\!\!\!\!\frac{\rho_{k}^{2}}{2}\gamma_{s_{1}}^{m}\gamma_{t_{1}}^{m}\gamma_{s_{2}}^{m}\gamma_{t_{2}}^{m}\!=\!\!\Big{(}\rho_{k}\sum_{s<t}^{N_{0}}\gamma_{s}^{m}\gamma_{t}^{m}\Big{)}^{2}\!\!\!\!-\!\!\rho_{k}^{2}\sum_{s<t}^{N_{0}}(\gamma_{s}^{m}\gamma_{t}^{m})^{2}.$$

$$\lambda ={\lambda }_{ij}^m:=\frac{\#|\{(u,v):cor(Y_{ij}^m,Y_{uv}^m)>0,Z_u=Z_v=k\}|}{N_k(N_k-1)/2-1}\text{for edge Yijm in community k}$$

$$\displaystyle P_{Z*}\Big{\{}\frac{L_{ind}(\bm{Y}|\bm{Z}=z;\Theta)}{L_{ind}(\bm{Y}|\bm{Z}=z^{*};\Theta)}>t\Big{\}}=\mathcal{O}(\exp\Big{\{}-C_{1}\frac{rNM}{1+\rho\kappa_{2}N\min(r,\kappa_{2}\lambda N)}\Big{\}}),$$

$$\displaystyle P_{Z*}\Big{\{}\frac{\tilde{L}(\bm{Y}|\bm{Z}=z;\Theta)}{\tilde{L}(\bm{Y}|\bm{Z}=z^{*};\Theta)}>t\Big{\}}=\mathcal{O}(\exp\Big{\{}-C\frac{r\lambda NM(1+\lambda N^{2})}{1+c\rho\kappa_{2}N\min(r,\kappa_{2}\lambda N)}\Big{\}}),$$

$$E\Vert {\mathbf{\alpha }}^{(\mathbf{s}+1)}-z^{\ast }{\Vert }_1\le c_{1}NK\exp (-c_2(1+\lambda )MN)+\frac{c_3{N}^{1+\frac{\eta }{4}}\Vert {\mathbf{\alpha }}^s-{\mathbf{z}}^{\ast }{\Vert }_1}{(1+\lambda N^{2+\frac{\eta }{4}})M},$$

$$x_{ij}^m\sim Unif(a_1,a_2)\text{if}{Z}_{iq}=Z_{jq}=1,$$

$$x_{ij}^m\sim Unif(b_1,b_2)\text{if}{Z}_{iq}\ne Z_{jq},q=1,\cdots ,K.$$

$$corr(Y_{i_1{j}_1}^m,Y_{i_2{j}_2}^m)=\tilde{\rho },\text{given}{z}_{i_1}=z_{j_1}=q,z_{i_2}=z_{j_2}=l,q\ne l,$$

$${\mu }_{qq}^m=\frac{exp({\beta }_{ql}{x}_{ij}^m)}{1+exp({\beta }_{ql}{x}_{ij}^m)}+{\gamma }_k,M\frac{k-1}{10}\le m\le M\frac{k}{10},$$

$$X_t^{+}=\{f_1(Y_{12}^t),f_1(Y_{13}^t),\cdots ,f_1(Y_{N-1,N}^t)\}$$

$$L=inf\{\alpha \ge 0:E(\exp (⟨z,X_t^{+}-E(X_t^{+})⟩))\le \exp \{{\alpha }^2⟨{\Sigma }_{1}z,z⟩\}/2,z\in R^{N(N-1)/2}\}\le C$$

$$\bm{P}\Big{(}\frac{1}{M}\sum_{t=1}^{M}\sum_{\begin{subarray}{c}i<j;k<g\\ (i,j)\neq(k,g)\end{subarray}}^{N}\alpha_{iq}\alpha_{jq}\alpha_{kq}\alpha_{gq}\hat{Y}^{t,q}_{ij}\hat{Y}^{t,q}_{kg}>0\Big{)}=1$$

Full text

60K35,

60K35,

60K35,

Bahadur Representation,

high-order approximation,

product trading network,

stochastic block model,

variational EM,

arXiv:0000.0000 \startlocaldefs

\endlocaldefs

T1This work is supported by NSF Grants DMS 1613190 and DMS 1821198.

and

1 Introduction

Network data has arisen as one of the most common forms of information collection. This is due to the fact that the scope of study not only focuses on subjects alone, but also on the relationships among subjects. Networks consist of two components: (1) nodes or vertices corresponding to basic units of a system, and (2) edges representing connections between nodes. These two main components can have various interpretations under different contexts of application. For example, nodes might be humans in social networks; molecules, genes, or neurons in biology networks, or web pages in information networks. Edges could be friendships, alliances, URLs, or citations. The combination of the nodes and the edges defines a network, which can be represented by an adjacency matrix to reflect direct connectivities among nodes.

In this paper, we are interested in identifying community structures, such as community detection of cluster nodes which have more concentrated connectivities in a subnetwork. Identifying communities is essential to provide deep understanding of relationships among nodes within a community and between communities to address scientific, social and political problems [64, 8, 22, 63, 48, 37, 41]. In terms of other applications, community detection plays an important role in decomposing original large-scale network structures [66, 61, 52] into several subnetworks with more simplified structures [12], and facilitates scalable computation for further analyses.

The major community detection methods can be summarized in the following three categories. One approach is to search a partition of nodes which optimizes a global criterion over all possible partitions. The corresponding criterion function measures the goodness of fit of a partition such as the modularity [45] or profile likelihood [9] on the observed networks [59, 38, 2] to capture densely connected communities. However, obtaining a global optimum based on this type of criterion is computationally infeasible. In addition, modularity also suffers from the resolution limit [21], which intrinsically ignores small communities. The second approach is the spectral method [56, 20, 4], which recovers dense connectivities through the eigenvectors of the adjacent matrix of the network. One critical drawback of the spectral method is that it lacks robustness in estimation, especially when networks are sparse or consist of high-degree nodes such as hubs [34]. The third approach is the maximum likelihood method for cluster networks. This includes the popular stochastic block model (SBM) [28] and its extensions to incorporate the heterogeneity of nodes’ degrees [31, 69], and latent distance modeling [25, 26] to handle overlapping communities [1, 7].

The SBM assumes that membership assignments for each node follow a multinomial distribution. Given the community memberships, edges in the same community are randomly generated from a specified distribution. The common key assumption for SBM algorithms is that connectivities are conditional independent given the membership of nodes. This assumption simplifies the complexity of the model, and the likelihood function can be explicitly formulated. However, the network data are likely dependent among connectivities, which are also considered in several random network modelings [27, 36, 33, 15]. For community detection, the conditional independency assumption typically does not hold in practice and therefore could lead to a misspecified model [53, 3, 62]. For example, friendships within a social community or functional connectivities in brain networks tend to be highly correlated.

In addition, under conditional independence, the community structure can only be identified based on the marginal mean discrepancy of connectivities between within-communities and across-communities. Specifically, as a fundamental assumption of the independent SBM, the marginal mean discrepancy is required to be greater than a sharp threshold to guarantee community detectability ([40, 42]). However, the marginal mean discrepancy assumption might not hold, while the correlations among edges could be non-negligible and highly informative in identifying community structures. We show that the proposed method is able to incorporate the correlation information to achieve consistent community detection when the marginal mean discrepancy is insignificant.

More recently, the SBM has been extended to address the heterogeneity feature of within-community for multiple network samples. For example, [60, 50] apply a fixed-effect model through an independent intercept without incorporating information from other networks. Alternatively, a random-effects model is proposed to incorporate heterogeneity [51, 67], which borrows information from multiple networks. However, both of these approaches require the specification of a distribution for the random effects. In addition, an EM-type algorithm is implemented to integrate out the random-effects, [51, 67] which could be computationally expensive when the size of the community or the network size is large.

In this paper, we propose a novel community detection method to jointly model community structures among multiple networks. The proposed method can simultaneously incorporate the marginal and correlation information to differentiate within-community and between-community connectivities. The key idea is to approximate the joint distribution of correlated within-community connectivities by using a truncated Bahadur representation [5]. Although the approximate likelihood function is not the true likelihood, it is able to maximize the true community memberships and serves as a tighter lower bound to the true likelihood compared with the independent SBM likelihood. Consequently, we identify communities via maximizing the approximate likelihood function, which also serves as a discriminative function for membership assignments of nodes. In particular, within-community correlations provide an additional community-concordance measurement to capture high-order discrepancy between within-community and across-community networks, and therefore increase discriminative power to identify communities.

The main advantages and contributions of the proposed method can be summarized as follows. The proposed method incorporates correlation information among connectivities to achieve more accurate community detection than the variational EM method using marginal information only. The improvement of the proposed community detection method is especially powerful when the marginal information is relatively weak in practice. In addition, compared to the existing random-effects model, the proposed method is more flexible in modeling the heterogeneity of communities for multiple networks and heterogeneity of correlations among edges. Furthermore, it does not require a distribution specification among within-community connectivities.

In addition, we establish the consistency of the community estimation for the proposed approximate likelihood under a general within-community edge correlation structure and show that the proposed method achieves a faster convergence rate of membership estimation compared to the independent likelihood. In terms of computational convergence, the proposed algorithm achieves a lower estimation bias and a faster convergence rate compared to the variational EM algorithm at each iteration via incorporating additional correlation information. The theoretical development in this paper is nontrivial, since establishing membership estimation consistency is more challenging under the framework of conditional dependency among edges compared to the existing ones assuming the conditional independent model. Furthermore, we show that the convergence of the variational EM algorithm [39] is a special case of our method under the conditional independent SBM.

Computationally, we develop a two-step iterative algorithm which is not sensitive to initial values as in the standard variational EM algorithm. In addition, compared to the existing fixed-effects SBM with independent intercepts or the random-effects SBM, the proposed method has lower computational complexity, as it does not involve integration of random effects as in [51], or estimating the fixed effects for each network as in [50]. Simulation studies and a real data application also confirm that the proposed method outperforms the existing variational EM significantly, especially when the marginal information of observed networks is weak.

This paper is organized as follows: Section 2 introduces the background of the proposed method. Section 3 introduces the proposed method to incorporate correlation information for community detection. Section 4 provides an algorithm and implementation strategies. Section 5 illustrates the theoretical properties of the proposed method. Section 6 demonstrates simulation studies, and Section 7 illustrates an application to world agricultural products trading data. The last section provides conclusions and some further discussion.

2 Background and Notation

In this section, we provide background and notation of the proposed community detection. The stochastic block model (SBM) [28] is a form of hierarchical modeling which captures the community structure for networks. Consider $M$ symmetric and unweighted sample networks $\bm{Y}=\{\bm{Y}^{m}\}_{m=1}^{M}=\{(Y_{ij}^{m})_{N\times N}\}_{m=1}^{M}$ with $N$ nodes for $K$ communities. Let $\{z_{i}\}_{i=1}^{N}$ be the membership for each node and $z_{i}\in\{1,2,\cdots,K\}$, and denote the membership assignment matrix $\bm{Z}=\{(Z_{iq})_{n\times K}\}\in\{0,1\}^{N\times K}$, where $Z_{iq}=\mathbbm{1}\{z_{i}=q\}$. Here $\bm{Z}$ has exactly one 1 in each row and at least one 1 in each column for no-null communities. The unknown membership $z_{i}\in\{1,2,\cdots,K\}$ can be modeled as a latent variable from a multinomial distribution:

[TABLE]

where $i=1,\cdots,N$, $\alpha_{i}=\{\alpha_{i1},\cdots,\alpha_{iK}\}$ and $\sum_{k=1}^{K}\alpha_{ik}=1$. Given the membership of nodes, the observed edges between two nodes $\{(Y_{ij}^{m})_{n\times n}\}_{m=1}^{M}$ typically follow a Bernoulli distribution:

[TABLE]

where $\mu_{ql}$ is the probability of nodes $i$ and $j$ being connected.

For the heterogeneous stochastic blocks model, the marginal mean $\mu_{ql}$ for each block in the $m$th network can be modeled as a logistic model to incorporate heterogeneity among edges:

[TABLE]

where $\{(x_{ij}^{m})_{N\times N}\}_{m=1}^{M}$ are edge-wise covariates, and edges within the same community preserve homogeneity by sharing a block-wise parameter $\beta_{ql}$. The joint likelihood function can be decomposed into a summation of edge-wise terms following the conditional independence assumption:

[TABLE]

The latent membership $\bm{Z}$ is estimated by $E(\bm{Z|Y})$ through the maximum likelihood estimator of model parameters $\Theta=\{\beta_{ql};q,l=1,\cdots,K;\;\alpha_{q};q=1,\cdots,K\}$ in (2.3). However, the classical EM algorithm is not applicable here, because the conditional distribution $\displaystyle P(\bm{Z|Y})=\frac{P(\bm{Y;Z})}{\sum_{\bm{Z}}P(\bm{Y;Z})}$ becomes intractable in the expectation step.

The variational EM algorithm [39, 29] is one of the most popular inference methods, and can be applied to approximate the likelihood $P(\bm{Z|Y})$ by a complete factorized distribution $R(\bm{Z},\bm{\tau})=\displaystyle\prod_{i=1}^{N}h(Z_{i};\tau_{i})$, where $h(\cdot)$ denotes a multinomial distribution, $\bm{\tau}\!=\!(\!\tau_{1},\!\cdots,\!\tau_{N}\!)$ and $\tau_{i}\!=\!(\!\tau_{i1},\!\cdots,\!\tau_{iK}\!)$ is a probability vector such that $\sum_{q=1}^{K}\tau_{iq}=1$. In the expectation step, the likelihood $logP(\bm{Y\!;\!Z})$ is averaged over $R(\bm{Z})$ such that for any $\bm{\tau}$, E_{R(\bm{Z},\bm{\tau})}\big{\{}logP(\bm{Y\!;\!Z})\big{\}}\leq E_{P(\bm{Z|Y})}\big{\{}logP(\bm{Y\!;\!Z})\big{\}} where,

[TABLE]

Instead of directly maximizing E_{P(\bm{Z|Y})}\big{\{}logP(\bm{Y;Z})\big{\}}, the variational EM approach alternatively maximizes its lower bound E_{R(\bm{Z},\bm{\tau})}\big{\{}logP(\bm{Y;Z})\big{\}} over model parameters $\Theta$ and variational parameters $\bm{\tau}$, and clusters nodes by $\bm{\tau}$ through $\hat{z}_{i}=\operatorname*{argmax}_{k}\{\hat{\tau}_{ik},k=1,\cdots,K\}$.

Throughout this paper, we consider the conditional version of SBM (CSBM) [9, 56, 16], where the true membership $\bm{Z^{*}}$ is fixed. The conditional stochastic block model framework assumes conditional independence among edges, i.e., $Y^{m}_{i_{1}j_{1}}$ and $Y^{m}_{i_{2}j_{2}}$ are independent given nodes’ membership $z_{i_{1}},z_{i_{2}},z_{j_{1}},z_{j_{2}}$, and the corresponding log-likelihood of observed sample networks is:

[TABLE]

The above log-likelihood can serve as a discriminant function in clustering membership $\bm{Z}$ in that if $logL_{ind}(\bm{Y|Z_{1}})>logL_{ind}(\bm{Y|Z_{2}})$ given two membership assignments $\bm{Z_{1}}$ and $\bm{Z_{2}}$, then $\bm{Z_{1}}$ is preferred over $\bm{Z_{2}}$, since the likelihood for the observed sample networks is higher. Naturally, $\bm{Z^{*}}$ can be estimated by

[TABLE]

The SBM in (2.4) allows one to differentiate within-community and between-community nodes via utilizing only the marginal information, in that the average connectivity rates within-communities are higher than those between-communities. However, the underlying conditional independence assumption among edges is too restrictive and practically infeasible. In most community detection problems it is common that edges within communities are more correlated. For example, social connections among friends are highly correlated in social networks. However, the dependency among edges is not captured by the traditional SBM, which could lead to significant information loss of the community structure.

3 Methodology

3.1 Community detection with dependent connectivity

In this paper, we incorporate within-community correlation to improve accuracy and efficiency in identifying communities, in addition to utilizing the edges’ marginal mean information, since within-community dependency contains additional information regarding the membership of nodes. This is especially effective when the marginal mean is not informative in differentiating between and within communities’ connectivity.

In this section, we propose an approximate likelihood function to capture the dependency among within-community edges. We assume that each observed sample network $Y^{m}_{n\times n}$ is generated from an underlying joint binary distribution $P(Y^{m})$ such that the correlation among within-community edges is nonnegative. Specifically, for the underlying distribution, the correlation among edges $Y^{m}_{i_{1}j_{1}},Y^{m}_{i_{2}j_{2}}$ within a community satisfies: $corr(Y^{m}_{i_{1}j_{1}},Y^{m}_{i_{2}j_{2}})=\rho_{q}(i_{1},i_{2},j_{1},j_{2})\in[0,1]$ given nodes $z_{i_{1}}$, $z_{i_{2}}$, $z_{j_{1}}$ and $z_{j_{2}}$ are in the same community $q$, where $1\leq i_{1}<j_{1}\leq N,1\leq i_{2}<j_{2}\leq N,(i_{1},j_{1})\neq(i_{2},j_{2})$ and $q=1,\cdots,K$. Note that correlations among each pair of edges could be different. The reason we consider nonnegative correlations among within-community edges is that this reflects the concordance among within-community edges in practice. Equivalently, the edges in community $k$ show concordance only when

[TABLE]

where $\hat{y}_{ij}^{m}$ and $\hat{y}_{uv}^{m}$ are normalized binary variables through the marginal mean. Note that positive correlation among edges has been considered for community detection on multiple networks. For example, [50, 51] utilize random effects to model the heterogeneity of the connectivity rate for an individual network, which infers a positive correlation among the edges within the same community. In practice, it has been demonstrated that assuming nonnegative correlations among edges is more sensible and interpretable. For example, the positive pairwise correlation among edges is more likely to produce star or triad relations which are widely observed in social networks [54, 55].

3.2 Approximate Likelihood

In this section, we propose an informative approximation of the true log-likelihood to cluster $\bm{Z}$ via incorporating interactions among edges within a community in addition to marginal mean information. This is because the exact joint likelihood function of correlated binary distribution $P(Y^{m})$ is computationally intractable. Specifically, we construct an approximate likeihood as a substitute of the true likelihood by facilitating the Bahadur representation [5]. That is, we retain the low-order dependency information among edges within-communities and discard the high-order dependency for computational efficiency. Although the approximate likelihood is not a true likelihood, it still serves the purpose of estimating the membership of nodes.

Consider $T$ dependent binary random variables, then the joint likelihood can be represented through the Bahadur representation:

[TABLE]

where

[TABLE]

and

[TABLE]

The idea of Bahadur representation is to approximate the joint distribution of dependent binary random variables as a function of moments with a sequential order. For the community detection problem, the binary random variables represent within-community edges, and the corresponding joint distribution can be explicitly decomposed into a marginal part and a correlation part. The marginal part consists of all the marginal mean $\mu_{ij}$ for each edge, which can be directly modeled through the dependency of the mean on covariates as in (2.2). The correlation part consists of interactions among all possible pairwise-associations of normalized edges, which add correlation information beyond a conditional independence likelihood model. Note that the conditional independence model is a special case of the proposed model when the correlation is zero, and the corresponding Bahadur representation collapses to a marginal part only, which is equivalent to the $\mbox{log}L_{ind}(\bm{Y|Z})$ in (2.4).

There are two major challenges in applying the Bahadur representation to model the interactions among within-community edges. First, the dimension of correlation parameters could be high if all the high-order interactions in (3.2) are incorporated, and this could lead to an increasing computational demand as the size of community grows. To solve this problem, we retain all the second-order interactions, but ignore interactions for higher orders beyond the second order, since the pairwise interactions among edges could be most important. In addition, we can further reduce the number of parameters via a homogeneous correlation structure such that all the pairwise correlations in each community are assumed to be the average within-community correlation, which can be simplified as an exchangeable correlation structure. The rationales of this simplification are based on the following. First, the pairwise correlation parameter $\rho_{q}(i_{1},i_{2},j_{1},j_{2})$ is a nuisance correlation parameter to enhance clustering. Second, both the numerical experiments and theoretical findings show that the density of pairwise correlation among within-community edges plays a more important role than the intensity of the correlation in affecting clustering performance.

The second challenge is that the range of the correlation coefficient could be constrained by the marginal means [19]. Consequently, the correlation parameter space is more restrictive if the variability of marginal means among edges is large. Nevertheless, our primary goal is to construct an objective function which can incorporate information from the marginal mean and correlations of edges within-community, and the objective function is not necessarily the true likelihood function. In the proposed method, we instead construct an approximate likelihood which is more flexible for incorporating highly dependent communities while still achieving computational efficiency.

Specifically, we construct an approximate likelihood $\tilde{L}(\bm{Y|Z})$ incorporating correlated within-community edges as follows:

[TABLE]

where $\mu_{ql}$ and $\hat{y}_{ij}^{m}$ are formulated in (2.2) and (3.2), and $\rho_{k}$ is the average of pairwise correlation in community $k$. Notice that the first term in (3.4) is the same as the marginal mean model, and the second term in (3.4) measures the concordance among edges within communities clustering $\bm{Z}$.

We denote the second term of (3.4) as

[TABLE]

Compared with $\mbox{log}L_{ind}(\bm{Y|Z})$ in (2.4), the proposed $\mbox{log}\tilde{L}(\bm{Y|Z})$ has more discriminative power over $\bm{Z}$, since it utilizes more information of the observed dependency within communities corresponding to clustering $\bm{Z}$. In addition, the nonnegativity of $\mbox{log}L_{cor}(\bm{Y|Z})$ ensure the fact that $\mbox{log}\tilde{L}(\mathbf{Y|Z})\geq\mbox{log}L_{ind}(\mathbf{Y|Z})$ is guaranteed, which implies that adding additional correlation information among edges can be more informative given within-community correlation exists. This leads to higher classification accuracy and estimation efficiency through maximizing (3.4).

The key part of the proposed method is to predict memberships of nodes through the Bayes factor constructed by the proposed $\mbox{log}\tilde{L}(\bm{Y|Z})$. Suppose the memberships of other nodes $\bm{Z_{-i}}$ are known, then we classify node $i$ based on the following Bayes factor:

[TABLE]

If the above Bayes factor $>1$, then the probability of node $i$ in community $q$ is larger than that of community $k$. The Bayes factor can be further decomposed as:

[TABLE]

which contains both the marginal ratio and the correlation ratio. It is clear that when the marginal information is weak in differentiating two communities, the marginal ratio is close to 1, and if the correlation ratio is informative, it can enhance the Bayes factor to improve community detection. In addition, the correlation ratio also serves as a correction to lower the estimation bias.

We illustrate the advantage of the proposed method in (3.4) over the conditional independent likelihood (2.4) using a simple numerical illustration. Specifically, we generate multiple networks based on the SBM with 30 nodes evenly split between two communities. The marginal means of within-community and between-community edges are the same at 0.5, implying that the marginal mean is not informative. We assume a true exchangeable correlation $\rho=0.6$ for within-community edges. Figure 1 illustrates that the likelihood function changes as memberships of nodes change with some misclassified nodes. The left graph is based on the conditional independent SBM utilizing only marginal information, which does not differentiate the two communities at all due to weak marginal information. However, the proposed approximate likelihood in the right graph has high differentiation power for the nodes’ memberships, and reaches maximum when the true memberships are selected.

4 Algorithm and Implementation

In this section, we propose a two-step algorithm to maximize the proposed approximate likelihood function. In addition, we provide implementation strategies to improve the stability and efficiency of the algorithm.

4.1 Algorithm

To estimate the true membership $\bm{Z^{*}}$ of nodes, we can ideally search through all the possible $\bm{Z}$ and choose the one with the largest $\mbox{log}\tilde{L}(\bm{Y|Z})$. However, this becomes infeasible when the number of nodes $N$ and the number of communities $K$ increases. In the following, we propose an iterative two-step algorithm to maximize $\mbox{log}\tilde{L}(\bm{Y|Z})$ in (3.2).

Here we directly maximize the approximate likelihood instead of a true likelihood as in the EM algorithm. In the expectation step, we alternatively update membership of each node while fixing other nodes, where $\tilde{L}(\mathbf{Y}|\alpha_{-i};Z_{ik})$ has the same formulation as $\tilde{L}(\bm{Y|Z})$ in (3.2) with $\{Z_{iq}\}_{N\times K}$ replaced by its expectation $\{\alpha_{iq}\}_{N\times K}$, except $Z_{ik}$. Note that $\alpha_{iq}$ is not the expectation under the true underlying joint distribution $P(Y,Z)=P(Y|Z)P(Z)$. Instead, it corresponds to the distribution defined by the approximate likelihood in (3.4). In the expectation step, the memberships are updated through the Bayes factor in (3.5) with the proposed $\tilde{L}(\mathbf{Y|Z})$. In the maximization step, we estimate the community-wise parameters $\beta_{ql}$ through the generalized estimating equation where the working correlation is exchangeable structure given the current membership of nodes and estimated average correlation $\rho_{q}$. Note that the variational EM is a special case of the proposed algorithm if the correlation information is ignored and the conditional independent model in (2.4) is assumed.

4.2 Computation and Implementation:

To ensure computational stability, the community-wise parameters $\beta_{ql}$ could be estimated through a simplified generalized estimation equation assuming an independent working correlation in algorithm 1. This is because the primary interest of community detection is classification accuracy, and the empirical studies show that correlation information plays a relatively minor role in parameter estimation.

We can achieve a better approximation to the true likelihood if higher-order moments are incorporated in the Bahadur representation in (3.2), which also increases its discrimination power. However, higher-order correlation could also increase the computational cost. Alternatively, we can recover partial higher-order interactions (e.g., the fourth order) derived from low order interactions (e.g., the second order). For example, consider four normalized edges $\hat{Y}^{m}_{i_{1}j_{1}},\hat{Y}^{m}_{i_{2}j_{2}},\hat{Y}^{m}_{i_{3}j_{3}}$ and $\hat{Y}^{m}_{i_{4}j_{4}}$ within the same community $k$ with a positive fourth order correlation among them, we have

[TABLE]

To simplify notation, denote $(Z_{1k}Z_{2k}\hat{Y}^{m}_{12},Z_{1k}Z_{3k}\hat{Y}^{m}_{13},\cdots,Z_{2k}Z_{3k}\hat{Y}^{m}_{23},\cdots,Z_{(N-1)k}Z_{Nk}\hat{Y}^{m}_{(N-1)N})$ as $(\gamma_{1}^{m},\gamma_{2}^{m},\cdots,\gamma_{N_{0}}^{m})$, where $N_{0}=\frac{N^{2}-N}{2}$. Then the second-order interaction term for the community $k$ in $L_{cor}(\bm{Y|Z})$ is

[TABLE]

Based on (4.2) and given $\bm{Z}$, we can approximate the fourth-order interaction for community $k$ under the exchangeable correlation structure by its lower bound:

[TABLE]

Note that the above lower bound of the fourth-order interaction can be calculated by the second-order interaction term in $L_{cor}(\bm{Y|Z})$. Therefore, we can still incorporate higher-order terms in $log\tilde{L}(\bm{Y|Z})$ without additional computational cost. For other types of non-exchangeable correlation structures, we can incorporate partial higher-order correlation similarly as above. The main difference is that each pair of edges is associated with a specific correlation given a dependency structure. Therefore, the simplified lower bound for higher-order correlations such as (5.1) does not hold in general, and could have a more complex form depending on the specific correlation structure.

In the following, we also provide some guidelines for determining the number of communities $K$ and initial membership of nodes. For a single network, the criterion-based methods choose $K$ to maximize a certain probabilistic criterion such as the integrated likelihood [23, 17, 35], composite likelihood BIC [58] or modularity criterion [10]. In addition, spectral methods estimate $K$ through the spectral property of the transformed adjacent matrix, such as a Laplacian matrix [44], non-backtracking matrix [11] or Bethe Hessian matrix [57]. In the hierarchical Bayesian framework, the number of communities is treated as a model parameter given a certain prior distribution and is jointly estimated with nodes’ memberships using the MCMC [23, 46, 47]. For multiple networks, we can extend the above techniques to estimate a consensus number of communities combining observed realizations of the SBM from each individual network.

In the context of the proposed within-community dependent modeling, we can first perform the modularity-maximizing method or spectral clustering on each individual network to obtain $K$, then take the average of these individual estimated $K$, which can be treated as a consensus number of communities. The above procedure is sensible under two considerations. First, each sample network is a realization of the SBM so that the individual estimation of $K$ is randomly distributed around the true underlying $K$. Thus the average of individual estimations provides an estimation of $K$ with low-bias and low-variance. Second, the spectral clustering or modularity methods are more favorable than other methods, due to their relatively low computational cost in estimating $K$. This is especially effective when the sample size of networks is large.

As an EM-type algorithm, the proposed optimization procedure can only guarantee the local maximum and requires multiple initializations to find the global maximum. In this paper, we obtain the membership initializations through spectral clustering on different sample networks, a benchmark algorithm for the traditional SBM. Spectral clustering is a model-free clustering algorithm and is able to provide a warm start for nodes’ memberships.

5 Theoretical Results

In this section, we establish the consistency of the estimated nodes’ membership based on the independent likelihood and the approximate likelihood approaches. In addition, we provide the computational convergence theorem for the proposed iterative algorithm in section 4. Compared to the independent likelihood approach, we show that the approximate likelihood approach leads to a computationally faster convergence rate regarding nodes’ membership estimation.

5.1 Consistency of nodes’ membership estimation

In this subsection, we study the consistency of the maximization likelihood estimator for both the independent likelihood and the approximate likelihood at the population level. With the independence assumption among within-community edges, the consistency and convergence rate of the MLE estimator can be obtained by [13, 68]. However, the convergence property of the MLE remains unknown if there exists a local dependence among edges.

One significant distinction using the independence assumption if the edges are correlated is that the increasing number of nodes and number of edges do not necessarily guarantee a lower misclassification rate and computationally faster convergence. This is because the discrepancy between marginal means from within-community and between-community is not accumulated due to the pairwise correlation, though it can be accumulated through increasing the number of sample networks. However, we show that the proposed approximated approach is able to benefit from the increasing number of nodes, and therefore achieves a faster computational convergence compared to the independent likelihood approach.

Without loss of generality, we assume that $x_{ij}^{m}=1,m=1,\cdots,M,i,j=1,\cdots,N$. That is, all the edges within the same block have the same marginal mean such that $\mu_{z_{i}z_{j}}:=E(Y_{ij}^{m}|i\in q,j\in l)=\frac{\exp(\beta_{ql})}{1+\exp(\beta_{ql})}$. We denote that the true marginal mean as $\Theta=\{\mu_{ql},1\leq q<l\leq K\}$. The following two regularity conditions regarding identifiability are standard:

(C1). Suppose for every $q\neq q^{\prime},1\leq q,q^{\prime}\leq K$ and $l\neq l^{\prime},1\leq l,l^{\prime}\leq K$, we have $\mu_{ql}\neq\mu_{q^{\prime}l^{\prime}}$. In addition, all the $\mu_{ql}$ are bounded such that $c_{2}\leq\mu_{ql}\leq c_{1},q,l=1,\cdots,K$ where $0<c_{2}\leq c_{1}<1$.

(C2). Community sizes from all sample networks are bounded above and below by $\kappa_{1}N\leq|\{i\in\{1,2,\cdots,N\}:Z^{*}_{iq}=1\}|\leq\kappa_{2}N,\;q=1,\cdots,K$, where $\kappa_{1}$ and $\kappa_{2}$ are constants such that $0<\kappa_{1}<\kappa_{2}<1$.

In the following, we establish the consistency of membership estimation for both the independent likelihood approach and the proposed approximate likelihood approach. For the within-community edges, we define the edgewise second-order pairwise correlation density as

[TABLE]

where $k=1,2,\cdots,K$ and $N_{k}(N_{k}-1)/2-1$ is the number of edges within community $k$ for the sample network $\bm{Y}^{m}$. For simplicity, we assume the homogeneous second-order correlation density such that $\lambda_{ij}^{m}=\lambda$ for all the within-community edges. Here $\lambda\in[0,1]$ determines the intensity of local dependency within a community. Specifically, $\lambda=0$ indicates that within-community edges are all independent, while $\lambda=1$ indicates that all the within-community edges are pair-wisely correlated. In addition, correlation density $\lambda$ is allowed to depend on the number of nodes, and increases such that it can model a more general class of correlation structure. For example, in a hub structure, an edge is only correlated with those sharing the same hub nodes and the density $\lambda=\frac{N_{k}-1}{N_{k}^{2}-1}=O_{N}(\frac{1}{N_{k}})$.

To establish asymptotic consistency for the proposed likelihood, we assume the sparsity of high-order correlation among within-community edges.

(C3). The number of third and fourth-order correlations defined in (3.2) among within-community edges do not exceed the order of the size of second-order correlations. Specifically, for edge $Y^{m}_{ij}$ in community $k$, $\#|\{(i,j),(u_{1},v_{1}),(u_{2},v_{2}):E(\hat{Y}_{ij}\hat{Y}_{u_{1}v_{1}}\hat{Y}_{u_{2}v_{2}})>0\}|\leq O_{N}(\lambda(N_{k}^{2}))$. In addition, $\#|\{(i,j),(u_{1},v_{1}),(u_{2},v_{2}),(u_{3},v_{3}):E(\hat{Y}_{ij}\hat{Y}_{u_{1}v_{1}}\hat{Y}_{u_{2}v_{2}}\hat{Y}_{u_{3}v_{3}})>0\}|\leq O_{N}(\lambda(N_{k}^{2})),k=1,2,\cdots,K$.

In general, assume that the pairwise correlations among the within-community edges are sufficient to cover a broad class of Markov dependence modeling under the general exponential random graph model. This includes the most commonly used edge dependence configurations such as a star, a triangular shape subnetwork [43] and the k-triangles shape [49]. Although considering that the additional higher-order edge correlation improves the model’s complexity, it could increase higher computational cost and instability. Empirically, it is sensible to assume that higher-order correlation only exists when second-order correlation already exists among edges, for the sake of identifiability and interpretability of the model. Otherwise, it could lead to the ’near degeneracy’ [24] when a higher-order dependency masks a lower-order dependency.

Let $P_{Z^{*}}:=\mathbb{P}(\cdot|Z=z^{*};\Theta)$ denote the conditional distribution of edges given the true membership of nodes and true parameters.

Theorem 5.1**.**

Under the regularity conditions (C1)-(C3), we establish the convergence rate of the membership estimator $z$ using the independent likelihood approach. That is, for every $t>0$ and $z\neq z^{*}$,

[TABLE]

where $r=\|z-z^{*}\|_{0}$ is the number of misclassified number of nodes up to the permutation of the label, and $\rho$ is the largest pairwise correlation among within-community edges. In addition, $C_{1}=\frac{c}{\max\{\log\frac{c_{1}}{c_{2}},\log\frac{1-c_{2}}{1-c_{1}}\}}$, where $c$ is a positive constant.

For the independent likelihood approach, the convergence rate depends on the number of sample network $M$ and the density of the pairwise correlation among within-community edges $\lambda$. The probability of true membership goes to 1 as $M$ or the node size $N$ increases given a relatively sparse pairwise correlation such that $\lambda N=o_{N}(1)$. If there is no pairwise correlation among edges, hence $\lambda=0$, then the convergence rate increases to $O_{N,M}(\exp(-C_{1}rNM))$, which degenerates in the conditional independent setting to the convergence rate established in [13].

In the case of the exchangeable correlation structure for within-community edges, hence $\lambda=1$, the convergence rate decreases to the order of $O_{N,M}(\exp(-C_{1}\frac{M}{\kappa_{2}\rho}))$, and therefore does not benefit from the increasing number of nodes. In this case, the consistency relies on accumulating independent sample networks. Theorem 5.1 also implies that the independent likelihood approach is unable to fully accumulate discriminative power from the increasing number of nodes when there exists dependency among within-community edges. Indeed, the convergence rate of the independent likelihood approach decreases in terms of network size N as the within-community correlation density $\lambda$ increases. However, we show that the proposed approximate likelihood approach still benefits from increasing nodes size even under the exchangeable correlation structure among edges.

Theorem 5.2**.**

Under the regularity conditions (C1)-(C3), we establish the convergence rate of the estimator $z$ using the proposed approximate likelihood approach. That is, for every $t>0$, $z\neq z^{*}$, and $\lambda>0$,

[TABLE]

where $M>O_{N}(\frac{1}{\lambda})$, $r=\|z-z^{*}\|_{0}$ is the number of misclassified nodes up to the permutation of the label, $C,c$ are positive constants, and $\rho$ is the largest within-community correlation.

Given the same number of network $M$ and node size $N$, the proposed approximate likelihood approach is able to achieve a constantly faster convergence rate compared with (5.2) since it has an additional term $\lambda^{2}N^{3}M$ on the numerator in (5.3) given $\frac{1}{\lambda}=o_{N}(N)$. Specifically, the proposed approach is most superior under the exchangeable correlation structure ($\lambda=1$), where the convergence rate of the independent likelihood is at the order of $O_{N,M}(\exp(-cM))$, in contrast to the proposed convergence rate of $O_{N,M}(\exp(-cNM))$. Intuitively, incorporating the correlation information increases the effective sample size of within-community edges. Under the sparsity assumption of higher-order correlation among edges, the proposed approach benefits from accumulating information on the second-order interactions among edges, while the independent likelihood approach only accumulates information from the first-order marginal mean of edges.

5.2 Computational convergence for the proposed algorithm

In this subsection, we provide the computational convergence property of the proposed algorithm in Section 4. The main difference between the proposed method and the variational EM lies in the Bayes factor of (3.5) in the expectation step from Algorithm 1. If we replace $\tilde{L}(\bm{Y|Z})$ by the conditional independent likelihood $L_{ind}(\bm{Y|Z})$ in (2.4) in the expectation step, the standard variational EM becomes a special case of Algorithm 1. Notice that [68] establishes computational convergence with the minimax rate of misclassification only when the within-community edges are independent. In addition, it assumes that the within-community marginal means are all the same, which is too restrictive in practice.

In the following, we establish the computational convergence for the proposed approximate likelihood. Specifically, we are able to show a faster convergence speed and a lower estimation bias compared to the existing one based on the independent likelihood in [68]. The following Theorem 5.3 also relaxes the homogeneous marginal mean assumption and allows the marginal means from within-community and between-community to be different. We denote the estimated memberships of nodes at the $s$th iteration as $\bm{\alpha^{(s)}}=(\alpha_{1}^{(s)},\cdots,\alpha_{N}^{(s)})$ from Algorithm 1. In addition to the assumptions (C1-C3) in Section 5.1, we require two regularity conditions for the following theorems:

(C4). Suppose the distance between initial membership $\bm{\alpha}^{(0)}$ and true membership $z^{*}$ is bounded: $\|\bm{\alpha}^{(0)}-z^{*}\|_{1}\leq cN^{1-\eta}$ where $0<\eta<1$.

A common issue for most EM-type algorithms including the one proposed is that they only guarantee convergence to a local optimum. If the likelihood function is unimodal, then the EM-type algorithm converges to the MLE as the unique global optimum. However, the proposed approximate likelihood is non-convex and multi-modal. Therefore, we assume that the initials are in the neighborhood of the MLE to ensure the convergence of the EM algorithm[6, 65]. Condition C4 is a common assumption to guarantee computational convergence for EM-type algorithms [68, 30, 32].

(C5). The estimated marginal mean $\hat{\mu}_{ql}$ has a bounded bias from the truth, i.e., $0<\gamma_{1}\leq\frac{\hat{\mu}_{ql}}{\mu_{ql}}\leq\gamma_{2},\;q,l=1\cdots,K,$.

Theorem 5.3**.**

Under the regularity conditions (C1)-(C5) and given $N$ is sufficiently large, we establish the convergence property of Algorithm 1 via incorporating correlation information. That is, with the correlation density $\frac{1}{\lambda}=o_{N}(N^{\frac{\eta}{2}})$, as $M$ and $N$ increases with $O_{N}(\frac{1}{\lambda})<M\leq o_{N}(N^{2-\frac{\eta}{2}})$, then

[TABLE]

where $c_{1},c_{2},c_{3}$ are positive constants.

In Theorem 5.3, the first term on the right side of the inequality represents the irreducible estimation bias which measures the discrepancy between the community structure and its realization. The second term provides a decreasing rate of misclassification along each iteration. Theorem 5.3 indicates that the estimated memberships are closer to the true memberships compared to the previous iteration step at a rate of $\frac{N^{1+\frac{\eta}{4}}}{(1+\lambda N^{2+\frac{\eta}{4}})M}$, where a larger sample size $M$ or node size $N$ contribute a faster convergence and a lower estimation bias. In general, Theorem 5.3 guarantees the convergence of the iterative algorithm even without incorporating correlation information, but improves the convergence rate and estimation bias when correlation information is incorporated.

Specifically, in contrast to the computational convergence rates in Theorem $3.1$ of [68], our Theorem 5.3 shows that incorporating the correlation information enables us to reduce the estimation bias and accelerate the convergence rate. Specifically, if we consider the $M$ sample networks with node size $N$ as a single network with $MN$ nodes, then the proposed approximate likelihood approach reduces the order of the estimation bias from $O_{N,M}(MN\exp(-cMN))$ in [68] to $O_{N,M}(N\exp\{-c^{\prime}(1+\lambda)MN\})$ in (5.4), and the order of the convergence rate from $O_{N,M}(\frac{1}{\sqrt{MN}})$ in [68] to $O_{N,M}(\frac{N^{1+\frac{\eta}{4}}}{(1+\lambda N^{2+\frac{\eta}{4}})M})$ in (5.4), respectively. Compared with the convergence rate of the membership estimator assuming conditional independence in (5.4) when the correlation density $\lambda=0$, incorporating within-community correlation accelerates the computational convergence when there is a sufficiently large number of within-community correlated edges corresponding to $\lambda>\frac{1}{\sqrt{MN}}$.

6 Numerical Studies

In this section, we conduct simulation studies to illustrate the performance of the proposed method on community detection in networks for dependent edges within-community. In particular, we compare our method to the existing variational EM method which assumes conditional independence among edges.

6.1 Study 1: Networks with dependent within-community connectivity

In the first simulation study, we consider networks where edges within the same community are correlated and compare the performance of various methods under different network sample sizes with various magnitudes of marginal means for within-community and between-community.

Suppose the memberships of nodes $\bm{Z^{*}}=\{\bm{Z}_{1},\cdots,\bm{Z}_{n}\}$ in the networks are given with $K$ communities, where $\bm{Z}_{i}$ is a binary indicator vector corresponding to the membership of nodes $i$. Conditional on $\bm{Z^{*}}$, edges in each sample network are generated following the Bernoulli marginal distribution as in (2.1), where within-community edges follow an exchangeable correlation structure as in (3.2). Here we assume that between-community edges are independent from each other. The block-wise marginal means $\mu_{ql}\;(q,l=1,\cdots,K)$ are associated with edgewise covariates through (2.2). In addition, the edgewise covariates follow a uniform distribution, where within-communities covariates

[TABLE]

and between-community covariates

[TABLE]

Although the probability of each edge is different, the edges within the same community share the same coefficient $\beta_{ql}$ in (2.2). In the following simulation studies, we generate correlated unweighted edges through the R package ”MultiOrd.”

Specifically, the sample networks consist of 40 nodes split into two communities. In a balanced community network, each community has 20 nodes. In an unbalanced case, two communities are comprised of 10 and 30 nodes, respectively. We compare the performance under different sample sizes of networks with $M=20,\,40$ and 60, and different intensities of within-community dependency with correlation coefficient $\rho=0$, $0.3$ and $0.6$.

To simulate a weak marginal signal case, we let the block-wise parameters be $\beta_{11}=1,\,\beta_{22}=1.5$ and $\beta_{12}=\beta_{21}=0$. The means of within-community and between-community covariates are 0 with $a_{1}=b_{1}=-0.2$ and $a_{2}=b_{2}=0.2$ in (6.1) and (6.2). Here, although the marginal mean of within-community edges is slightly larger than that of between-community edges on average due to the convexity of the logistic link function in (2.2), the marginal means of within-community edges and between-community edges are very close.

For a strong marginal signal case, the block-wise parameters are $\beta_{11}=0.3,\,\beta_{22}=0.6$ and $\beta_{12}=\beta_{21}=0.2$. The within-community covariates are generated via (6.1) with $a_{1}=0.9$ and $a_{2}=1.1$, and between-community covariates are generated from (6.2) with $b_{1}=-0.8$ and $b_{2}=-0.6$. Note that there is a distinct gap between within-community and between-community marginal means, thus the marginal signal is more dominant for nodes within communities.

We use the Adjusted Rand Index (ARI) to measure the performance of clustering. The ARI takes a value between $-1$ and $1$, where $1$ represents a perfect matching of true memberships and predicted memberships of clustering, [math] indicates a random clustering and a negative value indicates that the agreement is less than the expectation from a random result. In the following simulations, we choose five fixed initial memberships of nodes in both balanced and unbalanced communities. These initials can be obtained from spectral clustering on sample networks. The Adjusted Rand Indices based on these chosen initials range between $0.30$ to $0.34$ under the unbalanced community case and between $0.25$ to $0.29$ under the balanced community case, which are far from the true memberships.

We compare the performance of clustering and parameter estimation for the proposed method applying the second-order (Bahadur2nd) and the fourth-order (Bahadur4th) Bahadur approximation, and the variational EM (VEM) approach with only marginal information.

In Table 1 and Table 2, the proposed method with the second-order and fourth-order approximations outperform the variational EM in clustering. Specifically, under the weak marginal signal case in Table 1, the Adjusted Rand Index of the variational EM are 0.34 under different network sizes and correlation strengths, which are similar to the ones calculated by fixed initials. In addition, since the distributions of marginal means from within-community and between-community are similar, the variational EM marginal approach barely improves over the initial memberships as it only utilizes the marginal information. However, the proposed method with the second-order or fourth-order Bahadur representation improves on the ARI by about $280\%$, compared to the VEM when $\rho=0.3$ and $\rho=0.6$. In addition, the performance of the proposed method improves by $1\sim 5\%$ as the number of sample networks increases from 20 to 60. Furthermore, incorporating the fourth-order interaction can slightly improve the accuracy of clustering.

We notice that when the correlation is as moderate as 0.3, the proposed method still achieves significant improvement over the variational EM and almost fully recovers the true memberships of clustering. We consider this as an intrinsic advantage of the proposed method in capturing the relatively weak dependency among edges to improve the clustering. This is because the proposed method not only captures pairwise dependency but also reflects connectivities among nodes within a community. That is, even a weak dependency among pairwise connectivities can lead to an accumulative information recovery of clustering.

Table 2 illustrates the clustering performance when the marginal signal is strong. In contrast to Table 1, the variational EM significantly improves on clustering because of the large discrepancy between the within-community marginal mean and the between-community marginal mean. Nevertheless, incorporating the correlation among within-community edges still improves the clustering accuracy by $20\%$ to $26\%$ under various sample sizes of networks and intensities of correlation. The clustering accuracy of the proposed method improves when either the sample size or the correlation increases. In general, stronger correlation and a larger sample size lead to better performance when the marginal signal itself is strong.

In addition to clustering, we also provide estimation of the marginal parameters. Tables 3, 4 and 5 compare parameter estimation between the proposed method and the variational EM when the marginal signal is weak. For within-community parameters $\beta_{11}$ and $\beta_{22}$, the estimation of the proposed method consistently reduces bias $30\sim 99\%$ more than the variational method, except when $M=20$ and $\rho=0.6$. This is because the sample size $M=20$ is not sufficiently large to offset the high variance among highly-correlated within-community edges. For the between-community parameter $\beta_{12}$, the estimation bias of the proposed method consistently decreases more than $80\%$ compared to the VEM under all settings. Additionally, the standard errors of the proposed estimator decrease faster than the variational method as the sizes of networks increase.

We also investigate the clustering performance of the independent likelihood and the proposed approximate likelihood approach given different within-community second-order correlation density $\lambda$ in (5.1). The setting is similar to the weak marginal signal cases. Specifically, the sample networks contain two communities with identical pairwise within-community correlation $\rho=0.6$. The sizes of the sample networks and nodes are $M=40,N=40$. The density $\lambda$ increases from 0.01 to 1. The Adjusted Rand Index comparisons are illustrated in Figure 2. In general, the approximate likelihood approach has improving performance when the correlation connectivities among within-community edges increase, in contrast to the independent likelihood approach. Figure 2 shows that the true membership recovery using the approximate likelihood approach is high even when the second-order within-community correlation is relatively sparse ($\lambda=0.05$), while the independent likelihood approach performs poorly with a constant ARI regardless of $\lambda$ . This finding supports Theorem 5.1 and 5.2 in that the proposed method produces an accelerated decay in misclassification rate as $\lambda$ increases.

6.2 Study 2: Networks with additional dependence between different communities

In Study 2, we also investigate whether the proposed method holds for a more general dependency structure among edges from different communities, for example, correlation among edges between different communities

[TABLE]

where $\tilde{\rho}\leq\rho_{q}$ in (3.2) in general. While (3.2) characterizes the concordance of edges within a community, (6.3) also captures the heterogeneity of sample networks. The heterogeneity of multi-layer networks is common in community detection.

In this simulation, we demonstrate that the proposed method is still robust when there is heterogeneity of connectivities among sample networks. To simulate the dependency among inter-community connectivity, we split $M$ sample networks into 10 groups. Within each group, we add the random effects $\gamma_{k}$ to the within-community marginal means:

[TABLE]

where $\gamma_{k}\sim N(0,\sigma^{2})$, $k=1,\cdots,10$, $m=1,\cdots,M$, and $q=1,\dots,K$. The variance $\sigma$ of the random effect $\gamma_{k}$ captures the intensity of dependency among inter-community connectivities, which increases as $\sigma$ increases. We set $\sigma=0.5$ to represent a weak inter-community dependency and $\sigma=1.5$ for a strong inter-community dependency, while the other settings remain the same as in simulation Study 1. Our primary interest is to compare clustering performance between the proposed method and the variational method under the weak marginal signal case.

Tables 6 and 7 illustrate the clustering performance between the variational method and the proposed method under balanced and unbalanced community sizes, respectively. When the within-community correlation is moderate at 0.3, the proposed method improves the clustering accuracy by $170\%$ to $257\%$ for various network sizes and $\sigma$. For strong correlation $\rho=0.6$, the improvement is between $210\%$ to $257\%$. In particular, the proposed method has better performance when the networks have strong intra-community correlation and large sample sizes under both weak and strong inter-community correlation cases. In addition, using the fourth-order Bahadur representation improves the accuracy by $6\%$ and $14\%$ when $\sigma=0.5$ and $\sigma=1.5$ compared to the second-order Bahadur representation, indicating that the higher-order method still enhances the clustering outcome under the misspecified model. It is interesting to note that the performance of the proposed method decreases by $5\%$ to $15\%$ when the inter-community correlation is strong and the number of networks is small, compared to the same setting with weak inter-community correlation. However, the performances under both weak or strong inter-community correlation are similar when the sample size of networks increases. In conclusion, the proposed method is robust against misspecified dependency structure when the sample size increases.

7 Real Data Example

In this section, we apply the proposed method to the 2010 Worldwide Food Import/Export Network dataset [18] from the Food and Agriculture Organization of the United Nations (http://www.fao.org). We create 364 networks among 214 countries with a total of 318,346 edges, where each network captures the trading connections of a specific food product among countries.

The primary goal of the study is to identify the food and agricultural product trading network communities among different countries. One significant feature of these networks is that the average empirical correlation of the pairwise connection among trading countries is 0.29. Therefore, the SBM based on the conditional independent assumption among edges could possibly lead to a biased network clustering of countries.

We first preprocess the data to select nodes corresponding to the trading countries which are most relevant, the number of communities and the initial memberships of countries. Note that several major countries dominate the world economy and lead a high number of trading connectivities, while the other countries with limited agricultural product categories have fewer trading connections with other countries for specific product networks. Here we focus on the partial trading networks consisting of major countries whose corresponding degrees of nodes are larger than 9, which results in 51 countries with major economic impact in the world, such as the United States, mainland China, Japan and some European countries. The average empirical correlation of the trading connections among these countries is 0.22, indicating that the connectivity dependency should be considered in clustering these countries’ trading networks.

In general, there are two major procedures to select the number of communities. First, we can perform the Louvain method for community detection on each individual trading network to obtain the number of communities which maximizes the modularity and the size of the largest community. Next we take the average of the number of communities on networks whose number of communities is smaller than 10 and whose largest community size is larger than 14. This procedure removes the $18\%$ of the product trading networks whose countries are commercially isolated from other countries, as our goal is to detect the commercial communities among the countries which are more connected with other countries. After preprocessing, the average number of communities is 4.9 and we set it to be 4, and there are 296 sample networks remaining in the following analysis.

Table 8 and Figure 2 provide the estimated agricultural products trading communities among 51 countries based on the variational EM and the proposed method. For the proposed method, we implement the fourth-order Bahadur approximation since it can better capture high-order within-community connectivity dependency. Table 8 presents the clustering outcome among countries according to the variational method and the proposed method. The countries in the same community under the variational method are marked with the same color, while the newly formed communities based on the proposed method are illustrated on the right sides of Table 8 and Figure 2. In general, the Adjusted Rand Index for clustering between the variational method and the proposed method is 0.43, indicating that the communities detected by the two methods are quite different. The clustering results from the proposed method incorporating within-community dependency are more interpretable compared to the variational EM using only marginal information.

In particular, the proposed method identifies communities 1 and 2 (red and cyan color communities on the right panel of Figure 2) which are highly associated with their geographical and climate environments. However, these features are not detected by the variational method. For example, community 1 with the cyan color on the left of Figure 2 based on the variational method mainly consists of two types of countries: one group comprises Nordic and Eastern European countries, and the other group consists of countries in Latin American and Africa. In contrast, the proposed method clusters countries from geographically neighboring countries in east Europe, including Austria, Poland and Romania which are clustered with other communities by the variational method. Community 2 with blue color on the left of Figure 2 based on the variational method contains northern countries such as Canada as well as tropical countries. However, the proposed method identifies community 2 with tropical coastal countries and Arabian Peninsula countries, which provides more meaningful community clusters compared to the variational EM method.

The variational method and proposed method detect the same third community with orange color in Figure 2 which contains 7 major countries from the European Union: Belgium, France, Germany, Italy, Netherlands, Spain and the UK.

The fourth community from the variational method colored with red on the left of Figure 2 consists of 11 Eastern European countries, and all are categorized in community 1 from the proposed method. Community 4 with blue color on the right of Figure 2 in the proposed method includes countries with large populations or more developed agricultural product trading, such as mainland China, U.S.A, India and Japan.

In terms of parameter estimation, the average probability of having trading connections for communities 1 and 2 based on the variational method are 0.21 and 0.52, respectively. For the proposed method, the estimated correlations of connectivities within communities 1 and 2 are both 0.22, and the corresponding average within-communities connection rates are 0.28 and 0.22, respectively. The relatively low connection rates and correlations may be related to the low diversity and high overlaps of product categories due to more restrictive geographical and climate environments.

For community 3, the corresponding estimated marginal parameters $\beta_{33}$ from the proposed method and the variational method are 2.58 and 2.00 respectively, both of which indicate that the trading connection rate within European Union communities is greater than $88\%$ on average. This strong marginal signal of within-community connection explains that the additional correlation information is less influential in clustering. Additionally, the estimated correlation within the third community is 0.58, implying a high connection rate within-community. For community 4, the corresponding average connection rate is 0.49 based on the variational method, and the estimated within-community average connection rate and the correlation are 0.61 and 0.27, respectively. This is because community 4 involves large population countries with more frequent trading on product categories due to their higher food diversity than other countries.

8 Discussion

In this paper, we propose a new community detection method for networks incorporating the underlying dependency structure among connectivities. To model the correlation without specifying a joint likelihood for correlated edges, we construct an approximate likelihood based on the Bahadur representation which decomposes a joint distribution into a marginal term and high-order interaction terms. The proposed method provides flexible modeling on the correlation structure which can be specified through the interaction term in the approximate likelihood.

In theory, we establish the consistency of the nodes’ membership estimator based on the proposed approximate likelihood and show that it achieves a faster convergence than the independent method. In addition, we show that the proposed iterative algorithm possesses desirable convergence properties. In particular, we show that the proposed approximate approach can achieve a faster computational convergence and a lower clustering bias compared to the variational EM algorithm. Furthermore, we show that the variational EM algorithm is a special case of our algorithm under the conditional independent model, which confirms that incorporating correlation information improves the accuracy for community detection.

Our numeric studies indicate that incorporating the within-community correlation among edges can improve the clustering performance compared to the marginal model, even under a moderately misspecified model on inter-community dependency. The improvement of community detection is more significant when the marginal signal is weak, which is less informative for distinguishing between within-community and between-community networks. In addition, the proposed method enables us to achieve more accurate parameter estimation.

In this paper, we only consider incorporating the within-community dependency. It would be worthy of further research to investigate more generalized dependency structures to include between-community dependency as well.

APPENDIX: NOTATION AND PROOFS

.1 Notation

In the following, we denote the membership of node as random variable $z_{i},i=1,\cdots,N$. Then $\bm{Z}=\{z_{1},z_{2},\cdots,z_{N}\}$. Accordingly, we define the true membership of nodes as $z^{*}_{i}\in\{1,2,\cdots,K\},i=1,\cdots,N$ and $z^{*}=\{z_{1}^{*},z_{2}^{*},\cdots,z_{N}^{*}\}$. We denote $P^{*}(\cdot)=P(\cdot|\bm{Z}=z^{*})$ as the conditional probability of observed networks given the true nodes’ membership $z^{*}$. The number of misclassified nodes is denoted as $r$ such that $\|z-z^{*}\|_{0}=r$ for $z\neq z^{*}$. Define the t-$th$ sample network as $\bm{Y}^{t}=(Y^{t}_{ij})_{N\times N}$ and t-$th$ sample network standardized by $\hat{\mu}_{aa}$ as $\bm{\hat{Y}}^{t,a}=(\hat{Y}^{t,a}_{ij})_{N\times N}$ where $\hat{Y}_{ij}^{t,a}=\frac{Y^{t}_{ij}-\hat{\mu}_{aa}}{\sqrt{\hat{\mu}_{aa}(1-\hat{\mu}_{aa})}},a=1,\cdots,K,\;t=1,\cdots,M$. We further define the s-$th$ column of $\bm{\hat{Y}}^{t,a}$ as $\hat{Y}^{t,a}_{\cdot s}$.

Denote $\bm{\alpha}=(\alpha_{1},\cdots,\alpha_{N})$ as the estimated probability of nodes’ memberships. Specifically, let $\alpha_{i}=(\alpha_{i1},\cdots,\alpha_{iK})_{1\times K}$ be the probability of nodes $i$ belonging to each community where $\sum_{q=1}^{K}\alpha_{iq}=1,\;i=1,\cdots,N$. For simplicity of notation, if the subscripts indicate the community then $\alpha_{q}=(\alpha_{1q},\cdots,\alpha_{Nq})_{1\times N}$ represents the probability of each node belonging to community $q$, where $q=1,\cdots,K$. Similarly, $z_{q}^{*}=\{z^{*}_{1q},z^{*}_{2q},\cdots,z^{*}_{Nq}\}$ is a binary vector indicating nodes whose true membership belongs to community $q,q=1,\cdots,K$. Let $vec(\cdot)$ stand for the operation of vectorizing a matrix into a column.

.2 Lemmas

The following two lemmas are introduced as the technical steps in the proofs of Theorem 5.1, Theorem 5.2 and Theorem 5.3. The proofs of Lemma 1 and Lemma 2 are provided in the supplemental material.

Lemma 1**.**

Consider function f_{1}(x)=\sqrt{\big{\{}x\log\frac{\mu_{z_{i}z_{j}}}{\mu_{z_{i}^{*}z_{j}^{*}}}+(1-x)\log\frac{1-\mu_{z_{i}z_{j}}}{1-\mu_{z_{i}^{*}z_{j}^{*}}}\big{\}}_{+}} and denote

[TABLE]

where $\{Y_{ij}^{t}\}_{N\times N}$ are generated through the stochastic block model in section 3.1 and satisfy condition C1, C2 and C3. Define the covariance matrix of $X_{t}^{+}$ as $\Sigma_{1}$. Then $X_{t}^{+}$ is a subgaussian vector, i.e.,

[TABLE]

for some positive constant $C$.

Lemma 2**.**

Assume $\{Y_{ij}^{t}\}_{N\times N}$ are generated through the stochastic block model in section 3.1 and satisfy condition C1, C2 and C4. Then given $M>\mathcal{O}(\frac{1}{\lambda})$, we have

[TABLE]

as $M,N$ increase for $q=1,2,\cdots,K$.

.3 Proof of Theorem 5.1

Given the independent model in (2.4), we can simplify the likelihood ratio between a random membership $z$ and the true membership $z^{*}$ as

[TABLE]

We define two transformation functions $f_{1}(x)$ and $f_{2}(x)$ as:

[TABLE]

where $\{\}_{+}$ and $\{\}_{-}$ are positive part and negative part of a random variable. The previous summation can be decomposed as positive part and negative part:

[TABLE]

Define the vectorized edges in the $t$ th sample network as:

[TABLE]

Note that each element in $X_{t}^{+}$ or $X_{t}^{-}$ is a bounded binary random variable. In addition, as $f_{1}(Y_{ij}^{t})$ or $f_{2}(Y_{ij}^{t})$ only rescale $Y_{ij}^{t}$ then they preserve the within-community correlation among $Y_{ij}^{t}$. Then we consider the following quadratic forms

[TABLE]

such that

[TABLE]

For any $t>0$, we have

[TABLE]

Next, we estimate each of the term in (.3). Given the $\{Y_{ij}^{t}\}_{t=1}^{M}$ are binary random variables and the setting that any two within-community edges $Y_{i_{1}j_{1}}$ and $Y_{i_{2}j_{2}}$ have a nonnegative correlation $corr(Y_{i_{1}j_{1}},Y_{i_{2}j_{2}})\geq 0$. Notice that

[TABLE]

We denote the covariance matrix of $X_{t}^{+}$ and $X_{t}^{-}$ as $\Sigma_{1}$ and $\Sigma_{2}$. Notice that a term in (.1) is zero only when its corresponding node membership is misclassified. Define the the number of nonzero term in (1) as $N_{r}$ given $\|z-z^{*}\|_{0}=r$. Then we have $N_{r}=\frac{1}{2}rNM$. According to Lemma 1, $X_{t}^{+}$ is a subgaussian vector with a bounded subgaussian norm $L\leq C_{1}$ where $C_{1}$ is a positive constant and

[TABLE]

Next we estimate $\|\Sigma_{1}\|_{F},\|\Sigma_{1}\|_{op}$ and $\|\Sigma_{2}\|_{F},\|\Sigma_{2}\|_{op}$ where $\|\cdot\|_{F}$ is the matrix Frobenius norm and $\|\cdot\|_{op}$ is the matrix spectral norm. Denote

[TABLE]

Then $\|\Sigma_{1}\|_{op}=\|\Lambda R\Lambda\|_{op}\leq C_{2}\|R\|_{op}$ where $R$ is the correlation matrix of $X_{t}^{+}$ and

[TABLE]

Denote the largest eigenvalue of $R$ as $\lambda_{R}$. From the Gershgorin circle theorem, we have

[TABLE]

Denote the number of node in the largest community is $N_{k}$. Note that the misclassification number of node $\|z-z^{*}\|_{0}=r$ and edgewise correlation density $\lambda$ both affect the sparsity of $R$, we have for each row in $R$:

[TABLE]

where $\rho=\displaystyle\max_{i,j}R_{ij}$. Therefore, we have

[TABLE]

Similarly we have a same upper bound for $\|\Sigma_{2}\|_{op}$. Notice that the dimension of $R$ is $N_{r}\times N_{r}$ and $N_{r}\leq rN$. In each row of $R$, the number of non-zero elements is less than $1+N_{k}\min(r,\lambda N_{k})$. Therefore, we have

[TABLE]

Then we are able to estimate the upper bound for the first term in (.3). According to the generalized Hanson-Wright inequality in ([14]), we have:

[TABLE]

where $s=\frac{M\log t-E(Q_{1}-Q_{2})}{2}$, $A=\bm{I}_{M\times M}$ and $L$ is subgaussian norm of $X_{t}^{+}$ defined in (.4). Then we have $L\leq C_{1}$ and $\|A\|^{2}_{F}=M,\|A\|_{op}=1$. To estimate $s$, notice

[TABLE]

where there are total $N_{r}$ non-zero terms in the summation. We introduce the function

[TABLE]

Notice that $k(x,y)>0$ for every $x,y\in(0,1)$. Then we define:

[TABLE]

where the minimum are taken over $\left\{\left((q,l),\left(q^{\prime},l^{\prime}\right)\right)|\mu_{q,l}^{*}\neq\mu_{q^{\prime},l^{\prime}}^{*}\right\}$. Combined with $N_{r}=\frac{1}{2}rNM$, we have $-E(Q_{1}-Q_{2})>\frac{c^{*}}{2}rNM$. Then for any fixed $t>0$, $s>\mathcal{O}_{N}(\frac{c^{*}}{2}rNM)$. Therefore, we have

[TABLE]

where $C_{3}=\frac{\min(c^{*},c^{*2})}{2C_{1}^{2}\max\{\log\frac{c_{1}}{c_{2}},\log\frac{1-c_{2}}{1-c_{1}}\}}$. Hence for (.5) we have:

[TABLE]

where $C^{*}=\frac{C\min(c^{*},c^{*2})}{2C_{1}^{2}\max\{\log\frac{c_{1}}{c_{2}},\log\frac{1-c_{2}}{1-c_{1}}\}}=\frac{c}{\max\{\log\frac{c_{1}}{c_{2}},\log\frac{1-c_{2}}{1-c_{1}}\}}$ and $c:=\frac{C\min(c^{*},c^{*2})}{2C_{1}^{2}}>0$. Follow Lemma 1, $X_{t}^{-}$ is also subgaussian vector. Then we can obtain a same upper bound for

[TABLE]

in (.3) through the above procedure. Therefore,

[TABLE]

.4 Proof of Theorem 5.2

We continue use the notations in the previous proof of Theorem 5.1. First decompose the proposed approximate likelihood in two parts:

[TABLE]

Follow Lemma 2 without assuming condition C4, P\Big{(}\displaystyle\sum_{\begin{subarray}{c}i<j;u<v\\ (i,j)\neq(u,v)\end{subarray}}^{N}z^{*}_{ik}z^{*}_{jk}z^{*}_{uk}z^{*}_{vk}\hat{Y}_{ij}^{t,k}\hat{Y}_{uv}^{t,k}\geq 0\Big{)} goes to 1 as $M,N$ increase. Based on the mean value theorem, we have for some constant $C_{1}$ that

[TABLE]

Notice in summation (.4), the terms are non-zero only when $z_{ik}z_{jk}z_{uk}z_{vk}\neq z^{*}_{ik}z^{*}_{jk}z^{*}_{uk}z^{*}_{vk}$. We denote two node sets

[TABLE]

where $\#|\xi_{1}|=N_{1}$ and $\#|\xi_{2}|=N_{2}$. Given the number of misclassified nodes $\|z-z^{*}\|_{0}=r$, we have $N_{1}=\mathcal{O}(rN^{3})$ and $N_{2}=\mathcal{O}(rN^{3})$. In the following, we construct the augmented edge vectors for the $t$ th sample network by incorporating the vectorized pairwise edge interaction in (.4) such that:

[TABLE]

where $X_{t}^{+}$ and $X_{t}^{-}$ are defined in (.2). Denote the covariance matrix for $\tilde{X}_{t}^{+}$ and $\tilde{X}_{t}^{-}$ are $\tilde{\Sigma}_{1}$ and $\tilde{\Sigma}_{2}$ respectively. Since the second-order terms in $X_{t}^{+}$ and $X_{t}^{-}$ such as $\sqrt{\frac{C_{1}\rho_{k}}{2}\{\hat{Y}_{ij}^{t,k}\hat{Y}_{uv}^{t,k}\}_{+}}$ only rescale the original edgewise interaction $\hat{Y}_{ij}^{t,k}\hat{Y}_{uv}^{t,k}$ then they preserve the third-order and fourth-order correlation within communities such that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Notice that each element in $\tilde{X}_{t}^{+}$ or $\tilde{X}_{t}^{-}$ is a bounded binary random variable. Follow the same procedure in Lemma 1, we can show that both $\tilde{X}_{t}^{+}$ and $\tilde{X}_{t}^{-}$ are subgaussian random vectors such that $L_{1}\leq C_{2},L_{2}\leq C_{2}$ for some constant $C_{2}$ where $L_{1},L_{2}$ are subgaussian norm of $\tilde{X}_{t}^{+}$ and $\tilde{X}_{t}^{-}$. Then consider the following quadratic forms:

[TABLE]

Therefore, we have

[TABLE]

Similar to (.3), for any fixed $t>0$:

[TABLE]

Next we estimate $\|\tilde{\Sigma}_{1}\|_{F},\|\tilde{\Sigma_{1}}\|_{op}$ and $\|\tilde{\Sigma}_{2}\|_{F},\|\tilde{\Sigma}_{2}\|_{op}$. Denote

[TABLE]

then $\|\tilde{\Sigma}_{1}\|_{op}=\|\tilde{\Lambda}\tilde{R}\tilde{\Lambda}\|_{op}\leq C_{3}\|\tilde{R}\|_{op}$ where $\tilde{R}$ is the correlation matrix of $\tilde{X}_{t}^{+}$ and $C_{3}$ is the largest variance of elements in $\tilde{X}_{t}^{+}$. Denote the largest eigenvalue of $\tilde{R}$ as $\lambda_{\tilde{R}}$. From the Gershgorin circle theorem, we have

[TABLE]

Given that the misclassification number of node $\|z-z^{*}\|_{0}=r$, edgewise correlation density $\lambda$ and condition C3, for each row in $\tilde{R}$, there exists some constant $C_{4}$¿0 such that:

[TABLE]

where $\rho=\displaystyle\max_{i,j}\tilde{R}_{ij}$. Therefore, we have

[TABLE]

Similarly, $\|\tilde{\Sigma}_{2}\|_{op}$ follows a same upper bound. Notice that the dimension of $\tilde{R}$ is $(N_{r}+N_{1}+N_{2})\times(N_{r}+N_{1}+N_{2})$. Under the condition C3, in each row of $\tilde{R}$, the number of non-zero elements is less than $1+C_{4}N_{k}\min(r,\lambda N_{k})$. Therefore, we have for a constant $C^{\prime}>0$:

[TABLE]

According to the generalized Hanson-Wright inequality in ([14]):

[TABLE]

where $s=\frac{M\log t-E(\tilde{Q}_{1}-\tilde{Q}_{2})}{2}$, $A=\bm{I}_{M\times M}$ and $L_{1}$ is subgaussian norm of $\tilde{X}_{t}^{+}$. Notice $\|A\|^{2}_{F}=M,\|A\|_{op}=1$. Given (.4), we have

[TABLE]

Denote $\rho_{min}$ as the lower bound of all non-zero correlation among edges such that $E(\hat{Y}_{ij}^{t,k}\hat{Y}_{uv}^{t,k})=\rho_{ijuv}\geq\rho_{min}$. Given the edges from different communities are independent and within-community correlation density $\lambda$, we have for some positive constant $C_{5}$,

[TABLE]

[TABLE]

Assume that $r\leq cN$ for some constant $0<c<1$, we have for some constant $c_{0}>0$:

[TABLE]

Given any fixed $t>0$, $s>\mathcal{O}_{N}(c_{0}r(N+\lambda N^{3})M)$. For the first term in (.10),

[TABLE]

For the second term in (.10),

[TABLE]

Given $\lambda>0$, we have for some constant $C_{6}>0$

[TABLE]

Follow the same procedure we can show a upper bound for P^{*}\Big{\{}|\tilde{Q}_{2}-E\tilde{Q}_{2}|>s\Big{\}} with a same order to (.11). Combined with (.4) and (.10), we have for $\lambda>0$ and some constant $C>0$:

[TABLE]

.5 Proof of Theorem 5.3

Follow the notations introduced in Theorem 5.1 and Theorem 5.2, we further define that $\bm{w}=\max{\frac{P^{(s)}(Z_{i}=q)}{P^{(s)}(Z_{i}=l)}},\;i=1,\cdots,N,\;q,l=1,\cdots,K$. Let $\bm{E}$ stands for the operator of expectation step in Algorithm 1 in Section 4.

We first consider the misclassification of updated estimated membership for node $s$, e.g., $\bm{E}(z_{s})$ from the current estimation $\alpha_{s}$. We denote that $\bm{\alpha}_{-s}$ as the probability estimations of nodes’ memberships at current step except node $s$ and assume the true membership of node $s$ is $b$, i.e., $z^{*}_{s}=b$. If we use the marginal likelihood, then:

[TABLE]

Then given node $s$ belongs to different communities while the estimated membership for other nodes $\bm{\alpha}_{-s}$ are fixed. We decompose the proposed approximate likelihood into marginal part and correlation part in the following: $\log\tilde{L}(\bm{Y}|\bm{\alpha}_{-s};z_{s})=\log L_{mar}(\bm{Y}|\bm{\alpha}_{-s};z_{s})+\log L_{cor}(\bm{Y}|\bm{\alpha}_{-s};z_{s})$. The marginal likelihood $\log L_{mar}(\bm{Y}|\bm{\alpha}_{-s};z_{s})$,

[TABLE]

Therefore, the discrepancy among marginal likelihood is

[TABLE]

We can decompose the marginal discrepancy into four parts:

[TABLE]

For the correlation part, we consider the pairwise interaction terms in the $\log L_{cor}(\bm{Y}|\bm{\alpha})$. Notice that for $t=1,\cdots,M$

[TABLE]

where $A_{q}^{t}$ does not depend on $z_{s}$. Since the first term $(\sum_{i\neq s}^{N}\alpha_{sa}\alpha_{ia}\hat{Y}^{t,a}_{si})(\sum_{i<j}^{N}\alpha_{ia}\alpha_{ja}\hat{Y}^{t,a}_{ij})=o(N^{3})$ and the second term $\sum_{i\neq s}^{N}(\alpha_{ia}\hat{Y}^{t,a}_{si})^{2}=o(N)$, without loss of generality, we can keep the first dominating term when $N$ is large. For the correlation part $\log L_{cor}(\bm{Y}|\bm{\alpha}_{-s};z_{s})$, if $\alpha_{sq}=0,\;q\neq a$ and $\alpha_{sa}=1$, with probability approaching 1 from Lemma 2 as $M$ increases given $M>\mathcal{O}(\frac{1}{\lambda})$:

[TABLE]

Through the Taylor expansion, the discrepancy of correlation information when node $s$ belongs to different communities $a$ and $b$:

[TABLE]

where $\bm{C}_{A}$ is a constant relating to the gradient of function $\log(1+1/x)$ at $\bm{A}$. Then we set $\rho=\min{\rho_{q}},q=1,\cdots,K$

[TABLE]

For the simplicity of notation, we define and decompose the correlation discrepancy as followings:

[TABLE]

Notice that $\min\{1,\exp(x)\}\leq\exp(x_{0})+\sum_{l=0}^{m-1}\frac{1-\exp(x_{0})}{m}\mathbbm{1}\{x\geq(1-l/m)x_{0}\}$ and set $x_{0}=-\alpha^{\prime}MN$, where $\alpha^{\prime}>0$. Given (.5), we have:

[TABLE]

where E(\bm{L})=\mathbbm{P}\Big{(}\bm{A}+\bm{B}\geq\frac{m-l}{m}x_{0}\Big{)}. For some specific $t>0$,

[TABLE]

We then transfer $\bm{A_{3}}+\bm{B_{2}}$ into a quadratic form. For each community $q,q=1,\cdots,K$ define the transformations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Notice that the total number of non-zero terms in $X_{t}^{+}$ or $X_{t}^{-}$ is $N$. We define the node sets

[TABLE]

Note $\#|\tilde{\xi}_{a}|=o(N_{a}^{3})$ and $\#|\tilde{\xi}_{b}|=o(N_{b}^{3})$ where $N_{a}$ and $N_{b}$ are number of node in community $a$ and $b$. We further define augmented edges vectors:

[TABLE]

Denote the covariance of $\bar{X}_{t}^{+}$ and $\bar{X}_{t}^{-}$ as $\bar{\Sigma}_{1}$ and $\bar{\Sigma}_{2}$. Note that each element in $\bar{X}_{t}^{+}$ or $\bar{X}_{t}^{-}$ is a bounded binary random variable. Similarly, $\bar{X}_{t}^{+}$ and $\bar{X}_{t}^{-}$ are subgaussian vectors. Therefore,

[TABLE]

Denote $s=\frac{m-l}{m}x_{0}-t-\bm{A_{2}}-\bm{A_{4}}-E(\bm{A_{3}}+\bm{B_{2}})$, we estimate $E(\bm{A_{3}}+\bm{B_{2}})$, $\bm{A_{2}}$ and $\bm{A_{4}}$ in the following. Given $z^{*}_{s}=b$ and the result in (.6), we have for some constant $c>0$ and $q=1,\cdots,K$:

[TABLE]

Then

[TABLE]

Given edges from different communities are independent and correlation density $\lambda$, there exists a constant $C>0$ such that

[TABLE]

Therefore, $-E(\bm{A_{3}}+\bm{B_{2}})\geq c(N-1)+C\lambda(\kappa_{1}N)^{3}\geq c^{\prime}(N+\lambda N^{3})$ for some positive constant $c^{\prime}$. Based on condition C1 that $\mu_{ql},q,l=1,\cdots,K$ are bounded, it can be shown that $|EY_{is}^{t}\log\frac{\mu_{qa}}{\mu_{qb}}+(1-EY_{is}^{t})\log\frac{1-\mu_{qa}}{1-\mu_{qb}}|$ is bounded then $|\bm{A_{2}}|=\mathcal{O}_{N}(N)$.

From condition C5, we have

[TABLE]

Define $\gamma=\max\{-\log\frac{\gamma_{1}}{\gamma_{2}},\frac{\gamma_{2}}{\gamma_{1}},-\frac{1-\gamma_{2}}{1-\gamma_{1}},\frac{1-\gamma_{1}}{1-\gamma_{2}}\}$. Then we have

[TABLE]

Therefore we have $|\bm{A_{2}}+\bm{A_{4}}|=\mathcal{O}_{N}(N)$. We choose $t=-\frac{E(\bm{A_{3}}+\bm{B_{2}})}{2}$ and $x_{0}=-\alpha^{\prime}MN$ where $\alpha^{\prime}>0$. As the function of node size $N$, $M$ and $\lambda$ are constrained in the range $M\leq o(N^{2-\frac{\eta}{2}})$ and $\lambda N^{\frac{\eta}{2}}>1$, where $\eta$ is defined in condition C4. Then $\frac{m-l}{m}x_{0}=o_{N}(E(\bm{A_{3}}+\bm{B_{2}}))$. Therefore, $E(\bm{A_{3}}+\bm{B_{2}})$ is dominant term in $s$ such that $s\geq-C^{\prime}\lambda N^{3}$ where $C^{\prime}>0$ is a constant. Follow a similar discussion in (.9) and condition C3, we have the upper bound for $\|\bar{\Sigma}_{1}\|_{op}$:

[TABLE]

In addition, from $\#|X_{t}^{+}|=N$, $\#|\bar{\xi}_{a}|=o(N_{a}^{3})$, $\#|\bar{\xi}_{b}|=o(N_{b}^{3})$ and condition C3, we have the upper bound for $\|\bar{\Sigma}_{1}\|_{F}^{2}$:

[TABLE]

where $C_{1},C_{2},c,c^{*}$ are constants. Then we estimate the upper bound for the second term in (.5):

[TABLE]

According to the generalized Hanson-Wright inequality in ([14]):

[TABLE]

where $A=\bm{I}_{M\times M}$ and $\bar{L}_{1}$ is subgaussian norm of $\bar{X}_{t}^{+}$. Notice that

[TABLE]

Given $\lambda N^{\frac{\eta}{2}}>1$, we have for some constant $C^{*}>0$

[TABLE]

The upper bound for \mathbbm{P}\Big{\{}|\bar{Q}_{2}-E\bar{Q}_{2}|>\frac{Ms}{2}\Big{\}} can be similarly obtained. Therefore,

[TABLE]

Next, we estimate the term \mathbbm{P}\Big{(}\bm{A_{1}}+\bm{B_{1}}\geq t\Big{)}. Notice

[TABLE]

Given condition C4 such that $\|\bm{\alpha}-z^{*}\|_{1}=cN^{1-\eta},0<\eta<1$,

[TABLE]

Notice that for any community $a=1,\cdots,K$,

[TABLE]

Therefore, we have

[TABLE]

Since $\|\alpha_{a}-z_{a}^{*}\|_{2}=\sqrt{\|\alpha_{a}-z_{a}^{*}\|_{2}^{2}}\leq\sqrt{\|\bm{\alpha}-z^{*}\|_{1}}$ for any $a=1,\cdots,K$, then for some constant $C>0$,

[TABLE]

We define edge vectors $\tilde{Y}_{t},t=1,\cdots,M$ and membership vector $\bm{\theta}_{a,b}$ as:

[TABLE]

Notice for $a,b=1,\cdots,K$, we have

[TABLE]

where $\mu_{2}:=\max\{(\log\frac{\mu_{qa}}{\mu_{qb}}-\log\frac{1-\mu_{qa}}{1-\mu_{qb}})\},\;q=1,\cdots,K$ and $C_{1}>0$ is a constant. Then we can transform ${\rm Var}(\bm{A_{1}}+\bm{B_{1}})$ into

[TABLE]

From the condition C3 and same discussion in (.9), we have for some constant $C>0$ and $c>0$:

[TABLE]

Given $\frac{1}{\lambda}=o(N^{\frac{\eta}{2}})$, we have $E(\bm{A_{1}}+\bm{B_{1}})=o_{N}(E(\bm{A_{3}}+\bm{B_{2}}))$ then the $E(\bm{A_{3}}+\bm{B_{2}})$ is dominating in the term ${\{t-E(\bm{A_{1}}+\bm{B_{1}})\}^{2}}$. Based on the Markov inequality, for some constant $C_{2}>0$

[TABLE]

Combined upper bound of \mathbbm{P}\Big{(}\bm{A_{1}}+\bm{B_{1}}\geq t\Big{)} and \mathbbm{P}\Big{(}\bm{A_{3}}+\bm{B_{2}}\geq s\Big{)} with (.13), there exists positive constant $c_{1}>0,c_{2}>0,c_{3}>0$ such that:

[TABLE]

{supplement}\sdescription

Supplement to ”Community detection with dependent connectivity”. Due to space constraints, we relegate proofs of Lemma 1 and Lemma 2 to the supplement.