Fitting the Linear Preferential Attachment Model

TL;DR

This paper develops methods to fit a 5-parameter linear preferential attachment model to network data, providing estimators that are consistent and effective in both full history and snapshot scenarios, with applications to real data.

Contribution

It introduces novel estimation techniques for the linear preferential attachment model applicable to different data availability scenarios, including maximum likelihood and combined method of moments approaches.

Findings

MLE estimator is strongly consistent and asymptotically normal.

Proposed estimator performs well compared to MLE in snapshot data.

Methods are validated through simulations and real data examples.

Abstract

Preferential attachment is an appealing mechanism for modeling power-law behavior of the degree distributions in directed social networks. In this paper, we consider methods for fitting a 5-parameter linear preferential model to network data under two data scenarios. In the case where full history of the network formation is given, we derive the maximum likelihood estimator of the parameters and show that it is strongly consistent and asymptotically normal. In the case where only a single-time snapshot of the network is available, we propose an estimation method which combines method of moments with an approximation to the likelihood. The resulting estimator is also strongly consistent and performs quite well compared to the MLE estimator. We illustrate both estimation procedures through simulated data, and explore the usage of this model in a real data example.

Tables

Table 1. Table 5.1: Mean of 𝜽 ^ n M ​ L ​ E superscript subscript ^ 𝜽 𝑛 𝑀 𝐿 𝐸 \hat{\boldsymbol{\theta}}_{n}^{MLE} and 𝜽 ~ n subscript ~ 𝜽 𝑛 \tilde{\boldsymbol{\theta}}_{n} with ARE’s of 𝜽 ~ n subscript ~ 𝜽 𝑛 \tilde{\boldsymbol{\theta}}_{n} relative to 𝜽 ^ n M ​ L ​ E superscript subscript ^ 𝜽 𝑛 𝑀 𝐿 𝐸 \hat{\boldsymbol{\theta}}_{n}^{MLE} for 𝜽 = ( 0.3 , 0.5 , 2 , 1 ) 𝜽 0.3 0.5 2 1 \boldsymbol{\theta}=(0.3,0.5,2,1) under different choices of n 𝑛 n .

$n$	$Mean\left({\widehat{𝜽}}_n^{MLE}\right)$	$Mean\left({\widetilde{𝜽}}_n\right)$	$ARE\left({\widetilde{𝜽}}_n\right)$
$1000$	(0.300, 0.500, 2.076, 1.054)	(0.301, 0.500, 2.128, 1.066)	(0.408, 1.000, 0.397, 0.228)
$5000$	(0.300, 0.500, 2.022, 1.013)	(0.301, 0.500, 2.036, 1.010)	(0.414, 1.000, 0.386, 0.236)
$10000$	(0.300, 0.500, 2.011, 1.006)	(0.301, 0.500, 2.019, 1.006)	(0.408, 1.000, 0.388, 0.232)
$50000$	(0.300, 0.500, 2.003, 1.002)	(0.300, 0.500, 2.005, 1.002)	(0.399, 1.000, 0.393, 0.230)
$100000$	(0.300, 0.500, 2.001, 1.001)	(0.300, 0.500, 2.003, 1.000)	(0.392, 1.000, 0.382, 0.223)

Table 2. Table 5.2: Mean of 𝜽 ^ n M ​ L ​ E superscript subscript ^ 𝜽 𝑛 𝑀 𝐿 𝐸 \hat{\boldsymbol{\theta}}_{n}^{MLE} and 𝜽 ~ n subscript ~ 𝜽 𝑛 \tilde{\boldsymbol{\theta}}_{n} with ARE’s of 𝜽 ~ n subscript ~ 𝜽 𝑛 \tilde{\boldsymbol{\theta}}_{n} relative to 𝜽 ^ n M ​ L ​ E superscript subscript ^ 𝜽 𝑛 𝑀 𝐿 𝐸 \hat{\boldsymbol{\theta}}_{n}^{MLE} for ( n , δ in , δ out ) = ( 10 5 , 2 , 1 ) 𝑛 subscript 𝛿 in subscript 𝛿 out superscript 10 5 2 1 (n,\delta_{\text{in}},\delta_{\text{out}})=(10^{5},2,1) under different choices of ( α , β ) 𝛼 𝛽 (\alpha,\beta) .

$(\alpha ,\beta )$	$Mean\left({\widehat{𝜽}}_n^{MLE}\right)$	$Mean\left({\widetilde{𝜽}}_n\right)$	$ARE\left({\widetilde{𝜽}}_n\right)$
(0.001, 0.99)	(0.001, 0.990, 2.034, 1.016)	(0.001, 0.990, 2.071, 1.049)	(0.291, 1.000, 0.147, 0.316)
(0.01, 0.9)	(0.010, 0.900, 2.004, 1.001)	(0.010, 0.900, 2.008, 1.004)	(0.331, 1.000, 0.207, 0.381)
(0.1, 0.8)	(0.100, 0.800, 2.003, 1.001)	(0.100, 0.800, 2.004, 1.002)	(0.353, 1.000. 0.264, 0.216)
(0.2, 0.6)	(0.200, 0.600, 2.002, 1.001)	(0.200, 0.600, 2.003, 1.001)	(0.364, 1.000, 0.309, 0.236)
(0.5, 0.3)	(0.500, 0.300, 2.001, 1.001)	(0.500, 0.300, 2.002, 1.000)	(0.472, 1.000, 0.529, 0.202)
(0.7, 0.2)	(0.700, 0.200, 2.002, 1.000)	(0.700, 0.200, 2.002, 1.000)	(0.726, 1.000, 0.793, 0.217)
(0.1, 0.3)	(0.100, 0.300, 2.001, 1.001)	(0.100, 0.300, 2.002, 1.000)	(0.420, 1.000, 0.313, 0.629)

Equations

$$\text{P}[\text{choose w∈V(n−1)}]=\frac{D_{\mathrm{in}}^{(n-1)}(w)+{\delta }_{\text{i}n}}{n-1+{\delta }_{\text{i}n}N(n-1)}.$$

$$\textbf{P}[\text{choose $(v,w)$}]=\Bigl{(}\frac{D_{\rm out}^{(n-1)}(v)+\delta_{\text{o}ut}}{n-1+\delta_{\text{o}ut}N(n-1)}\Bigr{)}\Bigl{(}\frac{D_{\rm in}^{(n-1)}(w)+\delta_{\text{i}n}}{n-1+\delta_{\text{i}n}N(n-1)}\Bigr{)}.$$

$$\text{P}[\text{choose v∈V(n−1)}]=\frac{D_{\mathrm{out}}^{(n-1)}(v)+{\delta }_{\text{o}ut}}{n-1+{\delta }_{\text{o}ut}N(n-1)}.$$

$$\frac{N_{ij}(n)}{N(n)}\to f_{ij}=:\frac{p_{ij}}{1-\beta },n\to \infty .$$

$$p_0^{\text{in}}$$

$$p_i^{\text{in}}$$

$$a_1(\lambda ):=\frac{\alpha +\beta }{1+\lambda (1-\beta )},\lambda >0.$$

$$p_i^{\text{in}}$$

$${\iota }_{\text{in}}=1+\frac{1+{\delta }_{\text{in}}(\alpha +\gamma )}{\alpha +\beta }$$

$$p_j^{\text{out}}$$

$${\iota }_{\text{out}}=1+\frac{1+{\delta }_{\text{out}}(\alpha +\gamma )}{\beta +\gamma }.$$

$$v\leftarrow v_{\lceil W\rceil }^{(j)}{1}_{\{W\le n-1\}}+\lceil \frac{W-(n-1)}{{\delta }_{\text{i}n}}\rceil 1_{\{W>n-1\}}.$$

$$\text{P}\left(v=w\right)$$

$$L$$

$$=$$

$$\displaystyle\times\prod_{t=n_{0}+1}^{n}\left(\beta\Bigl{(}\frac{D_{\text{in}}^{(t-1)}(v^{(2)}_{t})+\delta_{\text{in}}}{t-1+\delta_{\text{in}}N(t-1)}\Bigr{)}\Bigl{(}\frac{D_{\text{out}}^{(t-1)}(v^{(1)}_{t})+\delta_{\text{out}}}{t-1+\delta_{\text{out}}N(t-1)}\Bigr{)}\right)^{\textbf{1}_{\{J_{t}=2\}}}$$

$$\times \prod _{t=n_0+1}^n{\left((1-\alpha -\beta )\frac{D_{\text{out}}^{(t-1)}(v_t^{(1)})+{\delta }_{\text{out}}}{t-1+{\delta }_{\text{out}}N(t-1)}\right)}^{{\text{1}}_{\{J_t=3\}}}$$

$$\log$$

$$=$$

