A network Poisson model for weighted directed networks with covariates

1 Introduction

1.1 Background

Many complex interactive behaviors can be conveniently represented as networks or graphs, where nodes denotes entities depending on different contexts and edges denote interactions between entities. Examples include friendships between people in social networks, emails between users in email networks, hyperlinks between internet webs in hyperlink networks, citations between papers and authors in citation networks, following behaviors between blogs in social media such as Twitter, chemical reactions between proteins in biological networks, to just name a few. Many statistical methodologies have been developed to analyze network data; see Goldenberg et al. (2010), Fienberg (2012), Robins et al. (2007a), Robins et al. (2007b), Bhattacharyya and Bickel (2016) and Kim et al. (2018) for reviews and references therein. The monograph by Kolaczyk (2009) provides a comprehensive introduction on statistical analysis of network data.

Networks could be undirected or directed, weighted or unweighted. Most realistic networks exhibit three typical features including sparsity, degree heterogeneity and homophily. Sparsity means that the density of networks is small, in which many nodes do not have direct connections. The degree heterogeneity describes such phenomenon that degrees of nodes vary greatly. Some nodes may have many connections while others may have relatively less connections. The homophily characterizes the tendency that individuals with same attributes such as age and sex are easier to form connections. For example, the directed friendship network in Lazega,E. (2001) displays a strong homophily effect as shown in Yan et al. (2019).

One of the most popular models to model the degree heterogeneity in undirected networks is the $β$ -model [Chatterjee et al. (2011)] that assigns one degree parameter to each node. It is an undirected version of the well-known $p_{1}$ model in Holland and Leinhardt (1981). Asymptotic theory is nonstandard because of an increasing dimension of parameter space. Exploring theoretical properties in the $β$ -model and its variants has received wide attentions in recent years [Chatterjee et al. (2011); Perry and Wolfe (2012); Hillar and Wibisono (2013); Yan and Xu (2013); Rinaldo et al. (2013); Yan et al. (2016b); Graham (2017); Yan (2018); Chen et al. (2021)]. In particular, Chatterjee et al. (2011) proved the uniform consistency of the maximum likelihood estimator (MLE); Yan and Xu (2013) derived its asymptotic normality by approximating the inverse of the Fisher information matrix. In the directed case, Yan et al. (2016a) proved the consistency and asymptotic normality of the MLE in the $p_{0}$ model that is an exponential random graph model with out-degree and in-degree sequences as sufficient statistics. In the framework of non-exponential random graph model, Wang et al. (2020)

proposed Probit model to model the degree heterogeneity of the directed networks and proved the consistency and asymptotic normality of the moment estimator. Besides,

Wang et al. (2021) use the probit distribution to model the degree heterogeneity of the affiliation networks and established the uniform consistency and the asymptotic normality of the moment estimator.

Yan et al. (2019) proposed the covariate- $p_{0}$ model to model the aforementioned three network features in unweighted directed networks. They established the consistency and asymptotic normality of the restricted MLE by using the restricted maximum likelihood method because of the challenge of exploring asymptotic theory [Fienberg (2012); Graham (2017)]. Wang (2021) further incorporated a sparsity parameter to the covariate- $p_{0}$ model to allow the sparsity and developed the unrestricted maximum likelihood theory including the consistency and asymptotic normality of the MLE. However, the covariate- $p_{0}$ model is only designed to unweighted directed networks. When we apply it to weighted networks, we need to neglect the weight information (i.e., treating all positive weight value as “1” and others as “0”). This may cause the information loss. As one example, all covariates are not significant when we apply the covariate- $p_{0}$ model to fit the well-known Enran email network [Cohen (2004)]; see Table 5. This motivates the present paper. We extend the covariate- $p_{0}$ model to weighted networks by using the Poisson distribution to model weighted edges.

1.2 The Model

We now introduce our model. Consider a weighted directed graph $G_{n}$ on $n$ ( $n ⩾ 2$ ) nodes labeled by $1, \dots, n$ . Let $A = (a_{i j})_{n \times n}$ be the adjacency matrix $G_{n}$ , where $a_{i j} \in {0, 1, \dots}$ is the weight of the directed edge from head node $i$ to tail node $j$ . We do not consider self-loops here, i.e., $a_{i i} = 0$ . Let $d_{i} = \sum_{j = 1, j \neq i}^{n} a_{i j}$ be the out-degree of vertex $i$ and $d = (d_{1}, \dots, d_{n})$ be the out-degree sequence of the graph G. Similarly, define $b_{j} = \sum_{i = 1, i \neq j}^{n} a_{i j}$ as the in-degree of vertex $j$ and $b = (b_{1}, \dots, b_{n})$ as the in-degree sequence. The pair ${b, d}$ or ${(b_{1}, d_{1}), \dots, (b_{n}, d_{n})}$

is the bi-degree sequence. We assume that all edges are independently distributed as Poisson random variables with the probability distributions:

P (a_{i j} = k) = \frac{e^{k (μ + α_{i} + β_{j} + Z_{i j}^{⊤} γ)}}{k!} exp (- e^{μ + α_{i} + β_{j} + Z_{i j}^{⊤} γ}), 1 \leq i \neq j \leq n .

