The method of moments and degree distributions for network models

Introduction

The analysis of network data has become an important component of doing research in many fields; examples include social and friendship networks, food webs, protein interaction and regulatory networks in genomics, the World Wide web and computer networks. On the algorithmic side, many algorithms for identifying important network structures such as communities have been proposed, mainly by computer scientists and physicists; on the mathematical side, various probability models for random graphs have been studied. However, there has only been a limited amount of research on statistical inference for networks, and on learning the network features by fitting models to data; to a large extent, this is due to the gap between the relatively simple models that are analytically tractable and the complex features of real networks not easily reproduced by these models.

Probability models on infinite graphs have a nice general representation based on results [Aldous 1981, Hoover 1979, Kallenberg 2005, Diaconis and Janson 2008], analogous to de Finetti’s theorem, for exchangeable matrices. Here, we give a brief summary closely following the notation of Bickel and Chen 2009. Graphs can be represented through their adjacency matrix $A$, where $A_{ij}=1$ if there is an edge from node $i$ to $j$ and 0 otherwise. We assume $A_{ii}=0$, that is, there are no self-loops. $A_{ij}$’s can also represent edge weights if the graph is weighted, and for undirected graphs, which is our focus here, $A_{ij}=A_{ji}$. For an unlabeled random graph, it is natural to require its probability distribution $P$ on the set of all matrices $\{[A_{ij}],i,j\geq 1\}$ to satisfy $[A_{\sigma_{i}\sigma_{j}}]\sim P$, where $\sigma$ is an arbitrary permutation of node indices. In that case, using the characterizations above one can write

where $\alpha$, ${\xi_{i}}$ and $\lambda_{ij}$ are i.i.d. random variables distributed uniformly on $(0,1)$, $\lambda_{ij}=\lambda_{ji}$ and $g$ is a function symmetric in its second and third arguments. $\alpha$ as in de Finetti’s theorem corresponds to the mixing distribution and is not identifiable. The equivalent of the i.i.d. sequences in de Finetti’s theorem here are distributions of the form $A_{ij}=g(\xi_{i},\xi_{j},\lambda_{ij})$. This representation is not unique, and $g$ is not identifiable. These distributions can be parametrized through the function

be the probability of an edge in the network. Then the density of $(\xi_{i},\xi_{j})$ conditional on $A_{ij}=1$ is given by

With this parametrization, it is natural to let $\rho=\rho_{n}$, make $w$ independent of $n$ and control the rate of the expected degree $\lambda_{n}=(n-1)\rho_{n}$ as $n\rightarrow\infty$. The case most studied in probability on random graphs is $\lambda_{n}=\Omega(1)$ [where $a_{n}=\Omega(b_{n})$ means $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$]. The case of $\lambda_{n}=1$ corresponds to the so-called phase transition, with the giant connected component emerging for $\lambda_{n}>1$.

Many previously studied probability models for networks fall into this class. It includes the block model [Holland, Laskey and Leinhardt 1983, Snijders and Nowicki 1997, Nowicki and Snijders 2001], the configuration model [Chung and Lu 2002] and many latent variable models, including the univariate [Hoff, Raftery and Handcock 2002] and multivariate [Handcock, Raftery and Tantrum 2007] latent variable models, and latent feature models [Hoff 2007]. In fact, dynamically defined models such as the “preferential attachment” model [which seems to have been first mentioned by Yule in the 1920s, formally described by de Solla Price 1965 and given its modern name by Barabási and Albert 1999] can also be thought of in this way if the dynamical construction process continues forever producing an infinite graph; see Section 16 of Bollobás, Janson and Riordan 2007.

The block model is very attractive from the analytical point of view and useful in a number of applications, but the class (2) is much richer than the block model itself. Moreover, the block model cannot deal with nonuniform edge distributions within blocks, such as the commonly encountered “hubs,” although a modification of the block model introducing extra node-specific parameters has been recently proposed by Karrer and Newman 2011 to address this shortcoming. It may also be difficult to obtain accurate results from fitting the block model by maximum likelihood when the graph is sparse.

In this paper, we develop an alternative approach to fitting models of type (2), via the classical tool of the method of moments. By moments, we mean empirical or theoretical frequencies of occurrences of particular patterns in a graph, such as commonly used triangles and stars, although the theory is for general patterns. While specific parametric models like the block model can be fitted by other methods, the method of moments applies much more generally, and leads to some general theoretical results on graph moments along the way. We note that related work on the method of moments was carried out for some specific parametric models in Picard et al. 2008.