$$+\sum _{t=n_0+1}^n\log \left(D_{\text{in}}^{(t-1)}(v_t^{(2)})+{\delta }_{\text{in}}\right){\text{1}}_{\{J_t\in \{1,2\}\}}+\sum _{t=n_0+1}^n\log \left(D_{\text{out}}^{(t-1)}(v_t^{(1)})+{\delta }_{\text{out}}\right){\text{1}}_{\{J_t\in \{2,3\}\}}$$

$$-\sum _{t=n_0+1}^n\log (t-1+{\delta }_{\text{in}}N(t-1)){\text{1}}_{\{J_t\in \{1,2\}\}}-\sum _{t=n_0+1}^n\log (t-1+{\delta }_{\text{out}}N(t-1)){\text{1}}_{\{J_t\in \{2,3\}\}}.$$

$$\frac{\partial }{\partial \alpha }\log L(\alpha ,\beta ,{\delta }_{\text{in}},{\delta }_{\text{out}}|G(n_0),(e_t{)}_{t=n_0+1}^n)=\frac{1}{\alpha }\sum _{t=n_0+1}^n{\text{1}}_{\{J_t=1\}}-\frac{1}{1-\alpha -\beta }\sum _{t=n_0+1}^n{\text{1}}_{\{J_t=3\}},$$

$$\frac{\partial }{\partial \beta }\log L(\alpha ,\beta ,{\delta }_{\text{in}},{\delta }_{\text{out}}|G(n_0),(e_t{)}_{t=n_0+1}^n)=\frac{1}{\beta }\sum _{t=n_0+1}^n{\text{1}}_{\{J_t=2\}}-\frac{1}{1-\alpha -\beta }\sum _{t=n_0+1}^n{\text{1}}_{\{J_t=3\}},$$

$$\frac{\partial }{\partial {\delta }_{\text{in}}}\log L(\alpha ,\beta ,{\delta }_{\text{in}},{\delta }_{\text{out}}|G(n_0),(e_t{)}_{t=n_0+1}^n)$$

$$=\sum _{t=n_0+1}^n\frac{1}{D_{\text{in}}^{(t-1)}(v_t^{(2)})+{\delta }_{\text{in}}}{\text{1}}_{\{J_t\in \{1,2\}\}}-\sum _{t=n_0+1}^n\frac{N(t-1)}{t-1+{\delta }_{\text{in}}N(t-1)}{\text{1}}_{\{J_t\in \{1,2\}\}},$$

$$\frac{\partial }{\partial {\delta }_{\text{out}}}\log L(\alpha ,\beta ,{\delta }_{\text{in}},{\delta }_{\text{out}}|G(n_0),(e_t{)}_{t=n_0+1}^n)$$

$$=\sum _{t=n_0+1}^n\frac{1}{D_{\text{out}}^{(t-1)}(v_t^{(1)})+{\delta }_{\text{out}}}{\text{1}}_{\{J_t\in \{2,3\}\}}-\sum _{t=n_0+1}^n\frac{N(t-1)}{t-1+{\delta }_{\text{out}}N(t-1)}{\text{1}}_{\{J_t\in \{2,3\}\}}.$$

$${\hat{\alpha }}^{MLE}$$

$${\hat{\beta }}^{MLE}$$

$$\sum _{t=n_0+1}^n\frac{1}{D_{\text{in}}^{(t-1)}(v_t^{(2)})+{\delta }_{\text{in}}}{\text{1}}_{\left\{J_t\in \{1,2\}\right\}}$$

$$N_i^{\text{in}}(n)=\sum _{j=0}^{\infty }{N}_{ij}(n),$$

$$\sum _{t=n_0+1}^n{\text{1}}_{\left\{D_{\text{in}}^{(t-1)}(v_t^{(2)})=i,J_t\in \{1,2\}\right\}}=N_{>i}^{\text{in}}(n)-N_{>i}^{\text{in}}(n_0),i\ge 1.$$

$$\sum _{t=n_0+1}^n{\text{1}}_{\left\{D_{\text{in}}^{(t-1)}(v_t^{(2)})=0,J_t\in \{1,2\}\right\}}=N_{>0}^{\text{in}}(n)-N_{>0}^{\text{in}}(n_0)-\sum _{t=n_0+1}^n{\text{1}}_{\{J_t=3\}},$$

$$\sum _{t=n_0+1}^n\frac{1}{D_{\text{in}}^{(t-1)}(v_t^{(2)})+{\delta }_{\text{in}}}{\text{1}}_{\left\{J_t\in \{1,2\}\right\}}$$

$$=\sum _{i=0}^{\infty }\frac{N_{>i}^{\text{in}}(n)-N_{>i}^{\text{in}}(n_0)}{i+{\delta }_{\text{in}}}-\frac{{\sum }_{t=n_0+1}^n{\text{1}}_{\{J_t=3\}}}{{\delta }_{\text{in}}}.$$

$$\frac{1}{n-n_0}\sum _{i=0}^{\infty }$$

$$-\frac{1}{{\delta }_{\text{in}}(n-n_0)}\sum _{t=n_0+1}^n{\text{1}}_{\{J_t=3\}}-\frac{1}{n-n_0}\sum _{t=n_0+1}^n\frac{N(t-1)}{t-1+{\delta }_{\text{in}}N(t-1)}{\text{1}}_{\{J_t\in \{1,2\}\}}=0,$$

$$\frac{1}{n-n_0}\sum _{j=0}^{\infty }$$

Full text

power laws,

multivariate heavy tail statistics,

preferential attachment,

estimation,

1 Introduction

The preferential attachment mechanism, in which edges and nodes are added to the network based on probabilistic rules, provides an appealing description for the evolution of a network. The rule for how edges connect nodes depends on node degree; large degree nodes attract more edges. The idea is applicable to both directed and undirected graphs and is often the basis for studying social networks, collaborator and citation networks, and recommender networks. Elementary descriptions of the preferential attachment model can be found in [5] while more mathematical treatments are available in [4, 21, 2]. Also see [10] for a statistical survey of methods for network data, [18] for consideration of statistics of an undirected network and [24] for asymptotics of a directed exponential random graph models. Limit theory for estimates of an undirected preferential attachment model was considered in [6].

For many networks, empirical evidence supports the hypothesis that in- and out-degree distributions follow a power law. This property has been shown to hold in linear preferential attachment models, which makes preferential attachment an attractive choice for network modeling [4, 21, 11, 12, 3]. While the marginal degree power laws in a simple linear preferential attachment model were established in [11, 12, 3], the joint regular variation (see [16, 15]) which is akin to a joint power law, was only recently established [19, 17]. In addition, it was shown in [22] that the joint probability mass function of the in- and out-degrees is multivariate regularly varying. This is a key result as the degrees of a network are integer-valued.

In this paper, we discuss methods of fitting a simple linear preferential attachment model, which is parametrized by $\boldsymbol{\theta}=(\alpha,\beta,\gamma,\delta_{\text{in}},\delta_{\text{out}})$. The first three parameters, $\alpha,\beta,\gamma$, correspond to probabilities of the 3 scenarios for adding an edge and hence sum to 1, i.e., $\alpha+\beta+\gamma=1$. The other two, $\delta_{\text{in}}$ and $\delta_{\text{out}}$, are tuning parameters related to growth rates. The tail indices of the marginal power laws for the in- and out-degrees can be expressed as explicit functions of $\boldsymbol{\theta}$ (see (2.6) and (2.7) below). The graph $G(n)=(V(n),E(n))$, where $V(n)$ is the set of nodes and $E(n)$ is the set of edges at the $n$th iteration, evolves based on postulates that describe how new edges and nodes are formed. This construction of the network is Markov in the sense that the probabilistic rules for obtaining $G(n+1)$ once $G(n)$ is known do not require prior knowledge of earlier stages of the construction.

The Markov structure of the model allows us to construct a likelihood function based on observing $G(n_{0}),G(n_{0}+1),\ldots,G(n_{0}+n)$. After deriving the likelihood function, we show that there exists a unique maximum at $\hat{\boldsymbol{\theta}}=(\hat{\alpha},\hat{\beta},\hat{\gamma},\hat{\delta}_{\text{in}},\hat{\delta}_{\text{out}})$ and that the resulting maximum likelihood estimator is strongly consistent and asymptotically normal. The normality is proved using a martingale central limit theorem applied to the score function. The limiting distribution also reveals that $(\hat{\alpha},\hat{\beta},\hat{\gamma})$, $\hat{\delta}_{\text{in}}$, and $\hat{\delta}_{\text{out}}$ are asymptotically independent. From these results, asymptotic properties of the MLE for the power law indices can be derived.

For some network data, only a snapshot of the nodes and edges is available at a single point in time, that is, only $G(n)$ is available for some $n$. In such cases, we propose an estimation procedure for the parameters of the network using an approximation to the likelihood and method of moments. This also produces strongly consistent estimators. These estimators perform reasonably well compared to the MLE where the entire evolution of the network is known but predictably there is some loss of efficiency.

