Detecting Latent Communities in Network Formation Models

Detecting Latent Communities in Network Formation Models

Introduction

In real world social and economic networks, individuals tend to form links with someones who are alike to themselves, resulting in assortative matching on observed individual characteristics (homophily). In addition, network data often exhibit natural communities such that individuals in the same community may share similar preferences for a certain type of homophily while those in different communities tend to have quite distinctive preferences. In many cases, such a community structure is latent and has to be identified from the data. The detection of such community structures is challenging yet crucial for network analyses. It prompts a couple of important questions that need to be addressed: How do we formulate a network formation model with individual characteristics, unobserved edge-wise fixed effects, and latent communities?\ When the model is formulated, how do we recover the community structure and estimate the community-specific parameters effectively in the model?

To address the first issue above, we propose a logistic undirected network formation model with observed measurements of homophily as regressors. We allow the regression coefficients to have a latent community structure such that the regression coefficient for covariate $l$ in the network formation model is $B_{l,k_{1}k_{2}}$ for any nodes $i$ and $j$ in communities $k_{1}$ and $k_{2}$, respectively. The edge-wise fixed effects are assumed to have a low-rank structure. This includes the commonly used discretized fixed effects and additive fixed effects as special cases. To address the second issue, we note that the estimation of this latent model is challenging, and it has to involve a multi-step procedure. In the first step, we estimate the coefficient matrices by a nuclear norm regularized logistic regression\ given their low-rank structures; we then obtain the estimators of their singular vectors which contain information about the community memberships via the singular value decomposition (SVD). Such singular vector estimates are only consistent in Frobenius norms but not in\ uniform row-wise Euclidean norm. A refined estimation is needed for accurate community detection. In the second step, we use the singular vector estimates from the first step as the initial values and iteratively run row-wise and column-wise logistic regressions to reestimate the singular vectors. \ The efficiency of the resulting estimator can be improved through this iterative procedure. In the third step, we apply the standard K-means algorithm to the singular vector estimates obtained in the second step. For technical reasons, we have to resort to sample-splitting techniques to estimate the singular vectors, and for numerical stability,\ both iterative procedures and multiple-splits are called upon. We establish the exact recovery of the latent community (strong consistency) under the condition that the expected degree of the network diverges to infinity at the rate $ \log n$ or higher order, where $n$ is the number of nodes. Under the exact recovery property, we can treat the estimated community memberships as the truth and further estimate the\ community-specific regression coefficients.

Our paper is closely related to three strands of literature in statistics and econometrics. First, our paper is closely tied to the large literature on the application of spectral clustering to detect communities in stochastic block models (SBMs). Since the pioneering work of HLL83, SBM has become the most popular model for community detection. The statistical properties of spectral clustering in such models have been studied by J15, JY16, LR15, PC20, QR13, RCY11, SB15, ss15, V18, ww87, YP14 , and YP16, among others. From an information theory perspective, AS15, ABH2016, MNS14, and V18 establish the phase transition threshold for the exact recovery of communities in SBMs, which requires the expected degree to diverge to infinity at a rate no slower than $\log n$. SWZ20 show that spectral clustering can achieve this information-theoretical minimum rate for the exact recovery. Nevertheless, most existing SBMs do not include covariates. A few exceptions include B17, B17, Weng_Feng2016, Yan_Sarkar2020 and Zhang_Levina_Zhu2016, who consider covariates-assisted community detection but not inferences on the underlying parameters. For more complicated models that can incorporate both covariates and community structures, people often resort to the variational EM algorithm, the performance of which highly hinges on the proper choice of initial values. In contrast, the\ network formation model proposed in this paper extends the SBM to a complex logistic regression model with both latent community structures and covariates, and our multi-step procedure provides an effective and reliable tool for the estimation of such a complex network model. Despite the fact that the regression coefficient matrices have to be estimated from the data in order to obtain the associated singular vectors for spectral clustering, we are able to obtain the exact recovery of the community structures at the minimal rate on the expected node degree and to conduct inferences on the underlying parameters in the model. \

Second, our paper is closely tied to the burgeoning literature on network formation models and panel structure models. For the former, see CDS11 , G17, Graham2019, Graham2020, Graham_de_Paula2019, HL83, Hoff_Raftery_Handcock2002, J19, L15, M17, RPF13, and YX13. We complement these works by allowing for community structures on the regression coefficients, which can capture a rich set of unobserved heterogeneity in the network data. In a working paper, M2017b also considers a network formation model with heterogeneous players and latent community structure. He\ assumes that the community structure follows an i.i.d. multinomial distribution and imposes a prior distribution over communities and parameters before conducting Bayesian estimation and inferences. In contrast, we treat the community memberships as\ fixed parameters and aim to recover them from a single observation of a large network. Our idea of introducing the community structure into the network formation model is mainly inspired by the recent works of BM15 and SSP16, who introduce latent group structures into panel data analyses. When the community structure is unobserved, it is analogous to the latent group structure in panel data models. For recent analyses of panel data models with latent group structures, see Ando_Bai2016, Chen2019, CSS2019, Dzemski_Okui2018, HJS2020, HJPS2021, LSZZ2020, Lu_Su2017, Su_Ju2018, SWJ2019, Vogt_Linton2020, Wang_Su2021, and Xu_Yue_Zhang2020, among others. In particular, Wang_Su2021 establish the connection between SBMs and panel data models with latent group structures and propose to adopt the spectral clustering techniques to recover the latent group structures in panel data models.

Last, our paper is related to the literature on the use of nuclear norm regularization in various contexts; see aal20, B19, BN19 , CHLZ18, FGZ2019, F19, KLT2011, M18, NW11, NRWY2012, and RT2011, among others. All these previous works focus on the error bounds (in Frobenius norm) for the nuclear norm regularized estimates, except M18 and CHLZ18 who study the inference problem in linear panel data models with a low-rank structure. Like M18 and CHLZ18, we simply use the nuclear norm regularization to obtain consistent initial estimates. Unlike M18 and CHLZ18, we study a logistic network formation model with a latent community structure and propose the iterative row- and column-wise logistic regressions to improve the error bounds (in row-wise Euclidean norm) for the singular vectors of the nuclear norm regularized estimates. Relying on such an improvement, we can fully recover the community memberships. Then, we can estimate the community specific parameters and make statistical inferences.

The rest of the paper is organized as follows. In Section (ref), we introduce the model and basic assumptions. In Section (ref) , we provide our multi-step estimation procedure. Section (ref) establishes the statistical properties of the proposed estimators of the singular vectors. Section (ref) studies the K-means estimation of the community memberships when the regression coefficient matrix is assumed to exhibit some community structure. Section (ref) studies the asymptotic properties of the regression coefficient estimates in the presence of latent community structures. Section (ref) discusses the determination of the ranks of the regression coefficient matrices. Section (ref) reports simulation results. In Section (ref), we apply the new methods to study the community structure of a Facebook friendship networks at one hundred American colleges and universities at a single time point. Section (ref) concludes. The online supplement provides the proofs of all theoretical results and the associated technical lemmas, and some additional technical details.

Notation. Throughout the paper, we write \textquotedblleft w.p.a.1" for \textquotedblleft with probability approaching one," $M=\{M_{ij}\}$ as a matrix with its $(i,j)$-th entry denoted as $M_{ij}$. We use $||\cdot ||_{op} $, $||\cdot ||_{F}$, and $||\cdot ||_{\ast }$ to denote matrix spectral, Frobenius, and nuclear norms, respectively. We use $[n]$ to denote $\{1,\cdots ,n\}$ for some positive integer $n$. For a vector $u$, $||u||$ and $u^{\top }$ denote its $L_{2}$ norm and transpose, respectively. For a vector $a=(a_{1},\cdots ,a_{n})$, let $\text{diag}(a)$ be the diagonal matrix whose diagonal is $a$. For a symmetric matrix $B\in \mathbb{R} ^{K\times K}$, we define

We define $\max (u,v)=u\vee v$ and $\min (u,v)=u\wedge v$ for two real numbers $u$ and $v$. We write $\mathbf{1}\{A\}$ to denote the usual indicator function that takes value 1 if event $A$ happens and 0 otherwise. Let $\odot $ denote Hadamard product.

The Model and Basic Assumptions

In this section, we introduce the model and basic assumptions.

The Model

For $i\neq j\in \left[ n\right] $, let $Y_{ij}$ denote the dummy variable for a link between nodes $i$ and $j$. It takes value 1 if nodes $i$ and $j$ are linked and 0 otherwise. Let $W_{ij}=(W_{1,ij},...,W_{p,ij})^{\top }$ denote a $p$-vector of measurements of homophily between nodes $i$ and $j$. Researchers observe the network adjacency matrix $\{Y_{ij}\}$ and covariates $\{W_{ij}\}$. We model the link formation between $i$ and $j$ is as

where $\{\zeta _{n}\}_{n\geq 1}$ is a deterministic sequence that may decay to zero and is used to control the expected degree in the network, $ W_{0,ij}=1,$ and $W_{l,ij}=W_{l,ji}$ for $j\neq i$ and $l\in [ p]$. For clarity, we consider undirected network so that $Y_{ij}=Y_{ji}$ and $\Theta _{l,ij}^{\ast }=\Theta _{l,ji}^{\ast }$ $\forall l$ if $i\neq j,$ $ \varepsilon _{ij}$ follows the standard logistic distribution for $i<j$, and $\varepsilon _{ij}=\varepsilon _{ji}$. Let $Y_{ii}=0$ for all $i\in \left[ n \right] $.

Apparently, without making any assumptions on $\Theta _{l}^{\ast }=\{\Theta _{l,ij}^{\ast }\}$ for $l\in \left[ p\right] \cup \{0\},$ one cannot estimate all the parameters in ((ref)) as the number of parameters can easily exceed the number of observations in the model. Specifically, we will follow the literature on reduced rank regressions and assume that each $ \Theta _{l}^{\ast }$ exhibits a certain low rank structure. Even so, it is easy to see that our model in (ref) is fairly general, and it includes a variety of network formation models as special cases.

Under the assumptions specified in the next subsection, it is easy to see that the expected degree of the network is of order $n\zeta _{n}.$ In the theory to be developed below, we allow $\zeta _{n}$ to shrink to zero at a rate as slow as $n^{-1}\log n$, so that the expected degree can be as small as $C\log n$ for some sufficiently large constant $C$ and the network is semi-dense.\footnote{ A network is dense if the expected degree grows at rate-$n$ and semi-dense if it diverges to infinity at a rate slower than $n$.} Of course, if $\zeta _{n}$ is fixed or convergent to a positive constant as $n\rightarrow \infty , $ the network becomes dense.