A well-studied class of random graph models where moments play a big role is the exponential random graph models (ERGMs). ERGMs are an exponential family of probability distributions on graphs of fixed size that use network moments such as number of edges, $p$-stars and triangles as sufficient statistics. ERGMs were first proposed by Holland and Leinhardt 1981 and Frank and Strauss 1986 and have then been generalized in various ways by including nodal covariates or forcing particular constraints on the parameter space; see Robins et al. 2007 and references therein. While the ERGMs are relatively tractable, fitting them is difficult since the partition function can be notoriously hard to estimate. Moreover, they often fail to provide a good fit to data. Recent research has shown that a wide range of ERGMs are asymptotically either too simplistic, that is, they become equivalent to Erdös–Renyi graphs, or nearly degenerate, that is, have no edges or are complete; see Handcock 2003 for empirical studies and Chatterjee and Diaconis 2011 and Shalizi and Rinaldo 2011 for theoretical analysis.

The rest of the paper is organized as follows. In Section 2, we set up the notation and problem formulation and study the distribution of empirical moments, proving a central limit theorem for acyclic patterns. We also work out examples for several specific patterns. In Section 3 we show how to use the method of moments to fit the block model, as well as identify a general nonparametric model of type (2). In Section 4, we focus on degree distributions, which characterize (asymptotically) the model (2). Section 5 discusses the relationship between normalized degrees and more complicated pattern counts that can be used to simplify computation of empirical moments. Section 6 concludes with a discussion. Proofs and additional lemmas are given in the Appendix.

The asymptotic distribution of moments

We start by setting up notation. Let $G_{n}$ be a random graph on vertices $1,\ldots,n$, generated by

where $w(u,v)\geq 0$, symmetric, $0\leq u,v\leq 1$, $\rho_{n}\rightarrow 0$. We cannot, unfortunately, treat $\rho_{n}$ and $w$ as two completely free parameters, as we need to ensure that $h\leq 1$. We can either assume that the sequence $\rho_{n}$ is such that $\rho_{n}w\leq 1$ for all $n$, or restrict our attention to classes where $w_{n}(u,v)=w(u,v)I(w(u,v)\leq\rho_{n}^{-1})\stackrel{{\scriptstyle L_{2}}}{{\rightarrow}}w(u,v)$. In either case, we can ignore the weak dependence of $w_{n}$ on $\rho_{n}$ and effectively replace $w_{n}$ with $w$.

Let $T\colon\;\mathcal{L}_{2}(0,1)\rightarrow\mathcal{L}_{2}(0,1)$ be the operator defined by

Thus $D_{i}$ is the degree of node $i$, $\bar{D}$ is the average degree and $L$ is the total number of edges in $G_{n}$.

Let $R$ be a subset of $\{(i,j)\colon\;1\leq i<j\leq n\}$. We identify $R$ with the vertex set $V(R)=\{i\colon\;(i,j)\mbox{ or }(j,i)\in R\mbox{ for some }j\}$ and the edge set $E(R)=R$. Let $G_{n}(R)$ be the subgraph of $G_{n}$ induced by $V(R)$. Recall that two graphs $R_{1}$ and $R_{2}$ are called isomorphic ($R_{1}\sim R_{2}$) if there exists a one-to-one map $\sigma$ of $V(R_{1})$ to $V(R_{2})$ such that the map $(i,j)\rightarrow(\sigma_{i},\sigma_{j})$ is one-to-one from $E(R_{1})$ to $E(R_{2})$.

Throughout the paper, we will be using two key quantities defined next:

Next, we give a proposition summarizing some simple relationships between $P$ and $Q$. The proof, which is elementary, is given in the Appendix. Similar results are implicit in Diaconis and Janson 2008.

If $G_{n}$ is a random graph, and $R$ a subset of $\{(i,j)\colon\;1\leq i<j\leq n\}$, then

where $\bar{R}=\{(i,j)\notin R,i\in V(R),j\in V(R)\}$. Further,

Here $R\subset S$ refers to $S\subset\{(i,j)\colon\;i,j\in V(R)\}$.

The quantities $P(R)$ and $Q(R)$ are unknown population quantities which we can estimate from data, that is, from the graph $G_{n}$. Define, for $R\subset\{(i,j)\colon\;1\leq i<j\leq n\}$ with $|V(R)|=p$,

where $N(R)$ is the number of graphs isomorphic to $R$ on vertices $1,\ldots,p$. For instance, if $R$ is a 2-star consisting of two edges $(1,2)$, $(1,3)$, then $N(R)=3$. Further, let

Here we use $R$ and $S$ to denote both a subset and a subgraph. Evidently,

The scaling here is controlled by the parameter $\rho_{n}$, the natural assumption for which is $\rho_{n}\rightarrow 0$. In that case, $P(R)\rightarrow 0$ for any fixed $R$ with a fixed number of vertices $p$. Therefore we consider the following rescaling of $P(R)$ and $Q(R)$: writing $|R|$ for $|E(R)|$, let