We illustrate the estimation procedure for both scenarios using simulated data. Simulation plays an important role in the process of modeling networks since it provides a way to assess the performance of model fitting procedures in the idealized setting of knowing the true model. Also, after fitting a model to real data, simulation provides a check on the quality of fit. Departures from model assumptions can often be detected via simulation of multiple realizations from the fitted network. Hence it is important to have efficient simulation algorithms for producing realizations of the preferential attachment network for a given set of parameter values. We adopt a simulation method, learned from Joyjit Roy, that was inspired by [1] and is similar to that of [20].

Our fitting methods are implemented in a real data setting using the Dutch Wiki talk network [14]. While one should not expect the simple 5-parameter (later extended to 7-parameter) linear preferential attachment model to fully explain a network with millions of edges, it does provide a reasonable fit to the tail behavior of the degree distributions. We are also able to detect important structural features in the network through fitting the model over separate time intervals.

Often it is difficult to believe in the existence of a true model, especially one whose parameters remain constant over time. Allowing, as we do, a preferential attachment model with only a few parameters and no possibility for node removal may seem simplistic and unrealistic for social network data. Of course, preferential attachment is only one mechanism for network formation and evidence for its use in fields outside data networks is mixed [8, 9] and we restrict attention to linear preferential attachment. Even imperfect models have the potential to capture salient properties in the data, such as heavy-tailedness of the in-degree and out-degree distributions, and to identify departures from model assumptions. While maximum likelihood estimation is essentially the gold standard for cases when the underlying model is a good representation of the data, it may perform poorly in case the model is far from being appropriate. In forthcoming work, we consider a semi-parametric estimation approach for network models that exhibit heavy-tailed degree distributions. This alternative estimation methodology borrows ideas from extreme value theory.

The rest of the paper is structured as follows. In Section 2, we formulate the linear preferential attachment network model and present an efficient simulation method for the network. Section 3 gives parameter estimators when either the full history is known or when only a single snapshot in time is available. We test these estimators against simulated data in Section 5 and then explore the Wiki talk network in Section 6.

2 Model specification and simulation

In this section, we present the linear preferential attachment model in detail and provide a fast simulation algorithm for the network.

2.1 The linear preferential attachment model

The directed edge preferential attachment model [3, 12] constructs a growing directed random graph $G(n)=(V(n),E(n))$ whose dynamics depend on five non-negative real numbers $\alpha,\beta,\gamma$, $\delta_{\text{i}n}$ and $\delta_{\text{o}ut}$, where $\alpha+\beta+\gamma=1$ and $\delta_{\text{in}},\delta_{\text{out}}>0$. To avoid degenerate situations, assume that each of the numbers $\alpha,\beta,\gamma$ is strictly smaller than 1. We obtain a new graph $G({n})$ by adding one edge to the existing graph $G({n-1})$ and index the constructed graphs by the number $n$ of edges in $E(n)$. We start with an arbitrary initial finite directed graph $G({n_{0}})$ with at least one node and $n_{0}$ edges. For $n>n_{0}$, $G(n)=(V(n),E(n))$ is a graph with $|E(n)|=n$ edges and a random number $|V(n)|=N(n)$ of nodes. If $u\in V(n)$, $D_{\rm in}^{(n)}(u)$ and $D_{\rm out}^{(n)}(u)$ denote the in- and out-degree of $u$ respectively in $G(n)$. There are three scenarios that we call the $\alpha$, $\beta$ and $\gamma$-schemes, which are activated by flipping a 3-sided coin whose outcomes are $1,2,3$ with probabilities $\alpha,\beta,\gamma$. More formally, we have an iid sequence of multinomial random variables $\{J_{n},n>n_{0}\}$ with cells labelled $1,2,3$ and cell probabilities $\alpha,\beta,\gamma$. Then the graph $G(n)$ is obtained from $G(n-1)$ as follows.

• •

If $J_{n}=1$ (with probability $\alpha$), append to $G(n-1)$ a new node $v\in V(n)\setminus V(n-1)$ and an edge $(v,w)$ leading from $v$ to an existing node $w\in V(n-1)$. Choose the existing node $w\in V(n-1)$ with probability depending on its in-degree in $G(n-1)$:

[TABLE]

• •

If $J_{n}=2$ (with probability $\beta$), add a directed edge $(v,w)$ to $E({n-1})$ with $v\in V(n-1)=V(n)$ and $w\in V(n-1)=V(n)$ and the existing nodes $v,w$ are chosen independently from the nodes of $G(n-1)$ with probabilities

[TABLE]

• •

If $J_{n}=3$ (with probability $\gamma$), append to $G(n-1)$ a new node $w\in V(n)\setminus V(n-1)$ and an edge $(v,w)$ leading from the existing node $v\in V(n-1)$ to the new node $w$. Choose the existing node $v\in V(n-1)$ with probability

[TABLE]

Note that this construction allows the possibility of having self loops in the case where $J_{n}=2$, but the proportion of edges that are self loops goes to 0 as $n\to\infty$. Also, multiple edges are allowed between two nodes.

2.2 Power law of degree distributions

Given an observed network with $n$ edges, let $N_{ij}(n)$ denote the number of nodes with in-degree $i$ and out-degree $j$. If the network is generated from the linear preferential attachment model described above, then from [3], there exists a proper probability distribution $\{f_{ij}\}$ such that almost surely

[TABLE]

Consider the limiting marginal in-degree distribution $p^{\text{in}}_{i}:=\sum_{j}p_{ij}$. It is calculated from [3, Equation (3.10)] that

[TABLE]

where

[TABLE]

Moreover, $p^{\text{in}}_{i}$ satisfies

[TABLE]

for some finite positive constant $C_{\text{in}}$, where the power index

[TABLE]

Similarly, the limiting marginal out-degree distribution has the same property:

[TABLE]

for $C_{\text{out}}$ positive and

[TABLE]

2.3 Simulation algorithm

We describe an efficient simulation procedure for the preferential attachment network given the parameter values $(\alpha,\beta,\gamma,\delta_{\text{i}n},\delta_{\text{o}ut})$, where $\alpha+\beta+\gamma=1$. The simulation cost of the algorithm is linear in time. This algorithm, which was provided by Joyjit Roy during his graduate work at Cornell University, is presented below for completeness. Note that this simulation algorithm is specifically designed for the case where the preferential attachment probabilities (2.1)–(2.3) are linear in the degrees. A similar idea for the simulation of the Yule-Simon process appeared in [20]. Efficient simulation methods for the case where the preferential attachment probabilities are non-linear are studied in [1], where their algorithm trades some efficiency for the flexibility to model non-linear preferential attachment.

Using the notation from the introduction, at time $t=0$, we initiate with an arbitrary graph $G(n_{0})=(V(n_{0}),E(n_{0}))$ of $n_{0}$ edges, where the elements of $E(n_{0})$ are represented in form of $(v_{i}^{(1)},v_{i}^{(2)})\in V(n_{0})\times V(n_{0})$, $i=1,\ldots,n_{0}$, with $v_{i}^{(1)},v_{i}^{(2)}$ denoting the outgoing and incoming vertices of the edge, respectively. To grow the network, we update the network at each stage from $G(n-1)$ to $G(n)$ by adding a new edge $(v_{n}^{(1)},v_{n}^{(2)})$. Assume that the nodes are labeled using positive integers starting from 1 according to the time order in which they are created, and let the random number $N(n)=|V(n)|$ denote the total number of nodes in $G(n)$.

Let us consider the situation where an existing node is to be chosen from $V(n)$ as the vertex of the new edge. Naively sampling from the multinomial distribution requires $O(N(n))$ evaluations, where $N(n)$ increases linearly with $n$. Therefore the total cost to simulate a network of $n$ edges is $O(n^{2})$. This is significantly burdensome when $n$ is large, which is usually the case for observed networks. Algorithm 1 describes a simulation algorithm which uses the alias method [13] for node sampling. Here sampling an existing node from $V(n)$ requires only constant execution time, regardless of $n$. Hence the cost to simulate $G(n)$ is only $O(n)$. This method allows generation of a graph with $10^{7}$ nodes on a personal laptop in less than 5 seconds.

To see that the algorithm indeed produces the intended network, it suffices to consider the case of sampling an existing node from $V(n-1)$ as the incoming vertex of the new edge. In the function Node_Sample in Algorithm 1, we generate $W\sim\text{Uniform}(0,n-1+N(n-1)\delta_{\text{i}n})$ and set

[TABLE]

Then

[TABLE]

which corresponds to the desired selection probability (2.1).

3 Parameter estimation: MLE based on the full network history

In this section, we estimate the preferential attachment parameter vector $(\alpha,\beta,\delta_{\text{in}},\delta_{\text{out}})$ under two assumptions about what data is available. In the first scenario, the full evolution of the network is observed, from which the likelihood function can be computed. The resulting MLE is strongly consistent and asymptotically normal. For the second scenario, the data only consist of one snapshot of the network with $n$ edges, without the knowledge of the network history that produced these edges. For this scenario we give an estimation approach through approximating the score function and moment matching, which produces parameter estimators that are also strongly consistent but less efficient than those based on the full evolution of the network. In both cases, the estimators are uniquely determined.

