Large-scale estimation of random graph models with local dependence
Sergii Babkin, Jonathan Stewart, Xiaochen Long, Michael Schweinberger
TL;DR
This paper introduces a scalable two-step estimation method for complex random graph models with local dependence, enabling efficient analysis of large networks through parallel computation and local structure estimation.
Contribution
It proposes a novel two-step estimation approach leveraging local graph structures, facilitating large-scale analysis of complex random graph models.
Findings
Effective estimation on networks with up to 10,000 nodes
Parallel computation significantly reduces estimation time
Successful application to large Amazon product network
Abstract
A class of random graph models is considered, combining features of exponential-family models and latent structure models, with the goal of retaining the strengths of both of them while reducing the weaknesses of each of them. An open problem is how to estimate such models from large networks. A novel approach to large-scale estimation is proposed, taking advantage of the local structure of such models for the purpose of local computing. The main idea is that random graphs with local dependence can be decomposed into subgraphs, which enables parallel computing on subgraphs and suggests a two-step estimation approach. The first step estimates the local structure underlying random graphs. The second step estimates parameters given the estimated local structure of random graphs. Both steps can be implemented in parallel, which enables large-scale estimation. The advantages of the two-step…
| Bayesian approach | Variational approach | Spectral clustering | |
|---|---|---|---|
| , , balanced | 42,579.6 | 23.83 | .18 |
| , , unbalanced | 32,480.9 | 24.00 | .18 |
| Term | Estimate | S.E. | Estimate | S.E. |
| Within-block edges | -1.410 | |||
| Within-block degrees (base parameter) | .742 | |||
| Within-block degrees (decay parameter) | .910 | |||
| Within-block shared partners (base parameter) | .303 | |||
| Within-block shared partners (decay parameter) | 1.106 | |||
| Between-block edges |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Large-scale estimation of random graph models with local dependence
Sergii Babkin
Microsoft
Jonathan Stewart
Department of Statistics, Rice University
Xiaochen Long
Department of Statistics, Rice University
Michael Schweinberger
Department of Statistics, Rice University
Abstract
A class of random graph models is considered, combining features of exponential-family models and latent structure models, with the goal of retaining the strengths of both of them while reducing the weaknesses of each of them. An open problem is how to estimate such models from large networks. A novel approach to large-scale estimation is proposed, taking advantage of the local structure of such models for the purpose of local computing. The main idea is that random graphs with local dependence can be decomposed into subgraphs, which enables parallel computing on subgraphs and suggests a two-step estimation approach. The first step estimates the local structure underlying random graphs. The second step estimates parameters given the estimated local structure of random graphs. Both steps can be implemented in parallel, which enables large-scale estimation. The advantages of the two-step estimation approach are demonstrated by simulation studies with up to 10,000 nodes and an application to a large Amazon product recommendation network with more than 10,000 products.
keywords:
exponential-family models \seplatent structure models \sepstochastic block models \sepvariational methods \sepEM algorithms \sepMM algorithms
††journal: Computational Statistics & Data Analysis
1 Introduction
The statistical analysis of network data is an emerging area in statistics (Kolaczyk, 2009). Network data arise in the study of insurgent and terrorist networks, contact networks facilitating the spread of infectious diseases (e.g., corona viruses such as SARS and MERS; Ebola; HIV), social networks, and the World Wide Web.
Many models of network data have been proposed, as described in recent review papers (Fienberg, 2012; Hunter et al., 2012; Salter-Townshend et al., 2012; Smith et al., 2019; Schweinberger et al., 2020; Hoff, 2020). Among the plethora of models, two broad classes of models can be distinguished: models with latent structure, including stochastic block models (e.g., Nowicki and Snijders, 2001; Bickel and Chen, 2009; Rohe et al., 2011) and latent space models (e.g., Hoff et al., 2002; Handcock et al., 2007; Sewell and Chen, 2015; Smith et al., 2019); and exponential-family models of random graphs (Lusher et al., 2013; Harris, 2013), including Erdős and Rényi random graphs, logistic regression models, and models resembling Markov random fields in spatial statistics (Besag, 1974). Both have advantages and disadvantages, as one of the pioneers of statistical network analysis observed:
“I expect that, especially for modelling larger networks (with, say, a few hundred or more nodes), the latent space models will not be able to represent network structures as expressed by subgraph counts sufficiently well and the exponential random-graph models will not be able to represent the cohesive structure sufficiently well. Models that combine important features of these two approaches may be the next generation of social network models” (Snijders, 2007, p. 324).
In other words, latent structure models capture who is close to whom, but are not flexible models of many network phenomena of interest (despite the fact that latent space models induce a weak tendency towards transitivity, see Hoff et al., 2002). Well-posed exponential-family models of random graphs can capture a vast range of network phenomena of interest (including, but not limited to, transitivity), but the underlying assumptions of many of those models make more sense in small networks than large networks. A case in point is the Markov random graphs of Frank and Strauss (1986). These models assume that possible edges and of pairs of nodes and are independent conditional on the rest of the graph provided and do not overlap, but are otherwise dependent. As a consequence, each possible edge depends on other possible edges, where is the number of nodes. If Markov random graphs were applied to online social networks such as Facebook, then each possible friendship would depend on billions of other possible friendships. Such dependence assumptions are implausible, and when coupled with homogeneity assumptions regarding parameters can give rise to the well-studied issue of model near-degeneracy (Handcock, 2003; Schweinberger, 2011; Chatterjee and Diaconis, 2013; Mele, 2017).
One class of next-generation random graph models was introduced in Schweinberger and Handcock (2015), which combines important features of stochastic block models and exponential-family models, with the goal of retaining the strengths of both of them while reducing the weaknesses of each of them. The basic idea is that a set of nodes is partitioned into blocks, and edges among nodes within and between blocks are governed by exponential-family models of random graphs with local dependence within blocks. A simple example is a model where edges between blocks are independent Bernoulli random variables, whereas edges within blocks are generated by an exponential-family model which encourages triangles within blocks but ensures that, for each pair of nodes within a block, the added value of additional edges and triangles decays. Such models induce local dependence within blocks and the overall dependence induced by the model is weak provided the blocks are not too large. We have shown elsewhere that such models are well-behaved—in contrast to the infamous triangle model first studied by Strauss (1986), Jonasson (1999), Häggström and Jonasson (1999), and others (e.g., Chatterjee and Diaconis, 2013)—and that statistical inference is possible and supported by statistical theory: e.g., we have established concentration and consistency results for canonical and curved exponential-family models of random graphs with local dependence under the assumption that the blocks are known and the sizes of the blocks are similar in a well-defined sense (Schweinberger and Stewart, 2020). In addition, when the blocks are unknown, the block memberships of most nodes can be recovered with high probability under weak dependence and smoothness conditions (Schweinberger, 2020).
While some progress has been made in terms of statistical theory, computing remains challenging. When the block structure is known (which is the case in multilevel networks, e.g., networks consisting of units of armed forces, Lazega and Snijders, 2016), statistical inference for parameters can rely on existing methods for estimating parameters of exponential-family models of random graphs (e.g., Strauss and Ikeda, 1990; Snijders, 2002; Hunter and Handcock, 2006; Caimo and Friel, 2011; Hummel et al., 2012; Okabayashi and Geyer, 2012; Atchade et al., 2013; Jin and Liang, 2013; Thiemichen and Kauermann, 2017; Krivitsky, 2017; Byshkin et al., 2018; Tan and Friel, 2018). If the block structure is unknown, however, it needs to be estimated based on the observed network. The recovery of unknown block structure resembles the recovery of unknown block structure in stochastic block models (e.g., Nowicki and Snijders, 2001; Bickel and Chen, 2009; Bickel et al., 2011; Rohe et al., 2011; Choi et al., 2012; Celisse et al., 2012; Priebe et al., 2012; Amini et al., 2013; Rohe et al., 2014; Zhang and Zhou, 2016; Binkiewicz et al., 2017; Gao et al., 2017). However, the recovery of unknown block structure is more challenging for exponential-family models with local dependence than for stochastic block models. The main challenge is the intractability of the complete-data likelihood function, i.e., the likelihood function given an observation of the network as well as the block structure. The intractability of the complete-data likelihood function is rooted in the intractability of the normalizing constants of within-block probability mass functions, which stems from the local dependence within blocks.
We present here a tractable approximation of the likelihood function, leveraging concentration results for random graphs with local dependence. Based on the approximation of the likelihood function, we propose a two-step estimation approach that exploits the local structure of random graphs with local dependence for the purpose of local computing. The first step estimates the block structure and decomposes the random graph into subgraphs with local dependence. The decomposition of the random graph relies on approximations of the likelihood function supported by theoretical results. The second step estimates parameters given the estimated block structure by using Monte Carlo maximum likelihood methods (Hunter and Handcock, 2006) or maximum pseudolikelihood methods (Strauss and Ikeda, 1990). Both steps can be implemented in parallel, which enables large-scale estimation on multi-core computers or computing clusters. We demonstrate the advantages of the two-step estimation approach by simulation studies with up to 10,000 nodes and an application to a large Amazon product recommendation network with more than 10,000 products.
The remainder of our paper is organized as follows. Section 2 introduces models. Section 3 discusses likelihood-based inference based on approximations of the likelihood function motivated by theoretical results, and Section 4 takes advantage of such approximations to estimate models. Section 5 presents simulation results and Section 6 presents an application.
2 Models
We consider random graphs with a set of nodes and a set of edges . Here, edges are regarded as random variables , where if nodes and are connected by an edge and otherwise. We focus on undirected random graphs without self-edges, although the methods introduced here can be extended to directed random graphs. Throughout, we denote by the set of possible edges of a random graph, and by the sample space of .
Since the pioneering research of Holland and Leinhardt (1970, 1972, 1976) in the 1970s, it is known that network data are dependent data. In small networks, it is very well possible that each possible edge depends on many other possible edges, and might in fact depend on all other possible edges. In large networks, however, it is not credible that each possible edge depends on all other possible edges: e.g., in networks with 10,000 nodes—including the network used in Section 6—it is implausible that each possible edge depends on all 50,000,000 other possible edges. Instead, it is tempting to believe that the dependence among edges in large networks is local, in the sense that each edge depends on a subset of other possible edges. We consider here a simple form of local dependence, following Schweinberger and Handcock (2015). Assume that is partitioned into subsets of nodes , called blocks, and let be the block memberships of nodes, where if node belongs to block and otherwise. We henceforth denote by the within-block edges among nodes in block () and by the between-block edges among nodes in block and nodes in block ().
Definition. Local dependence. A model of a random graph satisfies local dependence if the probability mass function of can be factorized as follows:
[TABLE]
where
- •
is the probability mass of graph , parameterized by ;
- •
is the probability mass of within-block subgraph , parameterized by ();
- •
is the probability mass of between-block edge , parameterized by ();
- •
the parameter vector consists of the subvectors () and ().
The map depends on the model. In general, may be a linear or non-linear function of a parameter vector of dimension . An example is given by the curved exponential-family parameterizations used in Section 6, where is a non-linear function of a parameter vector of dimension . Well-chosen curved exponential-family parameterizations help ensure that the added value of additional edges, triangles, and other subgraph configurations decays (Stewart et al., 2019), which is plausible in many applications and can improve the in-sample performance of models (Hunter et al., 2008) as well as the out-of-sample performance of models (Stewart et al., 2019). The factorization of the probability mass function of models with local dependence has at least three implications. First, edges among nodes in block can depend on other edges in block (). Second, the factorization implies that edges among nodes in block do not depend on edges that involve nodes in other blocks (). Third, edges between blocks are independent. In other words, the dependence is local in the sense that it is confined to within-block subgraphs.
We introduce here two examples of models with local dependence, to be used as motivating examples throughout Sections 3—5.
Example 1. Stochastic block model. An important special case of models with local dependence is given by stochastic block models (e.g., Nowicki and Snijders, 2001; Bickel and Chen, 2009; Rohe et al., 2011). Stochastic block models assume that possible edges between nodes and in blocks and are independent Bernoulli random variables, where the probability of an edge depends on the blocks and of nodes and , respectively. Stochastic block models are special cases of models with local dependence, having within-block probability mass functions of the form
[TABLE]
where
[TABLE]
The between-block probability mass functions are of the form
[TABLE]
where
[TABLE]
As a consequence, the joint probability mass function of random graph is proportional to
[TABLE]
However, while stochastic block models are popular, the assumption that edges within and between blocks are independent is restrictive, because edges among nodes that are close—in the sense of being members of the same block—may very well be dependent.
Example 2. Model with local dependence. To allow edges within blocks to be dependent, we build on the stochastic block model described above, but tilt the within-block probability mass function of the stochastic block model as follows:
[TABLE]
where and and is an indicator function, which is if nodes and are both connected to one or more other nodes in block , and is [math] otherwise. If and , the edge between nodes and is called transitive, because and form transitive triples with other nodes, which are known as triangles in the random graph literature. The first term of the within-block probability mass function shown above corresponds to the within-block probability mass function of the stochastic block model. The second term tilts the within-block probability mass function of the stochastic block model: If , the model rewards within-block subgraphs with transitive edges, whereas penalizes them, and neither rewards nor penalizes them—in which case the model reduces to the stochastic block model.
Using the same between-block probability mass functions as the stochastic block model, the joint probability mass function of random graph is proportional to
[TABLE]
The resulting model can capture an excess in the expected number of transitive edges within blocks, relative to the stochastic block model. To demonstrate, note that the resulting model is an exponential-family model and—by exponential-family theory (Brown, 1986, Corollary 2.5, p. 37)—the expected number of transitive edges in block satisfies
[TABLE]
where is the number of transitive edges in block and is the expectation of . In other words, the expected number of transitive edges in block is greater under the model with than under the stochastic block model with , assuming that both have the same edge parameters (). Therefore, the model with local dependence can capture an excess in the expected number of transitive edges within blocks, relative to the stochastic block model. In addition, models with local dependence can capture excesses in the expected number of other subgraph statistics within blocks by adding suitable model terms.
How models with within-block edge and transitive edge terms differ from the “triangle model”
It is worth noting that the model with local dependence induced by within-block edge and transitive edge terms differs in two important ways from the infamous triangle model with edge and triangle terms, which has been known to be ill-behaved since the 1980s (Strauss, 1986; Jonasson, 1999; Schweinberger, 2011; Chatterjee and Diaconis, 2013):
- •
The model with local dependence restricts dependence among edges to subsets of nodes, i.e., blocks. As long as the blocks are not too large, the overall dependence induced by the model is weak. By contrast, the triangle model does not restrict dependence to subsets of nodes, leading to undesirable behavior in large graphs.
- •
Within blocks, the model with local dependence and positive within-block transitive edge parameters assumes that, for each pair of nodes, the value added by the first triangle to the log odds of the conditional probability of an edge is positive, but additional triangles add nothing. By contrast, the triangle model assumes that, for each pair of nodes, each additional triangle has the same added value, which is unreasonable and results in undesirable behavior in large graphs.
These sensible assumptions ensure that models with local dependence have more desirable properties than the triangle model (Schweinberger and Stewart, 2020).
Exponential-family representations of models
It is worth noting that Models 1 and 2 can be represented in exponential-family form:
[TABLE]
where denotes the inner product of a vector of natural parameters and a vector of sufficient statistics , and ensures that sums to . While exponential-family representations are not needed to introduce Models 1 and 2, the properties of exponential families facilitate the theoretical results in Section 3.2, concerned with approximations of likelihood functions.
3 Likelihood-based inference
While it is tempting to base statistical inference for the unknown block structure and the unknown parameter vector on the likelihood function, likelihood-based inference for models with local dependence is challenging. The main reason is that the probability mass function is intractable, because the within-block probability mass functions are intractable (). The intractability of is rooted in the fact that its normalizing constant is a sum over all possible within-block subgraphs of block , which cannot be computed unless is small, i.e., unless ().
To facilitate likelihood-based inference, we introduce tractable approximations of the intractable probability mass function in Section 3.1 and support them by theoretical results in Section 3.2. A statistical algorithm that takes advantage of such approximations is introduced in Section 4.
3.1 Approximate likelihood functions: motivation
Suppose that we want to estimate both and . It is natural to estimate them by using an iterative algorithm that cycles through updates of and as follows:
- Step 1: Update given .
- Step 2: Update given .
The algorithm sketched above is generic and cannot be used in practice, but regardless of which specific algorithm is used—whether EM, Monte Carlo EM, variational EM, Bayesian Markov chain Monte Carlo, or other algorithms—most of them have in common that Step 1 is either infeasible or time-consuming, whereas Step 2 is less problematic than Step 1.
Step 1
Step 1 is either infeasible or time-consuming, because the probability mass function is intractable. To demonstrate, consider a Bayesian Markov chain Monte Carlo algorithm that updates given by Gibbs sampling. Gibbs sampling of turns out to be infeasible, because the full conditional distributions of depend on the intractable within-block probability mass functions (). One could approximate them by Monte Carlo samples of within-block subgraphs, but such approximations may not generate Markov chain Monte Carlo samples from the target distribution (Liang et al., 2016) and are problematic on computational grounds:
- •
Using Monte Carlo approximations of within-block probability mass functions is infeasible when the number of nodes is large, because such approximations are needed for each update of each of the block memberships .
- •
Worse, the block memberships cannot be updated in parallel, because the block membership of one node depends on the block memberships of other nodes.
As a consequence, Step 1 is infeasible when is large.
Step 2
Step 2 is less problematic than Step 1. While the probability mass function is intractable and may have to be approximated by Monte Carlo methods (Hunter and Handcock, 2006), such Monte Carlo approximations are needed once to update given , whereas Monte Carlo approximations are needed times to update given one by one. In addition, the probability mass function decomposes into between- and within-block probability mass functions () and hence within-block probability mass functions can be approximated in parallel, i.e., on multi-core computers or computing clusters.
Approximations
To enable feasible updates of given when is large, we are interested in approximating the intractable probability mass function by a tractable probability mass function. To do so, we focus on models with between-block edge terms and within-block edge terms along with additional within-block terms that induce local dependence within blocks. We denote the vector of between- and within-block edge parameters by and the vector of all other within-block parameters by . We assume that reduces the model to the stochastic block model. An example is given by the model with between-block edge terms and within-block edge and transitive edge terms described in Section 2: the parameter vector consists of the between-block edge parameters () and the within-block edge parameters (), whereas the parameter vector consists of the within-block transitive edge parameters (). If the parameter vector is set to , the model reduces to the stochastic block model, under which edges within and between blocks are independent.
Such models have two useful properties:
- •
The probability mass functions and impose the same probability law on between-block subgraphs.
- •
The probability mass function is tractable, because edges between and within blocks are independent under all .
We henceforth approximate by , which corresponds to the probability mass function of the stochastic block model.
The idea underlying the approximation is that when the blocks are not too large, most pairs of nodes are not members of the same block. Since and impose the same probability law on possible edges between pairs of nodes that are not members of the same block, and agree on most of the random graph. Therefore, can be approximated by for the purpose of updating given . Suppose, e.g., that we consider to update given by replacing by . We may decide to do so if the loglikelihood ratio
[TABLE]
is large: e.g., the acceptance probability of Metropolis-Hastings algorithms depends on the loglikelihood ratio above. If can be approximated by , we can base the decision on rather than , because
[TABLE]
Therefore, as long as
[TABLE]
is small, we have
[TABLE]
The advantage of approximating by is that there exist methods for stochastic block models to estimate the block structure from large random graphs (e.g., Daudin et al., 2008; Rohe et al., 2011; Amini et al., 2013; Vu et al., 2013). We take advantage of such methods in Section 4, but we first shed light on the conditions under which is small.
3.2 Approximate likelihood functions: theoretical results
We show that updates of given can be based on rather than by showing that
[TABLE]
is small with high probability when the blocks are not too large.
We start with a special case in Theorem 3.1 and then present more general results in Theorem 3.2. To prepare the ground for Theorems 3.1 and 3.2, let be the set of all block structures, and denote by the size of the largest block under . Let be any subset of block structures that includes the data-generating block structure , and denote by the size of the largest block among all block structures , so that for all . We allow the size of the largest block to increase as a function of the number of nodes : e.g., may be a constant multiple of or (). The size of the largest data-generating block is denoted by .
Theorem 3.1**.**
Consider the model with between-block edge terms and within-block edge and transitive edge terms described in Section 2. Let be the parameter space, be a compact subset of (), and be a compact subset of (). Choose as small as desired. Then there exists a universal constant such that, for all , with at least probability ,
[TABLE]
The proof of Theorem 3.1 can be found in the supplement. The basic idea underlying Theorem 3.1 is that the deviation cannot be too large when the blocks are not too large, because most pairs of nodes are not members of the same block and and impose the same probability law on possible edges between pairs of nodes that are not members of the same block.
Theorem 3.1 implies that the largest deviation of the loglikelihood function under the unrestricted model and the restricted model is smaller than with high probability, provided the blocks are not too large. To see that, observe that implies
[TABLE]
where is a function of but is not a function of the number of nodes . In other words, as long as the size of the largest block in satisfies , the maximum deviation is much smaller than with high probability. As a consequence, the largest deviation of the loglikelihood function under the unrestricted model and the restricted model is small with high probability—note that in dense random graphs many quantities are of order , so quantities of order less than may be considered small. For example, consider Bernoulli random graphs, under which possible edges are independent Bernoulli random variables (Erdős and Rényi, 1959, 1960). If a Bernoulli random graph is dense in the sense that does not decrease as a function of the number of nodes , then the expected number of edges is of order ,
[TABLE]
and so is the expected loglikelihood function of the natural parameter ,
[TABLE]
Other quantities are likewise of order in dense random graphs, so quantities of order less than may be considered small. Thus, the largest deviation of the loglikelihood function under the unrestricted model and the restricted model is small with high probability.
Last, but not least, note that implies that the number of blocks must satisfy . In other words, the number of blocks needs to grow as function of the number of nodes . It is worth noting that in the special case of stochastic block models, it is known that is allowed to grow as fast as in dense-graph settings (Zhang and Zhou, 2016).
It turns out that the result in Theorem 3.1 is not limited to the model with between-block edge terms and within-block edge and transitive edge terms, but is a special case of more general results. To introduce these more general results, we make the following assumptions. We assume that the probability mass function can be represented in exponential-family form,
[TABLE]
and satisfies local dependence as defined in Section 2; note that all models used in our paper can be represented in exponential-family form. We assume that and that is a subset of the interior of the natural parameter space of the exponential family (Brown, 1986). Let be the expectation under the data-generating parameter vector , where denotes the data-generating value of . We denote by the Hamming metric, which is defined by
[TABLE]
where is if and is [math] otherwise. The main assumptions can then be stated as follows.
- [C.1
] There exists such that, for all and all ,
[TABLE]
- [C.2
] There exists such that, for all and all ,
[TABLE]
where , , and denote the subvectors of , , and corresponding to the subgraph between blocks and () or the subgraph of block () and is a compact subset of ().
Conditions [C.1] and [C.2] are smoothness conditions: [C.1] states that the inner product of natural parameters and sufficient statistics must not be too sensitive to changes of edges, whereas [C.2] states that the inner product of between- and within-block natural parameters and expected sufficient statistics must not be too sensitive to changes of parameters.
An example of a model violating condition [C.1] is a model containing a sufficient statistic of the form
[TABLE]
Under such models, adding a single edge can change the inner product of natural parameters and sufficient statistics by a multiple of . That would violate [C.1], because [C.1] assumes that changing a single edge changes the inner product by at most .
An example of a model violating condition [C.2] is a model with within-block sufficient statistics that count the number of triangles within blocks,
[TABLE]
with within-block parameters of the form (). If all nodes belong to block , then the inner product of the block’s within-block natural parameter and expected sufficient statistic vector is a multiple of , where and are two possible values of the within-block triangle parameter of block . As a result, [C.2] would be violated, because [C.2] requires the inner product to be at most .
However, the fact that some model specifications violate conditions [C.1] and [C.2] is not necessarily a major concern, for two reasons. First, the specifications that violate conditions [C.1] and [C.2] are, more often than not, of limited interest in applications: e.g., we are not aware of any good reason for using the sufficient statistic in the first example, and the triangle term in the second example is known to induce model near-degeneracy in large networks and may therefore not be useful in practice (Strauss, 1986; Jonasson, 1999; Schweinberger, 2011; Chatterjee and Diaconis, 2013). Second, while some ill-posed specifications do not satisfy conditions [C.1] and [C.2], there are many well-posed specifications that do satisfy them: e.g., conditions [C.1] and [C.2] are satisfied by the model with between-block edge and within-block edge and transitive edge terms, which we verify in the proof of Theorem 3.1. In addition, conditions [C.1] and [C.2] cover the models with size-dependent parameterizations used in Sections 5 and 6.
The following result, Theorem 3.2, is a generalization of Theorem 3.1.
Theorem 3.2**.**
Consider a model with local dependence satisfying conditions [C.1] and [C.2]. Choose as small as desired. Then there exists a universal constant such that, for all , with at least probability ,
[TABLE]
The proof of Theorem 3.2 can be found in the supplement. An application of Theorem 3.2 to the model with between-block edge terms and within-block edge and transitive edge terms can be found in Theorem 3.1.
Trade-off between and the recovery of block structure
It is worth noting that the upper bound on the sizes of blocks cannot be too small, because it is not possible to recover the block structure with high probability when the blocks are too small. However, Zhang and Zhou (2016) showed that under stochastic block models the number of blocks is allowed to grow as fast as in dense-graph settings, which implies that the sizes of the blocks can be as small as . These important issues are studied in more depth in, e.g., Zhang and Zhou (2016) and Gao et al. (2017).
4 Two-step estimation approach
We propose a two-step estimation approach that takes advantage of the theoretical results of Section 3 and enables large-scale estimation of models with local dependence.
To describe the two-step estimation approach, assume that is the observed value of a random variable with distribution
[TABLE]
It is natural to base statistical inference on the observed-data likelihood function
[TABLE]
The problem is that \mathscr{L}(\bm{\theta},\mbox{\boldmath\pi}) is intractable, because is intractable and the set contains elements.
The first problem can be solved by taking advantage of the theoretical results of Section 3, which suggest that can be approximated by provided that the blocks are not too large. A complication is that and may not be close when the block structure is not contained in , in which case some of the blocks can be large and the within-block models of large blocks can induce strong dependence. In the worst case, all nodes belong to a single block, in which case can induce strong dependence throughout the random graph while induces no dependence at all. However, the basic inequality
[TABLE]
suggests that as long as the event is a rare event in the sense that \mathbb{P}_{\mbox{\boldmath\pi}}(\bm{Z}\in\mathbb{Z}\setminus\mathbb{S}) is close to [math], \mathscr{L}(\bm{\theta},\mbox{\boldmath\pi}) can be approximated by \mathscr{L}(\bm{\theta}_{1},\bm{\theta}_{2}=\bm{0},\mbox{\boldmath\pi}):
[TABLE]
The assumption that is a rare event—i.e., the probabilities are small—makes sense in a wide range of applications, because communities in real-world networks tend to be small (see, e.g., the discussion of Rohe et al., 2011). Therefore, as long as is a rare event, we can base statistical inference concerning the block structure on \mathscr{L}(\bm{\theta}_{1},\bm{\theta}_{2}=\bm{0},\mbox{\boldmath\pi}) rather than \mathscr{L}(\bm{\theta},\mbox{\boldmath\pi}). To simplify the notation, we write henceforth \mathscr{L}(\bm{\theta}_{1},\mbox{\boldmath\pi}) instead of \mathscr{L}(\bm{\theta}_{1},\bm{\theta}_{2}=\bm{0},\mbox{\boldmath\pi}).
The second problem can be solved by methods developed for stochastic block models, because \mathscr{L}(\bm{\theta}_{1},\mbox{\boldmath\pi}) is the observed-data likelihood function of a stochastic block model. There are many stochastic block model methods that could be used, such as profile likelihood (Bickel and Chen, 2009), pseudo-likelihood (Amini et al., 2013), spectral clustering (Rohe et al., 2011), and variational methods (Daudin et al., 2008; Vu et al., 2013). Among these methods, we found that the variational methods of Vu et al. (2013) work best in practice. In addition, the variational methods of Vu et al. (2013) have the advantage of being able to estimate stochastic block models from networks with hundreds of thousands of nodes due to a running time of for sparse random graphs and for dense random graphs (Vu et al., 2013). Some consistency and asymptotic normality results for variational methods for stochastic block models were established by Celisse et al. (2012) and Bickel et al. (2013).
Variational methods approximate \ell(\bm{\theta}_{1},\mbox{\boldmath\pi})=\log\mathscr{L}(\bm{\theta}_{1},\mbox{\boldmath\pi}) by introducing an auxiliary distribution with support and lower bound \ell(\bm{\theta}_{1},\mbox{\boldmath\pi}) by using Jensen’s ineqality:
[TABLE]
Each auxiliary distribution with support gives rise to a lower bound on \ell(\bm{\theta}_{1},\mbox{\boldmath\pi}). To choose the best auxiliary distribution—i.e., the auxiliary distribution that gives rise to the tightest lower bound on \ell(\bm{\theta}_{1},\mbox{\boldmath\pi})—we choose a family of auxiliary distributions and select the best member of the family. In practice, an important consideration is that the resulting lower bound is tractable. Therefore, we confine attention to a family of auxiliary distributions under which the resulting lower bounds are tractable. A natural choice is given by a family of auxiliary distributions under which the block memberships are independent:
[TABLE]
By the independence of block memberships under the auxiliary distribution, we obtain the following tractable lower bound on \ell(\bm{\theta}_{1},\mbox{\boldmath\pi}):
[TABLE]
where denotes the marginal probability mass function of and denotes the block memberships of all nodes excluding nodes and .
In practice, we obtain the best lower bound on \ell(\bm{\theta}_{1},\mbox{\boldmath\pi}) by maximizing \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) with respect to . Direct maximization of \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) with respect to is possible but inconvenient, because \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) contains products of and . As a consequence, a fixed-point update of depends on other terms and hence fixed-point updates tend to be time-consuming and get stuck in local maxima, as demonstrated by Vu et al. (2013). Vu et al. (2013) proposed an elegant approach to alleviating the problem by using minorization-maximization methods (Hunter and Lange, 2004). Such methods construct a minorizing function that approximates \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) but is easier to maximize than \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}). A function M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}) of minorizes \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) at point at iteration of an iterative algorithm for maximizing \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) if
[TABLE]
where \bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)} are fixed. In other words, M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}) is bounded above by \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) for all and touches \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) at . As a result, increasing M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}) with respect to increases \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}). Vu et al. (2013) showed that the following function minorizes \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) at point :
[TABLE]
The minorizing function M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}) is easier to maximize than \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}), because it replaces the products of and by sums of and . An additional advantage is that the maximization of M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}) amounts to quadratic programming problems, which can be solved in parallel.
We therefore propose a two-step estimation approach as described in Table 1. We discuss the two steps below and conclude with some comments on parallel computing.
Step 1
The first step estimates based on . We do so by increasing M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}) with respect to subject to the constraints and (). We increase rather than maximize M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}), because maximizing M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}) is more time-consuming and algorithms maximizing M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}) are more prone to end up in local maxima than algorithms increasing M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}). Since \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) and M(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi},\bm{\alpha}^{(t)}) depend on and and both are unknown, we iterate between updates of and updates of and . The updates of and are based on maximizing \hat{\ell}(\bm{\alpha};\bm{\theta}_{1},\mbox{\boldmath\pi}) with respect to and and are identical to the updates of Vu et al. (2013), because reduces the model to a stochastic block model. As a convergence criterion, we use
[TABLE]
where is a small constant. Upon convergence, we estimate the block memberhip indicators by if and otherwise (, ), where denotes the final value of .
Step 2
We estimate given by Monte Carlo maximum likelihood methods (Hunter and Handcock, 2006) or maximum pseudolikelihood methods (Strauss and Ikeda, 1990). Monte Carlo maximum likelihood methods exploit the fact that the loglikelihood function induced by , which is defined by
[TABLE]
can be written as
[TABLE]
where is a fixed parameter vector (e.g., may be an educated guess of ). In general, the expectation is intractable, but it can be estimated by a Monte Carlo sample average based on a Monte Carlo sample of graphs generated under . Therefore, we can approximate by
[TABLE]
where is a Monte Carlo approximation of based on a Monte Carlo sample of graphs generated by using . Hence given can be estimated by
[TABLE]
Additional details on Monte Carlo maximum likelihood methods can be found in Hunter and Handcock (2006). We note that the local dependence of the model facilitates parallel computing, which is discussed in the following paragraph. Standard errors of can be based on the estimated Fisher information matrix, although such standard errors are conditional on the estimated block structure and therefore do not reflect the uncertainty about . A parametric bootstrap approach would be an interesting alternative to capturing the additional uncertainty due to , but would be time-consuming.
An alternative is to estimate given by maximum pseudolikelihood estimators (Strauss and Ikeda, 1990). Maximum pseudolikelihood estimators were invented by Besag (1974) to sidestep intractable likelihood functions of Markov random fields in spatial statistics, and are known to be consistent estimators of Markov random fields with exponential-family parameterizations (e.g., Comets, 1992). Strauss and Ikeda (1990) adapted them to exponential-family random graph models. While believed to be inferior to maximum likelihood estimators when the dependence among edges is strong and propagates throughout the random graph (e.g., van Duijn et al., 2009), maximum pseudolikelihood estimators seem to perform well when the dependence is weak, e.g., when the dependence among edges is confined to non-overlapping or overlapping blocks (Stewart and Schweinberger, 2020). The maximum pseudolikelihood estimator is defined as
[TABLE]
where
[TABLE]
Here, denotes excluding , and denotes the conditional probability of given . Since the conditional distributions of given the rest of the random graph are Bernoulli distributions, computing maximum pseudolikelihood estimators consumes less time than computing Monte Carlo maximum likelihood estimators.
Parallel computing
In Step 1, the maximization of the minorizing function amounts to quadratic programming problems, which can be solved in parallel. In Step 2, the local dependence induced by the model implies that the contributions of the between- and within-block subgraphs to the loglikelihood function and its gradient and Hessian can be computed in parallel. Hence both steps can be implemented in parallel, which suggests that the two-step likelihood-based method can be used on a massive scale as long as the blocks are not too large and multi-core computers or computing clusters are available.
5 Simulation results
We assess the performance of the two-step estimation approach by conducting multiple simulation studies, with the number of nodes ranging from 30 to 10,000, the number of blocks ranging from 3 to 100, and the sizes of blocks ranging from 5 to 100:
- I.
A small-scale simulation study to compare the block recovery of the gold standard, the Bayesian approach of Schweinberger and Handcock (2015), to the two-step estimation approach. Since the Bayesian approach is too time-consuming to be applied to large networks, we use small networks with nodes and blocks. The blocks consist of 10 nodes (balanced case) or 5, 10, 15 nodes (unbalanced case).
- II.
A large-scale simulation study to assess the block recovery in Step 1 of the two-step estimation approach, using blocks of size .
- III.
A large-scale simulation study to assess how the block recovery in Step 1 of the two-step estimation approach is affected by the sparsity of between-block subgraphs, using blocks of size .
- IV.
A large-scale simulation study to assess the parameter recovery in Step 2 of the two-step estimation approach, using blocks of size .
In Step 1 of the two-step estimation approach, we use the variational approach described in Section 4. We compare the variational approach to the spectral clustering method described in Lei and Rinaldo (2015). Spectral clustering is an alternative to the variational approach. In Step 2 of the two-step estimation approach, we use maximum pseudolikelihood estimators, which are more scalable than Monte Carlo maximum likelihood estimators and facilitate simulation studies with up to 10,000 nodes. In each scenario, we generate 500 graphs from the model having between-block edge terms and within-block edge and transitive edge terms, as described in Section 2. To select sensible values of the parameter vector , note that the probabilistic behavior of random graphs with local dependence is sensitive to the choice of parameter values, and so is the block recovery. The same applies to stochastic block models: e.g., when the probabilities of edges within and between blocks are too low, we may not be able to recover the block structure with high probability (Zhang and Zhou, 2016; Gao et al., 2017). We have therefore selected the parameter vector based on the following considerations: We would like to ensure that, with high probability, the model generates graphs that resemble real-world networks in terms of sufficient statistics . In principle, we could select by inspecting the expectation of . The problem is that the expectation is not available in closed form. We have therefore selected based on simulating graphs, and checking whether the simulated graphs resemble real-world networks.
In simulation study I, to allow blocks of different sizes to have different parameters, we use size-dependent between-block edge parameters () and within-block edge and transitive edge parameters and (), where is the size of block under . The size-dependent parameterization is motivated by the sparsity of random graphs: e.g., if edges are independent Bernoulli random variables, it is tempting to believe that there exist constants and such that the expected number of edges of each node is bounded above by , because real-world networks are sparse. As a consequence, should be of order and should be of order , where . In more general models with edge terms as well as other model terms, all model terms should scale as the edge term, so that no model term can dominate any other model term. These considerations suggest that parameters should scale with the log number of nodes. In simulations studies II and IV, the sizes of blocks are identical, so we use between-block edge parameters () and within-block edge and transitive edge parameters and (). In simulation study III, we use the same within-block parameters, but between-block edge parameters (), with varying from and .
In all scenarios, we assess block recovery by using Yule’s -coefficient:
[TABLE]
where
[TABLE]
Here, is an indicator function, which is if the statement in parentheses is true and is [math] otherwise. The quantity is the number of pairs of nodes that are assigned to distinct blocks under both and , while the quantity is the number of pairs of nodes that are assigned to the same block under both and . The sum of the other two quantities, , is the number of pairs of nodes on which and disagree. If and agree on all pairs of nodes, Yule’s -coefficient is . In fact, Yule’s -coefficient is bounded above by , and is invariant to the labeling of the blocks.
The results of the small-scale simulation study I are shown in Figure 1. The results suggest that the two-step estimation approach is almost as good as the Bayesian approach in terms of block recovery in the balanced case ( blocks of size ), but is worse in the unbalanced case ( blocks of size ). The worse performance in the unbalanced case may be due to the fact that there are smaller blocks in the unbalanced case than in the balanced case, and recovering small blocks is more challenging than recovering large blocks. That said, the advantage of the Bayesian approach over the two-step estimation approach in terms of block recovery in the unbalanced case is limited, and comes at excessive costs: Table 2 reveals that the computing time of the Bayesian approach is thousands of times higher than the computing time of the two-step estimation approach.
The large-scale simulation studies II, III, and IV shed light on the performance of Steps 1 and 2 of the two-step estimation approach in large networks. The results of simulation study II, shown in Figure 2, reveal that Step 1 of the two-step estimation approach outperforms spectral clustering in terms of block recovery in most scenarios. Simulation study III shows how the block recovery in Step 1 of the two-step estimation approach is affected by between-block subgraph sparsity. According to Figure 3, the more sparse between-block subgraphs are, the more dense within-block subgraphs “stand out,” improving block recovery. Last, but not least, simulation study IV helps assess parameter recovery in Step 2 of the two-step estimation approach. Figure 4 shows that maximum pseudolikelihood estimates are close to the data-generating parameters, and more so when the blocks are small and the number of blocks is large.
6 Amazon product recommendation network
We use the two-step estimation approach to shed light on the structure of a large Amazon product recommendation network that is not captured by stochastic block models. The data on the Amazon product recommendation network were collected by Yang and Leskovec (2015) and can be downloaded from the website
http://snap.stanford.edu/data/com-Amazon.html
The network consists of products listed at www.amazon.com. Two products and are connected by an edge if and are frequently purchased together according to the “Customers Who Bought This Item Also Bought” feature at www.amazon.com. Amazon assigns all products to categories, which we consider to be ground-truth blocks. We use a subset of the network consisting of the top 500 non-overlapping categories with 10 to 80 products, where the ranking of categories is based on Yang and Leskovec (2015). The resulting network consists of 10,448 products and 33,537 edges.
To model the Amazon product recommendation network, we take advantage of curved exponential-family random graph models with within-block edge and geometrically weighted degree and edgewise shared partner terms (Snijders et al., 2006; Hunter and Handcock, 2006; Hunter et al., 2008). The resulting models are more general than the motivating example used in Sections 2 and 3—the model with between-block edge terms and within-block edge and transitive edge terms—and ensure that, for each pair of products, the added value of additional edges and triangles within blocks decays. In fact, transitive edge terms are special cases of geometrically weighted edgewise shared partner terms, and both of them are well-behaved alternatives to the ill-behaved triangle terms mentioned in Section 2. A full-fledged discussion of those models is beyond the scope of our paper. We refer the interested reader to the seminal papers of Snijders et al. (2006), Hunter and Handcock (2006), and Hunter et al. (2008).
The natural parameters of the within-block edge terms are given by
[TABLE]
where the logarithmic term arises from sparsity considerations and is a simple form of a size-dependent parameterization, as explained in Section 5. The within-block geometrically weighted degree terms are based on the number of products with edges in block . The natural parameters of within-block geometrically weighted degree terms are given by
[TABLE]
The within-block geometrically weighted edgewise shared partner terms are based on the number of connected pairs of products and in block such that and have shared partners in block . The natural parameters of the within-block geometrically weighted edgewise shared partner terms are given by
[TABLE]
To reduce computing time, it is convenient to truncate the two geometrically weighted model terms by setting () and (). The two thresholds and are motivated by the fact that no product has or more edges and less than 1% of all pairs of products has or more edgewise shared partners. Last, but not least, the natural parameters of the between-block edge terms are given by
[TABLE]
The resulting exponential family is a curved exponential family (Hunter and Handcock, 2006), because the natural parameter vector of the exponential family is a nonlinear function of given . The natural parameter vector depends on the sizes of blocks, because we do not want to force small and large blocks to have the same natural parameters, as explained in Section 5.
The curved exponential-family model specified above can capture an excess in the expected number of triangles within blocks relative to stochastic block models, while ensuring that, for each pair of products, the added value of additional edges and triangles within blocks decays (Stewart et al., 2019). An excess in the expected number of triangles within blocks relative to stochastic block models can arise when, e.g., (a) three products are similar (e.g., three books on the same topic); (b) three products are dissimilar but complement each other (e.g., a bicycle helmet, head light, and tail light); and (c) three products, either similar or dissimilar, were produced by the same source (e.g., three books written by the same author).
Since we know the number of ground-truth blocks, we set and estimate the block structure by using the two-step estimation approach, using the variational approach in Step 1 and Monte Carlo maximum likelihood estimates in Step 2. To assess the performance of the two-step estimation approach in terms of block recovery, we use Yule’s -coefficient. Yule’s -coefficient turns out to be , which indicates near-perfect recovery of the ground-truth block structure. An inspection of the estimated block structure reveals that, out of the 500 categories of products in the Amazon product recommendation network, products in 15 categories are misclassified. Some of the small categories are merged with large categories, while some unusual products of large blocks are misclassified as well. Some of the products are unusual in the sense of having few edges to other products in the same category, while others are unusual in the sense of having many edges to other products in the same category.
The Monte Carlo maximum likelihood estimates and standard errors of are shown in Table 3. The parameters of the geometrically weighted terms can be interpreted in terms of log odds of conditional probabilities of edges, given all other edges (Hunter, 2007). A work-out example can be found in Stewart et al. (2019). Table 3 suggests that there is evidence for transitivity. The observed tendency towards transitivity has advantages in practice: It enables Amazon to recommend co-purchases of brandnew products and existing products, even when there are no data on past co-purchases of these products. For example, when a brandnew product is introduced (e.g., a novel) and product is known to be related to existing product (e.g., a novel by the same author), and product tends to be co-purchased with product (e.g., a classic novel), then Amazon could recommend co-purchases of products and as well as products and , despite the fact that there are no data on past co-purchases of products and and there is no direct connection between products and (although there in indirect connection via product ).
To demonstrate that the curved exponential-family random graph model considered here can capture structural features of networks that simpler models—such as stochastic block models—cannot capture, we compare the goodness-of-fit of the curved exponential-family random graph model to the goodness-of-fit of stochastic block models. Since the two models impose the same probability law on between-block subgraphs, it is natural to compare the two models in terms of goodness-of-fit with respect to within-block subgraphs. We assess the goodness-of-fit of the two models in terms of the within-block geodesic distances of pairs of products, i.e., the length of the shortest path between pairs of products in the same block; the numbers of within-block dyadwise shared partners, i.e., the number of unconnected or connected pairs of products with shared partners in the same block; the numbers of within-block edgewise shared partners, i.e., the number of connected pairs of products with shared partners in the same block; and the number of within-block transitive edges, i.e., the number of pairs of products with at least one shared partner in the same block. Figures 6 and 6 compare the goodness-of-fit of the two models based on 1,000 graphs simulated from the estimated models. The figures suggest that the curved exponential-family random graph model considered here is superior to the stochastic block model in terms of both connectivity and transitivity.
7 Discussion
The two-step estimation approach proposed here enables large-scale estimation of models with local dependence and unknown block structure provided that the number of blocks is known. An important direction of future research are methods for selecting when is unknown. We note that even in the special case of stochastic block models, the issue of selecting has not received much attention, with the notable exception of recent work by Saldana et al. (2017), Wang and Bickel (2017), and others. Whether—and how—the developed methods can be extended to models with local dependence is an open question, but having scalable methods for selecting would doubtless be useful in practice.
The proposed methods are implemented in R packages hergm (Schweinberger and Luna, 2018). A stable and user-friendly version will be released in the near future.
Acknowledgements
We acknowledge support from the National Science Foundation (NSF) in the form of NSF awards DMS-1513644 (SB, JS, MS) and DMS-1812119 (JS, MS).
Proofs of theoretical results
We prove Theorems 3.1 and 3.2 of Section 3. To prove them, we need three additional results, Lemma 7.3 and Propositions 7.4 and 7.5. To state them, let
[TABLE]
where is considered as a function of for fixed . Observe that the expectation of the sufficient statistic vector exists as long as the data-generating natural parameter vector is contained in the interior of the natural parameter space of the exponential family (Brown, 1986, Theorem 2.2, pp. 34–35). As a consequence, the expectations and exist, because
[TABLE]
and
[TABLE]
We first state Lemma 7.3 and Propositions 7.4 and 7.5 and then prove Theorems 3.1 and 3.2.
Lemma 7.3**.**
Consider a model with local dependence. Let be a function of within-block edge variables that is Lipschitz with respect to the Hamming metric with Lipschitz coefficient and . Then there exists a universal constant such that, for all and all ,
[TABLE]
Proof of Lemma 7.3. The proof of Lemma 7.3 resembles the proof of Proposition 1 of Schweinberger and Stewart (2020) step by step and is therefore omitted.
Proposition 7.4**.**
Consider a model with local dependence satisfying conditions [C.1] and [C.2]. Choose as small as desired. Then there exists a universal constant such that
[TABLE]
Proof of Proposition 7.4. To show that the probability mass of concentrates around its expectation , observe that the Lipschitz coefficient of the function with respect to the Hamming metric is defined by
[TABLE]
Since the term of and cancels, we obtain
[TABLE]
By condition [C.1] and the fact that , there exists such that
[TABLE]
Thus, the Lipschitz coefficient of the function is bounded above by
[TABLE]
By applying Lemma 7.3 to the function of within-block edges with Lipschitz coefficient , we can conclude that there exists such that, for all ,
[TABLE]
A union bound over the block structures shows that
[TABLE]
Choose as small as desired and let
[TABLE]
where . Then
[TABLE]
Proposition 7.5**.**
Consider a model with local dependence satisfying conditions [C.1] and [C.2]. Then there exists a universal constant such that
[TABLE]
Proof of Proposition 7.5. By definition,
[TABLE]
where, for all ,
[TABLE]
By construction of and , the contributions of between-block subgraphs to the loglikelihood function are the same under both models, hence the expectation of the loglikelihood ratio reduces to the expectation of the loglikelihood ratio of within-block subgraphs:
[TABLE]
By the triangle inequality,
[TABLE]
The terms can be bounded above by invoking the triangle inequality:
[TABLE]
where , , and are the natural parameter vector, the sufficient statistics vector, and the log-normalizing constant of , and . We bound the two terms on the right-hand side of the inequality above one by one.
First term. By condition [C.2], there exist and such that
[TABLE]
where the last inequality follows from the assumption that is compact.
Second term. By the mean-value theorem along with classic exponential-family properties, there exists () such that
[TABLE]
It is worth noting that the existence of is ensured by the fact that the natural parameter space of exponential families is convex (Theorem 1.13, Brown, 1986, p. 19), which implies that an element in the interior of the natural parameter exists as long as and are contained in the interior of the natural parameter space. Therefore, the second term can be bounded along the same lines as the first term, which implies that there exist such that
[TABLE]
Conclusion. Collecting terms shows that that there exist such that
[TABLE]
Armed with Propositions 7.4 and 7.5, we can prove Theorem 3.1 and 3.2.
Proof of Theorem 3.1. Theorem 3.1 is an application of Theorem 3.2. To prove Theorem 3.1, all we need to do verify conditions [C.1] and [C.2] of Theorem 3.2.
Condition [C.1]. By the triangle inequality,
[TABLE]
The first term is bounded above by , whereas the other two terms are bounded above by , which implies that
[TABLE]
by the compactness of (). Thus, condition [C.1] is satisfied.
Condition [C.2]. Condition [C.2] is satisfied, because
[TABLE]
where denotes the -norm whereas denotes the -norm of the vector of interest.
Proof of Theorem 3.2.
Observe that, for all ,
[TABLE]
Choose as small as desired and let .
Then, by using
[TABLE]
along with Propositions 7.4 and 7.5, we can conclude that
[TABLE]
References
- Brown (1986)
Brown, L. (1986), Fundamentals of Statistical Exponential Families: With Applications in Statistical Decision Theory, Hayworth, CA, USA: Institute of Mathematical Statistics.
- Schweinberger and Stewart (2020)
Schweinberger, M., and Stewart, J. (2020), “Concentration and consistency results for canonical and curved exponential-family models of random graphs,” The Annals of Statistics, 48, 374–396.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Amini et al. (2013) Amini, A. A., Chen, A., Bickel, P. J., and Levina, E. (2013), “Pseudo-likelihood methods for community detection in large sparse networks,” The Annals of Statistics , 41, 2097–2122.
- 2Atchade et al. (2013) Atchade, Y. F., Lartillot, N., and Robert, C. (2013), “Bayesian computation for statistical models with intractable normalizing constants,” Brazilian Journal of Probability and Statistics , 27, 416–436.
- 3Besag (1974) Besag, J. (1974), “Spatial interaction and the statistical analysis of lattice systems,” Journal of the Royal Statistical Society, Series B , 36, 192–225.
- 4Bickel and Chen (2009) Bickel, P. J., and Chen, A. (2009), “A nonparametric view of network models and Newman-Girvan and other modularities,” in Proceedings of the National Academy of Sciences , Vol. 106, pp. 21068–21073.
- 5Bickel et al. (2011) Bickel, P. J., Chen, A., and Levina, E. (2011), “The method of moments and degree distributions for network models,” The Annals of Statistics , 39, 2280–2301.
- 6Bickel et al. (2013) Bickel, P. J., Choi, D., Chang, X., and Zhang, H. (2013), “Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels,” The Annals of Statistics , 41, 1922–1943.
- 7Binkiewicz et al. (2017) Binkiewicz, N., Vogelstein, J. T., and Rohe, K. (2017), “Covariate-assisted spectral clustering,” Biometrika , 104, 361–377.
- 8Brown (1986) Brown, L. (1986), Fundamentals of Statistical Exponential Families: With Applications in Statistical Decision Theory , Hayworth, CA, USA: Institute of Mathematical Statistics.