A sparse $p_{0}$ model with covariates for directed networks

There is a large body of work concerned with the degree heterogeneity. The $\beta$-model is one of the most popular models to model the degree heterogeneity, which is named by Chatterjee et al. (2011). Since the number of parameters increases with the number of nodes, asymptotic theory is nonstandard. Exploring the properties of the $\beta$-model and its generalizations has received wide attentions in recent years (Chatterjee et al., 2011; Perry and Wolfe, 2012; Hillar and Wibisono, 2013; Yan and Xu, 2013; Rinaldo and Fienberg, 2013). In particular, Chatterjee et al. (2011) proved the uniform consistency of the maximum likelihood estimator (MLE) and Yan and Xu (2013) derived the asymptotic normality of the MLE. Yan et al. (2016b) further established a unified theoretical framework under which the consistency and asymptotic normality of the moment estimator hold in a class of node-parameter network models. In the directed case, Yan et al. (2016a) proved the consistency and asymptotic normality of the MLE in the $p_{0}$ model, which is a special case of the $p_{1}$ model by Holland and Leinhardt (1981).

Graham (2017) incorporated covariates to the $\beta$-model to model the homophily in undirected networks. Dzemski (2019) and Yan et al. (2019) proposed directed graph models to model edges with covariates by using a bivariate normal distribution and a logistic distribution, respectively. The asymptotic normality of the restricted maximum likelihood estimator of the covariate parameter were derived under the assumption that the estimators for all parameters are taken in one compact set in these papers. Yan et al. (2019) further derived the asymptotic normality of the restricted MLE of the node parameter. However, the convergence rate of the covariate parameter is not investigated in their works. Yan (2021) worked with the unrestricted maximum likelihood estimation in a sparse $\beta$-model with covariates for undirected graphs and proved the uniformly consistency and asymptotic normality of the MLE, where the consistency is proved via the convergence rate of the Newton iterative sequence in a two-stage Newton method. We demonstrate that his method could be further extended here. We also note that extending his result to directed graphs is not trivial because of that the Fisher information matrix here is very different from that in the sparse $\beta$-model with covariates and that the dimension of parameter space here is much bigger.

Finally, we note that the model can be recast into a logistic regression model. There have been a considerable mount of literatures on asymptotic theory in generalized linear models with a growing number of parameters. Haberman (1977) establish the asymptotic theory of the MLE in exponential response models in a setting that allows the number of parameters goes to infinity. By assuming that a sequence of independent and identical distributed samples $\{X_{i}\}_{i=1}^{n}$ from a canonical exponential family $c(\theta)\exp(\theta^{\top}x)$ with an increasing dimension $p_{n}$, Portnoy (1988) shew $\|\widehat{\theta}-\theta^{*}\|_{\infty}=O_{p}(\sqrt{p_{n}/n})$ if $p_{n}/n=o(1)$ and $\widehat{\theta}$ is asymptotically normal if $p_{n}^{2}/n=o(1)$ under some additional conditions. For clustered binary data $\{(Y_{i},X_{i})\}_{i=1}^{n}$ with a $p_{n}$ dimensional covariate $X_{i}$ and regression parameter $\theta$, Wang (2011) derived the asymptotic theory of generalized estimating equations estimators if $p_{n}^{2}/n=o(1)$. Liang and Du (2012) establish the asymptotic normality of the maximum likelihood estimate when the number of covariates $p$ goes to infinity with the sample size $n$ in the order of $p=o(n)$, but the error bound between $\widehat{\theta}$ and $\theta$ is not investigated. We note that asymptotic theory in these papers needs the condition that $c_{1}N\leq\lambda_{\min}\leq\lambda_{\max}\leq c_{2}N$ for two positive constants $c_{1}$ and $c_{2}$, where $S_{N}=\sum_{i=1}^{N}x_{i}x_{i}^{\top}$ and $\lambda_{\min}$ and $\lambda_{\max}$ are the minimum and maximum eigenvalue of $S_{N}$, $x_{i}$ is a $p_{N}$-dimensional covariate vector. In our model, the first $n$ diagonal entries of $S_{N}$ is in the order of $n$ while the last $p$ diagonal entries of $S_{N}$ is in the order of $n^{2}$. Because of different orders of the diagonal elements of $S_{N}$, the ratio $\lambda_{\max}/\lambda_{\min}$ is in the order of $O(n^{2})$, far away from the constant order in Liang and Du (2012). It is also worth noting that the consistency and asymptotic normality of the MLE have been derived in exponential response model [Portnoy (1988)]. The asymptotic normality of $\widehat{\theta}$ does not contain a bias term in all mentioned these papers while the asymptotic distribution of $\widehat{\gamma}$ in our model contains a bias term. This is referred to as the incident bias problem in economic literature [e.g., Graham (2017)].