3.1 Likelihood calculation

Assume the network begins with the graph $G(n_{0})$ (consisting of $n_{0}$ edges) and then evolves according to the description in Section 2.1 with parameters $(\alpha,\beta,\delta_{\text{in}},\delta_{\text{out}})$, where $\delta_{\text{in}},\delta_{\text{out}}>0$ and $\alpha,\beta$ are non-negative probabilities. The $\gamma$ is implicitly defined by $\gamma=1-\alpha-\beta$. To avoid trivial cases, we will also assume $\alpha,\beta,\gamma<1$ for the rest of the paper. For MLE estimation we restrict the parameter space for $\delta_{\text{in}},\delta_{\text{out}}$ to be $[\epsilon,K]$, for some sufficiently small $\epsilon>0$ and large $K$. In particular, the true value of $\delta_{\text{in}},\delta_{\text{out}}$ is assumed to be contained in $(\epsilon,K)$. Let $e_{t}=(v^{(1)}_{t},v^{(2)}_{t})$ be the newly created edge when the random graph evolves from $G(t-1)$ to $G(t)$. We sometimes refer to $t$ as the time rather than the number of edges.

Assume we observe the initial graph $G(n_{0})$ and the edges $\{e_{t}\}_{t=n_{0}+1}^{n}$ in the order of their formation. For $t=n_{0}+1,\ldots,n$, the values of the following variables are known:

• •

$N(t)$, the number of nodes in graph $G(t)$;

• •

$D_{\text{in}}^{(t-1)}(v)$, $D_{\text{out}}^{(t-1)}(v)$, the in- and out-degree of node $v$ in $G(t-1)$, for all $v\in V(t-1)$;

• •

$J_{t}$, the scenario under which $e_{t}$ is created.

Then the likelihood function is

[TABLE]

and the log likelihood function is

[TABLE]

The score functions for $\alpha,\beta,\delta_{\text{in}},\delta_{\text{out}}$ are calculated as follows:

[TABLE]

Note that the score functions (3.3), (3.4) for $\alpha$ and $\beta$ do not depend on $\delta_{\text{in}}$ and $\delta_{\text{out}}$. One can show that the Hessian matrix of the log-likelihood for $(\alpha,\beta)$ is positive definite. Setting (3.3) and (3.4) to zero gives the unique MLE estimates for $\alpha$ and $\beta$,

[TABLE]

These estimates are strongly consistent by applying the strong law of large numbers for the $\{J_{t}\}$ sequence.

Next, consider the first term of the score function for $\delta_{\text{in}}$ in (3.5), and we have

[TABLE]

Observe that $\left\{D_{\text{in}}^{(t-1)}(v^{(2)}_{t})=i,J_{t}\in\{1,2\}\right\}$ describes the event that the in-degree of node $v^{(2)}_{t}{\in V(t-1)}$ is $i$ at time $t-1$ and is augmented to $i+1$ at time $t$. For each $i\geq 1$, such an event happens at some stage $t\in\{n_{0}+1,n_{0}+2,\ldots,n\}$ only for those nodes with in-degree $\leq i$ at time $n_{0}$ and in-degree $>i$ at time $n$. Let $N_{ij}(n)$ denote the number of nodes with in-degree $i$ and out-degree $j$ at time $n$, and $N^{\text{in}}_{i}(n)$ and $N^{\text{in}}_{>i}(n)$ to be the number of nodes with in-degree equal to $i$ and greater than $i$, respectively, i.e.,

[TABLE]

Then

[TABLE]

On the other hand, when $i=0$, $\left\{D_{\text{in}}^{(t-1)}(v^{(2)}_{t})=0,J_{t}\in\{1,2\}\right\}$ occurs for some $t$ if and only if all of the following three events happen:

• (i)

$v^{(2)}_{t}$ has in-degree $>0$ at time $n$;

• (ii)

$v^{(2)}_{t}$ does not have in-degree $>0$ at time $n_{0}$;

• (iii)

$v^{(2)}_{t}$ was not created under the $\gamma$-scheme (otherwise it would have been born with in-degree 1).

This implies:

[TABLE]

since there are, in total, $\sum_{t=n_{0}+1}^{n}\textbf{1}_{\{J_{t}=3\}}$ nodes created under the $\gamma$-scheme. Therefore,

[TABLE]

Setting the score function (3.5) for $\delta_{\text{in}}$ to 0 and dividing both sides by $n-n_{0}$ leads to

[TABLE]

where the only unknown parameter is $\delta_{\text{in}}$. In Section 3.2, we show that the solution to (3.1) actually maximizes the likelihood function in $\delta_{\text{in}}$. Similarly, the MLE for $\delta_{\text{out}}$ can be solved from

[TABLE]

where $N^{\text{out}}_{>j}(n)$ is defined in the same fashion as $N^{\text{in}}_{>i}(n)$.

Remark 3.1**.**

The arguments leading to (3.1) allow us to rewrite the likelihood function (3.1):

[TABLE]

Hence by the factorization theorem, $N(n_{0})$, $(J_{t})_{t=n_{0}+1}^{n}$, $(N^{\text{in}}_{>i}(n)-N^{\text{in}}_{>i}(n_{0}))_{i\geq 0}$, $(N^{\text{out}}_{>j}(n)-N^{\text{out}}_{>j}(n_{0}))_{j\geq 0}$ are sufficient statistics for $(\alpha,\beta,\delta_{\text{in}},\delta_{\text{out}})$.

3.2 Consistency of MLE

We remarked after (3.6) and (3.7) that $\hat{\alpha}^{MLE}$ and $\hat{\beta}^{MLE}$ converge almost surely to $\alpha$ and $\beta$. We now prove that the MLE of $(\delta_{\text{in}},\delta_{\text{out}})$ is also strongly consistent. Note that if we initiate the network with $G(n_{0})$ (for both $n_{0}$ and $N(n_{0})$ finite), then almost surely for all $i,j\geq 0$,

[TABLE]

and $(n-n_{0})/n\to 1$. In other words, $n_{0}$, $N^{\text{in}}_{>i}(n_{0})$, $N^{\text{out}}_{>j}(n_{0})$ are all $o(n)$. So for simplicity, we assume that the graph is initiated with finitely many nodes and no edges, that is, $n_{0}=0$ and $N(0)\geq 1$. In particular, these assumptions imply the sum of the in-degrees at time $n$ is equal to $n$.

Let $\Psi_{n}(\cdot),\Phi_{n}(\cdot)$ be the functional forms of the terms in the log-likelihood function (3.2) involving $\delta_{\text{in}}$ and $\delta_{\text{out}}$ respectively, normalized by $1/n$, i.e.,

[TABLE]

The following theorem gives the consistency of the MLE of $\delta_{\text{in}}$ and $\delta_{\text{out}}$.

Theorem 3.2**.**

Suppose $\delta_{\text{in}},\delta_{\text{out}}\in(\epsilon,K)\subset(0,\infty)$. Define

[TABLE]

Then these are the MLE estimators of $\delta_{\text{in}},\delta_{\text{out}}$ and they are strongly consistent; that is,

[TABLE]

Proof of Theorem 3.2.

We only verify the consistency of $\hat{\delta}_{\text{in}}^{MLE}$ since similar arguments apply to $\hat{\delta}_{\text{out}}^{MLE}$. Define

[TABLE]

Let us consider a limit version of $\psi_{n}$:

[TABLE]

where $p^{\text{in}}_{>i}(\delta_{\text{in}}):=\sum_{k>i}p^{\text{in}}_{k}(\delta_{\text{in}})$ with $p^{\text{in}}_{k}(\delta_{\text{in}}):=p^{\text{in}}_{k}$ as defined in (2.5), and

[TABLE]

Here we write $p^{\text{in}}_{i}(\delta_{\text{in}})$ to emphasize the dependence on $\delta_{\text{in}}$. In Lemmas A.1 and A.2, provided in the appendix, it is shown that $\psi(\cdot)$ has a unique zero at $\delta_{\text{in}}$, where $\psi(\lambda)>0$ when $\lambda<\delta_{\text{in}}$ and $\psi(\lambda)<0$ when $\lambda>\delta_{\text{in}}$, and

[TABLE]

Since $\psi$ is continuous, for any $\kappa>0$ arbitrarily small, there exists $\varepsilon_{\kappa}>0$ such that $\psi(\lambda)>\varepsilon_{\kappa}$ for $\lambda\in[\epsilon,\delta_{\text{in}}-\kappa]$ and $\psi(\lambda)<-\varepsilon_{\kappa}$ for $\lambda\in[\delta_{\text{in}}+\kappa,K]$. From (3.12),

[TABLE]

Note $\sup_{\lambda\in[\epsilon,K]}|\psi_{n}(\lambda)-\psi(\lambda)|<\varepsilon_{\kappa}/2$ implies