To proceed, let $\tau _{n}=\log (\zeta _{n})$, $\Gamma _{0,ij}^{\ast }=\tau _{n}+\Theta _{0,ij}^{\ast }$, $\Gamma _{ij}^{\ast }=(\Gamma _{0,ij}^{\ast },\Theta _{1,ij}^{\ast },...,\Theta _{p,ij}^{\ast })^{\top }$, and $ W_{ij}=(W_{0,ij},W_{1,ij},...,W_{p,ij})^{\top }$, where $W_{0,ij}=1$. Let $ \Gamma ^{\ast }=(\Gamma _{0}^{\ast },\Theta _{1}^{\ast },...,\Theta _{p}^{\ast }),$ where $\Gamma _{0}^{\ast }=\{\Gamma _{0,ij}^{\ast }\}$ and $ \Theta _{l}^{\ast }=\{\Theta _{l,ij}^{\ast }\}$ for $l\in \left[ p\right] .$ Then, we can rewrite the model in ((ref)) as

Below, we will let $\Gamma _{l}^{\ast }=\Theta _{l}^{\ast }$ for $l\in \left[ p\right] $ and impose some basic assumptions on the model in order to propose a multiple-step procedure to estimate the parameters of interest in the model.

Basic Assumptions

Now, we state a set of basic assumptions to characterize the model in (ref). The first assumption is about the data generating process (DGP).

Assumption (ref) specifies how the covariates and error terms are generated. In some applications, $e_{ij}$ is absent and $W_{l,ij}$ depend on $(X_{i},X_{j})$ only. For example, $W_{l,ij}=\left\Vert X_{i}-X_{j}\right\Vert $ for some $l$. We further assume that it is uniformly bounded to simplify the analysis. Assumption (ref).2 is standard.

The next assumption imposes some structures on $\left\{ \Theta _{l}^{\ast }\right\} _{0\leq l\leq p}.$

We note (ref) implies that $\Theta _{l,ij}^{\ast }=u_{i,l}^{\top }v_{j,l}.$ We view $\Theta _{0,ij}^{\ast }$ as the edge-wise fixed effects for the network formation model. We impose the normalization that $ \sum_{i,j\in \lbrack n]}\Theta _{0,ij}^{\ast }=0$ in the first part of Assumption (ref) because we have included the grand intercept term $ \tau _{n}(\equiv \log (\zeta _{n}))$ in ((ref)). The low-rank structure of $\Theta _{l}^{\ast }$ incorporates two special cases: (1) additive structure and (2) latent community structure, as illustrated in detail in Examples (ref) and (ref) below, respectively. When there are no covariates in regression and $\Theta _{0}$ belongs to the two cases in Examples (ref) and (ref), the model becomes the so-called Beta model and stochastic block model, respectively. We extend these models to the scenario with edge-wise characteristics and latent community structure for the slope coefficients.

For classification and inference, we need to impose the latent community structure as in Example (ref). This is summarized in the following assumption.

Two remarks are in order. First, Assumption (ref) implies that if $ \Theta_l^{\ast}$ has a latent community structure, the size of each community must be proportional to the number of nodes $n$. Such an assumption is common in the literature on network community detection and panel data latent structure detection. Second, it is possible to allow for $ \pi _{l,kn}$ and/or $\sigma _{k,l}$ to vary with $n$. In this case, one just needs to keep track of all these terms in the proofs.

To proceed, we state a lemma that shows Assumption (ref) is a special case of Assumption (ref) and lays down the foundation for our classification procedure Section (ref).

Lemma (ref) implies that, if $\Theta_l^{\ast}$ for some $l \in [p] \cup \{0\}$ has the community structure, its singular vectors $ \{v_{i,l}\}_{i \in [n]}$ contain information about the community structure. A similar result has been established in the community detection literature; see, e.g., RCY11 and SWZ20.

In Section (ref), we only require $\Theta _{l}^{\ast }$, $l\in \lbrack p]\cup \{0\}$ to be of low-rank and derive the uniform convergence rate of the estimators of $(u_{i,l},v_{i,l})$ across $i\in \lbrack n]$. In Section (ref), we further suppose that some coefficient $\Theta _{l}^{\ast }$ has a special community structure as in Assumption (ref) and apply the K-means algorithm to exactly recover their group identities. Last, for inference in Section (ref), we impose that all coefficients $\Theta _{l}^{\ast }$, $l\in \lbrack p]$ have (potentially different) community structures while $\Theta _{0}^{\ast }$ follows the structure in either Example (ref) or (ref).

The Estimation Algorithm

For notational simplicity, we will focus on the case of $p=1$. The general case with multiple covariates involves fundamentally no new ideas but more complicated notations.

First, we recognize that $\Gamma _{0}^{\ast }$ and $\Gamma _{1}^{\ast }$ are both low rank matrices with ranks bounded from above by $K_{0}+1$ and $K_{1}$ , respectively. We can obtain their preliminary estimates via the nuclear norm penalized logistic regression. Second, based on the normalization imposed in Assumption (ref).1, we can estimate $\tau _{n}$ and $ \Theta _{0}^{\ast }$ separately. We then apply the SVD to the preliminary estimates of $\Theta _{0}^{\ast }$ and $\Theta _{1}^{\ast }$ and obtain the estimates of $U_{l}$, $\Sigma _{l}$, and $V_{l}$, $l=0,1$. Third, we plug back the second step estimates of $\{V_{l}\}_{l=0,1}$ and re-estimate each row of $U_{l}$ by a row-wise logistic regression. We can further iterate this procedure and estimate $U_{l}$ and $V_{l}$ alternatively. Last, if we further impose $\Theta_1^*$ has a community structure, then we can apply the K-means algorithm to the final estimate of $V_{1}$ to recover the community memberships. We rely on a sample splitting technique along with the estimation. Throughout, we assume the ranks $K_{0}$ and $K_{1}$ are known. We will propose an singular-value-ratio-based criterion to select them in Section (ref).

Below is an overview of the multi-step estimation procedure that we propose.

Several remarks are in order. First, $\widehat{\tau }_{n}$ and $\{\widehat{ \Sigma }_{l}\}_{l=0,1}$ obtained in Step 1 are used in Steps 3-5 and to determine $\left\{ K_{l}\right\} _{l=0,1}$ in Section (ref) . Second, we employ the sample-splitting technique to create independence between the edges used for Steps 2 and 3. As $\widehat{V}_{l}^{(1)}$ in Step 2 is estimated by the nuclear-norm regularized logistic regression, we can only control the estimation error in Frobenius norm, as shown in Theorem (ref). On the other hand, to analyze the row-wise estimator, we need to control for the\ estimation error of $\widehat{V}_{l}^{(1)}$ in row-wise $L_{2}$ norm (denoted as $||\cdot ||_{2\rightarrow \infty }$). We overcome the discrepancy between $||\cdot ||_{F}$ and $||\cdot ||_{2\rightarrow \infty }$ by the independence structure. Third, one may propose to use each row of the full-sample lower-rank estimator $\widehat{V} _{l}$ as $\{\dot{v}_{j,l}^{0}\}_{j\in \lbrack n]}$, the initial estimates in Step 4. However, as $\widehat{V}_{l}$ is estimated using the full sample, it is not independent of, say, the $i$-th row of the edges if we want to estimate $(u_{i,0}^{\top },u_{i,1}^{\top })$. Fourth, in the literature, researchers overcome this difficulty by using the \textquotedblleft leave-one-out\textquotedblright\ technique. See, for example, abbe2017 , B13, JM15, SWZ20, and Z18, among others. Denote $\widehat{\Theta }_{l}^{(i)}$ as the low-rank estimator of $\Theta _{l}^{\ast }$ using all the edges except those on the $i$-th row and column and $\widehat{V}_{l}^{(i)}$ is obtained by applying the SVD on $\widehat{ \Theta }_{l}^{(i)}$. The key step for the \textquotedblleft leave-one-out\textquotedblright\ technique is to establish a perturbation theory to bound $\widehat{\Theta }_{l}^{(i)}-\widehat{\Theta }_{l}$, and thus, $\widehat{V}_{l}^{(i)}-\widehat{V}_{l}$. However, unlike the community detection literature, $\widehat{\Theta }_{l}$ and $\widehat{\Theta } _{l}^{(i)}$ are not directly observed but estimated by the nuclear-norm regularized logistic regression. It is interesting but very challenging, if possible, to establish such a perturbation theory. Fifth, although the sample-splitting can result in information loss, we compensate it in three aspects: (1) we just treat the sample-split estimator $\dot{v}_{j,l}^{(0,1)}$ as an initial value and in Step 4, we update it via an iterative algorithm which uses all the edges; (2) we can switch the roles of $I_{1}$ and $I_{2}$ and obtain $\dot{v}_{j,l}^{(H,1)}$ and $\dot{v}_{j,l}^{(H,2)}$ after $H$ iterations; (3) to mitigate the concern of the randomness caused by a single sample split, in Section (ref), we propose to repeat the sample-splitting $R$ times to obtain $R$ classifications, and select one of them based on the maximum-likelihood principle.

The Estimation of $(u_{i,l},v_{i,l})$

In the estimation of $(u_{i,l},v_{i,l})$ (see Steps 1--5 in the above procedure), we only require that $\Theta _{0}^{\ast }$ and $\Theta _{1}^{\ast }$ be of low-rank.

Full-Sample Low-Rank Estimation

Recall that $\Gamma _{0}^{\ast }=\tau _{n}+\Theta _{0}^{\ast }$ and $\Gamma _{1}^{\ast }=\Theta _{1}^{\ast }$. Let $\Gamma ^{\ast }=(\Gamma _{0}^{\ast },\Gamma _{1}^{\ast })$, $\Lambda \left( u\right) =\frac{1}{1+\exp \left( -u\right) }$ denote the standard logistic probability density function,

denote the conditional logistic log-likelihood function associated with nodes $i$ and $j$, and

We propose to estimate $\Gamma ^{\ast }$ by $\widetilde{\Gamma }=(\widetilde{ \Gamma }_{0},\widetilde{\Gamma }_{1})$ via minimizing the negative logistic log-likelihood function with the nuclear norm regularization:

where $Q_{n}(\Gamma )=\frac{-1}{n(n-1)}\sum_{i,j\in \lbrack n],i\neq j}\ell _{ij}\left( \Gamma _{ij}\right) $ and $\lambda _{n}>0$ is a regularization parameter. As mentioned above, we allow $\zeta _{n}$ to shrink to zero at a rate as slow as $n^{-1}\log n$ so that $\tau _{n}=\log \left( \zeta _{n}\right) $ is slightly smaller than $\log n$ in magnitude. So it is sufficient to consider a parameter space $\mathbb{T}\left( 0,\log n\right) $ that expands at rate-$\log n.$ Later on, we specify $\lambda _{n}=\frac{ C_{\lambda }(\sqrt{\zeta _{n}n}+\sqrt{\log n})}{n(n-1)}$ for some constant tuning parameter $C_{\lambda }$. Throughout the paper, we assume $W_{1,ij}$ has been rescaled so that its standard error is one. Therefore, we do not need to consider different penalty loads for $||\Gamma _{0}||_{\ast }$ and $ ||\Gamma _{1}||_{\ast }$. Many statistical softwares automatically normalize the regressors when estimating a generalized linear model. We recommend this normalization in practice before using our algorithm.