For the remainder of the paper, we proceed as follows. In Section 2, we give the details on the model considered in this paper. In section 3, we establish asymptotic results. Numerical studies are presented in Section 4. We provide further discussion and future work in Section 5. All proofs are relegated to the appendix.

In order to establish the existence and consistency of $(\widehat{\boldsymbol{\eta}},\widehat{\boldsymbol{\gamma}})$, we define a system of functions

	$\displaystyle F_{i}(\boldsymbol{\eta},\boldsymbol{\gamma})$	$\displaystyle=$	$\displaystyle\sum_{k=1;k\neq i}^{n}\frac{e^{\alpha_{i}+\beta_{k}+Z_{ij}^{\top}\boldsymbol{\gamma}}}{1+e^{\alpha_{i}+\beta_{k}+Z_{ij}^{\top}\boldsymbol{\gamma}}}-d_{i},~{}~{}~{}i=1,\ldots,n,$
	$\displaystyle F_{n+j}(\boldsymbol{\eta},\boldsymbol{\gamma})$	$\displaystyle=$	$\displaystyle\sum_{k=1;k\neq j}^{n}\frac{e^{\alpha_{k}+\beta_{j}+Z_{kj}^{\top}\boldsymbol{\gamma}}}{1+e^{\alpha_{k}+\beta_{j}++Z_{kj}^{\top}\boldsymbol{\gamma}}}-b_{j},~{}~{}~{}j=1,\ldots,n,$		(6)
	$\displaystyle F(\boldsymbol{\eta},\boldsymbol{\gamma})$	$\displaystyle=$	$\displaystyle(F_{1}(\boldsymbol{\eta},\boldsymbol{\gamma}),\ldots,F_{n}(\boldsymbol{\eta},\boldsymbol{\gamma}),F_{n+1}(\boldsymbol{\eta},\boldsymbol{\gamma}),\ldots,F_{2n-1}(\boldsymbol{\eta},\boldsymbol{\gamma}))^{\top},$

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

	$\displaystyle Q(\boldsymbol{\eta},\boldsymbol{\gamma})=\sum_{i,j=1;i\neq j}^{n}Z_{ij}(\mu_{ij}(\boldsymbol{\eta},\boldsymbol{\gamma})-a_{ij}),$		(7)
	$\displaystyle Q_{c}(\boldsymbol{\gamma})=\sum_{i,j=1;i\neq j}^{n}Z_{ij}(\mu_{ij}(\widehat{\boldsymbol{\eta}}_{\gamma},\boldsymbol{\gamma})-a_{ij}).$		(8)

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

$$F(\widehat{\boldsymbol{\eta}},\widehat{\boldsymbol{\gamma}})=0,~{}~{}F_{\gamma}(\widehat{\boldsymbol{\eta}}_{\gamma})=0,~{}~{}Q(\widehat{\boldsymbol{\eta}},\widehat{\boldsymbol{\gamma}})=0,~{}~{}Q_{c}(\widehat{\boldsymbol{\gamma}})=0.$$

We first present the error bound between $\widehat{\boldsymbol{\eta}}_{\gamma}$ with $\gamma\in B(\boldsymbol{\gamma}^{*},\epsilon_{n2})$ and $\boldsymbol{\eta}^{*}$. We construct the Newton iterative sequence $\{\boldsymbol{\eta}^{(k+1)}\}_{k=0}^{\infty}$ with initial value $\boldsymbol{\eta}^{*}$, where $\boldsymbol{\eta}^{(k+1)}=\boldsymbol{\eta}^{(k)}-[F_{\gamma}^{\prime}(\boldsymbol{\eta}^{(k)})]^{-1}F_{\gamma}(\boldsymbol{\eta}^{(k)})$ and $F_{\gamma}^{\prime}(\boldsymbol{\eta})=\partial F_{\gamma}(\boldsymbol{\eta})/\partial\boldsymbol{\gamma}^{\top}$. If the iterative converges, then the solution lies in the neighborhood of $\boldsymbol{\eta}^{*}$. This is done by establishing a geometrically fast convergence rate of the iterative sequence. Details are given in supplementary material. The error bound is stated below.