[TABLE]

and

[TABLE]

These jointly indicate that $\delta_{\text{in}}-\kappa\leq\hat{\delta}_{\text{in}}^{MLE}\leq\delta_{\text{in}}+\kappa$. Hence (3.13) implies

[TABLE]

for arbitrary $\kappa>0$. That is, $\hat{\delta}_{\text{in}}^{MLE}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\delta_{\text{in}}$.

∎

3.3 Asymptotic normality of MLE

In the following theorem, we establish the asymptotic normality for the MLE estimator

[TABLE]

Theorem 3.3**.**

Let $\hat{\boldsymbol{\theta}}_{n}^{MLE}$ be the MLE estimator for $\boldsymbol{\theta}$, the parameter vector of the preferential attachment model. Then

[TABLE]

where

[TABLE]

with

[TABLE]

In particular, $I(\boldsymbol{\theta})$ is the asymptotic Fisher information matrix for the parameters, and hence the MLE estimator is efficient.

Remark 3.4**.**

From Theorem 3.3, the estimators $(\hat{\alpha}^{MLE},\,\hat{\beta}^{MLE})$, $\hat{\delta}_{\text{in}}^{MLE}$, and $\hat{\delta}_{\text{out}}^{MLE}$ are asymptotically independent.

Proof of Theorem 3.3.

We first show the limiting distributions for $(\hat{\alpha}^{MLE},\,\hat{\beta}^{MLE})$, $\hat{\delta}_{\text{in}}^{MLE}$, and $\hat{\delta}_{\text{out}}^{MLE}$, respectively. From (3.6) and (3.7),

[TABLE]

where $\{J_{t}\}$ is a sequence of iid random variables. Hence the limiting distribution of the pair $\left(\hat{\alpha}^{MLE},\hat{\beta}^{MLE}\right)$ follows directly from standard central limit theorem for sums of independent random variables.

Next we show the asymptotic normality for $\hat{\delta}_{\text{in}}^{MLE}$; the argument for $\hat{\delta}_{\text{out}}^{MLE}$ is similar. Recall from (3.5) that the score function for $\delta_{\text{in}}$ can be written as

[TABLE]

where $u_{t}$ is defined by

[TABLE]

The MLE estimator $\hat{\delta}_{\text{in}}^{MLE}$ can be obtained by solving $\sum_{t=1}^{n}u_{t}(\delta)=0$. By a Taylor expansion of $\sum_{t=1}^{n}u_{t}(\delta)$,

[TABLE]

where $\dot{u}_{t}$ denotes the derivative of $u_{t}$ and $\hat{\delta}_{\text{in}}^{*}=\delta_{\text{in}}+\xi(\hat{\delta}_{\text{in}}^{MLE}-\delta_{\text{in}})$ for some $\xi\in[0,1]$. An elementary transformation of (3.18) gives

[TABLE]

To establish

[TABLE]

where $I_{\text{in}}$ is as defined in (3.15), it suffices to show the following two results:

1. (i)

$n^{-1/2}\sum_{t=1}^{n}u_{t}(\delta_{\text{in}})\stackrel{{\scriptstyle d}}{{\rightarrow}}N(0,I_{\text{in}})$, 2. (ii)

$n^{-1}\sum_{t=1}^{n}\dot{u}_{t}(\hat{\delta}_{\text{in}}^{*})\stackrel{{\scriptstyle p}}{{\rightarrow}}-I_{\text{in}}$.

These are proved in Lemmas A.3 and A.4 in the appendix, respectively.

To establish the joint asymptotic normality of the MLE estimator $\hat{\boldsymbol{\theta}}_{n}^{MLE}$, denote the joint score function vector for $\boldsymbol{\theta}$ by

[TABLE]

where $S_{n}(\alpha),S_{n}(\beta),S_{n}(\delta_{\text{in}}),S_{n}(\delta_{\text{out}})$ are the score functions for $\alpha,\beta,\delta_{\text{in}},\delta_{\text{out}}$, respectively. A multivariate Taylor expansion gives

[TABLE]

where $\dot{\mathbf{S}}_{n}$ denotes the Hessian matrix of the log-likelihood function $\log L(\boldsymbol{\theta})$, and $\hat{\boldsymbol{\theta}}_{n}^{*}=\boldsymbol{\theta}+\boldsymbol{\xi}\circ\left(\hat{\boldsymbol{\theta}}_{n}^{MLE}-\boldsymbol{\theta}\right)$ for some vector $\boldsymbol{\xi}\in[0,1]^{4}$, where “$\circ$” denotes the Hadamard product. From Remark 3.1, the likelihood function $L(\boldsymbol{\theta})$ can be factored into

[TABLE]

Hence

[TABLE]

as implied in the previous part of the proof, where $I(\boldsymbol{\theta})$ is as defined in (3.15) and is positive semi-definite.

Note that $(S_{n}(\alpha),S_{n}(\beta)),S_{n}(\delta_{\text{in}}),S_{n}(\delta_{\text{out}})$ are pairwise uncorrelated. As an example, observe that

[TABLE]

Using the Cramér-Wold device, the joint convergence of $\mathbf{S}_{n}(\boldsymbol{\theta})$ follows easily, i.e.,

[TABLE]

From here, the result of the theorem follows from (3.19) and (3.20).∎

4 Parameter estimation based on one snapshot

Based only on the single snapshot $G(n)$, we propose a parameter estimation procedure. We assume that the choice of the snapshot does not depend on any endogenous information related to the network. The snapshot merely represents a point in time where the data is available. Since no information on the initial graph $G(n_{0})$ is available, we merely assume $n_{0}$ and $N(n_{0})$ are fixed and $n\to\infty$.

Among the sufficient statistics for $(\alpha,\beta,\delta_{\text{in}},\delta_{\text{out}})$ derived in Remark 3.1, $\left(N^{\text{in}}_{>i}(n)\right)_{i\geq 0}$, $\left(N^{\text{out}}_{>j}(n)\right)_{j\geq 0}$ are computable from $G(n)$, but the $(J_{t})_{t=1}^{n}$ are not. However, when $n$ is large, we can use the following approximations according to the proof of Lemma A.2:

[TABLE]

and

[TABLE]

Substituting in (3.1), we estimate $\delta_{\text{in}}$ in terms of $\alpha$ and $\beta$ by solving

[TABLE]

Note that a strongly consistent estimator of $\beta$ can be obtained directly from $G(n)$:

[TABLE]

To obtain an estimate for $\alpha$, we make use of the recursive formula for $\{p^{\text{in}}_{i}\}$ in (A.1a):

[TABLE]

and replace $p^{\text{in}}_{0}$ by $N^{\text{in}}_{0}(n)/n$ for large $n$,

[TABLE]

Plug the strongly consistent estimator $\tilde{\beta}$ into (4.1) and (4.3), and we claim that solving the system of equations:

[TABLE]

gives the unique solution $(\tilde{\alpha},\tilde{\delta}_{\text{in}})$ which is strongly consistent for $(\alpha,\delta_{\text{in}})$.

Theorem 4.1**.**

The solution $(\tilde{\alpha},\tilde{\delta}_{\text{in}})$ to the system of equations in (4.4) is unique and strongly consistent for $(\alpha,\delta_{\text{in}})$, i.e.

[TABLE]

The proof of Theorem 4.1 is given in Section A.3.

The parameters $\tilde{\delta}_{\text{out}}$ and $\tilde{\gamma}$ can be estimated by a mirror argument. We summarize the estimation procedure for $(\alpha,\beta,\gamma,\delta_{\text{in}},\delta_{\text{out}})$ from the snapshot $G(n)$ as follows:

Estimate $\beta$ by $\tilde{\beta}=1-N(n)/n$. 2. 2.

Obtain $\tilde{\delta}_{\text{in}}^{0}$ by solving (i.e., matching (4.4a) and (4.4b))

[TABLE] 3. 3.

Estimate $\alpha$ by

[TABLE] 4. 4.

Obtain $\tilde{\delta}_{\text{out}}^{0}$ by solving

[TABLE] 5. 5.

Estimate $\gamma$ by

[TABLE]

Note that even though all three estimators $\tilde{\alpha}^{0},\tilde{\beta},\tilde{\gamma}^{0}$ are strongly consistent and hence $\tilde{\alpha}^{0}+\tilde{\beta}+\tilde{\gamma}^{0}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}1$, Step 1–5 do not necessarily imply the strict equality

[TABLE]

We recommend adding the following two steps for a re-normalization to overcome this defect.

Re-normalize the probabilities

[TABLE] 2. 7.

Plug $\tilde{\alpha}$ into (4.4a) to update the estimate of $\delta_{\text{in}}$, i.e., solve for $\tilde{\delta}_{\text{in}}$ from

[TABLE]

Similarly, solve for $\tilde{\delta}_{\text{out}}$ from

[TABLE]