(1)

The parameter $μ$ quantifies the network sparsity. $α_{i}$ quantifies the effect of establishing outbound edges from sender $i$ while $β_{j}$ quantifies the effect of attracting inbound edges from receiver $j$

. The vector parameter

γ

is a

p

-dimensional regression coefficient for the covariate

Z_{i j}

. We will call the above model the network Poisson model hereafter.

The covariate $Z_{i j}$ is either a vector associated with edges or a function of node-specific covariates. If $X_{i}$ denotes a $p$ -dimensional vector of node-level attributes, then these node-level attributes can be used to construct a vector $Z_{i j} = g (X_{i}, X_{j})$ , where $g (\cdot, \cdot)$ is a function of its arguments. As one example, if we let $g (X_{i}, X_{j})$ being equal to $∥ X_{i} - X_{j} ∥_{1}$ , then it measures the similarity between nodes $i$ and $j$ .

Motivated by techniques for the analysis of the unrestricted likelihood inference in the covariate- $p_{0}$ model for directed graphs in Wang (2021), we generalize their approaches to weighted directed graphs here. When the number of nodes $n$ goes to infinity, we derive the $ℓ_{\infty}$ -error between the MLE $(ˆ θ, ˆ γ)$ and its true value $(η, γ)$ . This is done by using a two-stage Newton process that first finds the error bound between ${ˆ θ}_{γ}$ and $θ$ with a fixed $γ$ and then derives the error bound between $ˆ γ$ and $γ$ . They are $O_{p} ((log n / n)^{1 / 2})$ for $ˆ η$ and $O_{p} (log n / n)$ for $ˆ γ$ , up to an additional factor on parameters. Further, we derive the asymptotic normality of the MLE. The asymptotic distribution of $ˆ γ$ has a bias term while $ˆ θ$ does not have such a bias, which collaborates the findings in Yan et al. (2019). This is because of different convergence rates for $ˆ γ$ and $ˆ θ$ . Wide simulations are carried out to demonstrate our theoretical findings. In our simulations, this bias is very small, even could be neglected, which is different from the significant bias effect in Yan et al. (2019) and Wang (2021). This may be due to that the weighted values attenuate the bias. The application to the Enran email data set illustrates the utility of the proposed model.

For the remainder of the paper, we proceed as follows. In Section 2, we give the maximum likelihood estimation. In section 3, we present theoretical properties of the MLE. Numerical studies are presented in Section 4. We provide further discussion in Section 5. The proofs of theorems are relegated to the Appendix. All supported lemmas and detailed calculations are in the Supplementary Material.

2 Maximum Likelihood Estimation

Let $α = (α_{1}, \dots, α_{n})^{⊤}$ and $β = (β_{1}, \dots, β_{n})^{⊤}$ . If one transforms $(μ, α, β)$ to $(μ + 2 c_{1}, α - c_{1} + c_{2}, β - c_{1} - c_{2})$

, then the probability in (

1) does not change. For the identifiability of the model, several possible restriction conditions immediately appear in our mind, including

α_{n} = 0

β_{n} = 0

, or

μ = 0

α_{n} = 0

, or

μ = 0

β_{n} = 0

. When we

α_{n} = 0

β_{n} = 0

, it will keep the density parameter

μ

The logarithm of the likelihood function is

ℓ (μ, α, β, γ) = \sum i \neq j (a_{i j} (μ + α_{i} + β_{j} + Z_{i j}^{⊤} γ) - e^{μ + α_{i} + β_{j} + Z_{i j}^{⊤} γ} - log (a_{i j}!))

(2)

The notation $\sum_{i, j = 1, i \neq j}^{n}$ is a shorthand for $\sum_{i = 1}^{n} \sum_{j = 1, j \neq i}^{n}$ . The score equations for the vector parameters $μ, α, β$ , $γ$ are easily seen as

\begin{matrix} n \sum i, j = 1, i \neq j a_{i j} & = & \sum i \neq j e^{μ + α_{i} + β_{j} + Z_{i j}^{⊤} γ}, n \sum i, j = 1, i \neq j a_{i j} Z_{i j} & = & \sum i \neq j Z_{i j} e^{μ + α_{i} + β_{j} + Z_{i j}^{⊤} γ}, d_{i} & = & n \sum j = 1, j \neq i e^{μ + α_{i} + β_{j} + Z_{i j}^{⊤} γ}, & i = 1, \dots, n, b_{j} & = & n \sum i = 1, i \neq j e^{μ + α_{i} + β_{j} + Z_{i j}^{⊤} γ}, & j = 1, \dots, n . \end{matrix}

(3)

Under the restriction $μ = 0$ and $β_{n} = 0$ , the first equation and the last equation with $j = n$ in the above system of equations will be excluded. Although there are a total of $2 n + 1 + p$ equations, the number of minimal equations is only $2 n - 1 + p$ .