Let $\widetilde{\tau }_{n}=\frac{1}{n(n-1)}\sum_{i\neq j}\widetilde{\Gamma } _{0,ij}$. We will show that $\widetilde{\tau }_{n}$ lies within $c_{\tau } \sqrt{\log n}$-neighborhood of the true value $\tau _{n},$ where $c_{\tau }$ can be made arbitrarily small provided that the expected degree is larger than $C\log n$ for some sufficiently large $C.$\footnote{ Let $\eta _{0n}=\sqrt{\frac{\log n}{n\zeta _{n}}}$ and $\eta _{n}=\eta _{0n}+\eta _{0n}^{2}.$ The proof of Theorem (ref).1 suggests that $\widetilde{\tau }_{n}-\tau _{n}=O_{p}(\eta _{n}\sqrt{\log n} ), $ which is\ $o_{p}(\sqrt{\log n})$ (resp. $o_{p}(1)$) if one assumes that the magnitude $n\zeta _{n}$ of the expected degree is of order higher than $ \log n$ (resp. $(\log n)^{2}$). But we will only assume that $\eta _{0n}\leq C_{F}\leq \frac{1}{4}$ for some sufficiently small constant $C_{F}$ below.} This rate is insufficient and remains to be refined. Given $\widetilde{\tau } _{n}$, we propose to reestimate $\Gamma ^{\ast }$ by $\widehat{\Gamma }=( \widehat{\Gamma }_{0},\widehat{\Gamma }_{1}),$ where

and $C_{M}$ is some constant to be specified later. Note that we now restrict the parameter space to expand at rate-$\sqrt{\log n}$ only.

Let $\widehat{\tau }_{n}=\frac{1}{n(n-1)}\sum_{i\neq j}\widehat{\Gamma } _{0,ij}$. Since $\Theta _{l}^{\ast }=\{\Theta _{l,ij}^{\ast }\}$ are symmetric, we define their preliminary low-rank estimators as $\widehat{ \Theta }_{l}=\{\widehat{\Theta }_{l,ij}\},$ where

$\delta _{l0}=\mathbf{1}\{l=0\},$ $f_{M}(u)=u\cdot \mathbf{1}\{|u|\leq M\}+M\cdot \mathbf{1}\{u>M\}-M\cdot \mathbf{1}\{u<-M\}$ is the round function, and $M$ is some positive constant. For $l=0,1$, we denote the SVD of $n^{-1}\widehat{\Theta }_{l}$ as

where $\widehat{\widetilde{\Sigma }}_{l}=\text{diag}(\widehat{\sigma } _{1,l},...,\widehat{\sigma }_{n,l})$, $\widehat{\sigma }_{1,l}\geq \cdots \geq \widehat{\sigma }_{n,l}\geq 0$, and both $\widehat{\widetilde{\mathcal{U }}}_{l}$ and $\widehat{\widetilde{\mathcal{V}}}_{l}$ are $n\times n$ unitary matrices. Let $\widehat{\mathcal{V}}_{l}$ consist of the first $K_{l}$ columns of $\widehat{\widetilde{\mathcal{V}}}_{l}$, such that $(\widehat{ \mathcal{V}}_{l})^{\top }\widehat{\mathcal{V}}_{l}=I_{K_{l}}$ and $\widehat{ \Sigma }_{l}=\text{diag}(\widehat{\sigma }_{1,l},\cdots ,\widehat{\sigma } _{K_{l},l})$. Then $\widehat{V}_{l}=\sqrt{n}\widehat{\mathcal{V}}_{l}.$

Split-Sample Low-Rank Estimation

We divide the $n$ nodes into two roughly equal-sized subsets $(I_{1},I_{2})$ . Let $n_{\ell }=\#I_{\ell }$ denote the cardinality of the set $I_{\ell }.$ If $n$ is even, one can simply set $n_{\ell }=n/2$ for $\ell =1,2.$

Now, we only use the pair of observations $(i,j)\in I_{1}\times [ n]$ to conduct the low-rank estimation. Let $\Gamma _{l}^{\ast }(I_{1})$ consist of\ the $i$-th row of $\Gamma _{l}^{\ast }$ for $i\in I_{1}$, $l=0,1.$ Let $ \Gamma ^{\ast }(I_{1})=(\Gamma _{0}^{\ast }(I_{1}),\Gamma _{1}^{\ast }(I_{1}))$. Define

We estimate $\Gamma ^{\ast }(I_{1})$ via the following nuclear-norm regularized estimation

where $Q_{n}^{(1)}(\Gamma )=\frac{-1}{n_{1}(n-1)}\sum_{i\in I_{1},j\in [ n],i\neq j}\ell _{ij}\left( \Gamma _{ij}\right) $, $\lambda _{n}^{(1)}=\frac{ C_{\lambda }(\sqrt{\zeta _{n}n}+\sqrt{\log n})}{n_{1}(n-1)}$, and the superscript $(1)$ means we use the first subsample ($I_1$) in this step.

Let $\widetilde{\tau }_{n}^{\left( 1\right) }=\frac{1}{n_{1}(n-1)}\sum_{i\in I_{1},j\in [ n],i\neq j}\widetilde{\Gamma }_{0,ij}^{(1)}$. As above, this estimate lies within $c_{\tau }\sqrt{\log n}$-neighborhood of the true value $\tau _{n}.$ To refine it, we can reestimate $\Gamma ^{\ast }(I_{1})$ by $ \widehat{\Gamma }^{(1)}=(\widehat{\Gamma }_{0}^{(1)},\widehat{\Gamma } _{1j}^{(1)}):$

Let $\widehat{\tau }_{n}^{(1)}=\frac{1}{n_{1}(n-1)}\sum_{i\in I_{1},j\in [ n],i\neq j}\widehat{\Gamma }_{0,ij}^{(1)}$. Noting that $\{\Gamma _{l}^{\ast }\}_{l=0,1}$ are symmetric, we define the preliminary low-rank estimates for the $n_{1}\times n$ matrices $\Theta _{l}^{\ast }(I_{1})$ by $\widehat{ \Theta }_{l}^{(1)}$ for $l=0,1$, where

and $\delta _{l0},$ $f_{M}(u)$ and $M$ are defined in\ Step 1. For $l=0,1$, we denote the SVD of $n^{-1}\widehat{\Theta }_{l}^{(1)}$ as

where $\widehat{\widetilde{\Sigma }}_{l}^{(1)}$ is a rectangular ($ n_{1}\times n$) diagonal matrix with $\widehat{\sigma }_{i,l}^{(1)}$ appearing in the $\left( i,i\right) $th position and zeros elsewhere, $ \widehat{\sigma }_{1,l}^{(1)}\geq \cdots \geq \widehat{\sigma } _{n_{1},l}^{(1)}\geq 0$, and $\widehat{\tilde{\mathcal{U}}}_{l}^{(1)}$ and $ \widehat{\tilde{\mathcal{V}}}_{l}^{(1)}$ are $n_{1}\times n_{1}$ and $ n\times n$ unitary matrices, respectively. Let $\widehat{\mathcal{V}} _{l}^{(1)}$ consist of the first $K_{l}$ columns of $\widehat{\widetilde{ \mathcal{V}}}_{l}^{(1)}$ such that $(\widehat{\mathcal{V}}_{l}^{(1)})^{\top } \widehat{\mathcal{V}}_{l}^{(1)}=I_{K_{l}}.$ Let $\widehat{\Sigma }_{l}^{(1)}= \text{diag}(\widehat{\sigma }_{1,l}^{(1)},\cdots ,\widehat{\sigma } _{K_{l},l}^{(1)})$. Then $\widehat{V}_{l}^{(1)}=\sqrt{n}\widehat{\mathcal{V}} _{l}^{(1)},$ and $(\widehat{v}_{j,l}^{(1)})^{\top }$ is the $j$-th row of $ \widehat{V}_{l}^{(1)}$ for $j\in \left[ n\right] $.

Split-Sample Row- and Column-Wise Logistic Regressions

We note that $\Theta_{l,ij}^{\ast} = u_{i,l}^{\top}v_{j,l}$ for $i \in I_2$ and $j \in [n]$. For the $i$-th row when $i \in I_2$, we can view $ \{v_{j,l}\}_{j \in [n]}$ and $u_{i,l}$ as regressors and the parameter, respectively, and estimate $u_{i,l}$ by the row-wise logistic regression. Although $\{v_{j,l}\}_{j \in [n]}$ are unobservable, we can replace them by their estimators obtained from the previous step.

Let $\mu =(\mu _{0}^{\top },\mu _{1}^{\top })^{\top }$ and $\Lambda _{ij}^{ \text{left}}(\mu )=\Lambda (\widehat{\tau }_{n}+\sum_{l=0}^{1}\mu _{l}^{\top }\widehat{v}_{j,l}^{(1)}W_{l,ij})$ and $\ell _{ij}^{\text{left}}\left( \mu \right) =Y_{ij}\log (\Lambda _{ij}^{\text{left}}(\mu ))$ $+(1-Y_{ij})\log (1-\Lambda _{ij}^{\text{left}}(\mu ))$, where the superscript ”left" means these functions are used to estimate the left singular vector $u_{i,l}$. Given the preliminary estimate $\{\widehat{v}_{j,l}^{(1)}\}$ obtained in Step 2, we can estimate the left singular vectors $\{u_{i,0},$ $u_{i,1}\}$ for each $i\in I_{2}$ by $\{\widehat{u}_{i,0}^{(1)},\widehat{u} _{i,1}^{(1)}\} $ via the row-wise logistic regression:

where $Q_{in,U}^{(0)}(\mu )=\frac{-1}{n_{2}}\sum_{j\in I_{2},j\neq i}\ell _{ij}^{\text{left}}\left( \mu \right)$ and the superscript $(0)$ means it is the initial step for the full sample iteration below. To keep the independence between $\{\widehat{v}_{j,l}^{(1)}\}_{j\in [n]}$ and the data in this regression, we only use $j \in I_2$ to run the regression.

