Affiliation networks with an increasing degree sequence

TL;DR

This paper investigates the asymptotic properties of estimators in affiliation networks with increasing size, establishing their consistency and normality, supported by simulations and real data.

Contribution

It introduces a new affiliation model based on degree sequences and proves the asymptotic behavior of its maximum likelihood estimator.

Findings

MLE is uniformly consistent and asymptotically normal

Simulation studies confirm theoretical results

Real data example supports model validity

Abstract

Affiliation network is one kind of two-mode social network with two different sets of nodes (namely, a set of actors and a set of social events) and edges representing the affiliation of the actors with the social events. Although a number of statistical models are proposed to analyze affiliation networks, the asymptotic behaviors of the estimator are still unknown or have not been properly explored. In this paper, we study an affiliation model with the degree sequence as the exclusively natural sufficient statistic in the exponential family distributions. We establish the uniform consistency and asymptotic normality of the maximum likelihood estimator when the numbers of actors and events both go to infinity. Simulation studies and a real data example demonstrate our theoretical results.

Tables

Table 1. Table 1: Estimated coverage probabilities of α i − α j subscript 𝛼 𝑖 subscript 𝛼 𝑗 \alpha_{i}-\alpha_{j} for pair ( i , j ) 𝑖 𝑗 (i,j) as well as the length of confidence intervals (in square brackets), and the probabilities that the MLE does not exist (in parentheses), multiplied by 100 100 100 .

m	$(i,j)$	$L=0$	$L=\log m$	$L={\left(\log m\right)}^{1/2}$	$L=\log \left(m\right)$
100	$(1,2)$	$94.93\left[0.40\right]\left(0\right)$	$95.17\left[0.55\right]\left(0\right)$	$94.68\left[0.70\right]\left(0.43\right)$	$\left(100\right)$
	$(50,51)$	$95.15\left[0.40\right]\left(0\right)$	$95.02\left[0.46\right]\left(0\right)$	$94.78\left[0.52\right]\left(0.43\right)$	$\left(100\right)$
	$(99,100)$	$95.66\left[0.40\right]\left(0\right)$	$95.17\left[0.42\right]\left(0\right)$	$94.75\left[0.44\right]\left(0.43\right)$	$\left(100\right)$

Table 2. Table 2: Estimated coverage probabilities of β i − β j subscript 𝛽 𝑖 subscript 𝛽 𝑗 \beta_{i}-\beta_{j} for pair ( i , j ) 𝑖 𝑗 (i,j) as well as the length of confidence intervals (in square brackets), and the probabilities that the MLE does not exist (in parentheses), multiplied by 100 100 100 .

n	$(i,j)$	$L=0$	$L=\log n$	$L={\left(\log n\right)}^{1/2}$	$L=\log \left(n\right)$
200	$(1,2)$	$94.77\left[0.60\right]\left(0\right)$	$94.13\left[0.83\right]\left(0\right)$	$94.75\left[1.17\right]\left(2.41\right)$	$\left(100\right)$
	$(100,101)$	$95.41\left[0.60\right]\left(0\right)$	$94.79\left[0.68\right]\left(0\right)$	$94.87\left[0.77\right]\left(2.41\right)$	$\left(100\right)$
	$(198,199)$	$96.60\left[0.60\right]\left(0\right)$	$95.04\left[0.60\right]\left(0\right)$	$94.86\left[0.63\right]\left(2.41\right)$	$\left(100\right)$

Table 3. Table 3: The Student Extracurricular network data: the estimated influence parameters 𝜽 ^ ^ 𝜽 \hat{\boldsymbol{\theta}} , 95 % percent 95 95\% confidence intervals (in square brackets) and their standard errors (in parentheses)

Student Organizations	Degree	${\hat{\alpha }}_i$
Spanish.Club	$199$	$-0.32[-0.49,-0.16]\left(0.08\right)$
Pep.Club	$157$	$-0.64[-0.82,-0.46]\left(0.09\right)$
NHS	$124$	$-0.94[-1.14,-0.75]\left(0.10\right)$
Latin.Club	$93$	$-1.29[-1.51,-1.07]\left(0.11\right)$
Orchestra..Symphonic	$89$	$-1.34[-1.57,-1.12]\left(0.11\right)$
Key.Club	$76$	$-1.52[-1.76,-1.28]\left(0.12\right)$
Spanish.Club..high.	$68$	$-1.65[-1.90,-1.40]\left(0.13\right)$
Drunk.Driving	$67$	$-1.67[-1.92,-1.41]\left(0.13\right)$
Forensics..National.Forensics.League.	$66$	$-1.68[-1.94,-1.43]\left(0.13\right)$
Choir..a.capella	$65$	$-1.70[-1.96,-1.44]\left(0.13\right)$
$⋮$	$⋮$	$⋮$
Chess.Club	$7$	$-4.04[-4.78,-3.29]\left(0.38\right)$
Volleyball..JV	$7$	$-4.04[-4.78,-3.29]\left(0.38\right)$
Teachers.of.Tomorrow	$5$	$-4.38[-5.26,-3.50]\left(0.45\right)$
Quiz.Bowl..all.	$5$	$-4.38[-5.26,-3.50]\left(0.45\right)$
Cheerleaders..Spirit.Squad	$5$	$-4.38[-5.26,-3.50]\left(0.45\right)$
Drunk.Driving.Officers	$4$	$-4.60[-5.58,-3.62]\left(0.50\right)$
Choir..vocal.ensemble..4.women.	$4$	$-4.60[-5.58,-3.62]\left(0.50\right)$
Choir..barbershop.quartet..4.men.	$4$	$-4.60[-5.58,-3.62]\left(0.50\right)$
Cross.Country..girls.8th	$3$	$-4.89[-6.02,-3.76]\left(0.58\right)$
Swim…Dive.Team..boys	$3$	$-4.89[-6.02,-3.76]\left(0.58\right)$

Table 4. Table 4: The Student Extracurricular network data: the estimated influence parameters 𝜽 ^ ^ 𝜽 \hat{\boldsymbol{\theta}} , 95 % percent 95 95\% confidence intervals (in square brackets) and their standard errors (in parentheses)