The MLE $(ˆ μ, ˆ α, ˆ β, ˆ γ)$ of the parameter vector $(μ, α, β, γ)$ is the solution of the above equations if it exist. The Newton-Raphson algorithm can be used to solve the above equations. We can also simply use the function “glm” in R language to calculate the solution.

3 Theoretical Properties

Let $˜ α = μ / 2 + α$ and $˜ β = μ / 2 + β$ , then $˜ α + ˜ β = μ + α + β$ . With this reparameterized technique, we could set $μ = 0$ for convenience. Further, we set $β_{n} = 0$ jointly as the identification condition in this section. The asymptotic properties of $ˆ μ$ in the restriction $α_{n} = β_{n} = 0$ is the same as those of ${ˆ α}_{n}$ in the restriction $μ = β_{n} = 0$ .

Notations. Let $R = (- \infty, \infty)$ be the real domain. For a subset $C \subset R^{n}$ , let $C^{0}$ and $¯ ¯¯ ¯ C$ denote the interior and closure of $C$ , respectively. For a vector $x = (x_{1}, \dots, x_{n})^{⊤} \in R^{n}$ , denote by $∥ x ∥$ for a general norm on vectors with the special cases $∥ x ∥_{\infty} = {max}_{1 \leq i \leq n} | x_{i} |$ and $∥ x ∥_{1} = \sum_{i} | x_{i} |$ for the $ℓ_{\infty}$ - and $ℓ_{1}$ -norm of $x$ respectively. When $n$ is fixed, all norms on vectors are equivalent. Let $B (x, ϵ) = {y : ∥ x - y ∥_{\infty} \leq ϵ}$ be an $ϵ$ -neighborhood of $x$ . For an $n \times n$ matrix $J = (J_{i, j})$ , let $∥ J ∥_{\infty}$ denote the matrix norm induced by the $ℓ_{\infty}$ -norm on vectors in $R^{n}$ , i.e.,

∥ J ∥_{\infty} = max x \neq 0 \frac{∥ J x ∥_{\infty}}{∥ x ∥_{\infty}} = max 1 \leq i \leq n n \sum j = 1 | J_{i, j} |,

and $∥ J ∥$ be a general matrix norm. Define the matrix maximum norm: $∥ J ∥_{max} = {max}_{i, j} | J_{i j} |$ . The notation $\sum_{i}$ denotes the summarization over all $i = 1, \dots, n$ and $\sum_{i \neq j}$ is a shorthand for $\sum_{i = 1}^{n} \sum_{j = 1, j \neq i}^{n}$ . The notation $f (x) ≲ g (x)$ or $f (x) = O (g (x))$ means that there exists a constant $c$ such that $| f (x) | \leq c | g (x) |$ .

For convenience of our theoretical analysis, define $θ = (α_{1}, \dots, α_{n}, β_{1}, \dots, β_{n - 1})^{⊤}$ . We use the superscript “*” to denote the true parameter under which the data are generated. When there is no ambiguity, we omit the super script “*”. Further, define

π_{i j} := α_{i} + β_{j} + Z_{i j}^{⊤} γ, λ (x) := e^{x} .

Write $λ^{'}$ , $λ^{''}$ and $λ^{'''}$ as the first, second and third derivative of $λ (x)$ on $x$ , respectively. Direct calculations give that $λ^{'} (x) = λ^{''} (x) = λ^{'''} (x) = e^{x}$ . Let $ϵ_{n 1}$ and $ϵ_{n 2}$ be two small positive number. When $θ \in B (θ^{*}, ϵ_{n 1})$ , $γ \in B (γ^{*}, ϵ_{n 2})$ , we have

| α_{i} + β_{j} + Z_{i j}^{T} γ | ⩽ max i \neq j ∣ ∣ α_{i}^{*} + β_{j}^{*} ∣ ∣ + 2 ϵ_{n 1} + p q ∥ γ^{*} ∥_{\infty} + p q ϵ_{n 2} := ρ_{n}

where $q := {max}_{i, j} ∥ Z_{i j} ∥_{\infty}$ , and we assume that $q$ is a constant. It is not difficult to verify, When $θ \in B (θ^{*}, ϵ_{n 1})$ and $γ \in B (γ^{*}, ϵ_{n 2})$ ,

e^{- ρ_{n}} ⩽ λ (π_{i j}) = λ^{'} (π_{i j}) = λ^{''} (π_{i j}) = λ^{'''} (π_{i j}) ⩽ e^{ρ_{n}} .

(4)

When causing no confusion, we will simply write $λ_{i j}$ stead of $λ_{i j} (θ, γ)$ for shorthand, where

λ_{i j} (θ, γ) = e^{α_{i} + β_{j} + Z_{i j}^{⊤} γ} = λ (π_{i j}) .

Write the partial derivative of a function vector $F (ˆ θ, γ)$ on $θ$ as

\frac{\partial F (ˆ θ, ˆ γ)}{} \partial θ^{⊤} = \frac{\partial F (θ, γ)}{\partial θ^{⊤}} {∣ ∣ ∣}_{θ = ˆ θ, γ = ˆ γ} .