Let $\nu =(\nu _{0}^{\top },\nu _{1}^{\top })^{\top }$ and $\Lambda _{ij}^{ \text{right}}(\nu )=\Lambda (\widehat{\tau }_{n}+\sum_{l=0}^{1}\nu _{l}^{\top }\widehat{u}_{i,l}^{(1)}W_{l,ij})$ and $\ell _{ij}^{\text{right} }\left( \nu \right) =Y_{ij}\log (\Lambda _{ij}^{\text{right}}(\nu))$ $ +(1-Y_{ij})\log (1-\Lambda _{ij}^{\text{right}}(\nu ))$, where the superscript ”right" means the functions are used to estimate the right singular vector $v_{j,l}$. Given $(\widehat{u}_{i,0}^{(1)},\widehat{u} _{i,1}^{(1)}),$ we update the estimate of the right singular vectors $ \{v_{i,0},v_{i,1}\}$ for each $j\in [ n]$ by $\{\dot{v}_{j,0}^{(0,1)},\dot{v} _{j,1}^{(0,1)}\}$ via the column-wise logistic regression:

where $Q_{jn,V}^{(0)}(\nu )=\frac{-1}{n_{2}}\sum_{i\in I_{2},i\neq j}\ell _{ij}^{\text{right}}\left( \nu \right) .$

Our final objective is to obtain accurate estimates of $\left\{ v_{j,l}\right\} _{j\in \left[ n\right] ,l=0,1}.$ To this end, we treat $\{ \dot{v}_{j,0}^{(0,1)},\dot{v}_{j,1}^{(0,1)}\}_{j\in \left[ n\right] }$ as the initial estimate in the following full-sample iteration procedure.

Full-Sample Iteration

Given the initial estimates, we use the full sample and iteratively run row- and column-wise logistic regressions to estimate $\{u_{i,l},v_{i,l}\}_{i \in [n]}$. For $h=1,2,...,H,$ let $\Lambda _{ij}^{\text{left,}h}(\mu )=\Lambda ( \widehat{\tau }_{n}+\sum_{l=0}^{1}\mu _{l}^{\top }\dot{v} _{j,l}^{(h-1,1)}W_{l,ij}))$ and $\ell _{ij}^{\text{left,}h}\left( \mu \right) =Y_{ij}\log (\Lambda _{ij}^{\text{left,}h}(\mu ))$ $+(1-Y_{ij})\log (1-\Lambda _{ij}^{\text{left,}h}(\mu )).$ Given $\{\dot{v}_{i,0}^{(h-1,1)}, \dot{v}_{i,1}^{(h-1,1)}\}$, we can compute $\{\dot{u}_{i,0}^{(h,1)},\dot{u} _{i,1}^{(h,1)}\}$ via the row-wise logistic regression

where\ $Q_{in,U}^{(h)}(\mu )=\frac{-1}{n}\sum_{j\in [ n],j\neq i}\ell _{ij}^{ \text{left,}h}\left( \mu \right) .$

Given $\{\dot{u}_{i,0}^{(h,1)},\dot{u}_{i,1}^{(h,1)}\}$, by letting $\Lambda _{ij}^{\text{right,}h}(\nu )=\Lambda (\widehat{\tau }_{n}+\sum_{l=0}^{1}\nu _{l}^{\top }\dot{u}_{i,l}^{(h,1)}W_{l,ij}))$ and $\ell _{ij}^{\text{right,} h}\left( \nu \right) =Y_{ij}\log (\Lambda _{ij}^{\text{right,}h}(\nu ))$ $ +(1-Y_{ij})\log (1-\Lambda _{ij}^{\text{right,}h}(\nu )),$ we compute $\{ \dot{v}_{j,0}^{(h,1)},\dot{v}_{j,1}^{(h,1)}\}$ via the column-wise logistic regression

where $Q_{jn,V}^{(h)}(\nu )=\frac{-1}{n}\sum_{i\in [ n],i\neq j}\ell _{ij}^{ \text{right,}h}\left( \nu \right) .$

We can stop iteration when certain convergence criterion is met for sufficiently large $H.$ Switching the roles of $I_{1}$ and $I_{2}$ and repeating the procedure in the last three steps, we can obtain the iterative estimates $\{\dot{u}_{i,0}^{(h,2)},\dot{u}_{i,1}^{(h,2)}\}_{i\in [ n]}$ and $ \{\dot{v}_{j,0}^{(h,2)},\dot{v}_{j,1}^{(h,2)}\}_{j\in [ n]}$ for $ h=1,2,\cdots ,H$.

K-means Classification

In this step, we further assume $\Theta _{1}^{\ast }$ has the latent community structure and $\Theta _{0}^{\ast }$ remains to be of low-rank. Recall that $\overline{v}_{j,1}=\left( \frac{(\dot{v}_{j,1}^{(H,1)})^{\top } }{||\dot{v}_{j,1}^{(H,1)}||},\frac{(\dot{v}_{j,1}^{(H,2)})^{\top }}{||\dot{v} _{j,1}^{(H,2)}||}\right) ^{\top }$, a $2K_{1}\times 1$ vector. We now apply the K-means algorithm to $\{\overline{v}_{j,1}\}_{j\in \lbrack n]}$. Let $ \mathcal{B}=\{\beta _{1},\ldots ,\beta _{K_{1}}\}$ be a set of $K_{1}$ arbitrary $2K_{1}\times 1$ vectors: $\beta _{1},\ldots ,\beta _{K_{1}}$. Define

and $\widehat{\mathcal{B}}_{n}=\{\widehat{\beta }_{1},\ldots ,\widehat{\beta }_{K_{1}}\}$, where $\widehat{\mathcal{B}}_{n}=\operatorname*{arg\,min}_{\mathcal{B}}\widehat{ Q}_{n}(\mathcal{B}).$ For each $j\in \lbrack n],$ we estimate the group identity by

where if there are multiple $k$'s that achieve the minimum, $\hat{g}_{j}$ takes value of the smallest one.

As mentioned previously, we can repeat Steps 2--6 $R$ times to obtain $R$ membership estimates, denoted as $\{\hat{g}_{j,r}\}_{{j\in \lbrack n],r\in \lbrack R]}}$. Recall that

which is a $K_{1}(K_{1}+1)/2$-vector. In addition, let $\chi _{1,ij}$ be the vectorization of the upper triangular part of the $K_{1}\times K_{1}$ matrix whose $(g_{i}^{0},g_{j}^{0})$ and $(g_{j}^{0},g_{i}^{0})$ entries are one and the rest entries are zero, i.e., $\chi _{1,ij}$ is a $K_{1}(K_{1}+1)/2$ vector such that the $((g_{i}^{0}\vee g_{j}^{0}-1)(g_{i}^{0}\vee g_{j}^{0})/2+g_{i}^{0}\wedge g_{j}^{0})$-th element is one and the rest are zeros, where $g_{i}^{0}\in \lbrack K_{1}]$ denotes the true group membership of the $i$-th node in $\Theta _{1}^{\ast }$. By construction,

Analogously, for the $r$-th split, denote $\hat{\chi}_{1r,ij}$ as a $ K_{1}(K_{1}+1)/2$ vector such that the $((\hat{g}_{i,r}\vee \hat{g}_{j,r}-1)( \hat{g}_{i,r}\vee \hat{g}_{j,r})/2+\hat{g}_{i,r}\wedge \hat{g}_{j,r})$-th element is one and the rest are zeros. We then estimate $B_{1}^{\ast }$ by $ \widehat{B}_{1,r}$, a symmetric matrix constructed from\ $\widehat{b}_{r}$ by reversing the vech operator:

where $\mathcal{L}_{n,r}(b)=\sum_{i<j}[Y_{ij}\log (\hat{\Lambda} _{ij}(b))+(1-Y_{ij})\log (1-\hat{\Lambda}_{ij}(b)))]$ with $\hat{\Lambda} _{ij}(b)=\Lambda (\widehat{\tau }_{n}+\widehat{\Theta }_{0,ij}+W_{1,ij}\hat{ \chi}_{1r,ij}^{\prime }b),$ $\widehat{\tau }_{n}$ is obtained in Step 1, $ \widehat{\Theta }_{0,ij}=[(\dot{u}_{i,0}^{(H,1)})^{\top }\dot{v} _{j,0}^{(H,1)}+(\dot{u}_{i,0}^{(H,2)})^{\top }\dot{v}_{j,0}^{(H,2)}]/2$, and $(\dot{u}_{i,0}^{(H,1)},\dot{v}_{j,0}^{(H,1)},$ $\dot{u}_{i,0}^{(H,2)},\dot{v }_{j,0}^{(H,2)})$ are obtained in Step 5.\footnote{ If we have multiple covariates $W_{l}$, $l\in \lbrack p]$, to compute $ \mathcal{L}_{n,r}(b)$, we let $\widehat{\Theta }_{l,ij}=[(\dot{u} _{i,l}^{(H,1)})^{\top }\dot{v}_{j,l}^{(H,1)}+(\dot{u}_{i,l}^{(H,2)})^{\top } \dot{v}_{j,l}^{(H,2)}]/2$ when $\Theta _{l}^{\ast }$ is only assumed to be of low-rank. For those $\Theta _{l}^{*}$'s that have latent communities, for the $r$-th split, we can estimate their memberships by $\hat{g}_{i,l,r}$ and construct $\hat{\chi}_{lr,ij}^{\prime }$ similarly. Then, we can define $ \mathcal{L}_{n,r}(b)$ and $\widehat{\mathcal{L}}(r)$ in the same manner.} Then, the likelihood of the $r$-th split is defined as $\widehat{\mathcal{L}} (r)=\mathcal{L}_{n,r}(\widehat{b}_{r}).$ Our final estimator $\{\hat{g} _{i,r^{\ast }}\}_{i\in \lbrack n]}$ of the membership corresponds to the $ r^{\ast }$-th split, where

Statistical Properties of the Estimators of $(u_{i,l},v_{j,l})$

In this section, we study the asymptotic properties of the estimators of $ (u_{i,l},v_{j,l})$ proposed in the last section.

Full- and Split-Sample Low-Rank Estimations

Suppose the singular value decomposition of $\Gamma_l^*$ is $\Gamma_l^* = \overline{U}_{l} \Sigma_l \overline{V}_l^{\top}$ for $l = 0,1$ and $ \overline{U}_{l,c}$ and $\overline{V}_{l,c}$ are the left and right singular matrices corresponding to the zero singular values. Let $\mathcal{P} _l(\Delta) = \overline{U}_{l,c}\overline{U}_{l,c}^{\top} \Delta \overline{V} _{l,c}\overline{V}_{l,c}^{\top}$ for some $n\times n$ matrix $\Delta$ and $ \mathcal{M}_j(\Delta) = \Delta - \mathcal{P}_j(\Delta)$. Define the restricted low-rank set as, for some $c_1>0$