5 Simulation study

We now apply the estimation procedures described in Sections 3 and 4 to simulated data, which allows us to compare the estimation results using the full history of the network with that using just one snapshot. Algorithm 1 is used to simulate realizations of the preferential attachment network.

5.1 MLE

For the scenario of observing the full history of the network, we simulated 5000 independent replications of the preferential attachment network with $10^{5}$ edges under the true parameter values

[TABLE]

For each realization, the MLE estimate of the parameters was computed and standardized as

[TABLE]

where $(\hat{\boldsymbol{\theta}}_{n})_{i}$ and $(\boldsymbol{\theta})_{i}$ denote the $i$-th components of $\hat{\boldsymbol{\theta}}_{n}^{MLE}$ and $\boldsymbol{\theta}$ respectively, and $\hat{\sigma}_{ii}^{2}$ is the $i$-th diagonal component of the matrix $\hat{\Sigma}:=\Sigma(\hat{\boldsymbol{\theta}}_{n}^{MLE})$. The explicit formula for the entries of $\hat{\Sigma}$ is

[TABLE]

where, see (3.15) and (3.16),

[TABLE]

By the consistency of the MLEs combined with the convergence of $\{N^{\text{in}}_{>i}(n)/n\}$ and $\{N^{\text{out}}_{>j}(n)/n\}$, see (A.8), we have that $\hat{\Sigma}_{n}\overset{a.s.}{\to}\Sigma$.

The QQ-plots of the normalized MLEs are shown in Figure 5.1, all of which line up quite well with the $y=x$ line (the red dashed line). This is consistent with the asymptotic theory described in Theorem 3.3. Confidence intervals for $\boldsymbol{\theta}$ can be obtained using this theorem. Given a single realization, an approximate $(1-\varepsilon)$-confidence interval for $(\boldsymbol{\theta})_{i}$ is

[TABLE]

where $z_{\varepsilon/2}$ is the upper $\varepsilon/2$ quantile of $N(0,1)$.

5.2 One snapshot

We used the same simulated data as in Section 5.1 to obtain parameter estimates $\tilde{\boldsymbol{\theta}}_{n}:=(\tilde{\alpha},\,\tilde{\beta},\,\tilde{\delta}_{\text{in}},\,\tilde{\delta}_{\text{out}})$ through only the final snapshot, i.e., the set of directed edges without timestamps, following the procedure described at the end of Section 4. For the purpose of comparison with MLE, Figure 5.2 gives the QQ-plots for the normalized estimates from the snapshots using the same standardizations for the MLEs, i.e.,

[TABLE]

where $(\tilde{\boldsymbol{\theta}}_{n})_{i}$ denotes the $i$-th components of $\tilde{\boldsymbol{\theta}}_{n}$. Again, the fitted lines in blue are the traditional QQ-lines and the red dashed lines are the $y=x$ line. The QQ-plot for $\tilde{\beta}$ exhibits the same shape as for $\hat{\beta}^{MLE}$, since the two estimates are identical.

From Figure 5.2, we see that the snapshot estimates of all four parameters are consistent and approximately normal, i.e., the QQ-plots are linear. However, the slopes of the QQ-lines for $\tilde{\alpha},\tilde{\delta}_{\text{in}},\tilde{\delta}_{\text{out}}$ are much steeper than the diagonal line, indicating a loss of efficiency for $\tilde{\boldsymbol{\theta}}_{n}$ compared with $\hat{\boldsymbol{\theta}}_{n}$. Indeed the estimator variance is inflated for all parameters except for $\beta$, where $\tilde{\beta}$ coincides with the true MLE. This is as expected since knowing only the final snapshot provides far less information than the whole network history.

Recall that for a consistent estimator $T_{n}$ of a one-dimensional parameter $\theta$ constructed from a random sample of size $n$, the asymptotic relative efficiencies (ARE) of $T_{n}$ is defined by

[TABLE]

where $T_{n}^{*}$ denotes the asymptotically efficient estimator. We may compute the ARE’s for the snapshot parameter estimates

[TABLE]

where $\widehat{\text{Var}}$ denotes the sample variance of the parameter estimate based on the 5000 replications. Note that $ARE(\tilde{\beta})=1$ since $\tilde{\beta}=\hat{\beta}^{MLE}$.

Given a single realization, the variances of the snapshot estimates can be estimated through resampling as follows.

Using the estimated parameter $\tilde{\boldsymbol{\theta}}_{n}$, simulate $10^{4}$ independent bootstrap replicates of the network with $n=10^{5}$ edges. For each simulated network, the snapshot estimate, $\tilde{\boldsymbol{\theta}}^{*}_{n}:=\left(\tilde{\alpha}^{*},\,\tilde{\beta}^{*},\,\tilde{\delta}_{\text{in}}^{*},\,\tilde{\delta}_{\text{out}}^{*}\right)$, is computed. The sample variance of these $10^{4}$ snapshot estimates can then be used as an approximation for the variance of $\tilde{\boldsymbol{\theta}}_{n}$ so that assuming asymptotic normality, a $(1-\varepsilon)$-confidence interval for $\boldsymbol{\theta}$ can be approximated by

[TABLE]

where $z_{\varepsilon/2}$ is the upper $\varepsilon/2$ quantile of $N(0,1)$.

5.3 Sensitivity test

Now we investigate the sensitivity of our estimates while values of the parameters $(n,\alpha,\beta,\delta_{\text{in}},\delta_{\text{out}})$ are allowed to vary. First consider the impact of $n$, the number of edges in the network. To do so we held the parameters fixed with values given by (5.1): $\left(\alpha,\beta,\delta_{\text{in}},\delta_{\text{out}}\right)=(0.3,\,0.5,\,2,\,1)$ and varied the value of $n$. The QQ-plots (not presented) for standardized estimates using both full MLE and one-snapshot methods were produced to check the asymptotic normality. When $n=500,1000$, diagnostics revealed departures from normality for both the MLE and the snapshot estimates. However, after increasing $n$ to $10000$, estimates obtained from both approaches appeared normally distributed as expected.

For each value of $n$ in Table 5.1, 5000 replicates of the network with $n$ edges and parameters $\boldsymbol{\theta}=(0.3,0.5,2,1)$ were generated. For each realization, the MLE’s $\hat{\boldsymbol{\theta}}_{n}^{MLE}$ were computed using the full history of the network and the one-snapshot estimates $\tilde{\boldsymbol{\theta}}_{n}$ were obtained using the 7-step snapshot method proposed in Section 4, pretending that only the last snapshot $G(n)$ was available. The mean for these two estimators were recorded in Table 5.1. There is little bias for both estimates of $\alpha$ and $\beta$, even for small values of $n$. On the other hand, there is some bias for estimated $\delta_{\text{in}}$ and $\delta_{\text{out}}$ for $n\leq 5000$. The magnitude of the biases for both types of estimates decrease as $n$ increases. Also the ARE’s of the snapshot estimator stay within a narrow band as $n$ increases.

Next we held $(n,\delta_{\text{in}},\delta_{\text{out}})=(10^{5},2,1)$ fixed and experimented with various values of $(\alpha,\beta)$ in Table 5.2. For each choice of $(\alpha,\beta)$, 5000 independent realizations of the network were generated and the means of the MLE $\hat{\boldsymbol{\theta}}_{n}^{MLE}$ and the one-snapshot estimates $\tilde{\boldsymbol{\theta}}_{n}$ were recorded. Overall, the biases for $\hat{\boldsymbol{\theta}}_{n}^{MLE}$ are remarkably small for virtually all combinations of parameter values, except for those parameter choices where one of $(\alpha,\beta)$ is extremely small. The biases for the snapshot estimates $\tilde{\boldsymbol{\theta}}_{n}$ exhibit a similar property, but the magnitudes of the biases are consistently larger than those in the MLE case.

In general, the snapshot estimators are able to achieve $20\%$–$50\%$ efficiency over the range of parameters considered. The loss of efficiency might be less than one would expect given the substantial reduction in the data available to produce the snapshot estimates. It is worth noting that in the case where $(\alpha,\beta)=(0.7,0.2)$, the efficiencies of the snapshot estimators for $\alpha$ and $\delta_{\text{in}}$ are much larger (.73 and .79, respectively). A heuristic explanation for this increase is that the parameter $\gamma=1-\alpha-\beta=0.1$ is relatively small. By the implicit constraints used for the snapshot estimates, we have

[TABLE]

that is, the snapshot estimate of the sum $\alpha+\gamma$ is the same as the MLE for the sum. Now if $\gamma$ is small, one would expect the resulting estimates to also be small so that $\tilde{\alpha}$ would be nearly the same as $\hat{\alpha}^{MLE}$. Hence the ARE would be close to 1. On the other hand, in the case of a larger $\gamma$, see the bottom row of Table 5.2 in which $\gamma=0.6$, the ARE for $\alpha$ is not as large (.42), but the ARE for $\tilde{\delta}_{\text{out}}$ is (.63).