Throughout the remainder of this paper, we make the following assumption.

Assumption 1.

Assume that $p$ , the dimension of $Z_{i j}$ , is fixed and that the support of $Z_{i j}$ is $Z^{p}$ , where $Z$ is a compact subset of $R$ .

The above assumption is made in Graham (2017) and Yan et al. (2019). In many real applications, the attributes of nodes have a fixed dimension and $Z_{i j}$ is bounded. As one example, if nodal variables are indicator such as sex, then the assumption holds. If $Z_{i j}$ is not bounded, we could make the transform ${~ Z}_{i j} = e^{Z_{i j}} / (1 + e^{Z_{i j}})$ .

3.1 Consistency

In order to establish asymptotic properties of $(ˆ θ, ˆ γ)$ , we define a system of functions

	$F_{i} (θ, γ)$	$= n \sum j = 1, j \neq i λ_{i j} (θ, γ) - d_{i}, i = 1, \dots, n,$
	$F_{n + j} (θ, γ)$	$= n \sum i = 1, i \neq j λ_{i j} (θ, γ) - b_{j}, j = 1, \dots, n,$
	$F (θ, γ)$	$= (F_{1} (θ, γ), \dots, F_{n} (θ, γ), F_{n + 1} (θ, γ), \dots, F_{2 n - 1} (θ, γ))^{⊤},$

which are based on the score equations for $ˆ θ$ . Define $F_{γ, i} (θ)$ be the value of $F_{i} (θ, γ)$ when $γ$ is fixed. Let ${ˆ θ}_{γ}$ be a solution to $F_{γ} (θ) = 0$ if it exists. Correspondingly, we define two functions for exploring the asymptotic behaviors of the estimator of $γ$ :

	$Q (θ, γ) = n \sum i, j = 1; i \neq j Z_{i j} (λ_{i j} (θ, γ) - a_{i j}),$		(5)
	$Q_{c} (γ) = n \sum i, j = 1; i \neq j Z_{i j} (λ_{i j} ({ˆ θ}_{γ}, γ) - a_{i j}) .$		(6)

$Q_{c} (γ)$ could be viewed as the concentrated or profile function of $Q (θ, γ)$ in which the degree parameter $θ$ is profiled out. It is clear that

F (ˆ θ, ˆ γ) = 0, F_{γ} ({ˆ θ}_{γ}) = 0, Q (ˆ θ, ˆ γ) = 0, Q_{c} (ˆ γ) = 0.

