Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
96,708 characters · 23 sections · 58 citation commands
Spectral estimation of large stochastic blockmodels with discrete nodal covariates
\onehalfspacing
The analysis and modeling of network data has far-reaching applications in economics, sociology, public health, computer science, neuroscience, and marketing, among other areas. Social networks have been shown to affect socioeconomic performance, such as education CalvoArmengolEtAl2009, DeGiorgiPellizzariRedaelli2009, CarrellEtAl2013, health and risky behaviors Nakajima2007, Badev2013, risk sharing arrangements FafchampsGubert2006, and employment opportunities Topa2001, Beaman2012, among others. Often, both observed and unobserved factors contribute to the global structure of networks and the processes that generate them. For example, in social networks, factors including gender, race, and personality affect the likelihood that two people interact. In applications, race and gender are typically observed, while personality is usually unobserved. It is therefore crucial to develop ways to disentangle the effect of observed and unobserved variables on link formation in networks.
The analysis of network models poses several econometric challenges, arising from the structure of correlation among links Chandrasekhar2016, DePaula2017, Graham2020, GrahamDePaula2020, Dzemski2017. In this paper, we focus on models with conditionally independent links, where the unobserved heterogeneity is modeled by latent positions in low-dimensional Euclidean space Graham2014, Auerbach2019. While this formulation rules out externalities and strategic considerations in link formation, which are dominant focal points in much of the econometric literature on strategic network formation games Mele2017, DePaula2017, Graham2020, GrahamDePaula2020, Menzel2017, there is a growing literature showing how models with conditionally independent links provide useful approximations of strategic models, at least in some special cases Mele2017, MeleZhu2020, DiaconisChatterjee2011, Graham2020.
In this paper, we analyze the stochastic blockmodel (SBM), a workhorse in the literature on community detection and clustering Abbe2018, NowickiSniders2011, AiroldiEtAl2008. In standard $K$-block SBMs, each node belongs to one of $K$ unobserved blocks (communities); conditional on the block assignments, links form independently as Bernoulli variables with probabilities that depend on the community memberships. Many existing works use SBMs to estimate or approximate unobserved block structure in networks. In contrast, applications involving SBMs that incorporate observed nodal attributes as covariates are comparatively few Sweet2015, ChoiEtAl2011, AtchadeEtAl2019. A possible reason is that estimation for stochastic blockmodels is computationally burdensome, and including covariates in the specification imposes significant additional challenges to modelling, estimation, and inference. Exact maximum likelihood estimation is infeasible, causing most estimation strategies to rely on approximations based on expectation--maximization algorithms and variational methods AiroldiEtAl2008,DaudinEtAl2008,BickelEtAl2013,LatoucheEtAl2012,Vu2013. However, these algorithms may converge slowly to the (approximate) solution and become impractical for networks with thousands of nodes.\footnote{Recent advances use further approximations and parallelization to improve computational efficiency AtchadeEtAl2019,Vu2013. We do not pursue such extensions in this paper.}
The goal of our paper is to develop an estimation method for SBMs with observed nodal covariates that is computationally feasible and scales well to large networks, while providing similar statistical guarantees as the other available methods. At a high level, our main strategy consists of formulating this goal as an inference problem in the context of the generalized random dot product graph (GRDPG) model AthreyaEtAl2018a,TangPriebe2018, TangEtAl2017, Rubin-DelanchyEtAl2018, AlidaeeEtAl2020, MuEtAl2021. In a GRDPG, each node is characterized by an unobserved latent position (vector), and each pair of nodes links with probability determined via a (possibly indefinite) inner product of the pair's latent positions; crucially, any SBM can be reformulated as a GRDPG where the latent positions are fixed within blocks. To address our computational goals, we turn to the use of spectral methods, which, in addition to enjoying theoretical guarantees TangEtAl2017, have been shown to be successful both in terms of feasibility and scalability in related settings PriebeFlashGraph, PriebeSSD2016,PriebeBillionNodesGraphs.
We provide several contributions to the literature on network econometrics. First, we present a GRDPG model framework that incorporates the effect of observed covariates on link probabilities. Second, we develop a spectral estimator for inference in stochastic blockmodels with covariates, adapting the spectral estimators developed for our new class of GRDPG models. Crucially, we obtain a new central limit theorem for the spectral estimator of the covariates' effect. Our estimator is asymptotically normal as long as the parameter(s) for covariate effects can be written as sufficiently well-behaved functions of the SBM block-specific probabilities. We provide explicit formulas for bias and variance properties of the estimator, and we show that the estimator is computationally fast, scaling well for large networks. Our method provides a statistical and algorithmic foundation for inference in a broad class of models for large network data, including networks that are relatively sparse in the sense that their average degree scales sub-linearly with network size.
Our exposition focuses on SBMs with a single binary (or discrete) observed covariate, though we emphasize that the theoretical and computational properties set forth in this work extend to settings involving multiple discrete covariates. Asymptotic normality continues to hold as long as the estimator for the effect of the covariate(s) can be expressed as a suitably well-behaved function of the SBM probabilities. We illustrate several examples with two binary covariates in a simulation study and showcase our method in an empirical application using Facebook data and three control variables. The case involving continuous covariates is more complicated and an active area of contemporaneous research. Current progress on this front is being facilitated by recently investigated Latent Structure Models (LSM) AthreyaLSM2018 and related ideas.
The development of our estimator depends crucially on several observations. First, a $K$-block stochastic blockmodel with one binary covariate can be reformulated as a (different yet related) $2K$-block stochastic blockmodel. Second, as discussed previously, a stochastic blockmodel graph can be viewed as a generalized random dot product graph whose latent positions are fixed within blocks AthreyaEtAl2018a,TangPriebe2018,TangEtAl2017. The behavior of our spectral estimation method is tied to the asymptotic behavior of spectral estimators for SBM block probability matrix entries recently studied in TangEtAl2017. Our asymptotic analysis provides explicit formulas for standard errors and establishes the existence of a bias term; however, this bias vanishes at a rate proportional to the size of the network.\footnote{Since we have a closed-form expression for the bias term, in principle we can naively correct for it in estimation, using a plug-in estimate. In our simulations we find that the bias term is usually so small that the correction is not necessary, at least for networks with a few thousand nodes. On the other hand, the bias is demonstrably substantial in the empirical application to Facebook data in Section (ref).}
The theoretical machinery used to perform inference extends methods developed for the analysis of latent positions network models AthreyaEtAl2018a,TangEtAl2017. In particular, we use Adjacency Spectral Embedding (ASE) for random graphs to embed the network in a low-dimensional space and to recover the latent positions of the nodes. Our method is motivated by the (verifiable) intuition that the adjacency matrix can be viewed as a (mild) perturbation of the probability matrix that generates the network data, and thus, that the eigenstructure of the adjacency matrix resembles that of the edge probability matrix TangPriebe2018,AthreyaEtAl2018a. In particular, spectrally decomposing the adjacency matrix provides accurate information about the structure of sufficiently large networks TangEtAl2017.
In addition to providing statistical guarantees, one of the advantages of our method is the speed of computation, obtained without sacrificing estimation accuracy. In our simulations (see Section (ref)) we compare our approach to the variational EM (VEM) algorithm DaudinEtAl2008, BickelEtAl2013, as implemented in the blockmodels package in R. Even for the simplest case of a stochastic blockmodel without covariates, our spectral method is faster by several orders of magnitude. For example, in a network with $n=5000$ nodes and $K=2$ blocks, we can estimate the model in few seconds using our spectral method, while it takes almost 10 minutes to estimate the model using the variational EM algorithm. When we add a binary covariate, our estimator converges in under 30 seconds, while in contrast it takes almost 10 hours when using a parallelized version of the VEM algorithm in blockmodels. Our methods are implemented in the R package grdpg available at \url{https://github.com/meleangelo/grdpg} and all replication files can be found at \url{https://github.com/meleangelo/grdpg_supplement}.
We perform a Monte Carlo study to examine the performance of our spectral estimator. We simulate networks in which the block structure is unbalanced and blocks have different sizes; we include scenarios in which the covariates are binary and independent, as well as cases in which the covariates are correlated. Our results confirm the existence of bias in small samples, but the bias decreases in magnitude when the size of the network increases. As expected, the bias is larger when the blocks are unbalanced and the covariates are correlated. These insights confirm that our estimator works best in very large networks, where the bias problem is attenuated, thus adding to the advantage of computational speed.
Finally, we apply our method to the study of Facebook friendship data using the Facebook 100 dataset, initially collected and analyzed in TraudEtAl2012. These data contain networks of friendships and node (person) information for 100 universities in the United States in the year 2005. We estimate a stochastic blockmodel for the Harvard University network, consisting of more than 13,000 nodes, using information on gender, off-campus residence, and university role (i.e., student, faculty, staff, etc.) of the users (see also AtchadeEtAl2019 for a related analysis). We find evidence of homophily, as suggested by the positive effect of gender, role, and off-campus residence on the probability of linking. This suggests that including information about observable covariates in the estimation may allow researchers to better recover unobservable block structure.
Some of the theoretical machinery presented herein will be useful to econometricians studying networks, as well as nonlinear panel data. Indeed, some of the GRDPG modeling and inference ideas are similar to the literature on interactive fixed effects ChenWeidnerFernandezVal2021, FERNANDEZVAL2016, MoonShumWeidner2018.\footnote{An application of SBM-related ideas is Bonhomme2017.} Another application of our model is to correct for the endogeneity of the network in empirical models of network effects ShaliziMcFowland2018,GoldsmithPinkhamImbens2013, BoucherFortin2016,Auerbach2019, JohnssonMoon2019. Currently most of these studies rely on an auxiliary model of network formation to capture unobserved heterogeneity that affects the outcome. Our model and computational method allow the researcher to perform this type of correction for large network data.
In a $K$-block Stochastic Blockmodel (SBM), nodes are randomly assigned to one of $K$ blocks; conditional on the blocks, nodes form links independently. A $K$-block SBM is characterized by the $K\times K$ matrix of probabilities $\bm{\theta}\in [0,1]^{K\times K}$, where the entry $\bm{\theta}_{k\ell}$ is the probability of a link occurring between nodes in blocks $k$ and $\ell$. The random variables comprising $\bm{\tau}=(\tau_1, \dots, \tau_n)$ describe the assignments of each node to a block, and they are i.i.d., such that the probability that node $i$ belongs to block $k$ is $\mathbb{P}(\tau_i=k)=\pi_k$, with $\bm{\pi} = (\pi_1, \dots, \pi_K)$. Conditional on the assignment to blocks $\bm{\tau}$, the probability that nodes $i$ and $j$ have a link is $\bm{P}_{ij}=\bm{\theta}_{\tau_i \tau_j}$. We use the $n\times n$ adjacency matrix $\bm{A}$ to describe the network, conditioning on the unobserved blocks. According to the SBM the entries of the adjacency matrix are generated as
and we write $ (\bm{A}, \bm{\tau}) \sim SBM(\bm{\theta},\bm{\pi})$ to denote the adjacency matrix drawn from a $K$-block SBM with probability matrix $\bm{\theta}$ and block assignment probabilities $\bm{\pi}$.\\
Conditioning on $\bm{\tau}$, the likelihood of the SBM with $K$ blocks is
The Generalized Random Dot Product Graph (GRPDG) model is an alternative model for network formation with conditionally independent links. In a GRDPG, each node $i$ is characterized by a $d$-dimensional vector (i.e., an unobserved latent position) $\bm{X}_{i}= (X_{i1},\dots,X_{id})\in \mathcal{X}_d\subseteq \mathbb{R}^d$. The latent positions are i.i.d. draws from a distribution $F$ with support $\mathcal{X}_d$, that is $\bm{X}_1, \bm{X}_2, \dots, \bm{X}_n \overset{iid}{\sim} F$. Let $\bm{X}=[\bm{X}_1, \dots, \bm{X}_n]^T$ denote the matrix formed by row-wise stacking all unobserved vectors $\bm{X}_{i}$.
Let $d_1 \ge 1$ and $d_2 \ge 0$ be integers, and define $d = d_1+d_2$. Let $\bm{I}_{d_1, d_2}$ be a $d\times d$ diagonal matrix containing $1$'s in $d_1$ diagonal entries and $-1$ in the remaining $d_2$ diagonal entries. For a GRDPG with signature $(d_1,d_2)$, the entries of the adjacency matrix $\bm{A}_{ij}$ are specified to be independent, after conditioning on the latent positions $\bm{X}_i$ and $\bm{X}_j$, namely
with link probability given by
For this setting, we write $(\bm{X},\bm{A}) \sim GRDPG_{d_1,d_2}(F)$.\footnote{ It must be noted that the support $\mathcal{\bm{X}}_d$ of $F$, is a subset of $\mathbb{R}^d$ such that $\bm{x}^T \bm{I}_{d_1,d_2} \bm{y}\in[0,1]$ for all $\bm{x},\bm{y}\in \mathcal{\bm{X}}_d$. }\\
A notable property of GRDPGs is that they encompass or approximate any conditionally independent low-rank network model. In particular, any SBM can be represented as a GRDPG with latent positions fixed within blocks. That is, the $K$ blocks are represented by a fixed location, so that each $\bm{X}_i$ can only take values $\bm{\nu} = \left[\bm{\nu}_1, \bm{\nu}_2, \dots, \bm{\nu}_K\right]$. Two nodes $i$ and $j$ belong to the same block $k$ if $\bm{X}_i=\bm{X}_j=\bm{\nu}_k$. The random variables $\bm{\tau}$ are such that $\tau_1,\dots,\tau_n \overset{iid}{\sim} Multinomial(1;\pi_1, \dots,\pi_K)$, with $\bm{\pi} \in (0,1)^K$ and $\sum_{k=1}^{K} \pi_k =1$.\footnote{Alternatively, we can think of a $K$-block stochastic blockmodel as a network where the $\bm{X}_i$'s are drawn from a mixture of degenerate distributions with mass centered at $\bm{\nu}$, i.e.,
} The GRDPG corresponding to model $(\bm{A}, \bm{\tau}) \sim SBM(\bm{\theta},\bm{\pi})$ can be obtained by an eigendecomposition of the matrix $\bm{\theta} = \bm{U}\bm{\Sigma} \bm{U}^T$ and by defining $\bm{\nu}_1,\bm{\nu}_2, \dots,\bm{\nu}_K$ as the rows of $\bm{U}\vert \bm{\Sigma}\vert^{1/2}$. The distribution $F$ is $F=\sum_{k=1}^{K}\pi_k \delta_{\bm{\nu}_k}$, where $\delta$ is the Dirac-delta; importantly, $d$ is the rank of the block-probabilities matrix $\bm{\theta}$, and $d_1, d_2$ are the number of positive and negative eigenvalues of matrix $\bm{\theta}$, respectively.\\
This paper utilizes inferential spectral methods for GRDPGs to estimate SBMs AthreyaEtAl2018a,TangEtAl2017. The same relationship between SBMs and GRDPGs holds for known link functions and $\bm{\theta}_{\tau_i\tau_j} = h\left(\bm{B}_{\tau_i\tau_j} \right)$, where $h$ is a known function that maps to $[0,1]$ and $\bm{B}$ is a $K\times K$ matrix of real numbers.\footnote{If $h$ is unknown we cannot in general expect to be able to accurately estimate the latent positions. See tang2013.} For example, $h$ could be the logistic function or the cumulative density function of the Gaussian distribution. Our stochastic blockmodel would have adjacency matrix $\bm{A}$ with elements
The stochastic blockmodel can be extended to include the effect of observed covariates ChoiEtAl2011, Sweet2015, AtchadeEtAl2019. Such models allow researchers to disentangle the effect of observed and unobserved nodal heterogeneity on the probability of linking. In particular, in social science, such models are used to estimate to what extent the network exhibits homophily or heterophily. Let node $i$ be characterized by an $r$-dimensional vector of observed covariates $\bm{Z}_{i}= (Z_{i}^{(1)}, \dots,Z_{i}^{(r)}) \in \mathcal{Z}\subseteq \mathbb{R}^r$ and let the stochastic blockmodel be
where $f$ is a known function, $\bm{\beta}$ is a vector of parameters, and where we allow $\bm{Z}_i$ to (possibly) depend on the latent blocks.\\
In this paper we will focus on the case of a single binary (or discrete) covariate $\bm{Z}_i$, and we will assume that the function $f$ is an indicator variables that indicates whether $i$ and $j$'s covariates have the same value, i.e.,
Here $\beta$ can be interpreted in terms of homophily. Namely, if $\beta>0$, the probability of a link between $i$ and $j$ is higher when their observables $\bm{Z}_{i}$ and $\bm{Z}_{j}$ are the same, and the network displays homophily in the observable variable. Viceversa, when $\beta<0$, we have heterophily. The extension to multiple discrete covariates has similar properties and will be discussed further below.\\
Our goal is to develop a general spectral method of inference for the parameter $\beta$ and for $\bm{B}_{\tau_i \tau_j}$ in the following stochastic blockmodel with a discrete nodal covariate:
We further wish to disentangle the effect of observed and unobserved heterogeneity on link probabilities. To achieve this, we need to extend results from previous work on GRDPGs and SBMs AthreyaEtAl2018a,TangPriebe2018,TangEtAl2017. In the following subsections we review some of the spectral methods we use in the paper, and we provide an example that highlights the core aspects of our method.
Estimation of SBMs for large networks, with or without observed covariates, is computationally challenging. The exact MLE problem is intractable because of the high-dimensional combinatorial problem of considering all possible partitions of the nodes in blocks BickelEtAl2013. Approximate methods are available, based on variational approximations DaudinEtAl2008,AiroldiEtAl2008,WainwrightJordan2008; however, even these methods are computationally prohibitive for large networks.\\ We make use of spectral methods, which have been shown in the literature to scale well with network size. Our spectral approach embeds the network into a low(er) dimensional space, thus reducing the dimensionality of the problem, while maintaining the geometric properties of the data. In particular we use the Adjacency Spectral Embedding (ASE) to estimate the latent positions of the GRDPG AthreyaEtAl2018a. In this sense, our method can be considered a dimension-reduction tool that decreases the complexity of the data by reducing the dimensionality of the space. The intuition about the spectral method is that if $\bm{P}$ is a low-rank matrix, then we can think of the adjacency matrix $\bm{A}$ as a perturbation of $\bm{P}$, that is $\bm{A}_{ij} = \bm{P}_{ij}+ \bm{E}_{ij}$, where $\bm{E}_{ij}$ is a matrix of independent stochastic perturbations.\footnote{In the Bernoulli case, $\bm{E}_{ij}$ is a shifted Bernoulli variable, with values $\bm{E}_{ij}=1-\bm{P}_{ij}$ with probability $\bm{P}_{ij}$ and $\bm{E}_{ij}=\bm{P}_{ij}$ with probability $1-\bm{P}_{ij}$.} If $\bm{A}$ and $\bm{P}$ are close enough, namely if $\bm{E}$ is small enough, then the leading eigenvalues and eigenvectors of $\bm{A}$ and $\bm{P}$ will be similar TangEtAl2017. As a consequence, the spectral decomposition of $\bm{A}$ will provide an estimate of the latent structure of the network, that is, the latent positions $\bm{X}$.
Consider first the case without observed covariates. Let $\bm{P}$ be positive semidefinite and let $h$ be the identity function $h(u) = u$. In this setting, we only have latent positions $\bm{X}$, that are unobserved. If we were able to observe $\bm{P}= \bm{X}\bm{X}^T $, estimation of $\bm{X}$ would be straightforward. Furthermore, we could use spectral embeddings for $\bm{P}$ by exploiting the fact that $\bm{P}$ is positive semidefinite of rank $d$ and has spectral decomposition $\bm{P} = \bm{U}_{\bm{P}} \bm{S}_{\bm{P}} \bm{U}_{\bm{P}}^{T}$, where $\bm{S}_{\bm{P}}$ is a diagonal matrix containing the largest $d$ eigenvalues (in absolute value) of $\bm{P}$ and $\bm{U}_{\bm{P}}$ is the matrix with the corresponding eigenvectors. This implies that a good estimate for $\bm{X}$ is $\widehat{\bm{X}} = \bm{U}_{\bm{P}} \vert \bm{S}_{\bm{P}} \vert^{1/2} $, where $\vert \cdot \vert$ denotes entrywise absolute values. The estimation problem arises because we only observe $\bm{A}$, a perturbed version of $\bm{P}$. The Adjacency Spectral Embedding of $\bm{A}$ into $\mathbb{R}^d$ is then $\widehat{\bm{X}} = \bm{U}_{\bm{A}} \vert \bm{S}_{\bm{A}}\vert^{1/2}$ where $\bm{S}_{\bm{A}}$ is a diagonal matrix containing the largest $d$ eigenvalues of $\bm{A}$ in absolute value and $\bm{U}_{\bm{A}}$ is the matrix with the corresponding eigenvectors.
In our asymptotic results for the above setup, we use the fact that ASE estimates of latent positions $\bm{X}$ asymptotically achieve perfect clustering (moreover, are asymptotically normal) and can be identified up to multiplication by an orthogonal matrix AthreyaEtAl2018a,TangEtAl2017. This implies that asymptotically the blocks are recovered exactly (up to relabeling). The same logic and results hold for non-positive definite matrices $\bm{P}$, allowing us to study more general stochastic blockmodels Rubin-DelanchyEtAl2018.
To illustrate the methodology and to develop intuition, we focus on the special case of a $K=2$ stochastic blockmodel with a single discrete covariate and with latent positions in the unit interval $[0,1]$, yielding $d_{1}=1, d_{2}=0,$ and $r=1$, where $Z_i\in \lbrace 0,1\rbrace$ is a binary variable (e.g., male/female, white/nonwhite, rich/poor, etc.) and the function $f(Z_i,Z_j;\beta)= \beta \mathbf{1}_{\lbrace Z_i = Z_j \rbrace}$ is an indicator for the equality of the covariates for $i$ and $j$, weighted by the parameter $\beta$. The main advantage of this approach is that we can illustrate the geometry of the method in a low-dimensional space. In our simple example, the matrix $\bm{B}$ is given by
where $p,q\in[0,1]$. We can conveniently re-write the matrix $\bm{B}$ as a dot-product of vector $\bm{\nu} = [p \ \ q]^T$, with $p,q\in[0,1]$, that is $\bm{B}=\bm{\nu}\bm{\nu}^T$, so that the SBM can be re-written as a random dot-product graph model with $\bm{X}_i = p$ if $i$ is in block 1, $\bm{X}_i = q$ if $i$ is in block 2. The probability of linking can then be written as
For ease of exposition the network blocks have the same probability, so $(\pi_1, \pi_2) = (0.5,0.5)$ and each community contains half males ($Z_i=1$) and half females ($Z_i=0$). However, we note that our algorithm and the theoretical results are valid when we allow the blocks to be of different size, and the observed covariates to be correlated with the unobserved blocks.
The model specified via ((ref)) corresponds to a 4-block stochastic blockmodel. Indeed, we have 2 unobserved blocks, that are split in two additional blocks by the observed binary variable. Therefore, the final result is a 4-block SBM. More generally, if there are $K$ latent blocks and one binary covariates, we will have a $\tilde{K}=2K$-block SBM.
The possible values of $\bm{X}_i^T \bm{X}_j$ are $\lbrace p^2, pq, q^2 \rbrace$. Therefore the 4-block model can be completely characterized by the $4\times 4$ matrix
The value $h(\bm{B}_{Z,11}) = h(p^2 + \beta)$ is the probability that two males in block 1 form a link; on the other hand, $h(\bm{B}_{Z,12})=h(p^2)$ is the probability that a male and a female in block 1 form a link; $h(\bm{B}_{Z,31})=h(pq+\beta)$ is the probability that two males, one in block 1 and one in block 2, form a link; and so on.
The above observations imply that, for this four block SBM, there exists a corresponding GRDPG with link probability matrix
for some $n\times d$ matrix of latent positions $\bm{Y}$ with $d_{1} \ge 1$, $d_{2} \ge 0$, and $d=d_1+d_2$.
To estimate the parameter $\beta$ and the latent positions $p$ and $q$ we use the following algorithmic approach.
In practice, we can estimate $\beta$ from multiple entries of the matrix $\bm{B}_Z$, for example $\beta = \bm{B}_{Z,11}-\bm{B}_{Z,12} = \bm{B}_{Z,33}-\bm{B}_{Z,34} $, and weight each estimate by the size of the blocks. This could improve the estimate, since some blocks are larger than others, so delivering more precise estimates. Our code implements this idea, which is more practical for empirical applications.
In this section, we derive a central limit theorem for the spectral estimator of $\beta$. For ease of exposition, we focus on the case of a single binary observed covariate and scalar $\beta$, though our method works for other specifications in which the effect of the observed covariates $\beta$ can be written as a function of the stochastic blockmodel's probability matrix $\bm{\theta}_Z$. Extensions to multiple binary or discrete observed covariates are straightforward albeit tedious.
We desire to estimate a stochastic blockmodel with observed covariates, where
We assume that the observed covariates are binary and can depend on the block assignment, that is $\bm{Z}_i\vert \tau_i \overset{ind}{\sim} Bernoulli(b_{\tau_i})$, where $b_{\tau_i}= P(\bm{Z}_i=1\vert \tau_i)$ . Our asymptotic results are easily extended to the case of discrete observed covariates with three or more possible outcomes.
As explained above in the simple example, our strategy consists of rewriting the SBM as a GRDPG. First, notice that the matrix $\bm{B}$ can be written as $\bm{B}_{\tau_i \tau_j} = \bm{X}_i^T \bm{X}_j$, where $\bm{X}_i$ is a $d \times 1$ vector of latent positions that has $K$ possible values $\bm{\nu}_1, \dots, \bm{\nu}_{K}$. In practice, $\bm{\nu}=(\bm{\nu}_1, \dots, \bm{\nu}_{K})$ are the centers of the $K$ blocks $\bm{X}$, such that $i$ and $j$ belong to unobserved block $k$ when $\bm{X}_i=\bm{X}_j=\bm{\nu}_k$. Let $\bm{\tau}$ be the function that assigns nodes to unobserved blocks; then $\tau_i=k$ if $\bm{X}_i=\bm{\nu}_k$. We can thus rewrite the stochastic blockmodel above as a random dot product graph with observed covariates as follows:
We first notice that both models are stochastic blockmodels with $\tilde{K}=2K$ blocks, because the indicator variable $\bm{1}_{\lbrace \bm{Z}_i=\bm{Z}_j\rbrace}$ splits each unobserved block in two blocks. The probabilities of belonging to a block $k$ for this $\widetilde{K}$-block SBM are denoted as $\bm{\eta} = (\eta_1, \dots, \eta_{\widetilde{K}}) = ( \pi_1 \cdot b_1, \pi_1 \cdot (1-b_1), \pi_2 \cdot b_2, \pi_2 \cdot (1-b_2), \dots , \pi_K \cdot b_K, \pi_K \cdot (1-b_K) )$; and the functions that assign nodes to blocks are $\bm{\xi} = (\xi_1,\dots,\xi_{n} )$, such that $\xi_i = 1$ if $ \tau_i = 1$ and $ \bm{Z}_i=0 $; $\xi_i =2$ if $\tau_i=1$ and $ \bm{Z}_i=1$; $\xi_i=3$ if $\tau_i=2$ and $ \bm{Z}_i=0 $; $\xi_i=4$ if $\tau_i=2$ and $ \bm{Z}_i=1$; and so on.\\
So we have a stochastic blockmodel $ (\bm{A}, \bm{\xi}, \bm{Z}) \sim SBM(\bm{\theta}_Z,\bm{\eta})$ with the $\widetilde{K}\times \widetilde{K}$ matrix of probabilities $\bm{\theta}_Z$
The stochastic blockmodel characterized by the matrix $\bm{\theta}_Z$ can be re-formulated as a GRDPG. Indeed, consider the eigendecomposition $\bm{\theta}_Z \equiv \bm{U}\bm{\Sigma} \bm{U}^T$, and define $\bm{\mu}=[\bm{\mu}_1,\bm{\mu}_2, \dots,\bm{\mu}_{\widetilde{K}}]$ as the rows of $\bm{U}\vert \bm{\Sigma}\vert^{1/2}$; then let $F=\sum_{k=1}^{\widetilde{K}}\eta_k \delta_{\bm{\mu}_k}$, where $\delta$ is the Dirac-delta; and $d_1$ and $d_2$ are the number of positive and negative eigenvalues of $\bm{\theta}_Z$, respectively. Then, the Generalized Random Dot Product Graph model $(\bm{Y},\bm{A}) \sim GRDPG_{d_1,d_2}(F)$ corresponding to our stochastic blockmodel $ (\bm{A}, \bm{\xi}, \bm{Z}) \sim SBM(\bm{\theta}_Z,\bm{\eta})$ is given by
where $d_1+d_2=\widetilde{d}=rank(\bm{\theta}_Z)$ and $\bm{Y}$ is the $n\times \widetilde{d}$ vector of latent positions with centroids $\bm{\mu}$.\\
We can now extend asymptotic results for estimation of RDPGs in AthreyaEtAl2018a,TangEtAl2017 to estimate block assignments and the effect of the covariates (see Rubin-DelanchyEtAl2018 for the corresponding generalization to GRDPGs). \\
Because the functions $\bm{\tau}$ that describe the assignments to blocks are unknown, the $\widetilde{K}$ SBM model assignment functions $\bm{\xi}$ are also unknown. Applying the Adjacency Spectral Embedding procedure, we recover an estimate $\widehat{\bm{\xi}}$.
We prove asymptotic normality for the parameter $\beta$, exploiting the fact that $\beta$ can be written as a function of the SBM probabilities, that is
If the blocks were known at the onset, we could use the estimator $\widehat{\beta} = h^{-1}(\widehat{\bm{\theta}}_{Z,11} )- h^{-1}(\widehat{\bm{\theta}}_{Z,12} )$. However, all that we have access to is the estimate $\widehat{\bm{\xi}}$, so it is crucial that this estimate be consistent. For RDPGs this is indeed the case, as one can prove that the latent blocks are recovered up to an orthogonal transformation matrix in the large $n$ limit (Lemma 4 in TangEtAl2017). Therefore we can recover the parameter $\beta$ up to relabeling of the blocks. This is summarized in the following theorem.
The previous theoretical result implicitly assumes a dense network. However, many social and economic networks of interest in applications display some degree of sparsity. This is an empirical regularity that social scientists have observed in many settings, as most people do not form many links. Economists rationalize sparsity with the fact that people have constraints on time to spend with their friends Jackson2008.
We follow the literature and assume that sparsity is an asymptotic feature of the data generating process. We multiply the probability $\bm{P}_{ij}$ by a scalar $\rho_n$ that governs the sparsity of the network, that is, the probability of a link between nodes $i$ and $j$ becomes
Our previous result in Theorem (ref) applies to dense networks; that is when $\rho_n\rightarrow c$ where $c\in (0,1]$ is a constant. For simplicity and without loss of generality, in Theorem (ref) we have assumed $c=1$.
In this section we consider formally the case of $\rho_n\rightarrow 0$ as $n\rightarrow\infty$. We have to limit the rate of convergence for $\rho_n$, because the network could become too sparse, not allowing estimation. We will describe this regime a semi-sparse, because we will allow $\rho_n\rightarrow 0$ but $n\rho_n = \omega(\sqrt{n})$, that is the average degree of the network grows sub-linearly in $n$.\footnote{The notation $n\rho_n = \omega(\sqrt{n})$ means that for any real constant $a>0$ there exists an $n_0\geq 1$ such that $\rho_n > a /\sqrt{n}\geq 0$ for every integer $n\geq n_0$.} The intuition for this restriction is that too much sparsity makes links “too rare” and therefore spectral estimation and inference are impeded by having too few observations.
Theorem (ref) says that as long as the network is not too sparse, the estimator of $\beta$ will be asymptotically normal. Notably, the bias term does not vanish asymptotically.
We note that our formulation of sparsity does not impose any restriction on the network data for a fixed $n$, as it is based on a large sample property.
The asymptotic results hold for discrete observed covariates and more general models, as long as the effect of the observed covariates on the probability of linking can be written as a function of the block probabilities. Let the observed variables $\bm{Z}_i = [\bm{Z}_{i}^{(1)},\bm{Z}_{i}^{(2)}]$ be two covariates. For simplicity, we consider the case of binary variables, and we assume $\bm{Z}_{i}^{(1)}\overset{ind}{\sim} Bernoulli(b_{\tau_i}^{(1)})$ and $\bm{Z}_{i}^{(2)}\overset{ind}{\sim} Bernoulli(b_{\tau_i}^{(2)})$. The results still hold for discrete variables. The model is
This stochastic blockmodel has $\tilde{K}=4K$ blocks, $ (\bm{A}, \bm{\xi}, \bm{Z}) \sim SBM(\bm{\theta}_{Z},\bm{\eta})$ with $\widetilde{K}\times \widetilde{K}$ matrix of probabilities $\bm{\theta}_{Z}$ given by
where each matrix $\bm{W}_{k\ell}$ is given by
The intuition is the same as the model with one covariate. The blocks can be inferred by clustering the diagonal elements of matrix $\bm{\theta}_{Z}$, and the parameters $\beta_1$ and $\beta_2$ are functions of the $\bm{\theta}_{Z}$ entries, namely
As such, the main characterization of the central limit theorem holds in this case with minimal modifications.
In many applications the researcher is interested in testing for differential homophily in observable characteristics. For example, homophily among males could be higher than homophily among females, other things being equal. This can be accomplished in our setting by a minor modification of the algorithm for estimation.
For ease of exposition, we will again consider the model with $K=2$ blocks and a binary covariate $\bm{Z}_i\in \lbrace 0,1 \rbrace$, where for concreteness we assume that $\bm{Z}_i = 0$ denotes “male”. The model with differential homophily is
where $\beta_1$ measures the impact of being both males on the probability of a link, while $\beta_2$ measures the effect of homophily among females. We can thus write the corresponding probability matrix $\bm{\theta}_Z$ for this model as follows
We notice that clustering the main diagonal cannot inform about the structure of the blocks as in the previous examples. However, if we were to know the block assignments we will be able to identify $\beta_1$ and $\beta_2$ a functions of the entries of $\bm{\theta}_Z$
Notice that both $\beta_1$ and $\beta_2$ are functions of the stochastic blockmodel probabilities, thus the main characterization of the central limit theorem still holds.
We can estimate the block structure in the following way. First, consider the subgraph among all nodes with $Z=0$. The corresponding probability matrix of this stochastic blockmodel is
The structure of the model ensures that this subgraph is a GRDPG and we can estimate the latent positions using our algorithm. We can then cluster the estimated latent positions to obtain the estimated block assignments.
Second, we repeat the same procedure for the subgraph among nodes with $Z=1$. Third, we estimate the GRDPG using our algorithm, and use the estimated block assignments in the first and second step. This allows us to estimate both $\beta_1$ and $\beta_2$ through the formula above.
For ease of exposition let's consider the model with one binary covariate in Theorem (ref).\footnote{The estimator for the other cases (multiple covariates, differential homophily, discrete covariates, etc. is obtained analogously.} Our central limit theorem focus on the differences of two entries of the matrix $\bm{\theta}_Z$. However, we can compute $\beta$ in several ways, using different entries of the matrix, e.g., $\beta = h^{-1}\left(\bm{\theta}_{Z,11}\right) - h^{-1}\left(\bm{\theta}_{Z,12}\right) = h^{-1}\left(\bm{\theta}_{Z,33}\right) - h^{-1}\left(\bm{\theta}_{Z,34}\right)$. Therefore, we rely on two ways to estimate the model. The first estimator consists of computing all the values of $\beta$ from the relevant pairs of entries of $h^{-1}\left(\bm{\theta}_Z\right)$ and then averaging out. The second estimator weights each estimated $\beta$ by the size of the corresponding blocks. Formally, for the first estimator, after estimating the block assignments $\widehat{\bm{\xi}}$, we compute the number of nodes in the block with a particular value of the covariate
and we assign each block to a value of the covariate $\bm{Z}_{\theta,k} = 1$ if $n_{1,k}>n_{0,k}$ and $\bm{Z}_{\theta,k} = 0$ otherwise. Let $\bm{\psi}=(\psi_1,\cdots,\psi_{2K} ) \in \lbrace 1, \cdots, K\rbrace^{2K}$ be the vector that assigns each element of the diagonal of $\theta_Z$ to the corresponding unobserved block; let $\widehat{\bm{\psi}}$ be the corresponding estimated assignment. We then consider the set of all pairs set $M = \lbrace(k\ell, k\ell^\prime), k, \ell, \ell^\prime \in \lbrace 1, 2, \cdots, 2\widehat{K} \rbrace \vert \widehat{\psi}_\ell = \widehat{\psi}_{\ell^{\prime}}, \bm{Z}_{\theta,k}=\bm{Z}_{\theta,\ell}, \bm{Z}_{\theta,k}\neq \bm{Z}_{\theta,\ell^{\prime}}\rbrace $. For each element $m = (k\ell, k\ell^\prime)\in M$ we compute the estimate $\widehat{\beta}_m$ as
We then average out the values of the $\widehat{\beta}_m$ to obtain the final estimate
where $\vert M \vert$ is the number of elements in set $M$, that is the number of paired entries in $h^{-1}\left(\widehat{\bm{\theta}}_Z\right)$ from which we can estimate $\beta$.
The second estimator weights each estimated $\widehat{\beta}_{k\ell,k\ell^\prime}=h^{-1}\left(\bm{\theta}_{Z,k\ell}\right) - h^{-1}\left(\bm{\theta}_{Z,k\ell^\prime}\right)$, by the size of the corresponding blocks used in its estimation. Formally we compute the weight
where $n_k =\sum_{i=1}^n \mathbf{1}_{\lbrace \widehat{\bm{\xi}}_i=k \rbrace} $. We then consider the pairs in the set $\Omega = \lbrace (\ell,\ell^\prime), \ell,\ell^\prime \in \lbrace 1,\cdots , 2\widehat{K}\rbrace \vert \widehat{\psi}_\ell = \widehat{\psi}_{\ell^\prime} \rbrace $ and estimate the weighted sum
Both estimators work well in practice, and we show some evidence in the examples and simulations.
We compare our spectral methods to a standard algorithm used in the literature, the variational EM algorithm, as implemented in the R package blockmodels. Our methods are implemented in the package grdpg, available on Github at \url{https://github.com/meleangelo/grdpg}. All the replication files for simulations and empirical application are at \url{https://github.com/meleangelo/grdpg_supplement}. We also note that the variational EM algorithm uses parallelization to increase computational efficiency, while our method is implemented without any parallelization and for networks with thousands of nodes.
In our first example, we do not include any covariates and we assume $h$ is the identity function, so that the link probabilities are defined by $\bm{P}_{ij}=\bm{X}_i^T \bm{X}_j$. We simulate networks with $n=2000, 5000, 10000,$ and $20000$ nodes, with latent space dimension $d=1$. In Table (ref) we report the results for $K=2$, with block centers $[p,q]=[0.1,0.7]$, and matrix of probabilities
For simplicity, we assume that blocks are equally likely, that is $(\pi_1,\pi_2 )= (0.5,.0.5)$. To evaluate the performance of the algorithms, we compare clustering accuracy and computational time. The assignment of nodes to the correct block is summarized by the Adjusted Rand Index (ARI) RandARI1971, and the computational time is given by the CPU time in seconds. Our point estimates are shown in Table (ref), and below we report the estimated block probabilities for $n=2000$.
The values $\hat{p}, \hat{q}$ shown in in Table (ref) are obtained by singular value decomposition of the estimated probability matrix (and rotation).
We notice that the VEM and GRDPG estimators produce similar point estimates and very precise clustering of the nodes, as indicated by the ARI. However, our GRDPG estimator converges much faster than the VEM. For networks with $n=10000$ nodes, our method provides estimates in approximately 30 seconds, while the VEM takes more than one hour to converge to the final approximation. When $n=20000$ and $n=30000$, our GRDPG approach converges in about 2 minutes and less than 7 minutes, respectively, while the VEM is impractical.
Here, a crucial choice is the number of dimensions for the spectral embedding. In our simulation we know that the rank of the matrix $\bm{\theta}$ is 1, therefore this is the optimal dimension (see AthreyaEtAl2018a). We choose $\hat{d}$ by profile likelihood methods as in ZhuGhodsi2006.\footnote{The screeplot, not shown, displays a huge step down in the (absolute) value of the eigenvalues of the adjacency matrix at the largest eigenvalue, which suggests that 1 dimension is sufficient to approximate the structure of the adjacency matrix.}
The clustering of the latent positions in blocks is performed using the MCLUST method implemented in the package Mclust in R FraleyRaftery1999.
In Table (ref) we report results from the same model with $K=5$ and latent positions $\bm{\nu}=(0.1,0.3,0.5,0.7,0.9)$. The results are comparable to the previous table, our estimator scales very well with the size of the network, while obtaining the same point estimates of the VEM algorithm. In this example, the difference in scaling for the two estimators is more pronounced. In particular, going from $K=2$ to $K=5$ blocks does not increase the computational burden too much for the GRDPG-based estimator.
We consider a model with a binary nodal covariate, $\bm{Z}_i \sim Bernoulli(0.5)$ and link probabilities
In this example we use $[p,q]=[-1.5, 1]$ and $\beta=1.5$, thus the matrix $\bm{\theta}$ is
while the full matrix $\bm{\theta}_Z$ that includes the effect of covariates is
We choose $\hat{d}$ by profile likelihood ZhuGhodsi2006. In Figure (ref) we show the screeplots. In the upper-left, we display the screeplot for the adjacency matrix, which suggests using $\hat{d}=4$ as an estimate of the dimension for $\widehat{\bm{Y}}$. We note that the fourth largest eigenvalue is negative, and the GRDPG model takes this into account. In the center-left plot, we show the screeplot of the adjacency matrix after netting out the effect of the covariates, which suggests the estimate $\hat{d}=1$ for determining the dimension of the unobserved latent positions $\bm{X}$.
The point estimates for $\bm{\theta}_Z$ (up to a permutation of the block labels) when the network has $n=2000$ nodes are respectively
According to our procedure, there are several ways to obtain an estimate of $\beta$. From matrix ((ref)), we group rows 1 and 2 in one block, and rows 3 and 4 in another block by clustering the diagonal entries. We can get an estimate of $\beta$ as $\widehat{\bm{B}}_{Z,11}-\widehat{\bm{B}}_{Z,12}$ or $\widehat{\bm{B}}_{Z,22}-\widehat{\bm{B}}_{Z,21}$ or $\widehat{\bm{B}}_{Z,33}-\widehat{\bm{B}}_{Z,34}$, etc. We know by our theorem that each of these estimators is asymptotically normal. Instead of choosing which entries to use to estimate $\beta$, we pool all possible estimates, weighting them by the proportion of observations that are assigned to each block. For example, the estimate $\widehat{\bm{B}}_{Z,11}-\widehat{\bm{B}}_{Z,12}$ is weighted by the proportion of links in the network that are used to estimate it.
The point estimates for $\beta$ reported in Table (ref) are $\widehat{\beta}_{GRDPG}=1.51201$ and $\widehat{\beta}_{VEM} = 1.50335$. The estimated latent positions are $\hat{p}=-1.49712$ and $\hat{q}=1.00067$ for the VEM; and $\hat{p}=-1.49454$ and $\hat{q}=0.99926$ for the GRDPG estimator. However, it takes almost 2 hours to obtain the VEM results, while it only takes 7 seconds with our estimator. The left plots in Figure (ref) show the latent positions $\widehat{\bm{Y}}$ of the GRDPG (including the effect of covariates) estimated by ASE. We plot the first coordinate against each of the other three. In the second and third plot from the top, we can notice that the latent positions nicely cluster into 4 blocks, as our theory predicts. In the bottom-left plot in Figure (ref) we display the estimated latent positions $\widehat{\bm{X}}$, estimated by netting out the effect of the covariates. The figure shows how the estimated latent positions $\widehat{\bm{X}}$ cluster around the true values $p$ and $q$ (the black vertical lines).
As explained above, a central advantage of our approach is computational speed. Indeed, in Table (ref) we show that when we increase the size of the network to $n=5000$, the estimated parameters are essentially the same for VEM and GRPDG. However, the GRDPG estimator take less than 30 seconds to converge; the VEM estimate takes almost 10 hours.
In Table (ref) we show estimates for models with latent positions $\bm{\nu}_1=(-1.5,-1.0)$ and $\bm{\nu}_2=(1.0,0.5)$. For the simulations in the first 3 rows we set $\beta=1.5$. It is quite remarkable that the computational time does not increase much, with respect to the case of $d=1$.
The second group of three rows shows the results of simulations with smaller $\beta=0.5$. This makes the estimation of the covariate effect more challenging. Indeed when $n=2000$ the classification in blocks and the point estimate are imprecise, as indicated by the adjusted Rand index (ARI). When we increase the network size to $n=5000$ and $n=10000$, the accuracy of the point estimates improves significantly. This example shows that our approach is extremely useful in very large networks, where VEM may become computationally impractical.
In summary, our simple examples and simulations show that our GRDPG-based estimator is quite fast and scales well to large networks. These good computational properties are obtained without sacrificing the accuracy of the estimates, as we prove that the algorithm produces the same point estimates as the variational EM in all the examples.
Given the computational times shown in the previous section, we run a simple Monte Carlo to understand the bias of the spectral estimator in networks of moderate size.\footnote{We do not compare the spectral estimator to the variational EM estimator, because the latter is too slow for a Monte Carlo with 1000 repetitions, even after parallelizing the execution. } We estimate a model with two binary observed covariates,
and vary the probabilities $b_z$ and $b_w$, as well as the correlation among the two variables. We estimate the following model in each Monte Carlo design
The Monte Carlo design considers networks of sizes $n=2000, 5000, 10000$ and we set the number of blocks to $K=2$. For all the simulations the parameter value that generates the data is $\bm{\beta} = (0.5,0.75)$ and the centers of the blocks are $\bm{\nu}=(-1.5,1.0)$. We summarize the designs in Table (ref). For each design and network size, we simulate 1000 networks and estimate the parameters of model (ref) with the simple mean estimator and the weighted mean estimator.
The first design corresponds to the examples in the previous section. Blocks are assumed to have same size and observables are independent Bernoulli variables with equal probability. The second design introduces correlation among the observables, as this is a realistic feature of many datasets. Designs 3 and 4 are intended to test the effect of unbalanced block size and unbalanced covariates, respectively, while maintaining the assumption of independence among observables. The final design assumes that we have unbalanced blocks, unbalanced covariates and correlated observables, allowing us to understand how the estimator behaves in a realistic setting. We expect that unbalancedness and correlation will increase bias, but this problem is less severe for larger networks, as our theory shows that the bias becomes vanishingly small as we increase the size of the network.
The results of our simulations are reported in Tables (ref) (simple mean estimator) and (ref) (weighted mean estimator). For each design, we report the absolute difference between the estimated parameter and the true value, the Monte Carlo standard error and the average time for estimation.\footnote{The time of estimation reported in the table includes the following steps: 1) compute the ASE from the adjacency matrix; 2) compute the matrix of latent positions; 3) cluster latent positions to recover blocks; 4) compute matrix $\widehat{\bm{B}}_Z$; 5) cluster diagonal entries of matrix $\widehat{\bm{B}}_Z$ to recover the unobservable block structure; 6) estimate $\bm{\beta}$ using the information on the block structure and the entries of matrix $\widehat{\bm{B}}_Z$; 7) compute simple mean and weighted mean of the estimated $\widehat{\bm{\beta}} $ according to (ref) and (ref). The simulation takes a little longer because we need to generate the data and the adjacency matrices for the Monte Carlo. Code is available in Github.} When using the simple mean estimator, the estimates are precise, while displaying a small bias. The most challenging design for our estimator is Design 5, where we impose different unobserved block size, different Bernoulli probabilities for the observables and correlation among observed characteristics. As expected, these features increases the bias in our estimates; however, this problem becomes less severe with larger networks.
Our weighted mean estimator has similar behavior.
In the empirical application we estimate the standard errors for the observed covariates effects according to the formula in Theorem (ref), using a plug-in estimate. To understand the behavior of this plug-in estimator, we ran several Monte Carlo experiments. In Table (ref) we report results using the model with independent covariates and balanced blocks (Design 1 in Table (ref)).
In columns 1 and 2 we show the true standard error computed using the formula in Theorem (ref) for the simple mean estimator; in columns 3 and 4, the standard error is computed using the weighted estimator with the true proportions of each block; in columns 5 and 6 we have our weighted estimator, using the estimated parameters and block proportions. Our plug-in estimator is very conservative. For a network with $5000$ nodes, the estimated standard error is extremely high, when compared to the true value. Even when the network size is increased to $10000$ nodes, the plug-in estimate for the standard error is quite large.
We conclude that using a plug-in estimator usually overestimates the standard errors. Our estimates are thus very conservative.
We apply our method to study a network of Facebook friendships, using the Facebook 100 dataset from TraudEtAl2012.\footnote{The entire dataset is available at https://archive.org/details/oxford-2005-facebook-matrix.} This network was extracted from the Facebook platform in September 2005, providing a snapshot of the friendship among students, faculty, staff and alumni at 100 U.S. universities.
We perform an analysis similar to AtchadeEtAl2019, using the Harvard University network data. The dataset consists of 15126 nodes and 7 nodal covariates: role, gender, major, minor, dorm, year, and high school. These are all discrete variables. We focus on dorm, gender, and role in the analysis. We make each of these variables binary. So rather than specify specific dorm information, our control variable indicates whether the student lives on or off-campus. The role is binarized to indicate whether the node is a student or not.\footnote{Roles include students, faculty, staff, alumni, etc. We focus on students because they are the ones mostly using the platform in 2005.}
The characteristics of the network are shown in Table (ref). In the first column we report the descriptive statistics for the original data, containing $n=15126$ nodes, with an average degree of $109.03$, an average clustering coefficient of $0.135$, with $46.6\%$ females, $50.8\%$ students and $22.9\%$ of people living off-campus. The average degree and clustering coefficient suggest that this is a moderately sparse network.
The second column contains the descriptive statistics for our processed sample. We keep all nodes with non-missing gender information; once we delete all missing gender information, there are no missing values for the other covariates; we then compute the largest connected component of the resulting network.\footnote{This is a standard procedure in the literature on SBMs AtchadeEtAl2019, AthreyaEtAl2018a, Abbe2018.} The final network contains a higher proportion of females ($53.9\%$) and students ($52\%$); and a smaller fraction of people living off-campus ($15.3\%$) than the original data. The average degree is slightly smaller ($105.02$), because we have deleted some nodes and edges; however, the clustering coefficients are of similar magnitude, $0.135$ and $0.137$, respectively.\footnote{Before we proceed to estimation, we regularize the adjacency matrix using the standard method proposed in LevinaRegularization2017. This regularization step avoids numerical issues with the spectral decomposition arising from significant node degree heterogeneity.}
We estimate the following model with three covariates
where $female_i=1$ if the node is female, $student_i=1$ if the node is a student and $off-campus_i=1$ if the node lives off-campus. We first estimate models with one binary covariate, using each of our control variables individually. Next we estimate the full model (ref) with three controls. For the Adjacency Spectral Embedding we choose the dimension of the latent space $\hat{d}=2$, using the profile likelihood method in ZhuGhodsi2006, as in our simulations.\footnote{Multiple methods exist for selecting the embedding dimension in practice, and this remains a topic of current research. In the context of networks, choosing a dimension smaller than the true $d$ will introduce bias in the estimated latent positions; on the other hand, using an embedding dimension larger than the true $d$ will increase the variance of the estimated latent positions. In this trade-off we prefer to err on the side of overestimating $d$. Specifically, we choose the value one plus the location of the first elbow in the screeplot.}
The clustering of the estimated latent positions is performed with a Gaussian mixture model, using the MCLUST implementation of FraleyRaftery1999 in R. We obtain latent positions estimates and $\widehat{\tilde{K}}=32$ blocks from the adjacency matrix. We then obtain the estimated matrix $\hat{\bm{B}}_Z$, and cluster its diagonal entries to recover the (unobserved) blocks $\hat{K}$ and estimate the vector of parameters $\bm{\beta}$.
The estimated parameters are shown in Table (ref). In the first three columns we report estimates for models with a single binary covariate. Each coefficient is precisely estimated, according to our naive plug-in standard error estimator; the estimated effects are all positive, which we interpret as evidence of homophily, a usual feature of many social networks. The point estimates are very similar when we estimate the full model (ref) (Column 4). These results are also consistent with the analysis in TraudEtAl2012 and AtchadeEtAl2019. Our method estimates $\widehat{K}=4$ unobserved blocks. If we choose to include only one covariate, the number of blocks estimated is $\widehat{K}=16$ (Columns 1--3).
The results indicate that there are additional unobserved characteristics that affect the network formation among Facebook users. In particular, the four blocks may capture shared interests, common preferences, similar schedules and additional information that is unobservable to the researcher.
In summary, our method is computationally faster than the traditional VEM, performs well in Monte Carlo simulations with different designs and data generating processes, and provides useful insights in real-data, empirical applications.
We have developed a spectral estimator for stochastic blockmodels with nodal covariates in large networks. The main theoretical contribution is an asymptotic normality result for the spectral estimator of the covariates' effect on the probability of linking. Our work leverages the relationship between generalized random dot product graphs and stochastic blockmodels, extending existing frameworks to include observed covariates and constructing an estimator that is fast and scalable for large networks. Our theoretical results also apply to moderately sparse graphs, which is important in a host of applications in economics and more generally in social sciences, public health, and computer science, where network data are often viewed as being sparse.
We have provided examples demonstrating that our method delivers the same accuracy as the variational EM algorithm, while converging much faster. Our Monte Carlo simulations and the empirical application show that this method works best in very large networks, when the variational EM becomes impractical.
We consider the present work a first step in the study of this class of models and the foundation for inference for SBMs and other latent position models for large networks with nodal covariates. While we have focused on binary and discrete covariates in this work, extensions to continuous covariates are currently being pursued via recently developed Latent Structure Models AthreyaLSM2018. In future work, similar ideas can also be applied to directed networks and bipartite networks AKM1999, BonhommeLamadonManresa2018, significantly expanding the realm of GRDPG applications in economics and social sciences.