if $\int w^{2(|R|+1)}(u,v)\,du\,dv<\infty$.

$\hat{\rho}_{n}=\frac{\bar{D}}{n-1}=\frac{2L}{n(n-1)}$ is the estimated probability of an edge. For these rescaled versions of $P$ and $Q$, we have the following theorem.

Suppose $\int_{0}^{1}\int_{0}^{1}w^{2}(u,v)\,dv\,du<\infty$.

for some $\sigma^{2}>0$. Suppose further $R$ is fixed, acyclic with $|V(R)|=p$ and $\int w^{2|R|}(u,v)\,du\,dv<\infty$. Then,

More generally, for any fixed $\{R_{1},\ldots,R_{k}\}$ as above with $|V(R_{j})|\leq p$,

Suppose $\lambda_{n}\rightarrow\lambda<\infty$. Conclusions (9)–(12) continue to hold save that $\sigma^{2}(R)$, $\Sigma(R)$ depend on $\lambda$ as well as $R$.

(1) Note that part (b) yields consistency and asymptotic normality of acyclic graph moment estimates across the phase transition to a giant component, that is, for $\lambda<1$ as well as $\lambda\geq 1$.

(2) Note that we are, throughout, estimating features of the canonical $w$. Unnormalized $P$ and $Q$ are trivially 0 if $\lambda_{n}$ is not of order $n$.

(4) Part (c) of the theorem shows that for graphs with $\lambda_{n}=\Omega(n)$, $\check{Q}$ always gives $\sqrt{n}$-consistent estimates of any pattern while $\check{P}$ is not consistent unless we assume acyclic graphs,

since the bias is of order $O(\lambda_{n}/n)=O(1)$. In the range $\lambda_{n}=o(n^{1/2})$ to $\Omega(n)$, what is possible depends on the pattern. For instance, if $\Delta=\{(1,2),(2,3),(3,1)\}$, a triangle, $\check{P}(\Delta)=\check{Q}(\Delta)$ (because there is no other graph on three nodes containing $\Delta$), and $\check{P}$ is $\sqrt{n}$-consistent if $\lambda_{n}\geq\varepsilon n^{1/3}$ by part (c) but otherwise only consistent if $\lambda_{n}\to\infty$.

2 Examples of specific patterns

Next we give explicit formulas for several specific $R$. Our main focus is on wheels (defined next), which, as we shall see, in principle can determine the canonical $w$.

A $(k,l)$-wheel is a graph with $kl+1$ vertices and$kl$ edges isomorphic to the graph with edges $\{(1,2),\ldots,(k,k+1);(1,k+2),\penalty\ldots,(2k,2k+1);\ldots,(1,(l-1)k+2),\ldots,(lk,lk+1)\}$.

In other words, a wheel consists of node $1$ at the center and $l$ “spokes” connected to the center, and each spoke is a chain of $k$ edges. We consider only $k\geq 2$. The number of isomorphic $(k,l)$-wheels on vertices $1,\ldots,p$ is $N(R)=(kl+1)!/l!$.

If the graph $R$ is a $(k,l)$-wheel, the theoretical moments have a simple form and can be expressed in terms of the operator $T$ as follows:

where the first equality holds by the definition of $Q$ and the second by the structure of a $(k,l)$-wheel.

order larger than $n^{1-2/(kl+1)}$. In the $\lambda_{n}$ range between $O(n^{1/2})$ and $O(n^{1-2/(kl+1)})$, we do not exhibit a $\sqrt{n}$-consistent estimate though we conjecture that by appropriate de-biasing of $\check{P}$ such an estimate may be constructed. However, $\lambda_{n}=o(n^{1/2})$ seems a reasonable assumption

for most graphs in practice, and then we can use the more easily computed $\check{P}$.

A $(\mathbf{k},\mathbf{l})$-wheel, where $\mathbf{k}=(k_{1},\ldots,k_{t})$, $\mathbf{l}=(l_{1},\ldots,l_{t})$ are vectors and the $k_{j}$’s are distinct integers, is the union $R_{1}\cup\cdots\cup R_{t}$, where $R_{j}$ is a $(k_{j},l_{j})$-wheel, $j=1,\ldots,t$, and the wheels $R_{1},\ldots,R_{t}$ share a common hub but all their spokes are disjoint.