6 Real network example

In this section, we explore fitting a preferential attachment model to a social network. As illustration, we chose the Dutch Wiki talk network dataset, available on KONECT [14] (http://konect.uni-koblenz.de/networks/wiki_talk_nl). The nodes represent users of Dutch Wikipedia, and an edge from node A to node B refers to user A writing a message on the talk page of user B at a certain time point. The network consists of 225,749 nodes (users) and 1,554,699 edges (messages). All edges are recorded with timestamps.

In order to accommodate all the edge formulation scenarios appeared in the dataset, we extend our model by appending the following two interaction schemes ($J_{n}=4,5$) in addition to the existing three ($J_{n}=1,2,3$) described in Section 2.1.

• •

If $J_{n}=4$ (with probability $\xi$), append to $G(n-1)$ two new nodes $v,w\in V(n)\setminus V(n-1)$ and an edge connecting them $(v,w)$.

• •

If $J_{n}=5$ (with probability $\rho$), append to $G(n-1)$ a new node $v\in V(n)\setminus V(n-1)$ with self loop $(v,v)$.

These scenarios have been observed in other social network data, such as the Facebook wall post network (http://konect.uni-koblenz.de/networks/facebook-wosn-wall), etc. They occur in small proportions and can be easily accommodated by a slight modification in the model fitting procedure. The new model has parameters $(\alpha,\beta,\gamma,\xi,\delta_{\text{in}},\delta_{\text{out}})$, and $\rho$ is implicitly defined through $\rho=1-(\alpha+\beta+\gamma+\xi)$. Similar to the derivations in Section 3, the MLE estimators for $\alpha,\beta,\gamma,\xi$ are

[TABLE]

and $\delta_{\text{in}},\delta_{\text{out}}$ can be obtained through solving

[TABLE]

We first naively fit the linear preferential attachment model to the full network using MLE. The MLE estimators are

[TABLE]

To evaluate the goodness-of-fit, 20 network realizations were simulated from the fitted model. We overlaid the empirical in- and out-degree frequencies of the original network with that of the simulations. If the model fits the data well, the degree frequencies of the data should lie within the range formed by that of the simulations, which gives an informal confidence region for the degree distributions. From Figure 6.1, we see that while the data roughly agrees with the simulations in the out-degree frequencies, the deviation in the in-degree frequencies is noticeable.

To better understand the discrepancy in the in-degree frequencies, we examined the link data and their timestamps and discovered bursts of messages originating from certain nodes over small time intervals. According to Wikipedia policy [23], certain administrating accounts are allowed to send group messages to multiple users simultaneously. These bursts presumably represent broadcast announcements generated from these accounts. These administrative broadcasts can also be detected if we apply the linear preferential attachment model to the network in local time intervals. We divided the total time frame down to sub-intervals of varying length each containing the formation of $10^{4}$ edges. The number $10^{4}$ is chosen to ensure good asymptotics as shown in Table 5.1. This process generated 155 networks,

[TABLE]

For each of the 155 datasets, we fit a preferential attachment model using MLE. The resulting estimates $(\hat{\delta}_{\text{in}},\hat{\delta}_{\text{out}})$ are plotted against the corresponding timeline on the upper left panel of Figure 6.2. Notice that $\hat{\delta}_{\text{in}}$ exhibits large spikes at various times. Recall from (2.1), a large value of $\delta_{\text{in}}$ indicates that the probability of an existing node $v$ receiving a new message becomes less dependent on its in-degree, i.e., previous popularity. These spikes appear to be directly related to the occurrences of group messages. This plot is truncated after the day 2016/3/16, on which a massive group message of size 48,957 was sent and the model can no longer be fit.

We identified 37 users who have sent, at least once, 40 or more consecutive messages in the message history. This is evidence that group messages were sent by this user. We presume these nodes are administrative accounts; they are responsible for about $30\%$ of the total messages sent. Since their behavior cannot be regarded as normal social interaction, we excluded messages from these accounts from the dataset in our analysis. We then also removed nodes with zero in- and out-degrees.

The re-estimated parameters after the data cleaning are displayed in the other three panels of Figure 6.2. Here all parameter estimates are quite stable through time.

The reduced network now contains 112,919 nodes and 1,086,982 edges, to which we fit the linear preferential attachment model. The fitted parameters based on MLE for our reduced dataset are

[TABLE]

Again the degree distributions of the data and 20 simulations from the fitted model are displayed in Figure 6.3. The out-degree distribution of the data agrees reasonably well with the simulations. For the in-degree distribution, the fit is better than that for the entire dataset (Figure 6.1). However, for smaller in-degrees, the fitted model over-estimates the in-degree frequencies. We speculate that in many social networks, the out-degree is in line with that predicted by the preferential attachment model. An individual node would be more likely to reach out to others if having done so many times previously. For in-degrees, the situation is complicated and may depend on a multitude of factors. For instance, the choice of recipient may depend on the community that the sender is in, the topic being discussed in the message, etc. As an example a group leader might send messages to his/her team on a regular basis. Such examples violate the base assumptions of the preferential attachment model and could result in the deviation between the data and the simulations.

Next we consider the estimation method of Section 4 applied to a single snapshot of the data. In order to implement this procedure, we donned blinders and assumed that our dataset consists only of the information of the wiki data at the last timestamp. That is, information about administrative broadcasts, and other aspects of the data learned by looking at the previous history of the data are unavailable. In particular, we would have no knowledge of the existence of the two additional scenarios corresponding to $J_{n}=4,5$. With this in mind, we fit the three scenario model using the methods in Section 4. The fitted parameters are

[TABLE]

The comparison of the degree distributions between the data and simulations from the fitted model is displayed in Figure 6.4 and is not too dissimilar to the plots in Figure 6.1 that are based on maximum likelihood estimation using the full network data. In particular, the out-degree distribution is matched reasonably well, but the fitted model does a poor job of capturing the in-degree distribution.

We see from this example that while the linear preferential attachment model is perhaps too simplistic for the Wiki talk network dataset, it has the ability to illuminate some gross features, such as the out-degrees, as well as to capture important structural changes such as the group message behavior. Consequently, despite its limitation, this model may be used as a building block for more flexible models. Modification to the existing model formulation and more careful analysis of change points in parameters is a direction for future research.

7 Acknowledgement

Research of the four authors was partially supported by Army MURI grant W911NF-12-1-0385. Don Towsley from University of Massachusetts introduced us to the model and within his group, James Atwood graciously supplied us with a simulation algorithm designed for a class of growth models broader than the one specified in Section 2.1; this later became [1]. Joyjit Roy, formerly of Cornell, created an efficient algorithm designed to capitalize on the linear growth structure. Finally, we appreciate the many helpful and sensible comments of the referees and editors.

Appendix A Proofs

A.1 For the proof of Theorem 3.2: Lemmas A.1 and A.2

Lemma A.1**.**

For ${\lambda>0}$, the function $\psi(\lambda)$ in (3.11) has a unique zero at $\delta_{\text{in}}$ and, $\psi(\lambda)>0$ when $\lambda<\delta_{\text{in}}$ and $\psi(\lambda)<0$ when $\lambda>\delta_{\text{in}}$.

Proof.

The probabilities $\{p^{\text{in}}_{i}(\lambda)\}$ satisfy the recursions in $i$ (cf. [3]):

[TABLE]

where $a_{1}(\lambda):=(\alpha+\beta)/(1+\lambda(1-\beta))$. Summing the recursions in (A.1) from [math] to $i$, we get (with the convention that $\sum_{i=0}^{-1}=0$)

[TABLE]

which can be simplified to

[TABLE]

From (2.4),

[TABLE]

Hence by rearranging (A.2), we have

[TABLE]

or equivalently,

[TABLE]

Now with the help of (A.3) and (A.4), we can rewrite $\psi(\lambda)$ in the following way:

[TABLE]

The series defining $C(\lambda)$ converges absolutely for any $\lambda>0$ since

[TABLE]

Summing over $i$ in $\eqref{sumb}$, we get by monotone convergence

[TABLE]

The infinite series converge because $p^{\text{in}}_{i}(\lambda)$ is a power law with index greater than 2; see (2.5) and (2.6). Solving for the infinite series we get

[TABLE]

Hence we have

[TABLE]

Now recall from (A.5) that $\psi(\lambda)$ is of the form

[TABLE]

where $C(\lambda)>0$ for all ${\lambda>0}$. Therefore $\psi(\cdot)$ has a unique zero at $\delta_{\text{in}}$ and $\psi(\lambda)>0$ when $\lambda<\delta_{\text{in}}$ and $\psi(\lambda)<0$ when $\lambda>\delta_{\text{in}}$.∎

We show the uniform convergence of $\psi_{n}$ to $\psi$ in the next lemma.

Lemma A.2**.**

As $n\to\infty$, for any $\epsilon>0$,

[TABLE]

Proof.

By the definition of $\psi$, $p^{\text{in}}_{>i}(\delta_{\text{in}})$ is a function of $\delta_{\text{in}}$ and is a constant with respect to $\lambda$. Hence we suppress the dependence on $\delta_{\text{in}}$ and simply write it as $p^{\text{in}}_{>i}$ when considering the difference $\psi_{n}-\psi$ as a function of $\lambda$:

[TABLE]

Thus,

[TABLE]

For the first term, note that for all $i\geq 0$,

[TABLE]

since the assumption on initial conditions implies the sum of in-degrees at $n$ is $n$. Therefore $N^{\text{in}}_{>i}(n)/n\leq i^{-1}$ for $i\geq 1$, and it then follows that

[TABLE]

Note that the last two terms on the right side can be made arbitrarily small uniformly on $[\epsilon,{\infty)}$ if we choose $M$ sufficiently large.

Recall the convergence of the degree distribution $\{N_{ij}(n)/N(n)\}$ to the probability distribution $\{f_{ij}\}$ in (2.4), we have

[TABLE]

Hence, for any fixed $M$,

[TABLE]

which implies further that choosing $M$ arbitrarily large gives

[TABLE]

The second term in (A.1) converges to 0 almost surely by strong law of large numbers, and the third term in (A.1) can be written as

[TABLE]

which is bounded by

[TABLE]

We have

[TABLE]

which converges to 0 almost surely by Cesàro convergence of random variables, since

[TABLE]

Further, by the strong law of large numbers,

[TABLE]

Hence the third term of (A.1) also goes to 0 almost surely as $n\to\infty$. The result of the lemma follows. ∎

A.2 For the proof of Theorem 3.3: Lemmas A.3 and A.4

Lemma A.3**.**

As $n\to\infty$,

[TABLE]

Proof.

Let $\mathcal{F}_{n}=\sigma(G(0),\ldots,G(n))$ be the $\sigma$-field generated by the information contained in the graphs. We first observe that $\{\sum_{t=1}^{n}u_{t}(\delta_{\text{in}}),\mathcal{F}_{n},n\geq 1\}$ is a martingale. To see this, note from (3.17) that $|u_{t}(\delta)|\leq 2/\delta$ and

[TABLE]

which satisfies the definition of a martingale difference. Hence $\left\{n^{-1/2}\sum_{r=1}^{t}u_{r}(\delta_{\text{in}})\right\}_{t=1,\ldots,n}$ is a zero-mean, square-integrable martingale array. The convergence (A.9) follows from the martingale central limit theory (cf. Theorem 3.2 of [7]) if the following three conditions can be verified:

1. (a)

$n^{-1/2}\max_{t}|u_{t}(\delta_{\text{in}})|\stackrel{{\scriptstyle p}}{{\rightarrow}}0$, 2. (b)

$n^{-1}\sum_{t}u_{t}^{2}(\delta_{\text{in}})\stackrel{{\scriptstyle p}}{{\rightarrow}}I_{\text{in}}$, 3. (c)

$\textbf{E}\left(n^{-1}\max_{t}u_{t}^{2}(\delta_{\text{in}})\right)$ is bounded in $n$.

Since $|u_{t}(\delta_{\text{in}})|\leq 2/\delta_{\text{in}}$, we have

[TABLE]

and

[TABLE]

Hence conditions (a) and (c) are straightforward.

To show (b), observe that

[TABLE]

Following the calculations in the proof of Lemma A.2, we have for $T_{1}$,

[TABLE]

We then rewrite $T_{2}$ as

[TABLE]

where

[TABLE]

by Cesàro’s convergence and

[TABLE]

where the equality follows from (A.4). For $T_{3}$, similar to $T_{1}$, we have

[TABLE]

Combining these results together,

[TABLE]

This completes the proof. ∎

Lemma A.4**.**

As $n\to\infty$,

[TABLE]

Proof.

The result of this lemma can be established by showing first

[TABLE]

and then

[TABLE]

We first observe that

[TABLE]

Recall the definition and convergence result for $T_{2}$ and $T_{3}$ in Lemma A.3, we have

[TABLE]

Also from (A.10),

[TABLE]

Hence

[TABLE]

and (A.11) is established.

By construction and definition, we have $\hat{\delta}_{\text{in}},\hat{\delta}_{\text{in}}^{*},\delta_{\text{in}}>0$. To prove (A.12), note that

[TABLE]

Then

[TABLE]

and

[TABLE]

From Theorem 3.2, $\hat{\delta}_{\text{in}}^{MLE}$ is consistent for $\delta_{\text{in}}$, hence

[TABLE]

We have

[TABLE]

This proves (A.12) and completes the proof of Lemma A.4. ∎

A.3 Proof of Theorem 4.1

Proof.

First observe that $\sum_{i}iN^{\text{in}}_{i}(n)$ sums up to the total number of edges $n$, so

[TABLE]

We can re-write (4.4a) as

[TABLE]

and (4.4b) as

[TABLE]

Then $\tilde{\delta}_{\text{in}}$ can be obtained by solving

[TABLE]

Similar to the proof of Theorem 3.2, we define the limit versions of $f_{n}$, and $g_{n}$ as follows:

[TABLE]

Now we apply the re-parametrization

[TABLE]

to $f$ and $g$, such that

[TABLE]

Note that for all $\eta\in\mathcal{I}$:

• •

Set $b_{i}(\eta):=(i^{-1}-(1-\beta))\eta$, then $1+b_{i}(\eta)>0$ for all $i\geq 1$. So $\tilde{f}(\eta)>0$ on $\mathcal{I}$;

• •

$\tilde{f}(\eta)\leq\frac{1}{1-(1-\beta)\eta}\sum_{i=0}^{\infty}p^{\text{in}}_{>i}\leq 1+(1-\beta)K<\infty$.

Meanwhile, $\tilde{g}$ is also well defined and strictly positive for $\eta\in\mathcal{I}$ because

[TABLE]

The first inequality holds since:

[TABLE]

We know $\alpha+\beta<1$ by our model assumption, thus verifying (A.15).

Define for $\eta\in\mathcal{I}$,

[TABLE]

then it follows that

[TABLE]

We now show that $\tilde{h}$ is concave and $\tilde{h}(\eta)\to 0$ as $\eta\to 0$, then the uniqueness of the solution follows.

First observe that

[TABLE]

We now claim that

[TABLE]

It suffices to check:

[TABLE]

Note that for $i\geq 1$,

[TABLE]

Recall (A.6), we then have

[TABLE]

Hence,

[TABLE]

which implies (A.17). Equation (A.18) then follows by a similar argument. Combining (A.16), (A.17) and (A.18) gives

[TABLE]

by the Cauchy-Schwarz inequality. Hence $\tilde{h}$ is concave on $\mathcal{I}$.

From Lemma A.1, $\psi(\delta_{\text{in}})=0$ where $\psi(\cdot)$ is as defined in (3.11). Hence we have $f(\delta_{\text{in}})=\alpha+\beta$ in a similar derivation to that of (A.13). Also from (4.2), we have $g(\delta_{\text{in}})=\alpha+\beta$. Hence, $\delta_{\text{in}}$ is a solution to $f(\delta)=g(\delta)$.

Under the $\delta\mapsto\eta$ reparametrization in (A.14), we have that $\tilde{f}(\eta_{\text{in}})=\tilde{g}(\eta_{\text{in}})$ where $\eta_{\text{in}}:=\delta_{\text{in}}/(1+\delta_{\text{in}}(1-\beta))$, and also

[TABLE]

This, along with the concavity of $\tilde{h}$, implies that $\eta_{\text{in}}$ is the unique solution to $\tilde{h}(\eta)=0$, or equivalently, to $\tilde{f}(\eta)=\tilde{g}(\eta)$ on $\mathcal{I}$.

Let $\tilde{f}_{n}(\eta):=f_{n}(\delta(\eta))$, $\tilde{g}_{n}(\eta):=g_{n}(\delta(\eta))$. We can show in a similar fashion that $\tilde{\eta}:=\tilde{\delta}_{\text{in}}/(1-\tilde{\delta}_{\text{in}}(1-\tilde{\beta}))$ is the unique solution to $\tilde{f}_{n}(\eta)=\tilde{g}_{n}(\eta)$. Using an analogue of the arguments in the proof of Theorem A.2, we have

[TABLE]

and therefore $\tilde{\eta}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\eta_{\text{in}}$. Since $\delta\mapsto\eta$ is a one-to-one transformation from $[\epsilon,K]$ to $\mathcal{I}$, we have that $\tilde{\delta}_{\text{in}}$ is the unique solution to $f_{n}(\delta)=g_{n}(\delta)$ and that $\tilde{\delta}_{\text{in}}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\delta_{\text{in}}$. On the other hand, $\tilde{\alpha}$ can be solved uniquely by plugging $\tilde{\delta}_{\text{in}}$ into (A.13) and is also strongly consistent, which completes the proof.

∎