Student ID	Degree	${\hat{\beta }}_j$
$122662$	$14$	$0.996[0.387,1.606]\left(0.311\right)$
$114850$	$12$	$0.792[0.146,1.438]\left(0.330\right)$
$888947$	$12$	$0.792[0.146,1.438]\left(0.330\right)$
$126259$	$11$	$0.680[0.011,1.348]\left(0.341\right)$
$139161$	$11$	$0.680[0.011,1.348]\left(0.341\right)$
$122638$	$10$	$0.559[-0.136,1.254]\left(0.355\right)$
$888981$	$10$	$0.559[-0.136,1.254]\left(0.355\right)$
$888988$	$10$	$0.559[-0.136,1.254]\left(0.355\right)$
$889059$	$10$	$0.559[-0.136,1.254]\left(0.355\right)$
$114037$	$9$	$0.428[-0.298,1.153]\left(0.370\right)$
$114671$	$9$	$0.428[-0.298,1.153]\left(0.370\right)$
$888892$	$9$	$0.428[-0.298,1.153]\left(0.370\right)$
$⋮$	$⋮$	$⋮$
$889102$	$1$	$-1.985[-3.974,0.004]\left(1.015\right)$
$889103$	$1$	$0.000[-0.844,0.844]\left(0.431\right)$
The full table is available by sending emails to zhang_yong@mails.ccnu.edu.cn

Equations

$$P(\mathcal{G}(m,n))=\exp ({\mathbf{\alpha }}^{\top }\mathbf{d}+{\mathbf{\beta }}^{\top }\mathbf{b}-\mathrm{Z}(\mathbf{\alpha },\mathbf{\beta })),$$

$$\exp ({\mathbf{\alpha }}^{\top }\mathbf{d}+{\mathbf{\beta }}^{\top }\mathbf{b})=\exp \left(\sum _{i=1}^m\sum _{j=1}^n({\alpha }_i+{\beta }_j)x_{i,j}\right)=\prod _{i=1}^m\prod _{j=1}^n\exp (({\alpha }_i+{\beta }_j)x_{i,j}),$$

$$\mathrm{Z}(\mathbf{\alpha },\mathbf{\beta })=\sum _{i=1}^m\sum _{j=1}^n\log (1+\exp ({\alpha }_i+{\beta }_j)).$$

$$\mathbb{P}(x_{ij}=1)=\frac{e^{{\alpha }_i+{\beta }_j}}{1+e^{{\alpha }_i+{\beta }_j}},$$

$$\ell (\mathbf{\theta })=\sum _{i=1}^m{\alpha }_i{d}_i+\sum _{j=1}^{n-1}{\beta }_j{b}_j-\sum _{i=1}^m\sum _{j=1}^n\log (1+e^{{\alpha }_i+{\beta }_j}).$$

$$\begin{array}[]{rcl}d_{i}&=&\sum_{j=1}^{n}\frac{e^{\alpha_{i}+\beta_{j}}}{1+e^{\alpha_{i}+\beta_{j}}},~{}~{}~{}i=1,\ldots,m,\\ b_{j}&=&\sum_{i=1}^{m}\frac{e^{\alpha_{i}+\beta_{j}}}{1+e^{\alpha_{i}+\beta_{j}}},~{}~{}~{}j=1,\ldots,n-1.\end{array}$$

$$v_{i,i}=\sum _{j=1}^n\frac{e^{{\alpha }_i+{\beta }_j}}{(1+e^{{\alpha }_i+{\beta }_j}{)}^2},i=1,\dots ,m,v_{m+j,m+j}=\sum _{i=1}^m\frac{e^{{\alpha }_i+{\beta }_j}}{(1+e^{{\alpha }_i+{\beta }_j}{)}^2},j=1,\dots ,n.$$