$(\mathbf{k},\mathbf{l})$-wheel has a total of $p=\sum_{j}l_{j}k_{j}+1$ vertices and $\sum_{j}l_{j}k_{j}$ edges. For example, a graph defined by $E=\{(1,2);(1,3),(3,4);(1,5),(5,6);(1,7),(7,8)$, $(8,9)\}$ is a $(\mathbf{k},\mathbf{l})$-wheel with $\mathbf{k}=(1,2,3)$ and $\mathbf{l}=(1,2,1)$. The number of distinct isomorphic $(\mathbf{k},\mathbf{l})$-wheels on $p$ vertices is $N(R)=p!(\prod_{j}l_{j}!)^{-1}$.

We can compute, defining $A(R)=\prod\{A_{ij}\colon\;(i,j)\in R\}$,

Thus $(\mathbf{k},\mathbf{l})$-wheels give us all cross moments of $T^{m}(\xi)$, $m\geq 1$. Note that all $(\mathbf{k},\mathbf{l})$-wheels are acyclic.

We are not aware of other patterns for which the moment formulas are as simple as those for wheels. For example, if $R$ is a triangle, then

where $h^{(2)}(u,w)=\int_{0}^{1}h(u,v)h(v,w)\,dv$ corresponds to $T^{2}f\equiv\int_{0}^{1}h^{(2)}(u,v)\times\penalty f(v)\,dv$.

Moments and model identifiability

We establish two results in this section: identifiability of block models with known $K$ using $\{\check{P}(R)\colon\;R$ a $(k,l)$-wheel, $1\leq l\leq 2K-1,2\leq k\leq K\}$, and

the general identifiability of the function $w$ from $\{\check{P}(R)\}$ using all $(\mathbf{k},\mathbf{l})$-wheels $R$.

Let $w$ correspond to a $K$-block model defined by parameters $\theta\equiv(\pi,\rho_{n},S)$, where $\pi_{a}$ is the probability of a node being assigned to block $a$ as before, and

Recall that the function $h$ in (2) is not unique, but a canonical $h$ can be defined. For the block model, we use the canonical $h$ given by Bickel and Chen 2009. Let $H_{ab}=S_{ab}\pi_{a}\pi_{b}$. Let the labeling of the communities $1,\ldots,K$ satisfy $H_{1}\leq\cdots\leq H_{K}$, where $H_{a}=\sum_{b}H_{ab}$ is proportional to the expected degree for a member of block $a$. The canonical function $h$ then takes the value $F_{ab}$ on the $(a,b)$ block of the product partition where each axis is divided into intervals of lengths $\pi_{1},\ldots,\pi_{K}$. Let $F\equiv\|F_{ab}\|$.

In view of (10), we will treat $\rho_{n}$ as known. Let $\{W_{kl}\colon\;1\leq l\leq 2K-1,2\leq k\leq K\}$ be the specified set of $(k,l)$-wheels, and let

Suppose $\theta=(\pi,S)$ defines a block model with known $K$, and the vectors $\pi,F\pi,\ldots,F^{K-1}\pi$ are linearly independent. Suppose $\varepsilon\leq\lambda_{n}=o(n^{1/2})$. Then:

$\{\tau_{kl}\colon\;l=1,\ldots,2K-1,k=2,\ldots,K\}$ identify the $K(K+3)/2-2$ parameters of the block model other than $\rho$ (i.e., the map $f$ is one to one).

If $f$ has a gradient which is of rank $\frac{K(K+3)}{2}-2$ at the true $(\pi_{0},S_{0})$, then ${f}^{-1}(P(\check{\tau}))$ is a $\sqrt{n}$-consistent estimate of $(\pi_{0},S_{0})$, where $\check{\tau}=\|\check{\tau}_{kl}\|$ and $P(\check{\tau})$ is the closest point in the range of $f$ to $\check{\tau}$.

Note that the linear independence condition rules out all matrices $F$ that have $1$ as an eigenvector. In particular, it rules out the case of $F_{aa}$ equal for all $a$, $F_{ab}$ equal for all $a\neq b$, which was studied in detail by Decelle et al. 2011. Using physics arguments, they showed that in that particular case, when $\lambda=O(1)$, there are regions of the parameter space where neither the parameters nor the block assignments can be estimated by any method.

2 The nonparametric model

In the general case, we express everything in terms of the operator $T_{w}\equiv T/\rho_{n}$ induced by the canonical $w$. We require that:

the joint distribution of $\{T_{w}^{l}(1)(\xi)\colon\;l\geq 1\}$ is determined by the cross moments of $(T_{w}^{l_{1}}(\xi),\ldots,T_{w}^{l_{k}}(\xi))$, for $l_{1},\ldots,l_{k}$ arbitrary.

A simple sufficient condition for (A) is $|w|\leq M<\infty$. A more elaborate one is the following:

Let $w$ characterize $T_{w}$, where $\int_{0}^{1}w^{2}(u,v)\,du\,dv<\infty$. By Mercer’s theorem,