We first present the error bound between ${ˆ θ}_{γ}$ with $γ \in B (γ^{*}, ϵ_{n 2})$ and $θ^{*}$ . This is proved by constructing the Newton iterative sequence ${θ^{(k + 1)}}_{k = 0}^{\infty}$ with initial value $θ^{*}$ and obtaining a geometrically fast convergence rate of the iterative sequence, where $θ^{(k + 1)} = θ^{(k)} - [F_{γ}^{'} (θ^{(k)})]^{- 1} F_{γ} (θ^{(k)})$ and $F_{γ}^{'} (θ) = \partial F_{γ} (θ) / \partial γ^{⊤}$ . Details are given in the Supplementary Material. The error bound is stated below.

Lemma 1.

If $γ \in B (γ^{*}, ϵ_{n 2})$ and $e^{ρ_{n}} = o ((n / log n)^{1 / 26})$ , then as $n$ goes to infinity, with probability at least $1 - O (1 / n)$ , ${ˆ θ}_{γ}$ exists and satisfies

∥ {ˆ θ}_{γ} - θ^{*} ∥_{\infty} = O_{p} (e^{7 ρ_{n}} \sqrt{\frac{log n}{n}}) = o_{p} (1) .

Further, if ${ˆ θ}_{γ}$ exists, then it is unique.

By the compound function derivation law, we have

	$0 = \frac{\partial F_{γ} ({ˆ θ}_{γ})}{\partial γ^{⊤}} = \frac{\partial F ({ˆ θ}_{γ}, γ)}{\partial θ^{⊤}} \frac{\partial {ˆ θ}_{γ}}{\partial γ^{⊤}} + \frac{\partial F ({ˆ θ}_{γ}, γ)}{\partial γ^{⊤}},$		(7)
	$\frac{\partial Q_{c} (γ)}{\partial γ^{⊤}} = \frac{\partial Q ({ˆ θ}_{γ}, γ)}{\partial θ^{⊤}} \frac{\partial {ˆ θ}_{γ}}{{\partial γ}^{⊤}} + \frac{\partial Q ({ˆ θ}_{γ}, γ)}{\partial γ^{⊤}} .$		(8)

By solving $\partial {ˆ θ}_{γ} / \partial γ^{⊤}$ in (7) and substituting it into (8), we get the Jacobian matrix $Q_{c}^{'} (γ)$ $(= \partial Q_{c} (γ) / \partial γ^{⊤})$ :

\frac{\partial Q_{c} (γ)}{\partial γ^{⊤}} = \frac{\partial Q ({ˆ θ}_{γ}, γ)}{\partial γ^{⊤}} - \frac{\partial Q ({ˆ θ}_{γ}, γ)}{\partial θ^{⊤}} {[\frac{\partial F ({ˆ θ}_{γ}, γ)}{\partial θ^{⊤}}]}^{- 1} \frac{\partial F ({ˆ θ}_{γ}, γ)}{\partial γ^{⊤}} .

(9)

The asymptotic behavior of $ˆ γ$ crucially depends on the Jacobian matrix $Q_{c}^{'} (γ)$ . Since ${ˆ θ}_{γ}$ does not have a closed form, conditions that are directly imposed on $Q_{c}^{'} (γ)$ are not easily checked. To derive feasible conditions, we define

H (θ, γ) = \frac{\partial Q (θ, γ)}{\partial γ^{⊤}} - \frac{\partial Q (θ, γ)}{\partial θ^{⊤}} {[\frac{\partial F (θ, γ)}{\partial θ^{⊤}}]}^{- 1} \frac{\partial F (θ, γ)}{\partial γ^{⊤}},

(10)

which is a general form of $\partial Q_{c} (γ) / \partial γ$ . Because $H (θ, γ)$ is the Fisher information matrix of the profiled log-likelihood $ℓ_{c} (γ)$ , we assume that it is positively definite. Otherwise, the network Poisson model will be ill-conditioned. When $θ \in B (θ^{*}, ϵ_{n 1})$ , we have the equation:

\frac{1}{n^{2}} H (θ, γ^{*}) = \frac{1}{n^{2}} H (θ^{*}, γ^{*}) + o (1),

(11)

whose proof is omitted. Note that the dimension of $H (θ, γ)$ is fixed and every its entry is a sum of $(n - 1) n$ terms. We assume that there exists a number $κ_{n}$ such that

sup θ \in B (θ^{*}, ϵ_{n 1}) ∥ H^{- 1} (θ, γ^{*}) ∥_{\infty} \leq \frac{κ_{n}}{n^{2}} .

(12)

If $n^{- 2} H (θ, γ^{*})$ converges to a constant matrix, then $κ_{n}$ is bounded. Because $H (θ, γ^{*})$ is positively definite,

κ_{n} = \sqrt{p} \times sup θ \in B (θ^{*}, ϵ_{n 1}) 1 / λ_{min} (θ),

where $λ_{min} (θ)$

is the smallest eigenvalue of

n^{- 2} H (θ, γ^{*})

. Now we formally state the consistency result.

Theorem 1.

If $κ_{n}^{2} e^{40 ρ_{n}} = o (n / log n)$ , then the MLE $(ˆ θ, ˆ γ)$ exists with probability at least $1 - O (1 / n)$ , and is consistent in the sense that

	$∥ ˆ γ - γ^{*} ∥_{\infty}$	$= O_{p} (\frac{κ_{n} e^{21 ρ_{n}} log n}{n}) = o_{p} (1)$
	$∥ ˆ θ - θ^{*} ∥_{\infty}$	$= O_{p} (e^{7 ρ_{n}} \sqrt{\frac{log n}{n}}) = o_{p} (1) .$

Further, if $ˆ θ$ exists, it is unique.

From the above theorem, we can see that $ˆ γ$ has a convergence rate $log n / n$ and $ˆ θ$ has a convergence rate $(log n / n)^{1 / 2}$ , up to an additional factor. If $∥ γ^{*} ∥_{\infty}$ and $∥ θ^{*} ∥_{\infty}$ are constants, then $ρ_{n}$ and $κ_{n}$ are constants such that the condition in Theorem 1 holds.

3.2 Asymptotic normality of $ˆ θ$ and $ˆ γ$

The asymptotic distribution of $ˆ θ$ depends crucially on the inverse of the Jacobian matrix $F_{γ}^{'} (θ)$ , which generally does not have a closed form. In order to characterize this matrix, we introduce a general class of matrices that encompass the Fisher matrix. Given two positive numbers $m$ and $M$ with $M \geq m > 0$ , we say the $(2 n - 1) \times (2 n - 1)$ matrix $V = (v_{i, j})$ belongs to the class $L_{n} (m, M)$ if the following holds:

\begin{matrix} 0 \leq v_{i, i} - \sum_{j = n + 1}^{2 n - 1} v_{i, j} \leq M, i = 1, \dots, 2 n - 1, v_{i, j} = 0, i, j = 1, \dots, n, i \neq j, v_{i, j} = 0, i, j = n + 1, \dots, 2 n - 1, i \neq j, m \leq v_{i, j} = v_{j, i} \leq M, i = 1, \dots, n, j = n + 1, \dots, 2 n, j \neq n + i, \end{matrix}

(13)

Clearly, if $V \in L_{n} (m, M)$ , then $V$ is a $(2 n - 1) \times (2 n - 1)$ diagonally dominant, symmetric nonnegative matrix. It must be positively definite. The definition of $L_{n} (m, M)$ here is wider than that in Yan et al. (2016a), where the diagonal elements are equal to the row sum excluding themselves for some rows. One can easily show that $F_{γ}^{'} (θ)$ belongs to this matrix class. With some abuse of notation, we use $V$ to denote the Fisher information matrix for the vector parameter $θ$ .

To describe the exact form of elements of $V$ , $v_{i j}$ for $i, j = 1, \dots, 2 n - 1$ , $i \neq j$ , we define

u_{i j} = λ (π_{i j}), u_{i i} = 0, u_{i \cdot} = n \sum j = 1 u_{i j}, u_{\cdot j} = n \sum i = 1 u_{i j} .

Note that $u_{i j}$

is the variance of

a_{i j}

. Then the elements of

V

are

v_{i j} = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} u_{i \cdot}, & i = j = 1, \dots, n, u_{i, j - n}, & i = 1, \dots, n, j = n + 1, \dots, 2 n - 1, j \neq i + n, u_{j, i - n}, & i = n + 1, \dots, 2 n, j = 1, \dots, n - 1, j \neq i - n, u_{\cdot j - n}, & i = j = n + 1, \dots, 2 n - 1, 0, & others . \end{matrix}

Yan et al. (2016a) proposed to approximate the inverse $V^{- 1}$ by the matrix $S = (s_{i, j})$ , which is defined as

s_{i, j} = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \frac{δ_{i, j}}{u_{i \cdot}} + \frac{1}{u_{\cdot n}}, & i, j = 1, \dots, n, - \frac{1}{u_{\cdot n}}, & i = 1, \dots, n, j = n + 1, \dots, 2 n - 1, - \frac{1}{u_{\cdot n}}, & i = n + 1, \dots, 2 n, j = 1, \dots, n - 1, \frac{δ_{i, j}}{u_{\cdot (j - n)}} + \frac{1}{u_{\cdot n}}, & i, j = n + 1, \dots, 2 n - 1, \end{matrix}

(14)

where $δ_{i, j} = 1$ when $i = j$ and $δ_{i, j} = 0$ when $i \neq j$ .

We derive the asymptotic normality of $ˆ θ$ by representing $ˆ θ$ as a function of $d$ and $b$ with an explicit expression and a remainder term. This is done via applying a second Taylor expansion to $F (ˆ θ, ˆ γ)$ and showing various remainder terms asymptotically neglect.

Theorem 2.

If $κ_{n} e^{23 ρ_{n}} = o (n^{1 / 2} / log n)$ , then for any fixed $k \geq 1$ , as $n \to \infty$ , the vector consisting of the first $k$ elements of $(ˆ θ - θ^{*})$ is asymptotically multivariate normal with mean $0$ and covariance matrix given by the upper left $k \times k$ block of $S$ defined in (14).

Remark 1.

By Theorem 2, for any fixed $i$ , as $n \to \infty$ , the asymptotic variance of ${^θ}_{i}$ is $1 / v_{i, i}^{1 / 2}$ , whose magnitude is between $O (n^{- 1 / 2} e^{ρ_{n}})$ and $O (n^{- 1 / 2})$ .

Remark 2.

Under the restriction $α_{n} = β_{n} = 0$

, the central limit theorem for the MLE

ˆ μ

is stated as:

(u_{n \cdot} + u_{\cdot n})^{- 1 / 2} (ˆ μ - μ)

, converges in distribution to the standard normality.

Now, we present the asymptotic normality of $ˆ γ$ . Let

V = \frac{\partial F (θ^{*}, γ^{*})}{\partial θ^{⊤}}, V_{γ γ} = \frac{\partial Q (θ^{*}, γ^{*})}{\partial γ^{⊤}}, V_{θ γ} = \frac{\partial F (θ^{*}, γ^{*})}{\partial γ^{⊤}} .

Following Amemiya (1985) (p. 126), the Hessian matrix of the concentrated log-likelihood function $ℓ^{c} (γ^{*})$ is $V_{γ γ} - V_{θ γ}^{⊤} V^{- 1} V_{θ γ}$ . To state the form of the limit distribution of $^γ$ , define

I_{n} (γ^{*}) = \frac{1}{(n - 1) n} (V_{γ γ} - V_{θ γ}^{⊤} V^{- 1} V_{θ γ}) .

(15)

Assume that the limit of $I_{n} (γ^{*})$ exists as $n$ goes to infinity, denoted by $I_{*} (γ^{*})$ . Then, we have the following theorem.

Theorem 3.

If $e^{9 ρ_{n}} = o (n^{1 / 2} / (log n)^{2})$ , then as $n$ goes to infinity, the $p$ -dimensional vector $N^{1 / 2} (ˆ γ - γ^{*})$

is asymptotically multivariate normal distribution with mean

I_{*}^{- 1} (γ^{*}) B_{*}

and covariance matrix

I_{*}^{- 1} (γ)

, where

N = n (n - 1)

and

B_{*}

is the bias term:

B_{*} = lim n \to \infty \frac{1}{2 \sqrt{N}} [n \sum i = 1 \frac{\sum_{j = 1, j \neq i}^{n} λ^{''} (π_{i j}^{*}) Z_{i j}}{\sum_{j = 1, j \neq i}^{n} λ^{'} (π_{i j}^{*})} + n \sum j = 1 \frac{\sum_{i = 1, i \neq j}^{n} λ^{''} (π_{i j}^{*}) Z_{i j}}{\sum_{i = 1, i \neq j}^{n} λ^{'} (π_{i j}^{*})}] .

Remark 3.

The limiting distribution of $ˆ γ$ is involved with a bias term

μ_{*} = \frac{I_{*}^{- 1} (γ) B_{*}}{\sqrt{n (n - 1)}} .

Since the MLE $ˆ γ$

is not centered at the true parameter value, the confidence intervals and the p-values of hypothesis testing for

ˆ γ

may not achieve the nominal level without bias-correction under the null:

γ^{*} = 0

. This is referred to as the so-called incidental parameter problem in econometric literature [Neyman and Scott (1948)]. The produced bias is due to different convergence rates of

ˆ γ

and

ˆ θ

. As discussed in Yan et al. (2019), we could use the analytical bias correction formula:

{ˆ γ}_{b c} = ˆ γ - {^I}^{- 1}^B / \sqrt{n (n - 1)}

, where

^I

and

^B

are the estimates of

I_{*}

and

B_{*}

by replacing

γ

and

θ

in their expressions with their MLEs

ˆ γ

and

ˆ θ

, respectively. In the simulation in next section, we can see that there is a little difference between uncorrected estimates and bias-corrected estimates, which is different from the covariate-

p_{0}

model for binary directed networks in Yan et al. (2019). This may be due to that the infinitely weighted values for edges attenuate the bias effect.

4 Numerical Studies

In this section, we evaluate the asymptotic results of the MLEs for model (1) through simulation studies and a real data example.

4.1 Simulation studies

Similar to Yan et al. (2019), the parameter values take a linear form. Specifically, we set $α_{i}^{*} = i c log n / n$ for $i = 0, \dots, n$ and let $β_{i}^{*} = α_{i}^{*}$ , $i = 0, \dots, n$ for simplicity. Note that there are $n + 1$ nodes in the simulations. We considered four different values for $c$ as $c \in {0, 0.2, 0.4, 0.6}$ . By allowing the true values of $α$ and $β$ to grow with $n$ , we intended to assess the asymptotic properties under different asymptotic regimes. Being different from Yan et al. (2019) with $p = 2$ , we set $p = 3$ here. The first element of $Z_{i j}$ was generated from standard normal distribution; the second element $Z_{i j 2}$ formed by letting $Z_{i j 2} = | X_{i 1} - X_{j 1} |$ , where $X_{i 1}$ is the first entry of the $2$ -dimensional node-specific covariate $X_{i}$ is independently generated from a $B e t a (2, 2)$ distribution; and the third element $Z_{i j 3} = X_{i 2} * X_{j 2}$ , where $X_{i 2}$ follows a discrete distribution taking values $1$ and $- 1$ with probabilities $0.3$ and $0.7$ . The density parameter is $μ = - log n / 4$ .

Note that by Theorems 2, ${^ξ}_{i, j} = [{^α}_{i} - {^α}_{j} - (α_{i}^{*} - α_{j}^{*})] / (1 / {^v}_{i, i} + 1 / {^v}_{j, j})^{1 / 2}$ , ${^ζ}_{i, j} = ({^α}_{i} + {^β}_{j} - α_{i}^{*} - β_{j}^{*}) / (1 / {^v}_{i, i} + 1 / {^v}_{n + j, n + j})^{1 / 2}$ , and ${^η}_{i, j} = [{^β}_{i} - {^β}_{j} - (β_{i}^{*} - β_{j}^{*})] / (1 / {^v}_{n + i, n + i} + 1 / {^v}_{n + j, n + j})^{1 / 2}$ are all asymptotically distributed as standard normal random variables, where ${^v}_{i, i}$ is the estimate of $v_{i, i}$ by replacing $(γ^{*}, θ^{*})$ with $(ˆ γ, ˆ θ)$ . We record the coverage probability of the 95% confidence interval, the length of the confidence interval, and the frequency that the MLE does not exist. The results for ${^ξ}_{i, j}$ , ${^ζ}_{i, j}$ and ${^η}_{i, j}$ are similar, thus only the results of ${^ξ}_{i, j}$ are reported. We simulated networks with $n = 100$ or $n = 200$ . Each simulation is repeated $5, 000$ times.

Table 1 reports the $95 %$ coverage frequencies for $α_{i} - α_{j}$ and the length of the confidence interval. As we can see, the length of the confidence interval increases with $c$ and decreases with $n$ , which qualitatively agrees with the theory. All simulated coverage frequencies are all close to the nominal level. This indicates that the conditions in the theorems may be relaxed greatly.

n	$(i, j)$	$c = 0$	$c = 0.2$	$c = 0.4$	$c = 0.6$
100	$(1, 2)$	$94.36 / 0.529$	$94.72 / 0.408$	$94.54 / 0.305$	$95.1 / 0.222$
	$(50, 51)$	$94.42 / 0.529$	$94.42 / 0.326$	$94.82 / 0.195$	$94.66 / 0.113$
	$(99, 100)$	$94.92 / 0.531$	$94.24 / 0.261$	$94.66 / 0.125$	$93.62 / 0.058$
	$(1, 100)$	$94.06 / 0.530$	$94.24 / 0.339$	$94.62 / 0.227$	$94.68 / 0.157$
	$(1, 50)$	$94.8 / 0.530$	$93.92 / 0.370$	$94.54 / 0.255$	$95.12 / 0.173$
	$(0, 0)^{*}$	$94.68 / 0.039$	$94.64 / 0.024$	$95.06 / 0.013$	$94.94 / 0.007$
200	$(1, 2)$	$94.34 / 0.401$	$94.4 / 0.301$	$94.64 / 0.215$	$95.32 / 0.148$
	$(50, 51)$	$94.34 / 0.402$	$94.64 / 0.233$	$94.98 / 0.126$	$94.56 / 0.068$
	$(99, 100)$	$94.66 / 0.406$	$94.6 / 0.177$	$93.66 / 0.075$	$94.26 / 0.031$
	$(1, 100)$	$94.5 / 0.404$	$94.68 / 0.242$	$94.86 / 0.155$	$95.26 / 0.102$
	$(1, 50)$	$94.04 / 0.402$	$94.16 / 0.267$	$94.58 / 0.173$	$94.98 / 0.112$
	$(0, 0)$	$95.1 / 0.020$	$95.26 / 0.011$	$95.2 / 0.006$	$94.9 / 0.003$

$(0, 0)$ indicates the simulated results for $ˆ μ$ .

Table 1: The reported values are the coverage frequency (

\times 100 %

) for

α_{i} - α_{j}

for a pair

(i, j)

/ the length of the confidence interval. The pair

(0, 0)

indicates those values for the density parameter

μ

Table 2 reports simulation results for the estimate $ˆ γ$ and bias correction estimate ${ˆ γ}_{b c} (= ˆ γ - {^I}^{- 1}^B / \sqrt{n (n - 1)})$ at the nominal level $95 %$

. The coverage frequencies for the uncorrected estimate are all close to those for corrected estimates. All simulated coverage frequencies achieves the nominal level. As expected, the standard error increases with

c

and decreases with

n

$n$	$c$	$γ_{1}$	$γ_{2}$	$γ_{3}$
$100$	$0$	$94.68 (94.66) / 0.039$	$94.7 (94.86) / 0.248$	$95.42 (95.46) / 0.066$
	$0.2$	$94.64 (94.64) / 0.024$	$94.04 (94.14) / 0.151$	$94.56 (94.56) / 0.040$
	$0.4$	$95.06 (95.06) / 0.013$	$94.5 (94.58) / 0.086$	$94.72 (94.74) / 0.023$
	$0.6$	$94.94 (94.94) / 0.007$	$95.44 (95.5) / 0.047$	$94.7 (94.7) / 0.012$
$200$	$0$	$95.1 (95.04) / 0.020$	$94.62 (94.8) / 0.132$	$95.18 (95.2) / 0.036$
	$0.2$	$95.26 (95.26) / 0.011$	$94.76 (94.96) / 0.074$

A network Poisson model for weighted directed networks with covariates

Authors

A sparse p_0 model with covariates for directed networks

Sparse Ising Models with Covariates

Undirected network models with degree heterogeneity and homophily

Incorporating Actor Heterogeneity into Large Network Models through Variational Approximations

Latent Poisson models for networks with heterogeneous density

Weighted directed networks with a differentially private bi-degree sequence

A Sparse Random Graph Model for Sparse Directed Networks

1 Introduction

1.1 Background

1.2 The Model

2 Maximum Likelihood Estimation

3 Theoretical Properties

Assumption 1.

3.1 Consistency

Lemma 1.

Theorem 1.

3.2 Asymptotic normality of $ˆ θ$ and $ˆ γ$

Theorem 2.

Remark 1.

Remark 2.

Theorem 3.

Remark 3.

4 Numerical Studies

4.1 Simulation studies

A network Poisson model for weighted directed networks with covariates

Authors

Related Research

A sparse p_0 model with covariates for directed networks

Sparse Ising Models with Covariates

Undirected network models with degree heterogeneity and homophily

Incorporating Actor Heterogeneity into Large Network Models through Variational Approximations

Latent Poisson models for networks with heterogeneous density

Weighted directed networks with a differentially private bi-degree sequence

A Sparse Random Graph Model for Sparse Directed Networks

This week in AI

1 Introduction

1.1 Background

1.2 The Model

2 Maximum Likelihood Estimation

3 Theoretical Properties

Assumption 1.

3.1 Consistency

Lemma 1.

Theorem 1.

3.2 Asymptotic normality of ˆθ and ˆγ

Theorem 2.

Remark 1.

Remark 2.

Theorem 3.

Remark 3.

4 Numerical Studies

4.1 Simulation studies

3.2 Asymptotic normality of $ˆ θ$ and $ˆ γ$