Several remarks are in order. First, Assumption (ref) is a slight generalization of CHLZ18 where, in terms of our notation, $\mathcal{C}_{1}(c_{2})$ and $\mathcal{C}_{2}(c_{3})$ take the forms:

Such a generalization is due to the fact that the network can be semi-dense, and thus, the convergence rates of our estimators of the singular vectors are slower than those of CHLZ18's estimators.

Second, suppose there are two sets of parameters $(\Gamma _{0}^{\ast },\Gamma _{1}^{\ast })$ and $(\Gamma _{0}^{\dagger },\Gamma _{1}^{\dagger })$ such that $\Gamma _{l}^{\ast }\neq \Gamma _{l}^{\dagger }$ for some $l\in \{0,1\}$. The singular value decomposition of $\Gamma _{l}^{\dagger }$ is $ \Gamma _{l}^{\dagger }=\tilde{U}_{l}\tilde{\Sigma}_{l}\tilde{V}_{l}^{\top }$ for $l=0,1$ and $\tilde{U}_{l,c}$ and $\tilde{V}_{l,c}$ are the left and right singular matrices corresponding to the zero singular values. Denote $ \widetilde{\mathcal{P}}_{l}(\Delta )=\tilde{U}_{l,c}\tilde{U} _{l,c}^{T}\Delta \tilde{V}_{l_{c}}\tilde{V}_{l,c}^{T}$ and $\widetilde{ \mathcal{M}}_{l}(\Delta )=\Delta -\widetilde{\mathcal{P}}_{l}(\Delta )$. Suppose that

Assumption (ref) holds for both $\mathcal{C}(c_{1})$ and $\tilde{ \mathcal{C}}(c_{1})$, and

Denote $\Delta _{l}=\Gamma _{l}^{\dagger }-\Gamma _{l}^{\ast }$, $l=0,1$. Then $\Delta _{0}+\Delta _{1}\odot W_{1}=0$ and it is possible to show that that $(\Delta _{0},\Delta _{1})$ belongs to either $\mathcal{C}(1)$ or $ \tilde{\mathcal{C}}(1)$.\footnote{ Without loss of generality, we assume that $||\Gamma _{0}^{\dagger }||_{\ast }+||\Gamma _{1}^{\dagger }||_{\ast }\leq ||\Gamma _{0}^{\ast }||_{\ast }+||\Gamma _{1}^{\ast }||_{\ast }$. Noting that

where the last equality holds due to CHLZ18, we have

which implies

i.e., $(\Delta _{0},\Delta _{1})\in \mathcal{C}(1)$.} If $(\Delta _{0},\Delta _{1})\notin \mathcal{C}_{1}(c_{2})$, then Assumption (ref) implies

and thus,

Therefore,

For any estimator $\hat{\Gamma}_{l}$ of $\Gamma _{l}^{\ast }$, $l=0,1$, we have, w.p.a.1,

Based on Assumption (ref) and other conditions in the paper, we can show that (see Theorem (ref) below)

where $C_{F,1}$ is some constant. This implies

and vise versa. The same conclusion holds if $(\Delta _{0},\Delta _{1})\in \mathcal{C}_{1}(c_{2})$. As a result, the ambiguity between $\Theta _{l}^{\ast }$ and $\tilde{\Theta}_{l}$ is asymptotically negligible and will not affect the convergence rates of their estimators.

Third, CHLZ18 provide a sufficient condition for Assumption (ref). Recall $W_{1,ij}=g_{1}(X_{i},X_{j},e_{ij})$. Following the same arguments in CHLZ18, it is possible to show that Assumption (ref) holds if $W_{1,ij}$ is bounded and $ Var(W_{1,ij}|X_{i},X_{j})>0$.\footnote{ In the general case with multiple covariates, they require $ \min_{i,j}\lambda _{\min }(\mathbb{E}W_{ij}W_{ij}^{\top }|X_{i},X_{j})\geq c>0$ where $\lambda _{\min }(A)$ is the minimum eigenvalue of matrix $A$ and $W_{ij}=(1,W_{1,ij},\cdots ,W_{p,ij})^{\top }$.} The sufficient condition basically requires the existence of $e_{ij}$ in $g_{l}(\cdot )$ which is a sequence of i.i.d. random variables across $i,j$.\footnote{ Note there are two key differences between the setups in our paper and CHLZ18. First, CHLZ18 consider the panel data with indexes $i\in \lbrack N]$ and $t\in \lbrack T]$ while we consider the network data with indexes $(i,j)\in \{1\leq i<j\leq n\}$. Second, CHLZ18 consider $ X_{it}=\mu _{it}+e_{it}$ such that given $\{\mu _{it}\}_{i\in \lbrack N],t\in \lbrack T]}$, $X_{it}$ is independent across both $t$ and $t$. Instead, we consider $W_{1,ij}=g_{1}(X_{i},X_{j},e_{ij})$ such that given $ \{X_{i}\}_{i\in \lbrack n]}$, $W_{1,ij}$ is independent across $1\leq i<j\leq n$. By examining the proofs of CHLZ18, we note that their argument does not rely on the special structure of $ X_{it}=\mu _{it}+e_{it}$ and works if $X_{it}=f(\mu _{it},e_{it})$ for some non-additive function $f$.} Note that the presence of $e_{ij}$ is sufficient, but may not be necessary. In our simulation, we generate $ W_{1,ij}=|X_{i}-X_{j}|$ with $\{X_{i}\}_{i\in \lbrack n]}$ being a sequence of i.i.d. standard normal random variables, and find that our method works well.

Fourth, Assumption (ref) rules out the case $ W_{1,ij}=g_{1}(X_{i},X_{j})$ when $X_{i}$ is discrete, which is equivalent to a community structure of $W_{1,ij}$. Suppose $W_{1,ij}=w_{k_{1}k_{2}}>0$ $ \forall k_{1},k_{2}$ where $i,j$ are in groups $k_{1}$ and $k_{2}$. Then, we can let $\Delta _{1}$ share the same community structure as $W_{1}$ and $ \Delta _{1,ij}=w_{k_{1}k_{2}}^{-1}$. Let $\Delta _{0}=-\iota _{n}\iota _{n}^{\top }$. Then we have

Because both $\Delta _{1}$ and $\Delta _{0}$ are of low-rank, we have

for some constant $C>0$. In addition, the singular value decomposition of $ \Delta _{0}$ is $\Delta _{0}=(-\iota _{n}/\sqrt{n})\times n\times (\iota _{n}/\sqrt{n})^{\top }$. It is possible to find some parameter $\Theta _{0}$ such that $||\mathcal{M}_{0}(\Delta _{0})||_{\ast }\geq cn$ for some $c>0$. \footnote{ This occurs, say, when $\iota _{n}/\sqrt{n}$ is in the spaces spanned by the left and right singular vectors of $\Theta _{0}$ that correspond to its nonzero singular values.} Then we can take $c_{1}=C/c$ so that

In this case, Assumption (ref) does not hold because $||\Delta _{0}||_{F}^{2}+||\Delta _{1}||_{F}^{2}\geq n^{2}>c_{2}\log (n)n/\zeta _{n}$ \footnote{ When $\zeta _{n}=C_{\varsigma }n^{-1}\log (n),$ we require that $ C_{\varsigma }$ be sufficiently large.} and

Assumption (ref) is a regularity condition. In particular, Assumptions (ref).2 implies the order of the average degree in the network is $n\zeta _{n}$. Assumption (ref).3 means that the average degree diverges to infinity at a rate that is not slower than $\log n$. Such a rate is the slowest for exact recovery in the SBM, as established by ABH2016, AS15, MNS14, and V18. As our model incorporates the SBM as a special case, the rate is also the minimal requirement for the exact recovery of $Z_{1}$, which is established in Theorem (ref) below. Assumption (ref).4 usually holds as the sample is split randomly and $\Theta _{0}^{\ast }$ satisfies the normalization condition in Assumption (ref).1. If $\Theta _{0}^{\ast }$ satisfies the additive structure as in Example (ref), then Assumption (ref).4 provided that $\frac{1}{n_{1}}\sum_{i\in I_{1}}\alpha _{i}=o(\sqrt{\frac{\log (n)}{n\zeta _{n}}})$. Such a requirement holds almost surely if $\alpha _{i}=a_{i}-\frac{1}{n}\sum_{i\in \lbrack n]}a_{i}$ and $\{a_{i}\}_{i=1}^{n}$ is a sequence of i.i.d. random variables with finite second moments. If $\Theta _{0}^{\ast }$ has the community structure as in Example (ref), then Assumption (ref).4 holds provided that $p_{0}^{\top }(I_{1})B_{0}^{\ast }p_{0}=o( \sqrt{\frac{\log (n)}{n\zeta _{n}}})$, where $p_{0}^{\top }(I_{1})=(\frac{ n_{1,0}(I_{1})}{n_{1}},\cdots ,\frac{n_{K_{0},0}(I_{1})}{n_{1}})$ and $ n_{k,0}(I_{1})$ denotes the size of $\Theta _{0}^{\ast }$'s $k$-th community for the subsample of nodes with index $i\in I_{1}$. As $p_{0}^{\top }B_{0}^{\ast }p_{0}=0$, the requirement holds almost surely if community memberships are generated from a multinomial distribution so that $ ||p_{0}-p_{0}(I_{1})||=o_{a.s.}(\sqrt{\frac{\log (n)}{n\zeta _{n}}})$.

Part 1 of Theorem (ref) indicates that despite the possible divergence of the grand intercept $\tau _{n},$ we can estimate it consistently up to rate $\eta _{n}.$ In the dense network, $\zeta _{n}\asymp 1$ where $a\asymp b$ denotes both $a/b$ and $b/a$ are stochastically bounded. In this case, $\tau _{n}$ $\asymp 1$ and it can be estimated consistently at rate-$\sqrt{(\log n)/n}.$ Note that the convergence rate of $ \widehat{\Theta }_{l}$ and $\widehat{\Theta }_{l}^{(1)}$ in terms of the Frobenius norm is also driven by $\eta _{n}.$ Similarly for $\widehat{\sigma }_{k,l}$ $\widehat{\sigma }_{k,l}^{(1)},$ $\widehat{V}_{l}/\sqrt{n}$ and $ \widehat{V}_{l}^{(1)}/\sqrt{n}.$ In part 4 of Theorem (ref) , the orthogonal matrices $\widehat{O}_l$ and $\widehat{O}_l^{(1)}$ are present because the singular values of $\Theta_l^*$ can be the same and its singular vectors can only be identified up to some rotation.

Split-Sample Row- and Column-Wise Logistic Regressions

Define two $\left( K_{0}+K_{1}\right) \times \left( K_{0}+K_{1}\right) $ matrices:

To study the asymptotic properties of the third step estimator, we assume that both matrices are well behaved uniformly in $i$ and $j$ in the following assumption.

Assumption (ref) assumes that $\Phi _{i}(I_{2})$ and $\Psi _{j}(I_{2})$ are positive definite (p.d.) uniformly in $i$ and $j$ asymptotically. Suppose $\Gamma _{1}$ follows the community structure as in Example (ref) with $K_{1}$ equal-sized communities and $ B_{1}^{\ast }=I_{K_{1}}$, then $\Pi _{1,n}=\text{diag}(1/K_{1},\cdots ,1/K_{1})$. By Lemma (ref) in the online supplement, if node $j$ is in community $k$, then $v_{j,1}=\sqrt{n}\sqrt{\frac{K_{1}}{n}}z_{j,1}=\sqrt{ K_{1}}e_{K_{1},k}$, where $e_{K_{1},k}$ denotes a $K_{1}\times 1$ vector with the $k$-th unit being 1 and all other units being 0. In addition, suppose $\Theta _{0}$ follows the specification in Example (ref). Then,

Suppose that $\alpha _{0,i}=a_{i}-\bar{a}$ for some i.i.d. sequence $ \{a_{i}\}_{i=1}^{n}$ with $\bar{a}=\frac{1}{n}\sum_{i=1}^{n}a_{i}$, and the group identities of $\Theta _{1}^{\ast }$ ($\{z_{i}\}_{i\in \lbrack n]}$) are independent of $\Theta _{0}^{\ast }$ and $\{X_{i}\}_{i\in \lbrack n]}$ and $\{e_{ij}\}_{i,j\in \lbrack n]}$. Further suppose $\mathbb{E} (W_{1,ij}a_{j}|X_{i})=0$, $\mathbb{E}(W_{1,ij}|X_{i})=0$, and $\mathbb{E} (W_{1,ij}^{2}|X_{i})\geq c>0$ for some constant $c$. Then, we can expect that, uniformly over $i\in I_{2}$,

which implies Assumption (ref) holds.

If $\Theta_0^{\ast}$ has the community structure as in Example (ref). Further suppose $\Theta _{0}^{\ast }$ and $\Theta _{1}^{\ast }$ share the same community structure $Z_{1}$, which is independent of $W_{1}$, $\mathbb{E}(W_{1,ij}|X_{i})=0$ and $\mathbb{E}(W_{1,ij}^{2}|X_{i})\geq c>0$ for some constant $c$, then one can expect that $\Phi _{i}(I_{2})$ has the same limit as above uniformly over $i\in I_{2}$.

The following theorem studies the asymptotic properties of $\widehat{u} _{i,l}^{(1)}$ and $\dot{v}_{j,l}^{(0,1)}$ defined in Step 3.

Theorem (ref) establishes the uniform bound for the estimation error of $\dot{v}_{j,l}^{(0,1)}$ up to some rotation. However, we only use half of the edges to estimate $\dot{v}_{j,l}^{(0,1)}$, which may result in information loss. In the next section, we treat $\dot{v} _{j,l}^{(0,1)}$ as an initial value and iteratively re-estimate $ \{u_{i,l}\}_{i\in [ n]}$ and $\{v_{j,l}\}_{i\in [ n]}$ using all the edges in the network. We will show that the iteration can preserve the error bound established in Theorem (ref).

Full-Sample Iteration

Define two $\left( K_{0}+K_{1}\right) \times \left( K_{0}+K_{1}\right) $ matrices:

To study the asymptotic properties of the fourth step estimators, we add an assumption.

The above assumption parallels Assumption (ref) and is now imposed for the full sample.

Theorem (ref) establishes the uniform bound for the estimation error in the iterated estimators $\{\dot{u}_{i,l}^{(h,1)}\}$ and $\{\dot{v} _{i,l}^{(h,1)}\}.$

By switching the roles of $I_{1}$ and $I_{2}$, we have, similar to Theorem (ref), that

where $\widehat{O}_{l}^{(2)}$ is a $K_{l}\times K_{l}$ rotation matrix that depends on $V_{l}$ and $\widehat{V}_{l}^{(2)}$. Then, following the same derivations of Theorems (ref) and (ref), we have, for $h=1,\cdots ,H$,

K-means Classification

If we further assume $\Theta_{1}^{\ast}$ has the community structure and satisfies Assumption (ref), then Lemma (ref) shows $ \{v_{j,1}\}_{j\in [ n]}$ contains information about the community memberships. It is intuitive to expect that we can use $\overline{v}_{j,l}$ defined in Section (ref) to recover the memberships as long as the estimation error is sufficiently small.

Let $g_{i}^{0}\in \left[ K_{1}\right] $ denote the true group identity for the $i$-th node in $\Theta_1^{\ast}$. To establish the strong consistency of the membership estimator $\hat{g}_{i}$ defined in (ref), we add the following side condition.

Apparently, Assumption (ref) is automatically satisfied in large samples if $\eta _{n}=o\left( 1\right) .$ The constant in the statement is not optimal.

Several remarks are in order. First, Theorem (ref) implies the K-means algorithm can exactly recover the latent community structure of $ \Theta _{1}^{\ast }$ $w.p.a.1$. Second, if we repeat the sample split $R$ times, we need to maintain Assumption (ref) for each split. Then, we can show the exact recovery of $\hat{g}_{{i,r}}$ for $r\in [ R]$ in the exact same manner, as long as $R$ is fixed. This implies $\hat{g}_{{i,r\ast } }$ for $r^{\ast }$ selected in (ref) also enjoys the property that

Third, if $\Theta _{0}^{\ast }$ also has the latent community structure as in Example (ref), we can apply the same K-means algorithm to $\{ \overline{v}_{j,0}\}_{j\in \left[ n\right] }$ with $\overline{v}_{j,0}\equiv (\dot{v}_{j,0}^{(H,1)\top }/||\dot{v}_{j,0}^{(H,1)}||,\dot{v} _{j,0}^{(H,2)\top }/||\dot{v}_{j,0}^{(H,2)}||)^{\top }$ to recover the group identities of $\Theta _{0}^{\ast }$. Last, if we further assume $ Z_{0}=Z_{1}=Z$ (which implies $K_{0}=K_{1})$, then we can catenate $ \overline{v}_{j,0}$ and $\overline{v}_{j,1}$ as a $4K_{1}\times 1$ vector and apply the same K-means algorithm to this vector to recover the group membership for each node.

Inference for $B_{1}^{\ast }$

In this section, we maintain the assumption that $\Theta _{1}^{\ast }$ has a latent community structure. In the general model with multiple covariates, we allow $\{\Theta _{l}^{\ast }\}_{l\in \lbrack p]}$ to have potentially different community structures. Note this includes the case that some of the $\Theta _{l}^{\ast }$'s are homogeneous. We can recover the community structures by applying the k-means algorithm in the previous section to each $\Theta _{l}^{\ast }$.

For the rest of the section, for notation simplicity, we continue to consider the case that there is only one covariate $W_1$ and $ \Theta_1^{\ast} $ has a latent community structure, which is estimated by $\{ \hat{g}_i\}_{i \in [n]}$ defined in the previous section. Given the exact recovery of the community memberships asymptotically, we can just treat $ \hat{g}_{i}$ as $g_{i}^{0}$.

We discuss the inference for $B_{1}^{\ast }$ for two specifications of $ \Theta_0^{\ast}$: (1) $\Theta_{0,ij}^{\ast}$ has an additive structure as in Example (ref) and (2) $\Theta_{0,ij}^{\ast}$ has a latent community structure as in Example (ref). In the first model, once the group membership of $\Theta_1^{\ast}$ is recovered, it boils to the one studied by G17. For the second model, when the memberships of both $ \Theta_0^{\ast}$ and $\Theta_1^{\ast}$ are recovered, it boils down to the standard logistic regression with finite-number of parameters.

Additive Fixed Effects

Suppose $\Gamma _{0,ij}^{\ast }=\tau _{n}+\alpha _{i}+\alpha _{j} $ and $ \Gamma _{1}^{\ast }=\Theta _{1}^{\ast }=Z_{1}B_{1}^{\ast }Z_{1}^{\top }$. Recall the definitions of $\chi _{1,ij},$ $\hat{\chi}_{1r,ij},$ and $\text{ vech}(B_{1}^{\ast })$ in Section (ref) such that $\chi _{1,ij}^{\top }$vech$(B_{1}^{\ast })=B_{1,g_{i}^{0}g_{j}^{0}}^{\ast }.$ We further denote $\hat{\chi}_{1,ij}$ as either $\hat{\chi}_{1,ij}$ if one single split is used or $\hat{\chi}_{1r^{\ast },ij}$ if $R$ splits are used and the $r^{\ast }$-th split is selected.

Corollary (ref) directly follows from Theorem (ref) and implies that we can treat $\chi _{1,ij}$ as observed. Then, (ref) can be written as

where $\omega _{1,ij}=W_{1,ij}\chi _{1,ij}$. This model has already been studied by G17. We can directly apply his Tetrad logit regression to estimate $\text{vec}(B_{1}^{\ast })$.

Let $S_{ij,i^{\prime }j^{\prime }}=Y_{ij}Y_{i^{\prime }j^{\prime }}(1-Y_{ii^{\prime }})(1-Y_{jj^{\prime }})-(1-Y_{ij})(1-Y_{i^{\prime }j^{\prime }})Y_{ii^{\prime }}Y_{jj^{\prime }}.$ Then, for an arbitrary $ K_{1}(K_{1}+1)/2$-vector $B$, the conditional likelihood of $S_{ij,i^{\prime }j^{\prime }}$ given $S_{ij,i^{\prime }j^{\prime }}\in \{-1,1\}$ is

where $\widetilde{\omega }_{1,ij,i^{\prime }j^{\prime }}=\omega _{1,ij}+\omega _{1,i^{\prime }j^{\prime }}-(\omega _{1,ii^{\prime }}+\omega _{1,jj^{\prime }})$. Further denote

Following G17, we define the tetrad regression estimator $\widehat{B}$ for $\text{vech}(B^{\ast })$ as

Let

be the indicator that the tetrad $\{i,j,i^{\prime },j^{\prime }\}$ take an identifying configuration, and thus, contributes to the tetrad logit regression. Further denote $t_{q,n}=\mathbb{P}(\top _{i_{1}i_{2}i_{3}i_{4}}=1,\top _{jj_{2}j_{3}j_{4}}=1)$ as the probability that tetrads $\{i_{1},i_{2},i_{3},i_{4}\}$ and $\{j,j_{2},j_{3},j_{4}\}$ both take an identifying configuration when sharing $q=0,1,2,3$, or 4 nodes in common. Then, we make the following assumption on the Hessian matrix.

The following theorem reports the asymptotic normality of $\widehat{B}.$