where the $\phi_{j}$ are orthonormal eigenfunctions and the $\lambda_{j}$ eigenvalues, $\sum\lambda_{j}^{2}<\infty$.

Suppose $\int_{0}^{1}\int_{0}^{1}w^{2}(u,v)\,du\,dv<\infty$. Assume the eigenvalues $\lambda_{1}>\lambda_{2}>\cdots$ of $T_{w}$ are each of multiplicity $1$ with corresponding eigenfunction $\phi_{j}$, and $\int_{0}^{1}\phi_{j}(u)\,du\neq 0$ for all $j$. The joint distribution of $(T_{w}(1)(\xi),\ldots,\penalty T_{w}^{m}(1)(\xi),\ldots)$ then determines, and is determined by, $w(\cdot,\cdot)$.

Note again that interesting cases are ruled out by the condition that all eigenfunctions of $T$ are not orthogonal to $1$. The general analogue to the block model case is that $P(A_{ij}=1|\xi_{i})$ cannot be constant for all $i$ and $j$. Constancy can be interpreted as saying that $A_{ij}$ and the latent variable $\xi_{i}$ associated with vertex $i$ are independent. The proof of Theorem 3 is given in the Appendix. The almost immediate application to wheels is stated next.

Since $\mathbf{T}_{l}\equiv(T(1)(\xi),\ldots,T^{l}(\xi))$ has a moment generating function converging on $0<|s|\leq\varepsilon_{l}$, the moments (including cross moments) determine the distribution of the vector. By (14), the $\tau_{\mathbf{kl}}$ give all moments of the vector $\mathbf{T}_{l}$ for all $l$. By Theorem 1, the $\check{\tau}_{\mathbf{kl}}$ are $\sqrt{n}$-consistent.

Degree distributions

The average degree $\bar{D}$ is, as we have seen in Theorem 1, a natural data dependent normalizer for moment statistics which eliminates the need to “know” $\rho_{n}$. In fact, as we show in this section, the joint empirical distribution of degrees and what we shall call $m$ degrees below can be used in estimating asymptotic approximations to $w(\cdot,\cdot)$ in a somewhat more direct way than moment statistics. They can also be used to approximate moment estimates based on $(\mathbf{k},\mathbf{l})$-wheels in a way that potentially simplifies computation.

The complexity of this computation is $O((n+m)\lambda_{n}^{m})$ (first term is for computing the row sums of $A^{m}$ and the second for eliminating the loops).

Define the empirical distribution of the vector of normalized degrees

Suppose $\lambda_{n}\to\infty$ and $|w_{2m}|<\infty$. Then ${\hat{F}}_{m}\stackrel{{\scriptstyle M_{2}}}{{\rightarrow}}F_{m}$ as $n\rightarrow\infty$, where $F_{m}$ is the distribution of $\bm{\theta}_{m}(\xi)=(\tau_{w}(\xi),\ldots,T_{w}^{m-1}(\tau_{w})(\xi))$, and $\tau_{w}(\xi)=\int_{0}^{1}w(\xi,v)\,dv$ is monotone increasing. Moreover, if $\hat{G}_{m}(\mathbf{x},\mathbf{y})$ is the empirical distribution of $(\mathbf{D}_{i}^{(m)},\bm{\theta}_{m}(\xi_{i}))$, then

There is an attractive interpretation of the last statement of Theorem 5. If $\lambda_{n}\to\infty$, $\lambda_{n}=o(n^{1/(m-1)})$, $m\geq 2$, then $D_{i}/\lambda_{n}$ can be identified with $\tau(\xi_{i})$ in the following sense: While $\xi_{i}$ is unobserved but $D_{i}/\bar{D}$ is, on average, $\tau(\xi_{i})$ and $D_{i}/\bar{D}$ are close. Since $\tau$ is monotone increasing in $\xi$, that is, is a measure of $\xi$ on another scale, we can treat $D_{i}/\lambda_{n}$ as the latent affinity of $i$ to form relationships.

Bollobás, Janson and Riordan 2007 show that if $m=1$, $\lambda_{n}=O(1)$, then the limit of the empirical distribution of the degrees can be described as follows: given $\xi\sim\mathcal{U}(0,1)$, the limit distribution is Poisson with mean $\tau_{w}(\xi)$. The limit of the joint degree distribution in this case can be determined but does not seem to give much insight.

Theorem 5 shows that the normalized degree distributions can be used for estimation of parameters only if $\lambda_{n}\to\infty$. If that is the case we can proceed as follows:

Let $\hat{\tau}_{1},\ldots,\hat{\tau}_{n}$ be the empirical quantiles of the normalized 1-degree distribution, and let $\hat{T}^{m}(\hat{\tau}_{k})$ be the $m$-degree of the vertex with normalized degree $\hat{\tau}_{k}$.