By the compound function derivation law, we have

	$\displaystyle 0=\frac{\partial F_{\gamma}(\widehat{\boldsymbol{\eta}}_{\gamma})}{\partial\boldsymbol{\gamma}^{\top}}=\frac{\partial F(\widehat{\boldsymbol{\eta}}_{\gamma},\gamma)}{\partial\boldsymbol{\eta}^{\top}}\frac{\partial\widehat{\boldsymbol{\eta}}_{\gamma}}{\boldsymbol{\gamma}^{\top}}+\frac{\partial F(\widehat{\boldsymbol{\eta}}_{\gamma},\boldsymbol{\gamma})}{\partial\boldsymbol{\gamma}^{\top}},$		(9)
	$\displaystyle\frac{\partial Q_{c}({\boldsymbol{\gamma}})}{\partial{\boldsymbol{\gamma}}^{\top}}=\frac{\partial Q(\widehat{\boldsymbol{\eta}}_{\gamma},{\boldsymbol{\gamma}})}{\partial\boldsymbol{\eta}^{\top}}\frac{\partial\widehat{\boldsymbol{\eta}}_{\gamma}}{{\boldsymbol{\gamma}}^{\top}}+\frac{\partial Q(\widehat{\boldsymbol{\eta}}_{\gamma},{\boldsymbol{\gamma}})}{\partial{\boldsymbol{\gamma}}^{\top}}.$		(10)

By solving $\partial\widehat{\boldsymbol{\eta}}_{\gamma}/\partial{\boldsymbol{\gamma}}^{\top}$ in (9) and substituting it into (10), we get the Jacobian matrix $Q_{c}^{\prime}(\boldsymbol{\gamma})$ $(=\partial Q_{c}(\boldsymbol{\gamma})/\partial\boldsymbol{\gamma}^{\top})$:

$\displaystyle\frac{\partial Q_{c}(\boldsymbol{\gamma})}{\partial\boldsymbol{\gamma}^{\top}}=\frac{\partial Q(\widehat{\boldsymbol{\eta}}_{\gamma},\boldsymbol{\gamma})}{\partial\boldsymbol{\gamma}^{\top}}-\frac{\partial Q(\widehat{\boldsymbol{\eta}}_{\gamma},\boldsymbol{\gamma})}{\partial\boldsymbol{\eta}^{\top}}\left[\frac{\partial F(\widehat{\boldsymbol{\eta}}_{\gamma},\boldsymbol{\gamma})}{\partial\boldsymbol{\eta}^{\top}}\right]^{-1}\frac{\partial F(\widehat{\boldsymbol{\eta}}_{\gamma},\boldsymbol{\gamma})}{\partial\boldsymbol{\gamma}^{\top}}.$

(11)

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

$$H(\boldsymbol{\eta},\boldsymbol{\gamma})=\frac{\partial Q(\boldsymbol{\eta},\boldsymbol{\gamma})}{\partial\boldsymbol{\gamma}}-\frac{\partial Q(\boldsymbol{\eta},\boldsymbol{\gamma})}{\partial\boldsymbol{\eta}}\left[\frac{\partial F(\boldsymbol{\eta},\boldsymbol{\gamma})}{\partial\boldsymbol{\eta}}\right]^{-1}\frac{\partial F(\boldsymbol{\eta},\boldsymbol{\gamma})}{\partial\boldsymbol{\gamma}},$$

(12)

which is a general form of $\partial Q_{c}(\boldsymbol{\gamma})/\partial\boldsymbol{\gamma}$. Because $H(\boldsymbol{\eta},\boldsymbol{\gamma})$ is the Fisher information matrix of the profiled log-likelihood $\ell_{c}(\boldsymbol{\gamma})$, we assume that it is positively definite. Otherwise, the covariate-$p_{0}$-model will be ill-conditioned. When $\boldsymbol{\eta}\in B(\boldsymbol{\eta}^{*},\epsilon_{n1})$, we have the equation:

$$\frac{1}{n^{2}}H(\boldsymbol{\eta},\boldsymbol{\gamma}^{*})=\frac{1}{n^{2}}H(\boldsymbol{\eta}^{*},\boldsymbol{\gamma}^{*})+o(1),$$

(13)