Theorem (ref) imposes two additional structures in order to make the inferences on $B^{\ast }$ by borrowing the asymptotic results from G17. One is that $\Gamma _{0}^{\ast }$ exhibits the usual additive fixed effects structure (with $K_{0}=2$) and the other is $\Gamma _{1}^{\ast }$ has a latent community structure. The model reduces to that of G17 in the special case of $K_{1}=1.$

Latent Community Structure in the Fixed Effects

Let $g_{i,0}^{0}$ be the true memberships of node $i$ for $\Theta _{0}^{\ast }$ and $\hat{g}_{i,0}$ be its estimator which can be computed by applying the K-means algorithm to $\{\overline{v}_{j,0}\}_{j\in \lbrack n]}$. Further note $Z_{0}\iota _{K_{0}}=\iota _{n}$ where recall that $\iota _{b}$ denotes a $b\times 1$ vector of ones. Therefore, $\Gamma _{0}^{\ast }=\tau _{n}\iota _{n}\iota _{n}^{\top }+Z_{0}B_{0}^{\ast }Z_{0}^{\top }=Z_{0}(B_{0}^{\ast }+\tau _{n}\iota _{K_{0}}\iota _{K_{0}}^{\top })Z_{0}^{\top }\equiv Z_{0}B_{0}^{\ast \ast }Z_{0}^{\top }$, i.e., $\Gamma _{0}^{\ast }$ shares the same community structure as $\Theta _{0}^{\ast }$. We then define $\chi _{0,ij}$ be a $K_{0}(K_{0}+1)/2\times 1$ vector whose $((g_{i,0}^{0}\vee g_{j,0}^{0}-1)(g_{i,0}^{0}\vee g_{j,0}^{0})/2+g_{i,0}^{0}\wedge g_{j,0}^{0})$ -th element is one and the rest are zeros and $\hat{\chi}_{0,ij}$ be a $ K_{0}(K_{0}+1)/2\times 1$ vector whose $((\hat{g}_{i,0}\vee \hat{g}_{j,0}-1)( \hat{g}_{i,0}\vee \hat{g}_{j,0})/2+\hat{g}_{i,0}\wedge \hat{g}_{j,0})$-th element is one and the rest are zeros. Similar to Corollary (ref), we have the following corollary.

We propose to estimate $\text{vech}(B^{\ast })\equiv (\text{vech} (B_{0}^{\ast \ast })^{\top },\text{vech}(B_{1}^{\ast })^{\top })^{\top }$ by

where

and

Let $\Lambda _{n,ij}(u)=\Lambda (\omega _{ij}^{\top }[$vech$(B^{\ast })+u(n^{2}\zeta _{n})^{-1/2}])$ and $\Lambda _{n,ij}\equiv \Lambda _{n,ij}(0),$ where\ $\omega _{ij}=(\chi _{0,ij}^{\top },\chi _{1,ij}^{\top }W_{1,ij})^{\top }$ is an\ $\mathcal{K}$-vector with $\mathcal{K} =\sum_{l=0}^{1}K_{l}(K_{l}+1)/2.$ Note that $\Lambda _{n,ij}=\Lambda (W_{ij}^{\top }\Gamma _{ij}^{\ast }).$

Although in theory, the inference for $B_{1}^{\ast }$ in the above two cases is straightforward, there are two finite-sample issues. First, the tetrad logistic regression does not scale with the number of nodes $n$ because the algorithm scans over all four-nodes figurations, which contains a total of $ O(n^{4})$ operations in a brutal force implementation. Although the Python code by G17 incorporates a number of computational speed-ups by keeping careful track of non-contributing configurations as the estimation proceeds, we still find in our simulations that the implementation turns extremely hard for networks with over 1000 nodes. One can, instead, use subsampling or divide-and-conquer algorithm for estimation. To establish the theoretical properties of such an estimator is an important and interesting topic for future research.

Second, for the specification in the second example, based on unreported simulation results, we find that $\widehat{B}_{1}$ has a small bias if there are some misclassified nodes. However, as the standard error of our estimator is even smaller, such a small bias may not be ignored in making inferences. If we further increase the sample size, then the classification indeed achieves exact recovery and such a bias vanishes quickly. However, in practice, researchers cannot know whether their sample size is sufficiently large. It is interesting to further investigate such a bias issue and make proper bias-corrections. This is, again, left as a topic for future research.

Determination of $K_{0}$ and $K_{1}$

In practice, $K_{0}$ and $K_{1}$ are unknown and need to be estimated from the data. In this case, we propose to replace them by a large but fixed integer $K_{\max }$ in the first step estimation to obtain the singular value estimates $\left\{ \hat{\sigma}_{k,l}\right\} _{k\in \left[ K_{\max } \right] ,l=0,1}.$ We propose a version of singular-value ratio (SVR) statistic in the spirit of the eigenvalue-ratio statistics of AH2013 and LY2012. That is, for $l=0,1,$ we estimate $K_{l}$ by

where $\bar{Y}=\frac{2}{n(n-1)}\sum_{1\leq i<j\leq n}Y_{ij},$ and $c_{l}$ is a tuning parameter to be specified. Without the indicator function in the above definition, $\widehat{K}_{l}$ is nothing but the SVR statistic. The use of the indicator function helps to avoid the overestimation of the ranks. Apparently, $n\bar{Y}$ consistently estimate the expected degree that is of order $n\zeta _{n}.$ By using Assumption (ref) and the results in Theorem (ref), we can readily establish the consistency of $\hat{K}_{l}$.

Monte Carlo Simulations

In this section, we conduct some simulations to evaluate the performance of our procedure.

Data generation mechanisms

We generate data from the following two models.

Model 1. We simulate the responses $Y_{ij}$ from the Bernoulli distribution with mean $\Lambda (\log (\zeta _{n})+\Theta _{0,ij}^{\ast }+W_{1,ij}\Theta _{1,ij}^{\ast })$ for $i<j$, where $\Theta _{0,ij}^{\ast }=\alpha _{i}+\alpha _{j}$ and $\Theta _{1}^{\ast }=ZB_{1}^{\ast }Z^{\top }$ . We generate $\alpha _{i}\overset{i.i.d}{\sim }$ $\mathcal{U}(-1/2,1/2)$ for $i=1,...,n$, and $W_{1,ij}=|X_{i}-X_{j}|$ for $i\neq j$, where $X_{i}$ $ \overset{i.i.d}{\sim }\mathcal{N}(0,1)$. For the $i^{\text{th}}$ row of the membership matrix $Z\in \mathbb{R}^{n\times K_{1}}$, the $C_{i}^{\text{th}}$ component is $1$ and other entries are $0$, where $C=(C_{1},...,C_{n})^{\top }\in \mathbb{R}^{n}$ is the membership vector with $C_{i}\in \left[ K_{1} \right] $.

Case 1. Let $K_{1}=2$ and $B_{1}^{\ast }=((0.6,0.2)^{\top },(0.2,0.7)^{\top })^{\top }$. The membership vector $C=(C_{1},...,C_{n})^{ \top }$ is generated by sampling each entry independently from $\{1,2\}$ with probabilities $\{0.4,0.6\}$. Let $\zeta _{n}=0.7n^{-1/2}\log n$.

Case 2. Let $K_{1}=3$ and $B_{1}^{\ast }=((0.8,0.4,0.3)^{\top },(0.4,0.7,0.4)^{\top },(0.3,0.4,0.8)^{\top })^{\top }$. The membership vector $C=(C_{1},...,C_{n})^{\top }$ is generated by sampling each entry independently from $\{1,2,3\}$ with probabilities $\{0.3,0.3,0.4\}$. Let $ \zeta _{n}=1.5n^{-1/2}\log n$.

Model 2. We simulate the responses $Y_{ij}$ from the Bernoulli distribution with mean $\Lambda (\log (\zeta _{n})+\Theta _{0,ij}^{\ast }+W_{1,ij}\Theta _{1,ij}^{\ast })$ for $i<j$, where $\Theta _{0}^{\ast }=ZB_{0}^{\ast }Z^{\top }$, $\Theta _{1}^{\ast }=ZB_{1}^{\ast }Z^{\top }$, and $W_{1,ij}$ is simulated in the same way as in Model 1. Note here we impose that the latent community structures for $\Theta _{0}^{\ast }$ and $ \Theta _{1}^{\ast }$ are the same. We then apply the K-means algorithm to the $4K_{1}\times 1$ vector $\{\overline{v}_{j,0}^{\top },\overline{v} _{j,1}^{\top }\}_{j\in [ n]}$ to recover the community membership, as described in Section (ref).

Case 1. Let $K_{0}=K_{1}=2$ and $B_{0}^{\ast }=((0.6,0.2)^{\top },(0.2,0.7)^{\top })^{\top }$, $B_{1}^{\ast }=((0.6,0.2)^{\top },(0.2,0.5)^{\top })^{\top }$. The membership vector $C=(C_{1},...,C_{n})^{ \top }$ is generated by sampling each entry independently from $\{1,2\}$ with probabilities $\{0.3,0.7\}$. Let $\zeta _{n}=0.5n^{-1/2}\log n$.

Case 2. Let $K_{0}=K_{1}=3$ and $B_{0}^{\ast }=((0.7,0.2,0.2)^{\top},(0.2,0.6,0.2)^{\top },(0.2,0.2,0.7)^{\top })^{\top }$ , $B_{1}^{\ast }=((0.7,0.3,0.2)^{\top },(0.3,0.7,0.2)^{\top },(0.2,0.2,0.6)^{\top })^{\top }$. The membership vector is generated in the same way as given in Case 2 of Model 1. Let $\zeta _{n}=1.5n^{-1/2}\log n$.

We consider $n=500,$ $1000,$ and $1500$. All simulation results are based on 200 realizations.

Simulation Results

We select the number of communities $K_{1}$ by an eigenvalue ratio method given as follows. Let $\widehat{\sigma }_{1,1}\geq \cdots \geq \widehat{ \sigma }_{K_{\max },1}$ be the first $K_{\max }$ singular values of the SVD decomposition of $\widehat{\Theta }_{1}$ from the nuclear norm penalization method given in Section (ref). We estimate $K_{1}$ by $\widehat{K} _{1}$ defined in ((ref)) by setting $c_{1}=0.1$ and $K_{\max }=10$. We set the tuning parameter $\lambda _{n}=C_{\lambda }\{\sqrt{n\overline{Y}}+ \sqrt{\log n}\}/\{n(n-1)\}$ with $C_{\lambda }=2$ and similarly for $\lambda _{n}^{\left( 1\right) }$. To require that the estimator of $\widehat{\Theta } _{l,ij}$ is bounded by finite constants, we let $M=2$ and $C_{M}=2$. The performance of the method is not sensitive to the choice of these finite constants. Define the mean squared error (MSE) of the nuclear norm estimator $\widehat{\Theta }_{l}$ for $\Theta _{l}$ as $\sum\nolimits_{i\neq j}( \widehat{\Theta }_{l,ij}-\Theta _{l,ij}^{\ast })^{2}/\{n(n-1)\}$ for $l=0,1$.

