Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
63,778 characters · 15 sections · 82 citation commands
Regression Modeling of the Count Relational Data with Exchangeable Dependencies
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Count data; Multiplicative model; Social network analysis; Weighted directed network; Weak exchangeability.
\spacingset{1.8}
The modeling of relational data has garnered profound interests across various domains as it leads to a more comprehensive understanding of relationships in complex systems. Examples include deciphering brain connectivity maps Zhao2014differentialNetworks, zhang2020mixed, understanding friendship relationships, Moody2011Popularity, and characterizing economic networks han2020individual. For such type of data, an important goal is to infer the mechanisms responsible for the observed relations, taking into account additional information on covariates within the relational structure.
Consider relational data measured on pairs of $n$ nodes, where the directed edges may be assigned specific weights. Efforts on modeling relational data have scattered in literature, with seminal examples including the social relations model Warner1979roundRobin,Wong1982roundRobin,kenny1984srm, snijders1999srm), which assumes normally distributed data and additive effects, and the row-column exchangeable model aldous1985exchangeability. In these modeling frameworks, the dependence structure within relational data are characterized by the latent variables. To further incorporate with the possibly additional covariates information, Hoff2002Latent developed latent space models, where they model the relational data as conditionally independent given the unobserved positions in social space of two nodes and the observed covariates that measure characteristics of the relational structure. Aligned with the latent space models, the latent factor models Hoff2005mixedEffects,westveld2011mem and additive and multiplicative effects (AME) models Hoff2013Likelihoods have been proposed as parametric models based on the latent factors. Consequently, in contrast to the aforementioned approach, a more generalized model is imposed, referred to as dyadic regression models graham2020dyadic, which operates without specifying a particular model form dependent on latent factors. Beyond modeling a single network, there also exists a rich body of literature on temporal modeling of dynamic relational data suening2017network,kim2018review and multiple relational data zhang2018network. Despite a handful of efforts on modeling relational data zhang2017estimating,Li2019predictNetwork,Hoff2021Additive,Le2022predictNetwork as well as relational arrays Harris2011Longitudinal,Banerjee2013diffusiion,Marrs2023exchangeableErrors, methods specifically designed for analyzing relational data with count observations remain limited.
Relational data characterized by weighted edges of count measurements are widely observed Krivitsky2012exponentialGraphs, which we referred to as “count relational data”. Examples can be found in various domains such as coauthorship and citation networks ji2016coauthorship, mobile phone communication networks Dong2012HMMs, email exchange networks diesner2005exploration, and transportation networks based on traffic counts Wang2009trafficCounts. A naive yet widely employed approach is to convert count edges to unweighted relational data with binary outcomes, which however may lead to information loss and spurious discoveries. Among limited studies focused on the count relational data, Krivitsky2012exponentialGraphs extended the exponential-family random graph models (ERGMs) to encompass the count outcomes. As a pioneering method for estimating covariate effects on network data, ERGMs Holland1981exponentialGraphs,frank1986markov relies on Markov chain Monte Carlo (MCMC) approximations to facilitate estimation that hinders its applicability for large relational data. Also, the maximum likelihood estimator for ERGMs could be time-consuming caimo2011bayesian,schmid2017exponential and has been found to be inconsistent under many network models Shalizi2013ConsistencyERGM. Apart from ERGMs, there exists some efforts by utilizing multivariate counting processes to model counts of interactions when incorporating continuous time. For instance, perry2013point proposed a continuous-time model for dynamic data featuring directed counts of interactions, building upon event history analysis and assuming that no two interactions take place simultaneously. We refer to chen2023degree for further discussions on such kind of dynamic models for continuous time relational data. Different from those works, we focus on a single network, employing a distinct approach to access the network dependencies.
Consider a motivating example of a food sharing network data collected by koster2014food, which contains the number of transferred gifts among households over a yearlong period. Figure (ref) illustrates the food sharing network between households together with their game harvests. Gift transfers, particularly those of larger volume, mainly flow from households with substantial harvests to those with smaller yields. Households with lower game harvests tend to receive more gifts than those with higher harvests. The authors adopt a Bayesian approach for analysis, noting that most existing methods are restricted to continuous response data, whereas the response variable in this case is a count. To address the limitations of current methods discussed above, we aim to propose a modeling framework for count relational data that accounts for network dependencies, which enjoys theoretical guarantees and computational efficiency.
To address the challenges for modeling count relational data with the incorporation of node and edge covariates, we propose a latent multiplicative Poisson model. The main advantage is that our model facilitates the characterization of the edge dependency in relational data $\{y_{ij}\}_{1\leq i \neq j \leq n}$ directly through the latent errors $\{e_{ij}\}_{1\leq i \neq j \leq n}$, which yields an extra level of flexibility in the network structure that cannot be handled by ERGMs. In addition, we avoid imposing any parametric specification on the latent errors. This sets our work apart from the latent space model Hoff2002Latent and other approaches Warner1979roundRobin,li2002unified which enforce a parametric model of the latent factors. The estimation of regression coefficients $\boldsymbol{\beta}$ therefore needs to be tailored for the absence of full knowledge of likelihood. To this end, we employ the pseudo-maximum likelihood (PML, Gourieroux1984pseudoMLE) approach. We establish the asymptotic properties of estimated coefficients under network dependencies. To lay the groundwork for the inference procedure, we further propose the estimation procedure of the covariance among $\{e_{ij}\}_{1\leq i \neq j \leq n}$. Specifically, edge dependency within relational data can be empirically estimated by the frequencies of small subgraphs between two or three nodes opsahl2009clustering. When incorporating the proposed regression model, after getting the consistent estimation of coefficients employing the PML method, we demonstrate the consistency of the estimated covariance parameters via function of network moments Zhang2022edgeworth. Our procedure is easily implemented and provides a computationally efficient approach for modeling count relational data compared to the existing methods that rely on MCMC approaches.
It is noteworthy that our proposed model differs from existing efforts that have utilized node-specific fixed effects to study the relational data graham2017econometric,dzemski2019empirical, chen2021nonlinear,zhang2023generalized. We leverage a more general formulation that takes advantage of introducing the weak exchangeable Silverman1976weaklyExch errors, which lends to us a concise representation of the covariance matrix. Furthermore, our model sets itself apart from the model in graham2020dyadic by relaxing the assumption that any pair of $(y_{ij},y_{kl})$ sharing at least one index in common are dependent, and accommodating edge covariates, which may not be encoded by the observable node information.
Below, we introduce the latent multiplicative Poisson model on count relational data in Section (ref) and describe the structure of latent errors. In Section (ref), we present the PML estimation procedure and analyze the asymptotic normality of proposed estimator. With the consistently estimated covariance parameters of latent errors proposed in Section (ref), we design valid inference procedures on the regression coefficients. Simulation studies are presented in Section (ref) to demonstrate the performance of our method. The application of our model to the food sharing network is given in Section (ref), where we analyze social activities among households in Nicaragua. We conclude with discussions on future directions in Section (ref). Technical proofs are given in the Supplementary Material.
Notation. The gradient and Hessian matrix of function $ g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\boldsymbol{\beta})$ with respect to $\boldsymbol{\beta}$ are represented by $\nabla g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\boldsymbol{\beta})$ and $\nabla^2 g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\boldsymbol{\beta})$, respectively. For ease of presentation, $\nabla g(\boldsymbol{\beta})$ and $\nabla^2 g(\boldsymbol{\beta})$ are used as their simplification. Denote the first and second derivatives of $g(\cdot) $ as $g'(\cdot) = \mathrm{d}g(z)/\mathrm{d}z $ and $g''(\cdot) = \mathrm{d}g'(z)/\mathrm{d}z$. The diagonal matrix with off-diagonals $\{a_{ij}\}$ is denoted as ${\rm diag}\{a_{ij}\}$. Unless specified otherwise, we write the $\ell_2$-norm of a vector as $\|{\mathbf a}\|$. Let $\xrightarrow{d}$, $\xrightarrow{p}$ and $\stackrel{\mathbb{P}_{\boldsymbol{\beta}}}{\rightarrow}$ represent the convergence in distribution, in probability, and in probability given $\boldsymbol{\beta}$, respectively. Throughout the rest of the paper, $\{\cdot_{ij}\}$ is used to denote the set of variables indexed by $i,j$, \sl i.e.\; $\{\cdot_{ij}\}_{1\leq i \neq j \leq n}$, for clarity.
Consider observing the count outcomes $\{y_{ij}\}$ among $n$ nodes indexed by $1,\ldots,n$ as well as the covariates $\{{\mathbf x}_{ij} : {\mathbf x}_{ij} \in \mathbb{R}^p\}$, where self-loops are excluded as node's interaction with itself is not of our interest. As discussed in Section (ref), existing models that focus on binary outcomes Hoff2002Latent or with additive errors minhas2019inferential,Le2022predictNetwork,Marrs2023exchangeableErrors do not naturally handle count data. Inspired by the conditional autoregressive models on counting process Zeger1988Count,Davis2000Count,Diggle2002Count, we introduce the following latent multiplicative Poisson model:
where $\boldsymbol{\beta} = (\beta_1, \beta_2, \ldots, \beta_p)^{{ \mathrm{\scriptscriptstyle T} }} \in \mathbb{R}^p$ is the vector of regression coefficients, $g: \mathbb{R} \rightarrow (0, \infty)$ serves as the link function, and the latent error $e_{ij}$ admits unit mean such that $\mathbb{E}(\lambda_{ij}) = g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\boldsymbol{\beta})$. The choice of $g(\cdot)$ is pre-specified, such as logistic function, arc-cotangent function, and exponential function, to name a few. In model (ref), the dependence across observations are modeled via latent errors $\{e_{ij}\}$, whose distribution does not need to be specified and therefore leads to great flexibility of our model.
The multiplicative nature of $\lambda_{ij}$ with respect to the regression component and the latent error benefits in two ways. First, it facilitates an analogous way to define the residual of count relational data without specifying a stringent parametric model. Second, it establishes a direct connection between the data dependence represented by the covariance of $\{y_{ij}\}$, and the dependence among latent errors characterized by the covariance of $\{e_{ij}\}$, as demonstrated in (ref). This cannot be achieved via an additive model for $\{\lambda_{ij}\}$ as the positivity of $\{\lambda_{ij}\}$ is not easily aligned with traditional assumption on $\mathbb{E}(e_{ij})=0$. Under model (ref), we have
and
Let $\xi_{ij}= y_{ij} \{g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\boldsymbol{\beta})\}^{-1}$, we have $\mathrm{Cov}(e_{ij},e_{i'j'}) = \mathrm{Cov}(\xi_{ij},\xi_{i'j'})$ by (ref), which hints a natural covariance estimator of $\{e_{ij}\}$ in Section (ref). Moreover, covariance terms of $\{e_{ij}\}$ in (ref), namely $\mathrm{Cov}(e_{ij},e_{ji}), \mathrm{Cov}(e_{ij},e_{ik}), \mathrm{Cov}(e_{ij},e_{kj})$, and $\mathrm{Cov}(e_{ij},e_{ki})$ represent the commonly-encountered network effects cranmer2014reciprocityEffect,du2024optimal, which we refer to as the reciprocity effect, same sender effect, same receiver effect, and sender-receiver effect, respectively. Such dependency structure is denoted as the social relations covariance model in Hoff2021Additive under parametric model assumptions.
We are now in position of modeling the latent errors $\{e_{ij}\}$, which lends a concise representation of the dependence among edges $\{y_{ij}\}$. Instead of imposing parametric assumptions, we only assume that the latent errors are weakly exchangeable Silverman1976weaklyExch. An array ${\mathbf z}=\{z_{i,j}\}$ is called weakly exchangeable if $\{z_{i,j}\}\stackrel{d}{=}\{z_{\pi(i), \pi(j)}\}$ for any simultaneous permutation $\pi(\cdot)$ of both the row and column labels. Such an assumption is desirable for relational data as it is of great interest to relate the outcomes involving node $i$ as a sender to that involving $i$ as a receiver. In this paper, we focus on the dissociated weakly exchangeable array to model latent variables, where any random variables within the array are independent whenever their indexing sets are disjoint.
A major appeal of the weak exchangeability is the concise parametrization of the covariance matrix of ${\mathbf e} =(e_{12}, e_{13}, \ldots, e_{n-1n})^{{ \mathrm{\scriptscriptstyle T} }}\in \mathbb{R}^{n^2-n}$, denoted by $\boldsymbol{\Omega}_e := \mathbb{E}\big\{\big({\mathbf e} - \mathbb{E}({\mathbf e})\big)\big({\mathbf e} - \mathbb{E}({\mathbf e})\big)^{{ \mathrm{\scriptscriptstyle T} }}\big\}$. In fact, six unique parameters are sufficient to parameterize the $O(n^4)$ entries of $\boldsymbol{\Omega}_e$, since there exist six distinguishable configurations of pairs drawn from $\{e_{ij}\}$ with unlabeled nodes Hoff2021Additive,Marrs2023exchangeableErrors. Specifically, the diagonal of $\boldsymbol{\Omega}_e$ admits $\eta_1=\mathrm{Var}(e_{ij})$, while off-diagonals $\mathrm{Cov}(e_{ij}, e_{ji}), \mathrm{Cov}(e_{ij}, e_{il}), \mathrm{Cov}(e_{ij}, e_{kj})$ and $\mathrm{Cov}(e_{ij}, e_{ki})$ are represented by $\eta_2, \eta_3, \eta_4$ and $\eta_5$, respectively. Additionally, we assume $\mathrm{Cov}(e_{ij}, e_{kl})=0$ for $\{i,j\}\cap\{k,\ell\}=\emptyset$ according to the dissociated array assumption Silverman1976weaklyExch. With such a parameterization, the multiplicities of $\eta_1$ to $\eta_5$ in each row/column of $\boldsymbol{\Omega}_e$ are $1, 1, n-2, n-2, 2n-4$, respectively, while remaining entries are zero. Denote $\boldsymbol{\eta}:=(\eta_1, \eta_2, \eta_3, \eta_4, \eta_5)^{{ \mathrm{\scriptscriptstyle T} }}$. For the non-negative definiteness of $\boldsymbol{\Omega}_e$, $\boldsymbol{\eta}$ should fall in the following parameter space:
where $\iota = (\eta_4^2+\eta_3^2)(n^2-2n+1) + 4\eta_5^2(n^2-6n+9) + 2\eta_3\eta_4(1-n^2+2n)$ and $\kappa = \eta_2\eta_5(8n- 24) + (\eta_3+\eta_4)\eta_5(12-4n) + 4\eta_2\{\eta_2-(\eta_3+\eta_4)\}$. Details are relegated to the Supplementary Material. In practice, this helps to establish an valid estimation of $\boldsymbol{\eta}$ to draw inference on $\boldsymbol{\beta}$, which will be discussed in Section (ref).
It is common to assume dependence between edges with sharing nodes, such as in the social relations model Warner1979roundRobin,Wong1982roundRobin,kenny1984srm,Gill2001srm, the conditionally independent dyad model chandrasekhar2016cid,graham2020cid, and the random-effects model gelman2006ref,westveld2011mem,aronow2015cluster,graham2021minimax. For relational data, those dependencies characterized by covariance terms as defined in Section (ref), representing the network effects between edges sharing common nodes. In our work, we assume at least one of $\{\eta_3,\eta_4,\eta_5\}$ is nonzero. In practice, the weakly exchangeable array can be easily generated from a variety of widely-used models, as exemplified in Example (ref). The proof of its weak exchangeability is given in the Supplementary Material.
As the primary task to analyze the relational data under model (ref), we estimate the regression coefficients $\boldsymbol{\beta} = (\beta_1, \beta_2, \ldots, \beta_p)^{{ \mathrm{\scriptscriptstyle T} }}$ by utilizing the pseudo-likelihood, from which the estimator's asymptotic distribution is carefully established to draw inference on $\boldsymbol{\beta}$.
The dependence among $\{e_{ij}\}$ imposes difficulty to estimate the regression coefficients, since integrating out the latent variables requires the knowledge of the joint distribution of $\{e_{ij}\}$. However, in our model, the weak exchangeability of $\{e_{ij}\}$ is a mild assumption, where the joint distribution of errors does not need to be specified. To address this, we employ the pseudo-maximum likelihood (PML, Gourieroux1984pseudoMLE) method. For data comes from a linear exponential family, PML necessitates only the accurate specification of the mean structure of data generating processes to ensure consistent estimation of coefficients. Though it can not easily work with dependence as pointed by Besag1975PML, the observations are independent conditional on the error terms. We therefore work on the conditional pseudo-likelihood since our target $\boldsymbol{\beta}$ is on the mean structure of model (ref). Enlightened by this approach, we maximize a likelihood function associated with a family of probability distributions. Conditional on latent errors $e_{ij}$, the distribution of $\{y_{ij}\}$ is $\prod^{n}_{i = 1} \prod^{n}_{\substack{j = 1; j \neq i}} \exp(-\lambda_{ij})(\lambda_{ij})^{y_{ij}}(y_{ij}!)^{-1}$ and therefore suggests an objective function under model (ref) as
Maximizer to $\ell_n (\boldsymbol{\beta})$, $\widehat \boldsymbol{\beta}_n$ is an estimator of $\boldsymbol{\beta}$. To facilitate our derivation, we introduce a few notations. For relational data with $n$ nodes, let $S_{1,n} := \{\{(i,j),(i,j)\}: \text{distinct~} i, j \in [n]\}$; $S_{2,n} := \{\{(i,j),(j,i)\}: \text{distinct~} i, j \in [n]\}$; $S_{3,n} := \{\{(i,j), (i,k)\}: \text{distinct~} i, j, k \in [n]\}$; $S_{4,n} := \{\{(i,j), (k,j)\}: \text{distinct~} i, j, k \in [n]\}$; and $S_{5,n} := \{\{(i,j), (k,i)\}: \text{distinct~} i, j, k \in [n]\} \cup \{\{(i,j), (j,k)\}: \text{distinct~} i, j, k \in [n]\}$. Now we impose the following conditions to draw inference on the coefficients.
Assumption (ref) establishes the regularity conditions for deriving the asymptotic behavior of $\widehat \boldsymbol{\beta}_n$. First, Conditions (ref) and (ref) impose mild assumptions on the gradient of the link function, extensively employed in Poisson model literature Davis2000Count,davis2016theory. Condition (ref) incorporates summability assumptions on the covariates and the integrability constraints on $g(\cdot)$, both of which are widely used in regression contexts phillips1999linear. It could be easily satisfied by linear functions and other functions discussed in Section (ref). Together, these three conditions ensure the consistency of $\widehat \boldsymbol{\beta}_n$ and regulate the limiting behavior of the asymptotic covariance of $\widehat \boldsymbol{\beta}_n$. Second, the asymptotic normality of $\widehat \boldsymbol{\beta}_n$ is guaranteed by conditions (ref) and (ref). Condition (ref) aligns with assumptions made by Lumley2003sparseCorrelations and Marrs2023exchangeableErrors. The moment condition in (ref) exhibits a higher degree of flexibility when compared to those found in high-dimensional generalized linear models tian2022transfer, which typically impose a light tail assumption on random noises. Indeed, there exists a broad spectrum of distributions for $\{e_{ij}\}$ satisfying Condition (ref), such as sub-Gaussian and sub-exponential families. Finally, Condition (ref) can be met by employing suitable encoding techniques to ensure that ${\mathbf X}$ remains within a compact domain and is trivially satisfied in the case of fixed design.
Now we are ready to formally provide the inference procedure of $\boldsymbol{\beta}$. The asymptotic property of $\widehat \boldsymbol{\beta}_n$ is summarized in Theorem (ref). Here we focus on the population version with respect to the true asymptotic covariance matrix.
Theorem (ref) paves a road for drawing inference on $\boldsymbol{\beta}$. As discussed in Section (ref), we focus on the setting with network dependence where at least one of $\{\eta_3,\eta_4,\eta_5\}$ is nonzero. Due to the complex network dependencies, the likelihood of $y_{ij}$ cannot be expressed as the sum of independent random variables. As a result, the standard central limit theorems and Lindeberg's condition are not directly applicable. To address this issue, one might consider leveraging U-statistics techniques to induce a summation of independent random variables graham2020dyadic,graham2020cid. However, such an approach requires additional assumptions on the distribution of $\{y_{ij}\}$. For example, graham2020dyadic assumes exchangeability of $\{y_{ij}\}$, and does not naturally accommodate edge covariates. In contrast, we assume exchangeability only on the error terms in the regression model. Moreover, graham2020dyadic assumes dependence between any pair $(y_{ij}, y_{kl})$ sharing at least one index, whereas our model does not impose this requirement. In fact, the unknown dependency structure can lead to indeterminate degenerate status of U-statistics under network setting du2024optimal. To overcome these challenges, we employ a key lemma from Bolthausen1982mixingCLT, which provides a sufficient condition for establishing the asymptotic normality of a sequence of measures based on the standard normal characteristic function. The proof of Theorem (ref) is given in the Supplementary Material.
The existence of $\mathbf{J}_n$ and $\mathbf{L}_n$ are guaranteed by conditions (ref) and (ref). By replacing $\boldsymbol{\beta}$ and $\boldsymbol{\eta}$ by $\widehat \boldsymbol{\beta}_n$ and $\widehat \boldsymbol{\eta}$ in $\mathbf{J}_n$ and $\mathbf{L}_n$ to get $\mathbf{\widehat J}_n$ and $\mathbf{\widehat L}_n$, the asymptotic covariance matrix can be estimated through $\mathbf{\widehat J}_n^{-1} \mathbf{\widehat L}_n \mathbf{\widehat J}_n^{-1}$. We will demonstrate its consistency in Section (ref). The following Example (ref) provides an illustration of Theorem (ref) with exponential link function, which yields a simple form of the asymptotic covariance matrix.
In this section, we construct a consistency estimator of the asymptotic covariance matrix in Theorem (ref). This involves the estimation of $\boldsymbol{\eta}$, which serves as the cornerstone for inferring $\boldsymbol{\beta}$. Recall that $\eta_1$ is the variance of error terms while $\eta_2$ through $\eta_5$ are the covariance terms in $\boldsymbol{\Omega}_e$, and $S_{1,n}$ through $S_{5,n}$ represent the sets of pairs of links corresponding to $\eta_1$ through $\eta_5$, as defined in Section (ref). First, we develop the estimation procedure of the covariance terms. By (ref), given the knowledge of $\{\xi_{ij}\}$, we consider the moment estimator of $\eta_2$ through $\eta_5$ as follows: \begingroup \allowdisplaybreaks
\endgroup for $\{(i,j),(k,i)\} \in S_{5,n}$ and $\{(i,j),(j,k)\} \in S_{5,n}$. Replacing $\{\xi_{ij}\}$ in $\overline \eta_2, \overline \eta_3, \overline \eta_4$, and $\overline \eta_5$ by $\{\widehat \xi_{ij}\}$ gives us $\widehat \eta_2, \widehat \eta_3, \widehat \eta_4$, and $\widehat \eta_5$.
For the variance term $\eta_1 = \text{Var}(\xi_{ij}) - \{g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\boldsymbol{\beta})\}^{-1}$, a natural estimator is $\overline \eta_{1,*} = |S_{1,n}|^{-1} \sum^n_{i\neq j} \xi^2_{ij} - |S_{1,n}|^{-2} (\sum^n_{i\neq j} \xi_{ij})^2 - |S_{1,n}|^{-1} \sum^n_{i\neq j} \{g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }} \boldsymbol{\beta})\}^{-1}$, given the knowledge of $\boldsymbol{\beta}$. In practice, replacing $\boldsymbol{\beta}$ by $\widehat \boldsymbol{\beta}_n$ obtained from (ref) gives:
However, the estimator in (ref) does not necessarily guarantee the positivity of the variance term, which may not lead to a legitimate $\widehat \boldsymbol{\Omega}_e$. To circumvent that difficulty, we refine the estimator using a hybrid procedure, which proceeds by first applying (ref), and then modify $\widehat \eta_{1,*}$ using a $k$-shorth estimator Pensia2019entangled. The $k$-shorth estimator outputs the center of the shortest interval containing at least $k$ points. While the traditional shorth estimator uses $k = N/2$ for sample size $N$, the estimator in Pensia2019entangled considered $k = c\log(N)$ for tunning parameter $c$, which provides the guaranteed optimality.
By (ref), we have $\eta_1 = \mathbb{E}(\xi^2_{ij}) - 1 - \{g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\boldsymbol{\beta})\}^{-1}$ for distinct $i, j \in [n]$. Let $\widehat \xi_{ij}= {y_{ij}} \{g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)\}^{-1}$ and $\widehat \zeta_{ij} = \widehat \xi^2_{ij} - 1 - \{g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)\}^{-1}$, then the estimator in (ref) can be rewritten as $(n^2 - n)^{-1} \sum^n_{i\neq j} \widehat \zeta_{ij}$. We consider the $n^2-n$ elements in $\{\widehat \zeta_{ij}\}$ to be the total points we apply the shorth estimation procedure on. Then we control the size of the $k$-shorth interval by letting $k = c\log(N)$ for $N=n^2 - n$, and take only the valid intervals with a positive center to constrain the estimator to give us a positive estimate of $\eta_1$, denoted as $\widehat \eta_{1, +}$. The proposed hybrid estimation procedure outputs the estimate in (ref) when it is greater than zero; otherwise, it outputs the positive $k$-shorth estimate. The hybrid estimator could be represented by $\widehat \eta_1 = \widehat \eta_{1,\mathrm{hybrid}} = \widehat \eta_{1, +} \cdot \mathbb{I}(\widehat \eta_{1,*} \leq 0) + \widehat \eta_{1,*}\cdot \mathbb{I}(\widehat \eta_{1,*} > 0).$
Let $\widehat \boldsymbol{\eta} = (\widehat \eta_1, \widehat \eta_2, \widehat \eta_3, \widehat \eta_4, \widehat \eta_5)^{{ \mathrm{\scriptscriptstyle T} }}$ denote the estimator of $\boldsymbol{\eta}$. Theorem (ref) below establishes the consistency of $\widehat \boldsymbol{\eta}$. When the number of nodes goes to infinity, $\widehat \boldsymbol{\eta}$ will fall in the parameter space discussed in Section (ref) by its consistency.
By Theorem (ref), we consider the covariance estimator of $\boldsymbol{\Omega}_0$ in Theorem (ref) takes the form $\widehat \boldsymbol{\Omega}_0 = \widehat{\text{Cov}}(y_{ij})$, where $\widehat{\text{Var}}(y_{ij}) = g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)^2 \widehat \eta_1 + g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)$; $\widehat{\text{Cov}}(y_{ij}, y_{ji}) = g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)g({\mathbf x}_{ji}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)\widehat \eta_2$; $\widehat{\text{Cov}}(y_{ij}, y_{il}) = g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)g({\mathbf x}_{il}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)\widehat \eta_3$; $\widehat{\text{Cov}}(y_{ij}, y_{kj}) = g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)g({\mathbf x}_{kj}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)\widehat \eta_4$; $\widehat{\text{Cov}}(y_{ij}, y_{jk}) = g({\mathbf x}_{ij}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)g({\mathbf x}_{jk}^{{ \mathrm{\scriptscriptstyle T} }}\widehat \boldsymbol{\beta}_n)\widehat \eta_5$, and ${\text{Cov}}(y_{ij}, y_{kl}) = 0$ by the dissociated assumption introduced in Section (ref). Note that to obtain the asymptotic properties of $\widehat \boldsymbol{\beta}_n$, $\boldsymbol{\Omega}_0$ needs to be invertible. The following proposition shows that the consistent estimators of $\eta_i$'s provide us a positive define $\widehat \boldsymbol{\Omega}_0$ as the estimator of covariance matrix of ${\mathbf Y}$.
Finally, we show the asymptotic covariance matrix in Theorem (ref) can be consistently estimated in practice.
Combining Theorem (ref) and Theorem (ref) leads to a formal inference procedure of $\boldsymbol{\beta}$. For example, for each $\ell \in [p]$, denote $\widehat \sigma^2_\ell$ the $\ell$th diagonal entry of $n^{-1}\mathbf{\widehat J}_n^{-1} \mathbf{\widehat L}_n \mathbf{\widehat J}_n^{-1}$, a $100(1-\alpha)\%$ confidence interval for the $\ell$th entry of $\boldsymbol{\beta}$, denoted by $\boldsymbol{\beta}_\ell$, is given by $\big[\widehat \boldsymbol{\beta}_\ell - \widehat \sigma_\ell \boldsymbol{\Phi}^{-1}(1-\alpha/2), \widehat \boldsymbol{\beta}_\ell + \widehat \sigma_\ell \boldsymbol{\Phi}^{-1}(1-\alpha/2)\big]$, where $\boldsymbol{\Phi}(\cdot)$ is the cumulative distribution function of standard normal distribution.
In this section, we evaluate the numerical performance of the proposed method for count relational data and compare it with other benchmarks. We illustrate the validity of our inference framework under the weakly exchangeable error setting in Section (ref). Since few models have been studied for count relational data and even fewer for the dependence structure introduced in this paper, we examine the performance of 95% confidence interval coverage probability among the three methods: (1) {\tt Our model}: the inference procedure proposed in our work; (2) {\tt Naive}: the inference procedure assuming no edge dependencies; (3) {\tt Oracle}: the inference procedure given the true covariance matrix of error terms. The oracle result with known value of $\boldsymbol{\eta}$ serves as a benchmark, while the naive method assumes observations ${\mathbf Y}$ are marginally independent. The results of the naive method will demonstrate the necessity of involving the dependency structure of relational data in the inference procedure. The following model is employed to generate count relational data:
We vary the number of nodes $n \in \{20 ,50, 100, 150\}$, fix $\boldsymbol{\beta} = (1, -0.5, -0.5, -1)^{{ \mathrm{\scriptscriptstyle T} }}$, and independently draw $x_{1ij}\sim N(2,1)$, $x_{2i}\sim \mathrm{Bernoulli}(0.6)$, $x_{3i}\sim N(1,1)$, and $x_{4ij}\sim N(1,1)$. We apply settings in Example (ref) to generate $\{e_{ij}\}$. Under each realization of ${\mathbf X}$, we simulate 1,000 error terms to calculate the empirical coverage probability, and repeat the experiment 15 times to evaluate the variation in the 95% confidence interval coverage of the three competing methods.
By construction, we have $\boldsymbol{\eta} = (1.1, 0, 0, 0, 0)$ under setting (ref), and $\boldsymbol{\eta} = (13, 2, 7, 2, 0.4) \times 10^{-2}$ under setting (ref). Although setting (ref) falls outside the scope of our primary interest due to the absence of edge dependencies, we can still estimate the asymptotic covariance matrix of $\widehat \boldsymbol{\beta}$ by $n^{-1}\mathbf{\widehat J}_n^{-1} \mathbf{\widehat L}_n \mathbf{\widehat J}_n^{-1}$. This serves as a configuration which admits the assumption of independence among edges in the naive approach. It is worth noting that the error generating procedures outlined in Example (ref) naturally ensure that $\boldsymbol{\eta}$ satisfies the constrains in Section (ref) for weakly exchangeable errors, thereby defining legitimate covariance matrices of errors on the parameter space.
In practice, to get the hybrid shorth estimate of $\eta_1$ from $\{\widehat \zeta_{ij}\}$, we apply cross validation to tune the parameter $c$ Pensia2019entangled as introduced in Section (ref). We set the possible range of $c$ to span from $2/\log(n^2-n)$ to $\sum_{i\neq j}^n\mathbb{I}\big[\widehat\zeta_{ij} > - \max(\{\widehat\zeta_{ij}\})\big]/\log(n^2-n)$ and denote the set of tunning parameters by $\mathcal{S}$. For each $c^* \in \mathcal{S}$, we randomly divide $\{\widehat \zeta_{ij}\}$ into $10$ folds of approximately equal size. After selecting a validation set, we apply $k$-shorth method on the remaining $9$ folds. The mean squared error, $\text{MSE}_\ell$, $\ell = 1,2,\ldots,10$, is computed using the observations in the held-out fold and the $k$-shorth estimate. The positive $k$-shorth estimator $\eta_{1, +}$ is calculated by setting $c = \operatorname*{arg\,min}_{c^*} \{c^* \in \mathcal{S}: \frac{1}{10}\sum^{10}_{\ell=1}\text{MSE}_\ell\}$. In practice, we enforce the positive semi-definiteness of $\widehat \Omega_e$ by applying an eigenvalue correction to $\widehat \eta_{1,\mathrm{hybrid}}$. Specifically, we adjust the smallest eigenvalue of $\widehat \boldsymbol{\Omega}_e$ to ensure it is nonnegative. Our numerical experiments indicate that such a minor perturbation has a negligible impact on the computational accuracy of the final results.
In this section, we present simulation results for inferring regression coefficients applying settings in Example (ref). Additional studies with $\{e_{ij}\}$ from Gamma distribution and experiments with different covariate configurations are given in the Supplementary Material.
As shown in Figure (ref), when error terms are independent and identically distributed, our method performs as good as the oracle results. The coverage probability of the naive method, however, is further from the nominal 95% level, and its variability across different realizations is larger than that of our method, especially when the number of nodes is close to or less than 50. For dependent error terms generated from setting (ref), our method performs extremely well as it recovers the dependence structure in the relational data. Specifically, our proposed inference procedure produces confidence intervals with coverage probability close to the nominal 95% level under all configurations. Its performance becomes better and closer to the oracle results as the size of relational data grows, whereas the coverage probability of naive method is far below the nominal level and becomes worse as the number of nodes increases.
In conclusion, our method outperforms the naive method across all settings, especially for weakly exchangeable dependent errors. Specifically, the empirical coverage probability of our method is approaching the nominal level and closely approximates the oracle benchmark as the number of nodes increases. Moreover, our method demonstrates robustness in terms of empirical coverage probability under different error generating procedure (heavy-tailed errors from Gamma distribution as well as light-tailed errors from truncated Normal distribution), and different configurations under model (ref). Comprehensive simulation results further supporting these findings are relegated to the Supplementary Material.
In this section, we apply the proposed model to investigate the food sharing network introduced in Section (ref). The data contains the number of transferred gifts over a yearlong period among 25 households of indigenous Mayangna and Miskito horticulturalists in Nicaragua, along with distance, relationship, and other covariates given in the Supplementary Material. The “association index" cairns1987comparison reflects the amount of time that households interact with one another, which characterizes the multi-faceted inter-household relationships. A complication which arises in the study is that not all households were present for the full duration of the yearlong study. For model interpretation, koster2014food accounts for the variation in the proportion of the year for which both members of each dyad were simultaneously present in the community by entering the natural logarithm of this exposure as an offset variable. This modification allows us to model the expected number of gifts per year. The social relations model (SRM) developed by kenny1984srm is applied in koster2014food to separate individual effects in the log mean structure from relationship effects in dyadic data. Their overall results indicate that food sharing networks largely correspond to kin-based networks of social interaction, suggesting that food sharing is embedded in broader social relationships between households.
Note that the marginal mean of $y_{ij}$ in our model does not depend on the network effects defined in Section (ref). The reciprocity effect, same sender effect, same receiver effect, and sender-receiver effect defined in our model are characterized by $\boldsymbol{\eta}$, which are different from the random effects defined in SRM. The analysis in koster2014food assumes random effects are normally distributed, but the authors find a noteworthy outlier (the number of gifts given between Household 1 and Household 25) in the relationship-level random effects and therefore include a dummy variable to represent this outlier as a fixed effect. Our model, however, benefits from the flexibility where we do not introduce strict distributional assumption on error terms except for weak exchangeability. Therefore, we exclude the artificial relationship covariate between Household 1 and Household 25 as introduced in koster2014food.
We choose the exponential link function and put all variables in model (ref) for preliminary analysis. it is worth noting that household dyads which spend considerable time together typically have close kinship ties hames1987garden, alvard2009kinship, koster2014food. This also agrees with the correlation of estimated coefficients given in the Supplementary Material, where we find strong correlation between the effect of association index and mother-offspring ties. Conceptually, association index is a proximity measure of close kin ties. To deal with the collinearity problem, we consider the model where mother-offspring, father-offspring, full sibling, or other close kin ties are omitted, whereas the association index ($\text{Association}_{ij}$) is reserved. We further consider a dummy variable $\text{Relatedness5}_{ij} = 1-\text{Relatedness1}_{ij}-\text{Relatedness2}_{ij}-\text{Relatedness3}_{ij}$, to denote weak ties as discussed in koster2014food, given the definition of those attributes in the Supplementary Material. Note that fishing is a common strategy for virtually all households koster2014food, and as mentioned in defrance2009zooarchaeology, meat circulated as a source of wealth and people generated wealth from the products that animals produced. These findings suggest a potential colinearity between meat harvesting and the wealth of a household. Together with the observations according to results guven in the Supplementary Material, we omit “Wealth" in the nodal covariates since there is notable correlation of its estimated coefficients with that of others (“Fish" and “Pig"). Lastly, we omit “Pastors" variable in our model since there is only 2 households with pastors among the 25. The structured sparsity introduced by it would make the inference procedure unstable.
We choose the exponential link function and put all variables in model (ref) for preliminary analysis. Pre-processing of the full sample is detailed in the Supplementary Material. Our final model takes the following form:
where $\{e_{ij}\}$ are weakly exchangeable with unit mean. In summary, after accounting for network effects, our model detects a statistically significant giver-game, giver-pigs, receiver-game, receiver-fish, weak kinship, distance, and association effects, but finds no evidence of effects for giver-game or receiver-pigs. It further suggests strong reciprocity effect and notable same sender/receiver effect in the relational data. Figure (ref) shows the estimation results, where $\widehat\boldsymbol{\eta} = (0.829, 0.427, 0.093, 0.111, 0.011)$.
Since the marginal mean of $y_{ij}$ in our model is different from that of SRM applied in koster2014food, making it difficult to compare the regression coefficients of our model and theirs. Nonetheless, we can compare the significance of estimated coefficients of our model to the existing results. We can also verify whether the sign of $\eta_2$ through $\eta_5$ agrees with the result obtained in koster2014food since the network effects can be presented by variance/covariance parameters in SRM model. Applying our definition of $\boldsymbol{\eta}$ to SRM gives $\eta_2 = 0.643$; $\eta_3 = 0.756$; $\eta_4 = 0.426$; and $\eta_5 = 0.041$, where the variance/covariance parameters are defined in koster2014food. Though not directly comparable, the positiveness of these four network effects defined in our model agree with the results in koster2014food.
{\it 6.3.1. Significant effects of distance and association index between households, together with giving and receiving behaviors.} As a result, our model finds 8 significant coefficients among the 10. The intercept is estimated to be 0.853, which is significantly different from 0, indicating the willingness of gift giving between the Households. Households who harvest more game ($\widehat\beta_1=0.357$) and own more pigs ($\widehat\beta_3=0.128$) are predicted to give significantly more gifts than households who harvest less game and own less pigs. However, there is no significant association between the amount of fish ($\widehat\beta_2=-0.974$) that households harvest and the number of gifts they tended to give to other households, which partly covariates to the relatively small proportion of fish that are sent as gifts. Additionally, the association between the amount of pigs ($\widehat\beta_6=-0.033$) the households own and gift receiving is not significant, having adjusted for the other factors in the model. Households located farther apart are predicted to exchange less gift ($\widehat\beta_8=-0.384$) than nearby households. Distance is entered as a log-transformed variable and so its coefficient has a partial elasticity interpretation: a 10% increase in the distance between two households is associated with a 3.8% decrease in the expected number of gifts exchanged between the two households. Regarding the association index, households who associated more frequently with one another were predicted to give more ($\widehat\beta_9=2.444$). Those results are similar to what is observed in koster2014food. Moreover, our model predicts less transfers between households with weaker ties ($\widehat\beta_7=-0.996$), after omitting the dummy variables denoting mother-offspring, father-offspring, full sibling, or other close kin ties. Similar observation on kin ties and food sharing could be found in helms1971asang and parsons1974between.
{\it 6.3.2. Different findings than prior research concerning the significant impact of receiving behavior on daily harvest of meat and fish.} As for the difference in results, our model suggests that households who harvest less game ($\widehat\beta_4=-0.175$) and less fish ($\widehat\beta_5=-1.020$) receive significantly more gifts than households who harvest more game and fish, while koster2014food finds no significant association between the amount of game/fish and the number of gifts they tended to receive from other households. This could result from that we omit “Wealth" in our analysis due to collinearity while koster2014food include this variable in their model.
{\it 6.3.3. Implication of strong reciprocity effect and notable same sender/receiver effect.} Now, we analyze the dependence structure in the food sharing network. As depicted in Figure (ref), the reciprocity effect dominates the same sender effect, same receiver effect, and sender-receiver effect. Specifically, it suggests the gift giving behavior is reciprocated. Though relatively smaller compared to $\eta_2$, the magnitude of $\eta_3$ and $\eta_4$ are larger than $\eta_5$, indicating notable dependencies between relations involving the same sender or receiver.
In this paper, we proposed a flexible multiplicative model on count edges for relational data. The model can handle different count distributions and is able to capture the underline pairwise dependence structure between edges given the observed data. For the regression setting, we proved that the proposed estimator is asymptotically normal and the estimate of covariance parameters is consistent, which delivers valid inference under the weak exchangeability assumption. Our work makes important progress toward the inference problem for modeling count edges in relational data. We also demonstrated the proposed model on a food sharing network.
Note that the estimation procedure of coefficients is robust to model misspecification if data comes from a linear exponential family with the same mean structure as in our model. The consistency of the estimator is guaranteed by the limit theory for the statistical agnostic, as discussed in Gourieroux1984pseudoMLE. Under the model assumption of weakly exchangeable errors in Section (ref), the asymptotic variance of $\widehat \boldsymbol{\beta}_n$ is dominated by $\eta_{3}, \eta_{4}$, and $\eta_{5}$, since in $\boldsymbol{\Omega}_e$, elements of $\boldsymbol{\eta}$ occur with multiplicity $n^2 - n$ (for both $\eta_{1}$ and $\eta_{2}$), $n^3 - 3n^2 + 2n$ (for both $\eta_{3}$ and $\eta_{4}$), and $2(n^3 - 3n^2 + 2n)$ for $\eta_{5}$. The variance-covariance estimation for Theorem (ref) naturally involves $\widehat\eta_1$ and $\widehat\eta_2$ which accounts for the asymptotically negligible bias mentioned in graham2020dyadic, though under a different modeling framework. As mentioned in graham2020cid, the non-zero covariance terms are crucial for understanding the sampling distribution of $\widehat \boldsymbol{\beta}_n$. It is not of our interest in this work when the same sender effect, same receiver effect, and sender-receiver effect do not exist, though it gives a faster convergence rate of $\widehat \boldsymbol{\beta}_n$. Under this scenario, the asymptotic normality of $\widehat \boldsymbol{\beta}_n$ does not necessarily hold except for i.i.d. errors, as less correlation does not imply less dependencies in the model. Related discussions are available in menzel2021bootstrap.
Our work can be extended in several directions. One is adapting the model to relational data where the pairwise edge dependencies can vary with the sample size, potentially affecting the convergence rate of model estimators. Another promising direction is to study the dependence structure within the relational data, offering insights for hypothesis testing to compare dependencies across different datasets.
The complete proofs, additional simulation and empirical illustration results, and further discussions on the parameter space of covariance of weakly exchangeable variables can be found in the Supplementary Material. Code for reproducing the experiments and real data analysis can be obtained from \url{https://github.com/WenqinDu/Count_Relational_Data_Modeling}.