whose proof is omitted. Note that the dimension of $H(\boldsymbol{\eta},\boldsymbol{\gamma})$ is fixed and every its entry is a sum of $(n-1)n$ terms. We assume that there exists a number $\kappa_{n}$ such that

$$\sup_{\boldsymbol{\eta}\in B(\boldsymbol{\eta}^{*},\epsilon_{n1})}\|H^{-1}(\boldsymbol{\eta},\boldsymbol{\gamma}^{*})\|_{\infty}\leq\frac{\kappa_{n}}{n^{2}}.$$

(14)

If $n^{-2}H(\boldsymbol{\eta},\boldsymbol{\gamma}^{*})$ converges to a constant matrix, then $\kappa_{n}$ is bounded. Because $H(\boldsymbol{\eta},\boldsymbol{\gamma}^{*})$ is positively definite,

$$\kappa_{n}=\sqrt{p}\times\sup_{\boldsymbol{\eta}\in B(\boldsymbol{\eta}^{*},\epsilon_{n1})}1/\lambda_{\min}(\boldsymbol{\eta}),$$

where $\lambda_{\min}(\boldsymbol{\eta})$ be the smallest eigenvalue of $n^{-2}H(\boldsymbol{\eta},\boldsymbol{\gamma}^{*})$. Now we formally state the consistency result.

From the above theorem, we can see that $\widehat{\boldsymbol{\gamma}}$ has a convergence rate $1/n$ and $\widehat{\boldsymbol{\eta}}$ has a convergence rate $1/n^{1/2}$, up to a logarithm factor, if $b_{n}$ and $\kappa_{n}$ is bounded above by a constant. If $\|\boldsymbol{\gamma}^{*}\|_{\infty}$ and $\|\boldsymbol{\eta}^{*}\|_{\infty}$ are constants, then $b_{n}$ and $\kappa_{n}$ are constants as well such that the condition in Theorem 1 holds.

As an application of the $\ell_{\infty}$-error bound, we use it to recover nonzero signals on node parameters. We need a separation measure $\Delta_{n}$ to distinguish between nonzero signals and zero signals for $\beta_{i}$s. Define $\Xi=\{i:\beta_{i}^{*}>\Delta_{n}\}$ as the set of signals. Because

$$\hat{\beta}_{i}\geq\beta_{i}^{*}-|\hat{\beta}_{i}-\beta_{i}^{*}|=\Delta_{n}-O\left(b_{n}^{3}\sqrt{\frac{\log n}{n}}\right),~{}i\in\Xi,$$

we immediately have the following corollary.

The asymptotic distribution of $\widehat{\boldsymbol{\eta}}$ depends crucially on the inverse of the Fisher information matrix $F^{\prime}_{\gamma}(\boldsymbol{\eta})$. It 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 $(2n-1)\times(2n-1)$ matrix $V=(v_{i,j})$ belongs to the class $\mathcal{L}_{n}(m,M)$ if the following holds:

$$\begin{array}[]{l}0\leq v_{i,i}-\sum_{j=1}^{2n-1}v_{i,j}\leq M,~{}~{}i=1,\ldots,2n-1,\\ v_{i,j}=0,~{}~{}i,j=1,\ldots,n,~{}i\neq j,\\ v_{i,j}=0,~{}~{}i,j=n+1,\ldots,2n-1,~{}i\neq j,\\ m\leq v_{i,j}=v_{j,i}\leq M,~{}~{}i=1,\ldots,n,~{}j=n+1,\ldots,2n,~{}j\neq n+i,\\ \end{array}$$

(15)