Fit smooth curves to $(\hat{\tau}_{k},\hat{T}^{m}(\hat{\tau}_{k}))$ viewed as observations of functions at $\hat{\tau}_{k}$, $k=1,\ldots,n$, for each $m$, and call these $\hat{T}^{m}(\cdot)$ (on $R$). By Theorem 5, $\hat{T}^{m}(t)\rightarrow T^{m-1}(\tau)(\tau^{-1}(t))$ for all $t$. If $T^{m-1}(\tau^{-1}(\cdot))$ are smooth, the convergence can be made uniform on compacts.

From the fitted functions $\hat{T}^{m}(\cdot)$, we can estimate the parameters of block models of any order consistently by replacing $\mathbf{v}_{m}$ in the proof of identifiability of block models by fitting the $\hat{T}^{m}(t)$ by $T^{m}(t)$ of the type specified by block models and then using the corresponding $\hat{\mathbf{v}}_{m}$. We only need the conditions of Theorem 5.

Computation of moment estimates and estimation of their variances

General acyclic graph moment estimates including those corresponding to patterns arising from $(\mathbf{k},\mathbf{l})$-wheels are computationally difficult. For $(k,l)$-wheels with small $k$ and $l$, we can use brute force counting, but unfortunately, the complexity of moment computation even for $(k,l)$-wheels appears to be $O(n\lambda_{n}^{k})$. Note that we need to count the sets of loopless paths of length $k$, $S_{i\mathbf{a}}$, for each $i$, where $S_{i\mathbf{a}}$ is the set of all paths of length $k$ originating at node $i$ which intersect another such path at $a_{1}<\cdots<a_{m}$, $1\leq m\leq k$, and $S_{i0}$ is the set of all paths of length $k$ from $i$ which do not intersect. The number of $(k,l)$-wheels with hub $i$ is then the number of $l$-tuples of such paths selected so that elements from $S_{i\mathbf{a}}$ appear at most once, with the remaining paths coming from $S_{i0}$. This is computationally nontrivial.

For very sparse graphs, however, intersecting paths can be ignored up to a certain order, and the wheel counts can be related to normalized $m$-degrees via a following approximation. If the conditions of Theorem 5 hold and $\lambda_{n}=o(n^{\alpha})$ for all $\alpha>0$, then

A similar formula holds for $\hat{\tau}_{\mathbf{k}\mathbf{l}}$.

The heuristic argument for (17) is that the expected number of paths of lengths $k$ from $i$ is $O(\lambda_{n}^{k})$. The expected number of pairs of such paths which intersect at least once is

if $\lambda_{n}=o(n^{\alpha})$ for all $\alpha>0$. Note that for $K$-block models this condition is not necessary for all $\alpha$, since we only need to count a finite number of $(k,l)$-wheels.

Estimation of variances of moment estimates even for $(\mathbf{k},\mathbf{l})$-wheels involve the counting of more complicated patterns. However, we propose the following bootstrap method:

Associate with each vertex $i$ the counts of $(\mathbf{k},\mathbf{l})$-wheels for which it is a hub, $S_{i}=\{n_{i\mathbf{k}\mathbf{l}}\mbox{: all }\mathbf{k},\mathbf{l}\}$, $i=1,\ldots,n$.

Sample without replacement $m$ vertices $\{i_{1},\ldots,i_{m}\}$, and let

For $R$ a $(\mathbf{k},\mathbf{l})$-wheel, define

Repeat this $B$ times to obtain $\check{P}_{1}^{*},\ldots,\check{P}_{B}^{*}$, and let

Then $\hat{\sigma}^{2}$ is an estimate of the variance of $\check{P}(R)$ if $\frac{m}{n}\rightarrow 0,m\rightarrow\infty$.

Discussion

Our Theorem 4 suggests that we might be able to construct consistent nonparametric estimates of $w_{\mathit{CAN}}$. That is, $\bm{\tau}_{M}=\{\tau_{\mathbf{k}\mathbf{l}}\colon\;|\mathbf{k}|\leq M,|\mathbf{l}|\leq M\}$ can be estimated at rate $n^{-1/2}$ for all $M<\infty$. But $\{\bm{\tau}_{M},M\geq 1\}$ determines $T_{w}$, and thus in principle we can estimate $T_{w}$ arbitrarily closely using $\{\hat{\tau}_{\mathbf{k}\mathbf{l}}\}$. This appears difficult both theoretically and practically. Theoretically, one difficulty seems to be that we would need to analyze the expectation of moments or degree distributions when the block model does not hold, which is doable. What is worse is that the passage to $w$ from moments is very ill-conditioned, involving first inversion via solution of the moment problem, and then estimation of eigenvectors and eigenvalues from a sequence of iterates $T_{w}(1),T_{w}^{2}(1)$, etc. If we assume $\lambda_{n}\rightarrow\infty$ so that we can use consistency of the degree distributions, we bypass the moment problem, but the eigenfunction estimation problem remains. A step in this direction is a result of Rohe, Chatterjee and Yu 2011 which shows that spectral clustering can be used to estimate the parameters of $k$ block models if $\lambda\rightarrow\infty$ sufficiently, even if $k\rightarrow\infty$ slowly. Unfortunately this does not deal with the problem we have just discussed, how to pick a block model which is a good approximation to the nonparametric model. For reasons which will appear in a future paper, smoothness assumptions on $w$ have to be treated with caution.