Table (ref) reports the MSEs for $\widehat{\Theta } _{l}$, the mean of $\widehat{K}_{1}$ and the percentage of correctly estimating $K_{1}$ based on the 200 realizations. We observe that the mean value of $\widehat{K} _{1}$ gets closer to the true number of communities $K_{1}$ and, the percentage of correctly estimating $K_{1}$ approaches to 1, as the samples size $n$ increases. When $n$ is large enough ($n=1500$), the mean value of $ \widehat{K}_{1}$ is the same as $K_{1}$ and the percentage of correctly estimating $K$ is exactly equal to 1.

Next, we use three commonly used criteria for evaluating the accuracy of membership estimation for our proposed method. These criteria include the Normalized Mutual Information (NMI), the Rand Index (RI) and the proportion (PROP) of nodes whose memberships are correctly identified. They all give a value between 0 and 1, where 1 means a perfect membership estimation. Table (ref) presents the mean of the NMI, RI and PROP values based on the 200 realizations for Models 1 and 2. The values of NMI, RI and PROP increase to 1 as the sample size increases for all cases. These results demonstrate that our method is quite effective for membership estimation in both models, and corroborate our large-sample theory.

\ \

Last, we estimate the parameters $B_{0}^{\ast }$ and $B_{1}^{\ast }$ by our proposed method given in Section (ref) for Model 2. Tables (ref) and (ref) show the empirical coverage rate (coverage) of the $95\%$ confidence intervals, the absolute value of bias (bias), the empirical standard deviation (emp_sd), and the average value of the estimated asymptotic standard deviation (asym_sd) of the estimates for $ B_{0}^{\ast }$ and $B_{1}^{\ast }$ in cases 1 and 2 of model 2, respectively, based on 200 realizations. We observe that the emp_sd and asym_sd decrease and the empirical coverage rate gets close to the nominal level $0.95$, as the sample size increases. Moreover, the value of emp_sd is similar to that of asym_sd for each parameter. This result confirms our established formula (in the online supplement) for the asymptotic variances of the estimators for the parameters. When the sample size is large enough $ (n=1500) $, the value of bias is very small compared to asym_sd, so that it can be negligible for constructing confidence intervals of the parameters.

Empirical applications

In this section, we apply the proposed method to study the community structure of social network datasets

Pokec social network

The dataset and model

Pokec is a popular on-line social network in Slovakia. The whole dataset has more than 1.6 million users, and it can be downloaded from https://snap.stanford.edu/data/soc-Pokec.html. In this social network, nodes are anonymized users of Pokec and edges represent friendships. Moreover, demographical features of the users are provided, including gender, age, hobbies, interest, education, etc. To illustrate our method, we select the first 10000 users. Each user is a node in the graph. After deleting the nodes with missing values in age and with degree less than 10, we have $1745$ nodes in our dataset. We use the continuous variable, age, as the covariate in our model, and use the friendship network to create an undirected adjacency matrix which has $1745$ nodes and $39650$ edges. The average degree in this dataset is $22.72$. The left panel of Figure (ref) shows the number of nodes in different age groups. We see that the age group of 25-29 is the largest group with 1175 users and the age groups of 20-24 and 30-34 have similar number of users. Around 98.8% of users are between the ages of 20 and 35 years old. Moreover, in the right panel of Figure (ref), we depict the boxplots of degrees (the number of users connected to each user) for the four age groups 20-24, 25-29, 30-34 and 35-39 that include most users. The plots of degrees vary across different age groups, indicating that age may play a role in the prediction of connections between users.

\

We consider fitting the model:

for $i=1,...,1745$, where $Y_{ij}$ is the observed value ($0$ or $1$) of the adjacency matrix in our dataset, and $W_{1,ij}=|X_{i}-X_{j}|/(\sqrt{ X_{i}^{2}+X_{j}^{2}})$, in which $X_{i}$ is the normalized age of the $i^{ \text{th}}$ customer.\footnote{ The variable $W_{1,ij}$ takes 1444 distinctive values. Given there are only 1745 nodes in our dataset, we can view $W_{ij}$ as continuous.} In this model, $(\tau _{n},\Theta _{0,ij}^{\ast },\Theta _{1,ij}^{\ast })$ are unknown parameters, and $\Theta _{0,ij}^{\ast }$ and $\Theta _{1,ij}^{\ast }$ have the latent group structures $\Theta _{0}^{\ast }=ZB_{0}^{\ast }Z^{\top } $ and $\Theta _{1}^{\ast }=ZB_{1}^{\ast }Z^{\top }$, respectively. Model ( (ref)) considered for this real application is similar to Model 2 in the simulation, and it allows for not only the main effect but also possible interaction effects of age and the latent community structure.

Estimation results

We first use the singular-value ratio method to obtain the estimated number of groups for $\Theta _{0}^{\ast }$ and $\Theta _{1}^{\ast }$: $\widehat{K} _{0}=2$ and $\widehat{K}_{1}=2$, i.e., we identify two subgroups in the friendship network.

Next, we use our proposed method to obtain the estimated membership for each node. As a result, we have identified $842$ nodes in one community and $903$ nodes in the other community. We reorganize the observed adjacency matrix according to the estimated memberships of the nodes, i.e., the nodes in the same estimated community are put together in the adjacency matrix. We use blue dots to represent the edges between nodes. The left panel of Figure (ref) displays the reorganized adjacency. We see that nodes within each community are generally more densely connected than nodes between communities. In the right panel of Figure (ref), we show the boxplots of age for the two identified subgroups. We can observe that in general, the values of age in group 1 are smaller than those in group 2.

Last, Table (ref) shows the estimates of $B_{0}^{\ast }$ and $B_{1}^{\ast }$ and their standard errors (s.e.). We obtain the p-value$ <0.01 $ for testing each coefficient in $B_{1}^{\ast }$ equal to zero, indicating that the covariate age has a significant effect on the prediction of the friendships between users.

Facebook friendship network

The dataset and model

The dataset contains Facebook friendship networks at one hundred American colleges and universities at a single point in time. It was provided and analyzed by TMP12, and can be downloaded from https://archive.org/details/oxford-2005-facebook-matrix. TMP12 used the dataset to illustrate the relative importance of different characteristics of individuals across different institutions, and showed that gender, dormitory residence and class year may play a role in network partitions by using assortativity coefficients. We, therefore, use these three user attributes as the covariates $X_{i}=(X_{i1},X_{i2},X_{i3})^{\top } $, where $X_{i1}=$binary indicator for gender, $X_{i2}=$multi-category variable for dorm number (e.g., \textquotedblleft 202\textquotedblright , \textquotedblleft 203\textquotedblright , etc.), and $X_{i3}=$integer valued variable for class year (e.g., \textquotedblleft 2004\textquotedblright , \textquotedblleft 2005\textquotedblright , etc.). We use the dataset of Rice University to identify the latent community structure interacted with the covariates by our proposed method. \

We use the dataset to fit the model:

where $Y_{ij}$ is the observed value ($0$ or $1$) of the adjacency matrix in the dataset, and $W_{1,ij}=\{\sum\nolimits_{k=1}^{3}(2D_{ij,k}/\Delta _{k})^{2}\}^{1/2}$, where $\Delta _{k}=\max (D_{ij,k})-\min (D_{ij,k})$ and $ D_{ij,k}=X_{ik}-X_{jk}$ for $k=1,2,3$.\footnote{ We note that $W_{1,ij}$ takes 1512 distinctive values. Given there are just 3073 nodes in the dataset, we can view $W_{1,ij}$ as continuous.} In this model, $(\tau _{n},\Theta _{0,ij}^{\ast },\Theta _{1,ij}^{\ast })$ are unknown parameters, and $\Theta _{0,ij}^{\ast }$ and $\Theta _{1,ij}^{\ast }$ have the latent group structures $\Theta _{0}^{\ast }=ZB_{0}^{\ast }Z^{\top } $ and $\Theta _{1}^{\ast }=ZB_{1}^{\ast }Z^{\top }$, respectively. Following model 2 in the simulation, we impose that $\Theta _{0}^{\ast }$ and $\Theta _{1}^{\ast } $ share the same community structure. It is worth noting that RAM19 fit a similar regression model as ((ref)) but let the coefficient of the pairwise covariate be an unknown constant with respect to $(i,j)$ such that $\Theta _{1,ij}^{\ast }=\Theta _{1}^{\ast }$. Although RAM19's RAM19 model can take into account the covariate effect for community detection, it does not consider possible interaction effects of the observed covariates and the latent community structure. As a result, it may cause the number of estimated groups to be inflated. In the dataset of Rice University, we delete the nodes with missing values and with degree less than 10, and consider the class year from 2004 to 2009. After the cleanup, there are $n=3073$ nodes and 279916 edges in the dataset for our analysis.

Estimation results

We first use the eigenvalue ratio method to obtain the estimated number of groups for $\Theta^{\ast } _{0}$ and $\Theta^{\ast } _{1}$: $\widehat{K} _{0}=4$ and $\widehat{K}_{1}=4.$

Next, we use our proposed method to obtain the estimated membership for each node. Table (ref) presents the number of students in each estimated group for female and male, for different class years, and for different dorm numbers. It is interesting to observe that most female students belong to either group 2 or group 4, and most male students belong to either group 1 or group 3. There is a clear community division between female and male; within each gender category, the students are further separated into two large groups. Moreover, most students in the class years of 2004 and 2005 are in either group 1 or group 2, while most students in the class years of 2008 and 2009 are in either group 3 or group 4. Students in the class years of 2006 and 2007 are almost evenly distributed across the four groups, with a tendency that more students will join groups 3 and group 4 when they are in later class years. This result indicates that students tend to be in different groups as the gap between their class years becomes larger. Last, Table (ref) shows the estimates of $B^{\ast }_{0}$ and $B^{\ast }_{1} $ and their standard errors (s.e.). We obtain the p-value$ <0.01$ for testing each coefficient in $B^{\ast }_{1}$ equal to zero, indicating that the three covariates are useful for identifying the community structure.

Conclusion

In this paper, we proposed a network formation model which can capture heterogeneous effects of homophily via a latent community structure. When the expected degree diverges at a rate no slower than rate-$\log n$, we established that the proposed method can exactly recover the latent community memberships almost surely. By treating the estimated community memberships as the truth, we can then estimate the regression coefficients in the model by existing methods in the literature.