Clearly, if $V\in\mathcal{L}_{n}(m,M)$, then $V$ is a $(2n-1)\times(2n-1)$ diagonally dominant, symmetric nonnegative matrix. It must be positively definite. The definition of $\mathcal{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 the Fisher information matrix for the vector parameter $\boldsymbol{\eta}$ belongs to $\mathcal{L}_{n}(m,M)$. With some abuse of notation, we use $V$ to denote the Fisher information matrix for the vector parameter $\boldsymbol{\eta}$.

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

$$u_{ij}=\mu^{\prime}(\pi_{ij}),~{}~{}u_{ii}=0,~{}~{}u_{i\cdot}=\sum_{j=1}^{n}u_{ij},~{}~{}u_{\cdot j}=\sum_{i=1}^{n}u_{ij}.$$

Note that $u_{ij}$ is the variance of $a_{ij}$. Then the elements of $V$ are

$$v_{ij}=\begin{cases}u_{i\cdot},&i=j=1,\ldots,n,\\ u_{i,j-n},&i=1,\ldots,n,j=n+1,\ldots,2n-1,j\neq i+n,\\ u_{i-n,j},&i=n+1,\ldots,2n,j=1,\ldots,n-1,j\neq i-n,\\ u_{\cdot j-n},&i=j=n+1,\ldots,2n-1,\\ 0,&\mbox{others}.\end{cases}$$

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}=\left\{\begin{array}[]{ll}\frac{\delta_{i,j}}{u_{i\cdot}}+\frac{1}{u_{\cdot n}},&i,j=1,\ldots,n,\\ -\frac{1}{u_{\cdot n}},&i=1,\ldots,n,~{}~{}j=n+1,\ldots,2n-1,\\ -\frac{1}{u_{\cdot n}},&i=n+1,\ldots,2n,~{}~{}j=1,\ldots,n-1,\\ \frac{\delta_{i,j}}{u_{\cdot(j-n)}}+\frac{1}{u_{\cdot n}},&i,j=n+1,\ldots,2n-1,\end{array}\right.$$

(16)

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

We derive the asymptotic normality of $\widehat{\boldsymbol{\boldsymbol{\eta}}}$ by approximately representing $\boldsymbol{\widehat{\boldsymbol{\eta}}}$ as a function of $\mathbf{d}$ and $\mathbf{b}$ with an explicit expression. This is done via applying a second Taylor expansion to $F(\widehat{\boldsymbol{\eta}},\widehat{\boldsymbol{\gamma}})$ and showing various remainder terms asymptotically neglect. Details are given in the supplementary mateiral.

With similar arguments for proving the asymptotic normality of the restricted MLE for ${\boldsymbol{\gamma}}$ in Yan et al. (2019), we have the same asymptotic normality of the unrestricted MLE. We only present the result here and the proof is omitted. Let

$$V_{\eta\gamma}=\frac{\partial F(\boldsymbol{\eta}^{*},\boldsymbol{\gamma}^{*})}{\partial\boldsymbol{\eta}^{\top}},~{}~{}V_{\gamma\gamma}=\frac{\partial Q(\boldsymbol{\eta}^{*},\boldsymbol{\gamma}^{*})}{\partial\boldsymbol{\gamma}^{\top}}.$$

Following Amemiya (1985) (p. 126), the Hessian matrix of the concentrated log-likelihood function $\ell^{c}(\boldsymbol{\gamma}^{*})$ is $V_{\gamma\gamma}-V_{\eta\gamma}^{\top}V^{-1}V_{\eta\gamma}$. To state the form of the limit distribution of $\hat{\boldsymbol{\gamma}}$, define

$$I_{n}(\boldsymbol{\gamma}^{*})=\frac{1}{(n-1)n}(V_{\gamma\gamma}-V_{\eta\gamma}^{\top}V^{-1}V_{\eta\gamma}).$$

(17)

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

The limiting distribution of $\boldsymbol{\widehat{\gamma}}$ is involved with a bias term

$$\mu_{*}=\frac{I_{*}^{-1}(\boldsymbol{\gamma})B_{*}}{\sqrt{n(n-1)}}.$$

If all parameters $\boldsymbol{\gamma}$ and $\boldsymbol{\theta}$ are bounded, then $\mu_{*}=O(n^{-1/2})$. It follows that $B_{*}=O(1)$ and $(I_{*})_{i,j}=O(1)$ according to their expressions. Since the MLE $\boldsymbol{\widehat{\gamma}}$ is not centered at the true parameter value, the confidence intervals and the p-values of hypothesis testing constructed from $\boldsymbol{\widehat{\gamma}}$ cannot achieve the nominal level without bias-correction under the null: $\boldsymbol{\gamma}^{*}=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 the appearance of additional parameters. As discussed in Yan et al. (2019), we could use the analytical bias correction formula: $\boldsymbol{\widehat{\gamma}}_{bc}=\boldsymbol{\widehat{\gamma}}-\hat{I}^{-1}\hat{B}/\sqrt{n(n-1)}$, where $\hat{I}$ and $\hat{B}$ are the estimates of $I_{*}$ and $B_{*}$ by replacing $\boldsymbol{\gamma}$ and $\boldsymbol{\theta}$ in their expressions with their MLEs $\boldsymbol{\widehat{\gamma}}$ and $\boldsymbol{\widehat{\theta}}$, respectively.