While $\lambda_{n}\rightarrow\infty$ has not occurred in practice in the past, networks with high average degrees are now appearing routinely. In particular, university Facebook networks have $\lambda$ of 15 or more with $n$ in the low thousands. In any case $\lambda_{n}\rightarrow\infty$ can still be useful as an asymptotic regime that can help us understand some general patterns, in the same way that the sample size going to infinity does in ordinary statistics. Note that most of the time we do not specify the rate of growth of $\lambda_{n}$, which can be very slow.

2 Adding covariates and directed graphs

In principle, adding covariates $X_{i}$ at each vertex or $X_{ij}$ at each edge simply converts our latent variable model, $w(\cdot,\cdot)$ into a mixed model

which can be turned into a logistic mixed model. Special cases of such models have been considered in the literature; see Hoff 2007 and references therein. We do not pursue this here. The extension of this model to directed graphs is also straightforward.

3 Dynamic models

Many models in the literature have been specified dynamically; see Newman 2010. For instance, the “preferential attachment” model constructs an $n$ graph by adding 1 vertex at a time, with edges of that vertex to previous vertices formed with probabilities which are functions of the degree of the candidate “old” vertex. If we let $n\rightarrow\infty$, we obtain models of the type we have considered whose $w$ function can be based on an integral equation for $\tau(\xi)$, our proxy for the degree of the vertex with latent variable $\xi$. We shall pursue this elsewhere also.

Appendix: Additional lemmas and proofs

[Proof of Proposition 1] The first line of (6) is immediate, conditioning on $\{\xi_{1},\ldots,\xi_{n}\}$. The second line in (6) follows by expanding the second product. Finally, (6) follows directly from the definitions of $P$ and $Q$.

The following standard result is used in the proof of Theorem 1.

Suppose $(U_{n},V_{n})$ are random elements such that,

in probability. Then $U_{n}$, $V_{n}$ are asymptotically independent,

Since $\lambda_{n}=(n-1)\rho_{n}$, the first term is

The second term is a $U$-statistic of order 2, which is well known to be $O(n^{-1})$. Thus, (9) follows in case (a).

To establish (10) and (b), we note that the conditional distribution of $\sqrt{n\lambda_{n}}T_{1}$ given $\bm{\xi}$ is that of a sum of independent random variables with conditional variance

Applying Lemma 1, we see that if $\lambda_{n}=O(1)$, (b) follows. On the other hand, if $\lambda_{n}\rightarrow\infty$, $\sqrt{n}T_{1}$ is negligible, and the Gaussian limit is determined by $T_{2}$.

The proof of (1) and (12) is similar. We shall decompose $\check{P}(R)$ as $U_{1}+U_{2}$ as we did $\frac{L}{n\lambda_{n}}$. If $\lambda_{n}\rightarrow\infty$, it is enough to prove that

replacing $\bar{D}$ by $n\rho_{n}=\lambda_{n}$ gives a perturbation of order $(n\lambda_{n})^{-{1/2}}=o(n^{-{1/2}})$.

In case (b), it is enough to show that the joint distribution of $\sqrt{n}((\hat{P}(R)-P(R))\rho_{n}^{-|R|},T_{1},T_{2})$ is Gaussian

in the limit, since in view of (9) and (10) we can apply the delta method to $\check{P}(R)$. Let $p\equiv|V(R)|$, $q\equiv|R|$. Each term in $\check{P}(R)$ is of the form

Condition on $\bm{\xi}=\{\xi_{1},\ldots,\xi_{n}\}$. Then terms $T(S)$, as above, yield

We begin by considering $\operatorname{Var}(U_{1}|\xi)$ which we can write as

where the sum ranges over all $S_{1}\sim R$, $S_{2}\sim R$.