$$s_{i,j}=\left\{\begin{array}[]{ll}\frac{\delta_{i,j}}{v_{i,i}}+\frac{1}{v_{m+n,m+n}},&{i,j=1,\dots,m,}\\ -\frac{1}{v_{m+n,m+n}},&{i=1,\dots,m,~{}j=m+1,\dots,m+n-1,}\\ -\frac{1}{v_{m+n,m+n}},&{i=m+1,\dots,m+n-1,j=1,\dots,m,}\\ \frac{\delta_{i,j}}{v_{i,i}}+\frac{1}{v_{m+n,m+n}},&{i,j=m+1,\dots,m+n-1,}\end{array}\right.$$

$$\Vert \hat{\mathbf{\theta }}-{\mathbf{\theta }}^{\ast }{\Vert }_{\infty }=O_p\left(\frac{(\log n{)}^{1/2}{e}^{6\parallel {\mathbf{\theta }}^{\ast }{\parallel }_{\infty }}}{n^{1/2}}\right)=o_p(1).$$

$$\Vert J{\Vert }_{\infty }=\underset{\mathbf{x}\ne 0}{\max }\frac{\Vert J\mathbf{x}{\Vert }_{\infty }}{\Vert \mathbf{x}{\Vert }_{\infty }}=\underset{1\le i\le n}{\max }\sum _{j=1}^n|J_{i,j}|.$$

$$\Vert F(\mathbf{x})(\mathbf{v})-F(\mathbf{y})(\mathbf{v}){\Vert }_{\infty }\le \lambda \Vert \mathbf{x}-\mathbf{y}{\Vert }_{\infty }\Vert \mathbf{v}{\Vert }_{\infty },$$

$$\begin{array}[]{l}q\leq v_{i,i}-\sum_{j=m+1}^{m+n-1}v_{i,j}\leq Q,i=1,\dots,m;v_{m,m}=\sum_{j=m+1}^{m+n-1}v_{m,j},\\ v_{i,j}=0,i,j=1,\dots,m,i\neq j,\\ v_{i,j}=0,i,j=m+1,\dots,m+n-1,i\neq j,\\ q\leq v_{i,j}=v_{j,i}\leq Q,i=1,\dots,m,j=m+1,\dots,m+n-1,\\ v_{i,i}=\sum_{k=1}^{m}v_{k,i}=\sum_{k=1}^{m}v_{i,k},i=m+1,\dots,m+n-1.\end{array}$$

$$\Vert V^{-1}-S\Vert \le \frac{c_1{Q}^2}{q^{3}mn},$$

$$v_{i,i}$$

$$q\le v_{i,j}\le Q.$$

$$F=T^{-1}-S=(T^{-1}-S)(I-TS)+S(I-TS)=FU+W.$$

$$u_{i,j}$$

$$u_{i,j}=\left\{\begin{array}[]{ll}(\delta_{i,j}-1)\frac{v_{i,j}}{v_{j,j}}+\frac{v_{i,m+n}}{v_{m+n,m+n}}&i=1,\ldots,m;j=m+1,\ldots,m+n-1,\\ (\delta_{i,j}-1)\frac{v_{i,j}}{v_{j,j}}&i=m+1,\ldots,m+n-1;j=1,\ldots,m,m+1,\ldots,m+n-1.\end{array}\right.$$

$$u_{i,j}=({\delta }_{i,j}-1)\frac{v_{i,j}}{v_{j,j}}+(1_{\{i\le m,j>m\}}-1_{\{i\le m,j\le m\}})\frac{v_{i,m+n}}{v_{m+n,m+n}},$$

$$w_{i,j}$$

$$w_{i,j}=\left\{\begin{array}[]{ll}(\delta_{i,j}-1)\frac{v_{i,j}}{v_{i,i}v_{j,j}}+\frac{v_{i,m+n}}{v_{i,i}v_{m+n,m+n}}&i=1,\ldots,m;j=m+1,\ldots,m+n-1,\\ (\delta_{i,j}-1)\frac{v_{i,j}}{v_{i,i}v_{j,j}}+\frac{v_{m+n,j}}{v_{j,j}v_{m+n,m+n}}&i=m+1,\ldots,m+n-1;j=1,\ldots,m,\\ (\delta_{i,j}-1)\frac{v_{i,j}}{v_{i,i}v_{j,j}}&i=m+1,\ldots,m+n-1;j=m+1,\ldots,m+n-1.\\ \end{array}\right.$$

$$0\le \frac{v_{i,j}}{v_{i,i}{v}_{j,j}}\le \frac{Q}{q^{2}mn},$$

$$|w_{i,i}|$$

$$|w_{i,j}|$$

$$|w_{i,j}-w_{i,k}|$$

$$|w_{i,i}-w_{i,k}|$$

$$\max(|w_{i,j}|,|w_{i,j}-w_{i,k}|)\leq\frac{3Q}{q^{2}mn}~{}~{}~{}~{}~{}\mbox{for all $i,j,k$}.$$

$$f_{i,j}=\sum _{k=1}^{m+n-1}{f}_{i,k}[({\delta }_{k,j}-1)\frac{v_{k,j}}{v_{j,j}}+(1_{\{k\le m,j>m\}}-1_{\{k\le m,j\le m\}})\frac{v_{k,m+n}}{v_{m+n,m+n}}]+w_{i,j}.$$

$$\Vert V^{-1}\mathbf{x}{\Vert }_{\infty }\le \frac{2c_1{Q}^2}{q^{3}mn}+\frac{|x_{m+n}|}{v_{m+n,m+n}}+\underset{i=1,\dots ,m+n-1}{\max }\frac{|x_i|}{v_{i,i}},$$

$$\begin{array}[]{rcl}F_{i}({\boldsymbol{\theta}})&=&d_{i}-\sum_{k=1}^{n}f(\alpha_{i}+\beta_{k}),i=1,\dots,m,\\ F_{m+j}({\boldsymbol{\theta}})&=&b_{j}-\sum_{k=1}^{m}f(\alpha_{k}+\beta_{j}),j=1,\dots,n-1,\\ F({\boldsymbol{\theta}})&=&(F_{1}({\boldsymbol{\theta}}),\dots,F_{m}({\boldsymbol{\theta}}),F_{m+1}({\boldsymbol{\theta}}),\dots,F_{m+n-1}({\boldsymbol{\theta}}))^{T},\end{array}$$

$$\parallel [F^{{}^{\prime }}(\mathbf{x})-F^{{}^{\prime }}(\mathbf{y})]\mathbf{v}{\parallel }_{\infty }\le K_1\parallel \mathbf{x}-\mathbf{y}{\parallel }_{\infty }\parallel \mathbf{v}{\parallel }_{\infty },$$

$$\underset{i=1,\dots ,m+n-1}{\max }\parallel F_i^{{}^{\prime }}(\mathbf{x})-F_i^{{}^{\prime }}(\mathbf{y}){\parallel }_{\infty }\le K_2\parallel \mathbf{x}-\mathbf{y}{\parallel }_{\infty },$$

Full text

Key words: Affiliation networks; Asymptotic normality; Consistency; Maximum likelihood estimators

Mathematics Subject Classification: 62E20, 62F12.

Running title: Affiliation networks

1 Introduction

Affiliation network is one kind of two-mode social network that consists of two different types of sets of nodes, namely, a set of actors and a set of social events. The network edges indicate the affiliation of actors with social events. Such network data are commonly used to represent memberships between social organizations and their members, for example, the affiliation of the researchers to the academic institutions or interlocking directors to companies or actors to movies. Other scenarios include ceremonial events attended by members faculty sit, social events people attend, trade partners of major oil exporting nations and so on.

In affiliation networks, the actors are brought together to jointly participate in social events. Joint participation in events not only provides the opportunity for actors to interact, but also increases the probability that links (e.g., friendship) between actors form. For example, belonging to the same organizations (boards of directors, political party, labor union, and so on) provides the opportunity for people to meet and interact, and thus a link between individuals is more easily to form in these circumstances. Similarly, when actors participate in more than one event, two events are connected through these actors. There has been increasing interest in analyzing affiliated network data in recent years. A number of approaches are proposed to analyze or model affiliation network data [e.g., Conyon and Muldoon (2004), Robins and Alexander (2004), Snijders et al. (2013)]. Iacobucci and Wasserman (1990) proposed a $p_{2}$ exponential family distribution using the degree sequence as the sufficient statistics to model the weighted affiliation network, which is a close relative of $p_{1}$ model introduced by Holl and Leinhardt (1981). Latapy et al. (2008) extended the basic network statistics used to analyze one-mode networks to give a description of analysis for two-mode networks systematically. Snijders et al. (2013) proposed a stochastic actor-oriented model for the co-evolution of two-mode and one mode networks. By extending exponential random graph models for the one-mode networks, Wang et al. (2009) proposes a number of two-mode specifications as the sufficient statistics in exponential family graph models for two-mode affiliation networks and compared the goodness of fit results obtained using the maximum likelihood and pseudo-likelihood approaches by simulation.

At present, little theoretical results are obtained in affiliation network models although many properties of statistical models for one-mode networks are derived [e.g., Chatterjee and Diaconis (2013), Shalizi and Rinaldo (2013), Bhattacharyya and Bickel (2016)]. Even in the aforementioned simple Iacobucci and Wasserman (1990)’s model by assuming that all edges are independent, the asymptotic properties of the maximum likelihood estimator (MLE) are still not addressed due to a growing dimension of parameter space. In this paper, we study the asymptotic properties of the MLE in an affiliation model with the degree sequence as the exclusively natural sufficient statistic in the exponential family distributions. This model is identical to the Iacobucci and Wasserman (1990) model for unweighed edges (i.e., binary edges). We establish the uniform consistency and asymptotic normality of the maximum likelihood estimator when the number of actors and events both go to infinity. A key step to the proof is that we make use of the approximate inverse of the Fisher information matrix with small approximation errors, which is the extension of that used in Yan et al. (2016a).

The rest of this paper is organized as follows. We introduce the model in Section 2. In Section 3, we present the asymptotic results including the uniform consistency and asymptotic normality of the MLE. Simulation studies and a real application to the student extracurricular affiliation data are given in Section 4. Some further discussions are devoted to Section 5. All the proofs are putted in the appendix.

2 Model

2.1 Notations

Let $A$ be an event set with $m$ events denoted by $\{1,\ldots,m\}$, and $P$ be an actor set with $n$ actors denoted by $\{1,\ldots,n\}$. An affiliation network ${\cal{G}}(m,n)$ records the affiliation of each actor with each event in an affiliation matrix $X=(x_{i,j})_{m\times n}$; $x_{i,j}=1$ if actor $i$ is affiliated with event $j$ and $x_{i,j}=0$ otherwise. Each row of $X$ describes an actor’s affiliation with the events and each column describes the memberships of the event. In practice, $n$ is usually large and $m$ relatively small. Therefore, we assume $m\leq n$ hereafter. The affiliation network $\mathcal{G}(m,n)$ can also be represented in a bipartite graph, in which the nodes are partitioned into two subsets for the actors and the events and the edges exist only between pairs of nodes belonging to different subsets. In bipartite graphs, no two actors are adjacent and no two events are adjacent. If pairs of actors are reachable, it is only via paths containing one or more events. Similarly, if pairs of events are reachable, it is only via paths containing one or more actors.

Define $d_{i}=\sum_{j=1}^{n}x_{i,j}$ as the degree of vertex $i\in A$ and $\mathbf{d}=(d_{1},\dots,d_{m})^{\top}$. Similarly, define $b_{j}=\sum_{i=1}^{m}x_{i,j}$ as the degree of vertex $j\in P$ and $\mathbf{b}=(b_{1},\dots,b_{n})^{\top}$. The pair $\{\mathbf{d,b}\}$ is the degree sequence of the affiliation network ${\cal{G}}(m,n)$.

2.2 An affiliation network model

In this subsection, we present an exponential random bipartite graph model for affiliation networks with the degree sequence as the exclusively natural sufficient statistic. The probability mass function on the affiliation network $\mathcal{G}(m,n)$ is assumed to be of exponential form:

[TABLE]

where $\mathrm{Z}(\boldsymbol{\alpha},\boldsymbol{\beta})$ is the normalizing constant, $\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{m})^{\top}$ and $\boldsymbol{\beta}=(\beta_{1},\dots,\beta_{n})^{\top}$ are parameter vectors. Each affiliation network with the same degree sequence is equally judged. The parameter $\alpha_{i}$ quantifies the popularity of the event $i$ and $\beta_{j}$ quantifies the activity of the actor $j$ to participate in events. Note that

[TABLE]

which implies that the $mn$ random variables $x_{i,j}$ are mutually independent and $\mathrm{Z}(\boldsymbol{\alpha},\boldsymbol{\beta})$ can be expressed as

[TABLE]

Therefore, $x_{i,j}$, $1\leq i\leq m,1\leq j\leq n$ are mutually independent Bernoulli random variables with the success probability:

[TABLE]

which is the Iacobucci and Wasserman (1990) model for affiliated networks with binary edges.

Since the sample is just one realization of the bipartite random graph, the density or probability mass function (2.1) is also the likelihood function. If one transforms $(\boldsymbol{\alpha},\boldsymbol{\beta})$ to $(\boldsymbol{\alpha}-c,\boldsymbol{\beta}+c)$, the likelihood does not change. Following Yan et al. (2016a), we set $\beta_{n}=0$ for the identifiability of the parameter.

The consistency and asymptotic normality of the MLE in the exponential random graph model with the degree sequence for undirected one-model networks or the bi-degree sequence for directed one-mode networks have been established recently [Chatterjee et al. (2011), Hillar and Wibisono (2013), Yan and Xu (2013), Yan et al. (2016a), Yan et al. (2016b)]. For the bipartite networks, the asymptotic theory for the MLE in the model (2.1) has not been explored. The model is closely related to the Rasch model [Rasch (1960)] for dichotomous item response experiments, which assumes that item $i$ correctly gives a response to subject $j$ with probability $\exp(\alpha_{i}-\beta_{j})/(1+\exp(\alpha_{i}-\beta_{j}))$. By assuming that all the parameters in the Rasch model are bounded, Haberman (1977) proved the consistency and asymptotic normality of the MLE when the number of items and subjects goes to infinity simultaneously.

3 Asymptotic results

Let ${\boldsymbol{\theta}}=(\alpha_{1},\ldots,\alpha_{m},\beta_{1},\ldots,\beta_{n-1})^{\top}$ and $\mathbf{g}=(d_{1},\ldots,d_{m},b_{1},\ldots,b_{n-1})^{\top}$. The log-likelihood function is

[TABLE]

The likelihood equations are:

[TABLE]

Let $\widehat{{\boldsymbol{\theta}}}=(\widehat{\alpha}_{1},\dots,\widehat{\alpha}_{m},\widehat{\beta}_{1},\cdots,\widehat{\beta}_{n-1})^{\top}$ be the MLE of ${\boldsymbol{\theta}}$ and $\widehat{\beta}_{n}=0$. If $\widehat{{\boldsymbol{\theta}}}$ exists, then it is the solution to the system of equation (3.1).

Let $V=(v_{i,j})$ be the Fisher information matrix of the parameter vector $\boldsymbol{\theta}$, which is a diagonal dominant matrix with nonnegative entries. The diagonal elements of $V$ are

[TABLE]

Motivated by the approximate inverse proposed by Yan et al. (2016a) for the Fisher information matrix in the directed one-mode network model involved with the bi-degree sequence, we proposed a generalized simple matrix $S=(s_{i,j})$ to approximate the $V^{-1}$, which is defined as

[TABLE]

where $\delta_{i,j}=1$ when $i=j$ and [math] otherwise. For a vector $\mathbf{x}=(x_{1},\dots,x_{n})^{\top}\in\mathbb{R}^{n}$, denote the $\ell_{\infty}$ norm of $\mathbf{x}$ by $\|\mathbf{x}\|_{\infty}=\max_{1\leq i\leq n}|x_{i}|$. We present the consistency of $\widehat{{\boldsymbol{\theta}}}$ here, whose proof is given in Appendix A.

Theorem 1**.**

Assume that ${{\boldsymbol{\theta}}}^{*}\in\mathbb{R}^{m+n-1}$ with $\|{{\boldsymbol{\theta}}}^{*}\|_{\infty}\leq\tau\log n$, where $\tau\in(0,1/24)$ is a constant, and that $X\sim\mathbb{P}_{{{\boldsymbol{\theta}}}^{*}}$, where $\mathbb{P}_{{{\boldsymbol{\theta}}}^{*}}$ denotes the probability distribution (2.1) on $X$ under the parameter ${{\boldsymbol{\theta}}}^{*}$. If $m/n=O(1)$, then as $n$ goes to infinity, with probability approaching one, the MLE $\widehat{{\boldsymbol{\theta}}}$ exists and satisfies

[TABLE]

Further, if the MLE exists, it is unique.

Next, we present the central limit theorem of $\widehat{{\boldsymbol{\theta}}}$, whose proof is given in Appendix B.

Theorem 2**.**

Assume that $X\sim\mathbb{P}_{{{\boldsymbol{\theta}}}^{*}}$. If $m/n=O(1)$ and $\|{{\boldsymbol{\theta}}}^{*}\|_{\infty}\leq\tau\log n$, where $\tau\in(0,1/36)$ then for any fixed $k\geq 1$, as $n\rightarrow\infty$, the vector consisting of the first $k$ elements of $(\widehat{{\boldsymbol{\theta}}}-{{\boldsymbol{\theta}}}^{*})$ is asymptotically multivariate normal with mean zero and covariance matrix given by the upper left $k\times k$ block of $S^{*}$. where $S^{*}$ is the matrix by replacing ${{\boldsymbol{\theta}}}$ in $S$ given in (3.2) with its true value ${{\boldsymbol{\theta}}}^{*}$.

Remark 1**.**

By Theorem 2, for any fixed $i$, as $n\rightarrow\infty$ , the convergence rate of ${\widehat{{{\theta}}}}_{i}$ is $1/v_{i,i}^{1/2}$. Since $me^{-2\parallel{{\boldsymbol{\theta}}^{*}}\parallel_{\infty}}\leq v_{i,i}\leq n/4$, the rate of convergence is between $O(m^{-1/2}e^{\|{{\boldsymbol{\theta}}}^{*}\|_{\infty}})$ and $O(n^{-1/2})$.

4 Simulation studies

We carry out the numerical simulations to evaluate Theorem 2. Following Yan et al. (2016a), the parameter values take a linear form. Specifically, we set $\alpha_{i+1}^{*}=(m-1-i)L/(m-1)$ for $i=0,1,\dots,m-1$, and considered four different values [math], $\log(\log(m))$, $\log(m)^{1/2}$ and $\log(m)$ for $L$. For the parameter vector $\boldsymbol{\beta}$, let $\beta_{j+1}^{*}=(n-1-i)L/(n-1),j=0,1,\dots,n-2$ for simplicity and $\beta_{n}^{*}=0$ by default.

By Theorem 2, $\hat{\xi}_{i,j}=[\hat{\alpha}_{i}-\hat{\alpha}_{j}-(\alpha^{*}_{i}-\alpha^{*}_{j})]/(1/\hat{v}_{i,i}+1/\hat{v}_{j,j})^{1/2}$, $\hat{\eta}_{i,j}=[\hat{\beta}_{i}-\hat{\beta}_{j}-(\beta^{*}_{i}-\beta^{*}_{j})]/(1/\hat{v}_{m+i,m+i}+1/\hat{v}_{m+j,m+j})^{1/2}$, are asymptotically distributed as standard normal distribution, where $\hat{v}_{i,i}$ is the estimate of ${v}_{i,i}$ by replacing $\boldsymbol{\theta}$ with $\boldsymbol{\widehat{\theta}}$. We assess the asymptotic normality of $\hat{\xi}_{i,j},~{}\hat{\eta}_{i,j}$ using the quantile-quantile (QQ) plot. Further, 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. Each simulation is repeated $10,000$ times.

We only simulate a single combination for $(m,n)$ with $m=100,n=200$, and present the QQ-plots of $\hat{v}_{ii}^{-1/2}(\hat{\alpha}_{i}-\alpha_{i})$ and $\hat{v}_{ii}^{-1/2}(\hat{\beta}_{i}-\beta_{i})$ in Figure 1 and Figure 2, respectively. The horizontal and vertical axes are the theoretical and empirical quantiles, respectively, and the straight lines correspond to the reference line $y=x$. In Figure 1, we can see that the empirical quantiles coincide with the theoretical ones very well. In Figure 2, there are slight deviations when $L=(\log m)^{1/2}$. When $L=(\log n)^{1/2}$, there are a little derivations on both tails of plots. When $L=\log n$, the MLE does not exist in all repetitions, the QQ plots are not available in this case.

The coverage probability of the $95\%$ confidence interval for $\alpha_{i}-\alpha_{j}$ and $\beta_{i}-\beta_{j}$, the length of the confidence interval, and the frequency that the MLE did not exist, which are reported in Table 1 and Table 2, respectively. There are similar results for the two tables. We can see that the length of estimated confidence interval increases as $L$ increases for fixed $n$, and decreases as $n$ increases for fixed $L$. The coverage frequencies are all close to the nominal level $95\%$. When $L=(\log n)^{1/2}$ (conditions in Theorem 2 no longer hold), the MLE does not exist with a positive probability; when $L=\log(n)$, the MLE did not exist with $100\%$ frequencies.

A data example. We analyze a student extracurricular affiliation network data collected by Dan McFarland in 1996, which can be downloaded from http://dl.dropbox.com/u/25710348/snaimages/mag_act96.txt. It consists of $1295$ students (anonymized) and $91$ student organizations in which they are members (e.g. National Honor Society, wrestling team, cheerleading squad, etc.). In order to guarantee the existence of the MLE, we remove those $438$ individuals that don’t belong to any organizations. The MLEs of the parameters for remaining students and organizations and their standard errors as well as the $95\%$ confidence intervals are reported in Tables 3 and 4. The value of estimated parameters reflect the size of degrees. For example, the largest five degrees in student organizations are $199,157,124,93,89$ for organizations “Spanish.Club, Pep.Club, NHS, Latin.Club, Orchestra.Symphonic”, which also have the top five influence parameters $-0.32,-0.64,-0.94,-1.92,-1.34$. On the other hand, the organizations with the five smallest influence parameters $-4.60,-4.60,-4.60,-4.89,-4.89$ have degrees $4,4,4,3,3$.

5 Summary and discussion

Statistical models for affiliation networks provide insight into the formulation of complex social affiliation between actors and events. They also indirectly reflect how events create ties among actors and the actors create ties among events. Meanwhile, the asymptotic inference in these models are challenge like other network models due to that the structure of the network data is non-standard. In this paper, we derive the uniform consistency and asymptotic normality of the MLE in the exponential random bipartite graph models for affiliation networks with the degree sequence as the exclusively sufficient statistic. The conditions imposed on $\tau$ that guarantee the good asymptotic properties of the MLE may not be best possible. In particular, the conditions guaranteeing the asymptotic normality are stronger than those guaranteeing the consistency. Simulation studies suggest that the conditions on $\tau$ might be relaxed. The asymptotic behavior of the MLE depends not only on $\tau$, but also on the configuration of all the parameters. We will investigate this problem in the future work.

We only consider dyadic independence assumption. Like the model specifications in the exponential random graph models for the one-mode network, one can add the counts of $k$-stars of various sizes as the model terms in the model (2.1) to model the dependence structure. However, such model incurs the degeneracy problem [e.g., Chatterjee and Diaconis (2013)] in which the generated graphs are almost full or empty. To overcome this disadvantage, Wang et al. (2009) proposes a number of new two-mode specifications such as the k-two-paths and three-path statistics. Although Wang et al. (2009) show that these model specifications have good performance by simulations, the theoretical properties of the model are still unknown.

Appendix

In this appendix, we will present the proofs for Theorems 1 and 2. We start with some preliminaries. For an $n\times n$ matrix $J=(J_{i,j})$, $\|J\|_{\infty}$ denotes the matrix norm induced by the $\|\cdot\|_{\infty}$-norm on vectors in $\mathbb{R}^{n}$:

[TABLE]

Let $D$ be an open convex subset of $\mathbb{R}^{n}$. We say an $n\times n$ function matrix $F(\mathbf{x})$ whose elements $F_{ij}(\mathbf{x})$ are functions on vectors $\mathbf{x}$, is Lipschitz continuous on $D$ if there exists a real number $\lambda$ such that for any $\mathbf{v}\in R^{n}$ and any $\mathbf{x,y}\in D,$

[TABLE]

where $\lambda$ may depend on $n$ but independent of $\mathbf{x}$ and $\mathbf{y}$. For fixed $n$, $\lambda$ is a constant.

We introduce a class of matrices. Given two positive numbers $q,Q$, we say the $(m+n-1)\times(m+n-1)$ matrix $V=(v_{i,j})$ belongs to the class ${\cal L}_{m,n}(q,Q)$ if the following holds:

[TABLE]

If $V\in{\cal L}_{m,n}(q,Q)$, then $V$ is a $(m+n-1)\times(m+n-1)$ diagonally dominant, symmetric nonnegative matrix. Define $v_{m+n,i}=v_{i,m+n}:=v_{i,i}-\sum_{j=1}^{m+n-1}v_{i,j}$ for $i=1,\dots,m+n-1$ and $v_{m+n,m+n}=\sum_{i=1}^{m+n-1}v_{m+n,i}.$ Then $q\leq v_{m+n,i}\leq Q$ for $i=1,\dots,m,v_{m+n,i}=0$ for $i=m,m+1,\dots,m+n-1$ and $v_{m+n,m+n}=\sum_{i=1}^{m}v_{i,m+n}=\sum_{i=1}^{m}v_{m+n,i}.$ The Fisher information matrix of the parameter vector $\boldsymbol{\theta}$, $V$, belongs to the matrix class $V\in{\cal L}_{m,n}(q,Q)$. The approximate error using $S$ in (3.2) to approximate the inverse of $V$ is given in the lemma below, whose proof is the extension of that for Proposition 1 in Yan et al. (2016a).

Lemma 1**.**

If $V\in{\cal L}_{m,n}(q,Q)$ with $Q/q=o(n)$ and $m/n=O(1)$, then for large enough $n$,

[TABLE]

where $c_{1}$ is a constant that dose not depend on $q,Q,m$ and $n$, and $\|A\|:=\max_{i,j}|a_{i,j}|$ for a general matrix $A=(a_{i,j})$.

Proof.

Recall that if $V\in\mathcal{L}_{m,n}(q,Q)$, then for $i=1,\ldots,m+n$,

[TABLE]

and if $v_{i,j}>0$ for $1\leq i\leq m,1\leq j\leq n,i\neq j$, then

[TABLE]

The above equation and inequality will be repeatedly used in this proof.

Let $I$ denote the $(m+n-1)\times(m+n-1)$ identity matrix. Define $F=(f_{i,j})=V^{-1}-S$, $U=(u_{i,j})=I-VS$ and $W=(w_{i,j})=SU$. Then we have the recursion

[TABLE]

Note that for $i=1,\ldots,m,~{}j=1,\ldots,m$, we have

[TABLE]

Similarly, we also have

[TABLE]

In all, $u_{i,j}$ can be written in a unified form:

[TABLE]

where $1_{\{\cdot\}}$ is an indicator function. Similarly, for $i=1,\ldots,m,~{}j=1,\ldots,m$, we have

[TABLE]

and

[TABLE]

Further, when $1\leq i\neq j\leq m+n$,

[TABLE]

and it is not difficult to show that, when $i,j,k$ are different from each other, we have

[TABLE]

It follows that

[TABLE]

Next we use the recursion (5.3) to obtain an upper bound of the approximate error $\|F\|$. By (5.3) and (5.4), for any $i\in\{1,\ldots,m+n-1\}$, we have that for $j=1,\ldots,m+n-1$

[TABLE]

Thus, to prove Lemma 1, it is sufficient to show that $|f_{i,j}|\leq c_{1}Q^{2}/(q^{3}mn)$ for any $i,j$. The condition $m/n=O(1)$ guarantees that $Q/q=o(m)$ if $Q/q=o(n)$. The remainder of the proof of Lemma 1 is similar to the proof of Proposition 1 in Yan et al. (2016a), and we omit the details here.

∎

Note that if $Q$ and $q$ are bounded constants, then the upper bound of the above approximation error is on the order of $(mn)^{-1},$ indicating that $S$ is a high-accuracy approximation to $V^{-1}$. Further, based on the above proposition, we immediately have the following lemma.

Lemma 2**.**

If $V\in{\cal L}_{n}(q,Q)$ with $Q/q=o(n)$ and $m/n=O(1)$, then for a vector $\mathbf{x}\in R^{m+n-1}$,

[TABLE]

where $x_{m+n}:=\sum_{i=1}^{m}x_{i}-\sum_{i=m+1}^{m+n-1}x_{i}$.

Similar to Theorem 8 in Yan et al. (2016a), we have the following lemma for the rate of convergence for the Newton’s iterative sequence to solve the system of the likelihood equations, whose proof is similar to that in Yan et al. (2016a) and we omit it here.

Lemma 3**.**

Define a system of equations:

[TABLE]

where $f(\cdot)$ is a continuous function with the third derivative. Let $D\subset\mathbb{R}^{m+n-1}$ be a convex set and assume for any $\mathbf{x,y,v}\in D$, we have

[TABLE]

[TABLE]

where $F^{{}^{\prime}}({\boldsymbol{\theta}})$ is the Jacobian matrix of $F$ on ${\boldsymbol{\theta}}$ and $F^{{}^{\prime}}_{i}({\boldsymbol{\theta}})$ is the gradient function of $F_{i}$ on ${\boldsymbol{\theta}}$. Consider ${\boldsymbol{\theta}}^{(0)}\in D$ with $\Omega({\boldsymbol{\theta}}^{(0)},2r)\subset D$ where $r=\parallel[F^{{}^{\prime}}({\boldsymbol{\theta}}^{(0)})]^{-1}F({\boldsymbol{\theta}}^{(0)})\parallel_{\infty}$ for any ${\boldsymbol{\theta}}\in\Omega({\boldsymbol{\theta}}^{(0)},2r)$. We assume that $m/n=O(1)$ and

[TABLE]

For $k=1,2,\ldots$, define the Newton iterates ${\boldsymbol{\theta}}^{(k+1)}={\boldsymbol{\theta}}^{(k)}-[F^{{}^{\prime}}({\boldsymbol{\theta}}^{(k)})]^{-1}F({\boldsymbol{\theta}}^{(k)})$. Let

[TABLE]

If $\rho<1/2,$ then ${\boldsymbol{\theta}}^{(k)}\in\Omega({\boldsymbol{\theta}}^{(0)},2r),k=1,2,\ldots$, are well defined and satisfy

[TABLE]

Further, $\lim_{k\rightarrow\infty}{\boldsymbol{\theta}}^{(k)}$ exists and the limiting point is precisely the solution of $F({\boldsymbol{\theta}})=0$ in the rage of ${\boldsymbol{\theta}}\in\Omega({\boldsymbol{\theta}}^{(0)},2r)$.

Appendix A: Proofs for Theorem 1

We define a system of functions:

[TABLE]

Note that the solution to the equation $F({\boldsymbol{\theta}})=0$ is precisely the MLE. Then the Jacobin matrix $F^{{}^{\prime}}({\boldsymbol{\theta}})$ of $F({\boldsymbol{\theta}})$ can be calculated as follows. For $i=1,\dots,m,$

[TABLE]

[TABLE]

and for $j=1,\dots,n-1$,

[TABLE]

[TABLE]

Since $e^{x}/(1+e^{x})^{2}$ is a decreasing function on $x$ when $x\geq 0$ and an increasing function when $x\leq 0$. Consequently, for any $i,j$, we have

[TABLE]

According to the definition of ${\cal L}_{m,n}(q,Q)$, we have that $-F^{{}^{\prime}}({\boldsymbol{\theta}})\in{\cal L}_{m,n}(q,Q)$, where

[TABLE]

Therefore, Lemma 3 and Lemma 1 can be applied. Let ${{\boldsymbol{\theta}}}^{*}$ denote the true parameter vector. The constants $K_{1},K_{2}$ and $r$ in the upper bounds of Lemma 3 are given in the following lemma.

Lemma 4**.**

Take $D=R^{m+n-1}$ and ${{\boldsymbol{\theta}}}^{0}={{\boldsymbol{\theta}}}^{*}$ in Lemma 3. Assume

[TABLE]

If $m/n=O(1)$, then we can choose the constants $K_{1},K_{2}$ and $r$ in Lemma 3 as

[TABLE]

where $c_{11},c_{12}$ are constants.

Proof.

For fixed $m,n$, we first derive $K_{1}$ and $K_{2}$ in the inequalities of Lemma 3. Let $\mathbf{x},\mathbf{y}\in R^{m+n-1}$ and

[TABLE]

Then, for $i=1,\dots,m$, we have

[TABLE]

[TABLE]

[TABLE]

Note that

[TABLE]

By the mean value theorem for vector-valued functions (Lang (1993), p.341), we have

[TABLE]

where

[TABLE]

Therefore,

[TABLE]

Similarly, for $i=m+1,\dots,m+n-1,$ we also have $F^{{}^{\prime}}_{i}(\mathbf{x})-F^{{}^{\prime}}_{i}(\mathbf{y})=J^{(i)}(\mathbf{x}-\mathbf{y})$ and $\sum_{s,l}|J^{(i)}_{(s,l)}|\leq m$. Consequently,

[TABLE]

and for $\forall~{}\mathbf{v}\in R^{m+n-1}$,

[TABLE]

So we can choose $K_{1}=n$ and $K_{2}=n/2$ in Lemma 3.

It is obvious that $-F^{{}^{\prime}}({{\boldsymbol{\theta}}}^{*})\in{\cal L}_{n}(q_{*},Q_{*})$, where

[TABLE]

Note that

[TABLE]

By the assumption of A1 and Lemma 2, if $m/n=O(1)$, then we have

[TABLE]

where $c_{11},c_{12}$ are constants. ∎

The following lemma assures that condition (A1) holds with a large probability.

Lemma 5**.**

With probability at least $1-4/n$, we have

[TABLE]

Proof.

Note that $m<n$. By Hoeffding (1963)’s inequality, we have

[TABLE]

Therefore,

[TABLE]

Similarly, we have

[TABLE]

Consequently,

[TABLE]

This is equivalent to Lemma 5. ∎

Proof of Theorem 1.

Assume that condition (A1) holds. Recall that the Newton’s iterates in Lemma 3, ${\boldsymbol{\theta}}^{(k+1)}=[F^{{}^{\prime}}({\boldsymbol{\theta}}^{(k)})]^{-1}F({\boldsymbol{\theta}}^{(k)})$ with ${\boldsymbol{\theta}}^{(0)}={\boldsymbol{\theta}}^{*}.$ If ${\boldsymbol{\theta}}\in\Omega({\boldsymbol{\theta}}^{*},2r)$, then $-F^{{}^{\prime}}({{\boldsymbol{\theta}}}^{*})\in{\cal L}_{m,n}(q,Q)$ with

[TABLE]

By Lemma 4 and condition(A1), for sufficient small $r$,

[TABLE]

If ${{\boldsymbol{\theta}}}^{*}\in\mathbb{R}^{m+n-1}$ with $\|{{\boldsymbol{\theta}}}^{*}\|_{\infty}\leq\tau\log n$, where $0<\tau<1/24$ is a constant, then as $n\rightarrow\infty$,

[TABLE]

Therefore, $\rho r\rightarrow 0$ as $n\rightarrow\infty$. By Lemma 3, $\lim_{n\rightarrow\infty}{\widehat{{\boldsymbol{\theta}}}}^{(n)}$ exists. Denote the limit as $\widehat{\boldsymbol{\theta}}$. Then it satisfies

[TABLE]

By Lemma 5, condition (A1) holds with probability approaching to one, thus the above inequality also holds with probability approaching to one. The uniqueness of the MLE comes from Proposition 5 in Yan et al. (2016a). ∎

Appendix B: Proofs for Theorem 2

We first present one proposition. Since $d_{i}=\sum_{k}a_{i,k}$ and $b_{j}=\sum_{k}a_{k,j}$ are sums of $m$ and $n$ independent random variables, by the central limit theorem for the bounded case in Loève (Loeve (1977), page 289), we know that ${v_{i,i}}^{-1/2}(d_{i}-\mathbb{E}(d_{i}))$ and ${v_{m+j,m+j}}^{-1/2}(b_{j}-\mathbb{E}(b_{j}))$ are asymptotically standard normal if $v_{i,i}$ and $v_{m+j,m+j}$ diverge, respectively. Note that

[TABLE]

Then we have the following proposition.

Proposition 1**.**

Assume that $X\sim\mathbb{P}_{{{\boldsymbol{\theta}}}^{*}}$. If $e^{\parallel{{\boldsymbol{\theta}}}^{*}\parallel_{\infty}}=o(n^{1/2})$, then for any fixed $k\geq 1$, as $n\rightarrow\infty$, the vector consisting of the first $k$ elements of $S\{\mathbf{g}-\mathbb{E}(\mathbf{g})\}$ is asymptotically multivariate normal with mean zero and covariance matrix given by the upper left $k\times k$ block of $S$.

To complete the proof of Theorem 2, we need two lemmas as follows.

Lemma 6**.**

Let $R=V^{-1}-S$ and $U=Cov[R\{\mathbf{g}-\mathbb{E}\mathbf{g}\}]$. Then

[TABLE]

Proof.

Note that

[TABLE]

where $I$ is a $(m+n-1)\times(m+n-1)$ diagonal matrix, and by (5.5), we have

[TABLE]

Thus,

[TABLE]

∎

Lemma 7**.**

Assume that the conditions in Theorem 1 hold. If $\|{{\boldsymbol{\theta}}}^{*}\|_{\infty}\leq\tau\log n$ with $\tau<1/24$ and $m/n=O(1)$, then for any $i$,

[TABLE]

Proof.

By Theorem 1, if $m/n=O(1)$, then we have

[TABLE]

Let $\widehat{\gamma}_{i,j}=\widehat{\alpha}_{i}+\widehat{\beta}_{j}-\alpha_{i}^{*}-\beta_{j}^{*}$. By the Taylor expansion, for any $1\leq i\leq m,1\leq j\leq n,i\neq j$,

[TABLE]

where

[TABLE]

and $0\leq\phi_{i,j}\leq 1$. By the likehood equations (3.1), it is not difficult to verify that

[TABLE]

where $\mathbf{h}=(h_{1},\dots,h_{m+n-1})^{T}$ and

[TABLE]

Equivalently,

[TABLE]

By (A2), it is easy to show

[TABLE]

Because

[TABLE]

Therefore

[TABLE]

Note that $(S\mathbf{h})_{i}=h_{i}/v_{i,i}+(-1)^{1_{\{i>m\}}}h_{m+n}/v_{m+n,m+n},$ and $(V^{-1}\mathbf{h})_{i}=(S\mathbf{h})_{i}+(R\mathbf{h})_{i}$. Then we have

[TABLE]

by Lemma 1, we have

[TABLE]

If $\parallel{{\boldsymbol{\theta}}}^{*}\parallel_{\infty}\leq\tau\log n$, and $\tau<1/36$, then $\mid(V^{-1}\mathbf{h})_{i}\mid\leq\mid(S\mathbf{h})_{i}\mid+\mid(R\mathbf{h})_{i}\mid=o(n^{-1/2})$.

∎

Proof of Theorem 2.

By (C1), we have

[TABLE]

By Lemma 6, we have

[TABLE]

By Chebyshev’s inequality, if $\parallel{{\boldsymbol{\theta}}}^{*}\parallel_{\infty}\leq\tau\log n$, $m/n=O(1)$, and $\tau<1/36$, then

[TABLE]

Therefore, we have

[TABLE]

By Lemma 6 and Lemma 7, we have

[TABLE]

Theorem 2 follows directly from Proposition 1. ∎

Acknowledgements

We are very grateful to one anonymous referee and the Editor for their valuable comments that have greatly improved the manuscript. Qin’s research is partially supported by the National Natural Science Foundation of China (No.11271147, 11471135). Yan’s research is partially supported by the National Natural Science Foundation of China (No.11401239) and the self-determined research funds of CCNU from the colleges’s basic research and operation of MOE (CCNU15A02032, CCNU15ZD011) and a fund from KLAS (130026507). Zhang’s research is supported by a fund from CCNU (2016CXZZ157).