If $E(S_{1})\cap E(S_{2})=\phi$ the covariance is 0. In general, suppose the graph $S_{1}\cap S_{2}$ has $c$ vertices and $d$ edges. Since $R$ is acyclic any subgraph is acyclic. By Corollary 3.2 of Chartrand, Lesniak and Behzad 1986 for every acyclic graph, $|V(S)|\geq|E(S)|+1$. Now,

There are $O(n^{2p-c})$ terms in (18) which have $c$ vertices in common. Therefore by (19) the total contribution of all such terms to $\operatorname{Var}(U_{1})$ is

if $\lambda_{n}\to\infty$. On the other hand

Thus, $\sqrt{n}(U_{1},U_{2})$ are jointly asymptotically Gaussian; see, for instance, Serfling 1980.

Since if $\lambda_{n}\to\infty$, $T_{1},U_{1}=o_{P}(n^{-{1/2}})$, the result follows if $\lambda_{n}\to\infty$. If $\lambda_{n}=O(1)$, we note that $\sqrt{n}(T_{1},U_{1})$ are sums of $q$ dependent random variables in the sense of Bulinski [see Doukhan 1994] and hence, given $\bm{\xi}$, are jointly asymptotically Gaussian. It is not hard to see that the limiting conditional covariance matrix is independent of $\xi$, as it was for $T_{1}$ marginally. By Lemma 1 again $(T_{1},U_{1})$ and $(T_{2},U_{2})$ are asymptotically independent and (a) and (b) follow.

Since $\rho=\frac{\lambda_{n}}{n}$ we obtain

For fixed $c\geq 1$ this is maximized by $d=\frac{c(c-1)}{2}$ and $n^{1-{2/c}}$ is maximized for $c\leq p$ by $c=p$.

[Proof of Theorem 2] Since $T$ corresponds to the canonical $h$,

Continuing we see that the $(K-1)(2K-1)$ moments $\{\tau_{kl}\colon\;2\leq k\leq K,1\leq l\leq 2K-1\}$ yield

for $j=1,\ldots,K$ where $v^{(0)}\equiv\pi$.

Given $\pi,v^{(1)},\ldots,v^{(K)}$ linearly independent, we can compute $F$ since by (21), we can write

where $V^{(1)}=(v^{(0)},\ldots,v^{(K-1)})^{T}$ and $V^{(2)}=(v^{(1)},\ldots,v^{(K)})^{T}$ and hence

Consistency and $\sqrt{n}$-consistency follow from Theorem 1 and the delta method. {proof}[Proof of Proposition 2] Note that

by the arithmetic/geometric mean and Minkowski inequalities. By Hölder’s inequality (23) is bounded by

is computable since we know $T_{w}1(\cdot)$ and the eigenfunction $\phi_{1}$ and eigenvalue $\lambda_{1}$. More generally, $T_{w}^{k}g_{1}^{(2)}$, $|g_{k-1}^{(2)}|$ can be similarly determined. Then, by the same argument as before, using 1 not orthogonal to $\phi_{2}$, we obtain $g_{k}^{(1)}\rightarrow_{L_{2}}\lambda_{2}\phi_{2}$ and $g_{k}^{(1)}/|g_{k}^{(1)}|\rightarrow_{L_{2}}\phi_{2}$. Now form $g_{0}^{(3)}\equiv 1-\lambda_{1}(1,\phi_{1})\phi_{1}-\lambda_{2}(1,\phi_{2})\phi_{2}$ and proceed as before, and continue to determine $\lambda_{k},\phi_{k}$ for all $k$. This and (15) complete the proof. {proof}[Proof of Theorem 5] Note first that (16) implies that the $M_{2}$ distance between $\hat{F}_{m}$ and the empirical distribution of $\{\bm{\theta}_{m}(\xi_{i})\}$ tends to 0. The first conclusion of the theorem now follows by the Glivenko–Cantelli theorem and the Law of Large Numbers.

where $w_{E(R)}=\prod_{(a,b)\in E(R)}w(\xi_{a},\xi_{b})$. Further, (Appendix: Additional lemmas and proofs) is a $U$-statistic of order $m$ under $|w_{2m}|<\infty$ and

Note that $R=\{(i,i_{1}),(i_{1},i_{2}),\ldots,(i_{m-1},j)\}$ is acyclic if all vertices are distinct. As in the proof of Theorem 1, all nonzero covariance terms in (Appendix: Additional lemmas and proofs) are of order $\rho^{2m-d}n^{2m-c}$ where $c\geq d$ since the intersection graphs all have $i$ in common but are otherwise acyclic. The largest order term corresponds to $c=d=m$, so that

where $C$ depends on $|w_{2m}|$ only. Thus (25) holds if $\lambda_{n}\rightarrow\infty$.

Acknowledgment

Thanks to Allan Sly for a helpful discussion.

References