# A Generalized Estimating Equation Approach to Network Regression

Riddhi Pratim Ghosh  Department of Mathematics and Statistics, Bowling Green
State University  Jukka-Pekka Onnela  Department of Biostatistics, Harvard
University  Ian Barnett  Department of Biostatistics, University of
Pennsylvania

###### Abstract

Regression models applied to network data where node attributes are the
dependent variables poses a methodological challenge. As has been well
studied, naive regression neither properly accounts for community structure,
nor does it account for the dependent variable acting as both model outcome
and covariate. To address this methodological gap, we propose a network
regression model motivated by the important observation that controlling for
community structure can, when a network is modular, significantly account for
meaningful correlation between observations induced by network connections. We
propose a generalized estimating equation (GEE) approach to learn model
parameters based on clusters defined through any single-membership community
detection algorithm applied to the observed network. We provide a necessary
condition on the network size and edge formation probabilities to establish
the asymptotic normality of the model parameters under the assumption that the
graph structure is a stochastic block model. We evaluate the performance of
our approach through simulations and apply it to estimate the joint impact of
baseline covariates and network effects on COVID-19 incidence rate among
countries connected by a network of commercial airline traffic. We find that
during the beginning of the pandemic the network effect has some influence,
the percentage of urban population has more influence on the incidence rate
compared to the network effect after the travel ban was in effect.

Keywords: network regression; transportation networks; generalized estimating
equations; COVID-19

##  1 Introduction

Network data provide quantitative information to study and unveil the pattern
of interactions among various objects from individuals to countries. The
inferential task in scientific applications requires learning about the
influence that an individual has on others, which is further complicated as
the individuals are often nested in communities. Two important methodological
challenges of network regression are (1) accounting for correlation induced by
network structure, and (2) allowing for node attributes to appear in the data
both as outcome as well as covariates in the design matrix. For example,
network regression was applied in the study of the nature and extent of the
person-to-person spread of obesity, Christakis Fowler (2007) found that
obesity appears to spread through social ties, a phenomenon that may in part
be explained by the notion of homophily – that birds of a feather flock
together (Shrum ., 1988; Igarashi ., 2005; Tifferet, 2019). Sensitivity
analyses suggest that contagion effects for obesity and smoking cessation are
reasonably robust to possible latent homophily or environmental confounding;
those for happiness and loneliness are somewhat less so (VanderWeele, 2011).
To further investigate the causal relationship, Shalizi Thomas (2011) provided
three factors underlying such interaction: homophily, or the formation of
social ties due to matching individual traits; social contagion; the causal
effect of an individual’s covariates on his or her measurable responses; and
an individual’s response. This has led to the development of two types of
models: (1) community detection models that aim to find distinct communities
or clusters of similar individuals (Holland ., 1983; Newman, 2018), and (2) a
more general framework of network regression that can model such phenomenon
using a direct link between observed individual attributes or covariates and
the network interactions (Holland Leinhardt, 1981; PD. Hoff ., 2002; P. Hoff,
2021). Some approaches that aim at merging the above two models exploit
distributional assumptions to incorporate covariate information in detecting
communities in the network (Binkiewicz ., 2017; Mu ., 2022). These examples
highlight the need for regression techniques that are robust to non-trivial
network structure because in all these studies community labels are known
which makes the evaluation of covariate effects easy. However, in a more
realistic situation, community labels are unknown (unobserved) with
differences in covariates existing through the unobserved communities making
the inference challenging.

In this article, we aim to develop a network regression model motivated by the
important observation that the differences between the communities in a
network can be attributed to the differences in the influences of covariates
responsible for an edge formation. Our model is based on a flexible network
effect assumption, and it allows us to perform tests for the model parameters
in addition to the estimation. Aside from detecting communities solely based
on degree, the degree corrected stochastic block model was one of the first
model that allows within community heterogeniety (Qin Rohe, 2013). There has
been a surge of interests among psychologists and social scientists to study
such behavior- Aukett . (1988) shows that gender difference plays an important
role in friendship patterns: women shows a preference for a few, closer,
intimate same-sex friendships based on sharing emotions whereas men build up
friendship based on the activities they do together. Staber (1993) studies how
women and men form entrepreneurial relationships, concluding that women’s
networks are wider with more strangers and higher proportion of cross-sex
ties. In each of these network regression examples, community labels are known
which alleviates the inferential task of finding the effects of covariates and
learning community structures. However, in many other scenarios, community
labels are unknown which makes the inferential task more challenging. For
example, here we will consider the network regression problem of seeing the
impact of air travel network flow between countries on COVID-19 incidence
rates, where network structure is unknown and must be estimated. The current
literature lacks powerful methods that can address this important issue. Our
focus is to bridge this gap by leveraging covariate-assisted information to
propose an effective tool which also accounts for the network modeling through
adjacency matrix.

Generalized estimating equations (GEE) (Liang Zeger, 1986; Zeger, 1986)
provide a popular approach that is often used to analyze longitudinal and
other types of correlated data (Burton ., 1998; Diggle ., 2002; Mandel .,
2021) We propose first performing community detection in order to estimate
community membership before then using a GEE regression model to account for
the resulting estimated community structure. This approach is general and
agnostic to the community detection algorithm used so long as each node can
belong to at most one community. Given that network regression needs to
account for correlation between attributes from nodes that are connected, a
GEE allows for arbitrary correlation between nodes within the same community.
Of note, this approach is best suited for highly modular networks so that the
GEE assumption of independence between different communities is more accurate,
though we explore its performance in less modular networks with more between-
community mixing.

The rest of this article is organized as follows. In Section 2 we introduce
our network regression model along with our GEE extension, and we describe the
transportation network we use to model country-specific COVID-19 incidence
rates. In Section 3 we present the theoretical results followed by extensive
simulation results and real data analysis. Finally, Section 4 concludes the
article with a discussion.

##  2 Methods

In this section, we introduce our network regression model and describe the
air travel newtorks for each month of the first quarter of 2020 (the beginning
of the COVID-19 pandemic) with country-specific information such as COVID-19
incidence rate, GDP, and population size. Next, we present a generalized
estimating equation (GEE) approach to network regression that accounts for
network-induced correlation between observations.

###  2.1 Model and notation

We consider a directed network of $n$ nodes and an $n\times n$ adjacency
matrix $\bm{A}=(a_{ij})$ with $0$’s on the diagonal (i.e. no self-loops). We
denote the feature variable of the $i$th node by $y_{i}$, the $l$-dimensional
vector of covariates by $\bm{\alpha}$, and the corresponding design matrix by
$\bm{x}_{i}$. Denoting the coefficient of interest by $\beta$, the network
regression model for node attribute $y_{i}$ we consider:

$\displaystyle y_{i}=\bm{\alpha}^{\top}\bm{x}_{i}+\beta\cdot\sum_{j\neq
i}A_{ji}y_{j}/(n-1)+\epsilon_{i},$ (1)

where $\epsilon_{i}\overset{iid}{\sim}N(0,\sigma^{2})$ and
$\bm{\alpha}\in\mathbb{R}^{l}$. One can note that the above model is
reminiscent of the first-order autoregressive spatial model of Kelejian Prucha
(1998) which frequently contains a spatial lag of the dependent variable as a
covariate that is spatially autoregressive.

Using vector and matrix notation, the above model can also be written as

$\displaystyle\bm{y}=\bm{X}^{\top}\bm{\alpha}+\beta\cdot\bm{A}\bm{y}/(n-1)+\bm{\epsilon},$

where $\bm{y}$ is the concatenation of the $y_{i}$s in a vector of length $n$,
and $\bm{X}=[\bm{x}_{1}:\bm{x}_{2}:...:\bm{x}_{n}]$, is the matrix of
covariates. Note, model (1) can be trivially adapted to generalized linear
models with non-linear link functions and non-Gaussian data in the scope of
GEE models.

###  2.2 Data description

With 622 million confirmed cases and 6.5 million deaths globally as of October
01, 2022, the COVID-19 pandemic has had a tremendous impact on the world,
shrinking the global economy by 5.2%, the largest recession in the history
post World War II (Bank, 2020). The travel bans in places worldwide have
severely affected the tourism industry, with estimated losses of 900 billion
to 1.2 trillion USD and tourism down 58%-78% (Le ., 2022). The airline
industry has also suffered heavily, with 43 airlines declaring bankruptcy and
193 of 740 European airlines at risk of closing. Here we focus on the start of
the pandemic covering the transition to travel bans across the world to study
the relative effectiveness of travel bans for controlling and contributing to
COVID-19 incidence rates.

We use pandemic data from the Johns Hopkins University coronavirus data
repository through April 30, 2020, (COVID, 19). Flight data are from the
Official Airline Guide (OAG) (Strohmeier ., 2021). Because only data for
January and February 2020 are available from OAG, we used the estimated fight
data for other time periods using the OpenSky Network database (Schäfer .,
2014; Strohmeier ., 2021). This database tracks the number of fights from one
country to another over time, which we use to estimate country-to-country
flight data for other months. We include as covariates the GDP, total
population, and percentage of the urban population for each country in our
network. Constructing a network based on by the flight data and incorporating
the above country-specific attributes such as GDP, population, etc. as
covariates through $\bm{\alpha}$ in (1), we aim to estimate the effectiveness
of travel bans through $\beta$ in model (1).

###  2.3 Generalized estimation equation (GEE) approach

Out network contains $K$ communities where community $k$ is defined by

$E_{k}=\\{i:g_{i}=k\\},$

where index $g_{i}$ represents the community membership of node $i$,
$|E_{k}|=n_{k}$ (number of nodes in community $k$) and $\sum n_{k}=n$. Let
$\bm{y}_{k}$ and $\bm{\epsilon}_{k}$ denote the concatenated vectors of
$y_{ij}$s and $\epsilon_{ij}$s respectively, and
$\bm{X}_{k}=[\bm{x}_{1},\bm{x}_{2},...,\bm{x}_{n_{k}}]$ denote the $l\times
n_{k}$ sub-matrix of covariates corresponding to cluster $k$.

To fit our network regression model in the GEE framework, we use communities
as clusters and use the following equation to model the node attributes of
members of the $k$th cluster:

$\displaystyle\bm{y}_{k}=\bm{X}_{k}^{\top}\bm{\alpha}+\beta\bm{Z}_{k}+\bm{\epsilon}_{k},$
(2)

where $\bm{A}_{k}$ is the $n_{k}\times n_{k}$ sub-matrix of $\bm{A}$
pertaining to the cluster $k$, and $\bm{Z}_{k}=\bm{A}_{k}\bm{y}_{k}/(n-1)$.

The marginal mean $\bm{\mu}_{k}$ of $\bm{y}_{k}$ has the form:

$\bm{\mu}_{k}=E(\bm{y}_{k}|\bm{X}_{k})=(\bm{I}_{n_{k}}-\beta\bm{A}_{k}/(n-1))^{-1}\bm{X}_{k}^{\top}\bm{\alpha},\;\;\;\
k=1,2,...,K,$

where $\bm{I}_{n_{k}}$ is the identity matrix of order $n_{k}$.

Adopting a GEE approach (Liang Zeger, 1986), the resulting estimating equation
is given by

$\displaystyle\sum_{k=1}^{K}\bm{D}_{k}^{\top}\bm{V}_{k}^{-1}(\bm{y}_{k}-\bm{\mu}_{k})=\bm{0},\;\;\
k=1,2,...,K,$ (3)

where
$\bm{D}_{k}=\frac{\partial\bm{\mu}_{k}}{\partial(\beta,\bm{\alpha})^{\top}}$
is of dimension $n_{k}\times(l+1)$ and $\bm{V}_{k}$ is the $n_{k}\times n_{k}$
working covariance matrix of $\bm{y}_{k}$. The explicit form for $\bm{D}_{k}$
is

$\bm{D}_{k}=[-\underbrace{(\bm{I}_{n_{k}}-(\beta/(n-1))\bm{A}_{k})^{-1}\bm{A}_{k}(\bm{I}_{n_{k}}-(\beta/(n-1))\bm{A}_{k})^{-1}\bm{X}^{\top}\bm{\alpha}}_{n_{k}\times
1}\;\;\ :\;\;\
\underbrace{(\bm{I}_{n_{k}}-(\beta/(n-1))\bm{A}_{k})^{-1}\bm{X}_{k}^{\top}}_{n_{k}\times
l}].$

One can note that $\bm{D}_{k}$ consists of two partitioned matrices where the
first one corresponds to the network parameter $\beta$ and the second one is
due to the covariate $\bm{\alpha}$. We can solve for $\hat{\bm{\alpha}}$ and
$\hat{\beta}$ in equation (3) through iterative reweighted least squares, and
use the robust sandwich covariance estimator to perform inference on
$\bm{\alpha}$ and $\beta$.

##  3 Results

###  3.1 Theoretical results

In this section, we prove the asymptotic normality of the resulting GEE
estimator for $\beta$ and $\bm{\alpha}$ jointly. Towards this, we assume
constant probabilities of edge formation between and within communities where
we denote these quantities by $p$ and $q$ respectively as in a stochastic
block model (Holland ., 1983). Our proof of asymptotic normality hinges on
Theorem 2 of Liang Zeger (1986) which establishes the asymptotic normality of
the regression parameter in the classical GEE approach under the assumption
that the correlation parameter appropriately scaled by the number of
communitites is consistent. Our primary distinction from this approach is that
we must account for probabilities of edge formation instead of correlation
between observations. Therefore, we first establish the consistency of $p$ and
$q$ following Proposition 1 of Chen . (2021) and subsequently show asymptotic
normality of $(\hat{\beta},\hat{\bm{\alpha}})$ from the estimating equation
(3).

####  3.1.1 Consistency of $\hat{p}$ and $\hat{q}$

######  Proposition 3.1.

Consider a network generated from a stochastic block model of $K$ communities
of size $m$ so that total number of nodes is $n=Km$. Assume that
$m^{\gamma}p\rightarrow p^{*}$ and $m^{\gamma}q\rightarrow q^{*}$ as
$n\rightarrow\infty$, where $p^{*}$ and $q^{*}$ are positive fixed constants,
and $\gamma\in[0,2)$, then $m^{1+\gamma/2}(\hat{p}-p)$ and
$m^{1+\gamma/2}(\hat{q}-q)$ are both $o_{P}(1)$.

###### Proof.

See the appendix in Section 5.1. ∎

Guided by the above proposition, we establish the asymptotic normality of the
model parameters in the following theorem. The key difference with the
classical formula is the inclusion of the community size in the covariance
coming from the consistency of $p$ and $q$.

####  3.1.2 Asymptotic normality of $\hat{\beta}$

######  Theorem 3.2.

Under the conditions of Proposition 3.1
$K^{1/2}m^{1+\gamma/2}\big{(}(\hat{\beta},\hat{\bm{\alpha}})-(\beta,\bm{\alpha})\big{)}^{\top}$
is asymptotically multivariate normal with zero mean and covariance given by

$V=\underset{K\rightarrow\infty}{\text{lim}}Km^{2+\gamma}\Big{(}\sum_{k=1}^{K}\bm{D}^{\top}_{k}\bm{V}_{k}^{-1}\bm{D}_{k}\Big{)}^{-1}\Big{(}\sum_{k=1}^{K}\bm{D}^{\top}_{k}\bm{V}_{k}^{-1}cov(\bm{Y}_{k})\bm{V}_{k}^{-1}\bm{D}_{k}\Big{)}\Big{(}\sum_{k=1}^{K}\bm{D}^{\top}_{k}\bm{V}_{k}^{-1}\bm{D}_{k}\Big{)}^{-1}$

.

###### Proof.

See the appendix in Section 5.2. ∎

Remark: The variance formula involves the term $m^{2+\gamma}$, where $m$ is
the community size. If $m$ is large, then the covariance will increase at a
rate $m^{2+\gamma}$. This is reminiscent of the well-known fact that sandwich
estimator $\hat{V}$ of the covariance of $(\beta,\bm{\alpha})^{\top}$ is not
stable if $m$ is large relative to $K$. In essence, if $m$ grows at a similar
rate to $K$ the sandwich estimator becomes and unstable estimator of
covariance. In practice this implies that the GEE approach works best when the
network contains many smaller communities rather than only a few larger
communities.

###  3.2 Simulations

We simulate a network of $n$ nodes having balanced communities of size $m=10$
via a stochastic block model and vary $n$ in $\\{200,400\\}$ with the number
of communities ($K$) being $20$ and $40$ for both values of $n$, respectively.
Let $p$ and $q$ denote the within community and between community
probabilities of an edge formation as in a stochastic block model. In each
setting, we vary $(p,q)\in\\{(0.8,0),(0.7,0.1),(0.6,0.2),(0.5,0.3)\\}$.

####  3.2.1 Estimation of $\beta$ and $\alpha$

The choice of true model parameters and the corresponding data-generating
process are detailed to span networks with varying degrees of modularity. We
set $\beta_{0}=0.5$ and
$\bm{\alpha}_{0}=(1,1,1,1,0.5,0.5,0.5,-0.5,-0.5,2)^{\top}$ ($l=10$ in equation
(1)). We simulate the $l\times n$ matrix $\bm{X}$ by a multivariate normal
distribution in the following manner. The $j$th column of $\bm{X}$, if it
belongs to community $k$ ($k=1,2,...,K$), follows
MVN$((k/10)\bm{1}_{l},0.0001\bm{I}_{l})$, where $\bm{1}_{l}$ and $\bm{I}_{l}$
are the vector with 1s and identity matrix of dimension $l$, respectively. In
each setting, the adjacency matrix $\bm{A}$ of dimension $n$ is simulated by
the stochastic block model which has $K$ communities of the aforementioned
size such that within and between community edge probabilities are $p$ and
$q$, respectively. Finally, the response variable $y_{i}$ is generated
according to equation (1) with $\sigma$ = 0.01.

To estimate $\beta$ and $\bm{\alpha}$, we first perform a community detection
on the directed graph obtained from the adjacency matrix $\bm{A}$ as in
Rosvall Bergstrom (2007). Next, with the resulting communities we fit a GEE
using the geepack R package (Halekoh ., 2006) and report the estimates of bias
and variance by taking the average of over $B=1000$ replications.

The squared bias and variance of the estimated $\beta$ increase as the
networks become less modular (i.e. $q$ increases) for both our GEE approach as
well as naive least squares (see Figure 1). The naive least squares method
does not assume any community structure and reports the parameter estimates
from assuming independence between observations using equation (1) directly.

From Figure 1 we see that less modular networks are more difficult to fit in
general based on the increase in bias for both methods. Despite this our GEE
approach uniformly is less biased that naive least squares, which demonstrates
that controlling for correlation within communities is effectively done by
GEE. This also exposes the weakness of our GEE approach: if the network is not
modular with high degrees of mixing between communities (i.e. large $q$) the
GEE framework cannot accommodate this due its assumption of independence
between communities. However, even when $q$ is large we still observe that GEE
somewhat mitigates the impact of correlation induced by the network at least
partially, which explains why its bias is smaller than that of naive least
squares. Essentially, even when the network structure suggests the GEE
assumptions are incorrect, one is still generally better off partially
accounting for network structure using our GEE approach rather than ignore
community structure altogether.

The smaller standard errors demonstrated by least squares in Figure 1 reflect
the inaccuracy of the method. By ignoring community structure and assuming
independence between observations in network regression settings, we expect to
see anti-conservative hypothesis tests, too-narrow confidence intervals, and
general overconfidence. This overconfidence is reflected in the too-small
standard errors for naive least squares as well as in the anti-conservative
Type I errors we demonstrate next.

Figure 1: Bias squared and standard error of estimates of $\beta$ with varying
degrees of network modularity. In order to estimate $\beta$ $B=1000$
replications were performed for $n=200,400$ with the average degree 7 and 16
respectively. The values of $(p,q)$ are varied over
$\\{(.8,0),(.7,.1),(.6,.2),(.5,.3)\\}$.

####  3.2.2 Hypothesis testing of $\beta$

We consider an approach to testingthe hypothesis of $H_{0}:\beta=0$ against
the alternative $H_{A}:\beta\neq 0$. We perform a simulation study to obtain
Type I error in a variety of network structures. We consider two working
correlation structures for our GEE model: independence and exchangeable.
First, we simulate our data as in Section 3.2 with $\beta=0$. We obtain an
empirical null distribution by replicating this procedure $B=1000$ times to
obtain $\hat{\beta}^{(1)},\hat{\beta}^{(2)},...,\hat{\beta}^{(B)}$. We
estimate the P-value by $\sum_{b=1}^{B}I(|\hat{\beta}^{(b)}|<\hat{\beta})/B$.
We summarize our Type I error results at the $0.05$ significance level in
Table 1 which demonstrates that as network modularity decreases ($q$ gets
larger), the hypothesis test for $H_{0}$ becomes anti-conservative. This is
true for both least squares as well as GEE, although least squares is more
adversely impacted. In addition, while GEE excels in the highly modular
setting, as expected, LS is still anti-conservative. This demonstrates the
efficacy of GEE as a better inferential tool over naive least squares in
network regression settings, even when GEE assumptions of between-community
independence do not hold. It is also instructive to note that although the
P-values for independent correlation structure are generally less than those
for exchangeable structure, the differences are not overwhelming.

Table 1:  Type I error for the test of $H_{0}:\beta=0$. Comparison of Type-I error at the $0.05$ significance level for a network of $n$ nodes with $K$ balanced communities for different choices of within community edge probability ($p$) and between community edge probability ($q$) among GEE with independent and exchangeable correlation structure and naive least squares. $(n,K)$ | $(p,q)$ | GEE (independent) | GEE (exchangeable) | LS
---|---|---|---|---
| (0.8, 0) | 0.049 | 0.055 | 0.062
$(200,20)$ | (0.7, 0.1) | 0.050 | 0.055 | 0.079
| (0.6, 0.2) | 0.058 | 0.061 | 0.091
| (0.5, 0.3) | 0.072 | 0.075 | 0.095
| (0.8, 0) | 0.050 | 0.050 | 0.060
$(400,40)$ | (0.7, 0.1) | 0.051 | 0.056 | 0.075
| (0.6, 0.2) | 0.059 | 0.062 | 0.090
| (0.5, 0.3) | 0.070 | 0.075 | 0.096

###  3.3 Real data analysis

Air travel networks were constructed from flight data retrieved from the
Official Airline Guide (OAG), with node attribute outcomes being the country-
specific COVID-19 incidence rates by month available from the Johns Hopkins
University coronavirus data repository through April 30, 2020, (COVID, 19). In
addition we included country-specific GDP (WBG, 20201),total population (WBG,
20202), and percentage of the urban population (WBG, 20203)of the countries
from the website of the World Bank as covariates. We study the effectiveness
of travel bans on the incidence rate of the COVID-19 using our network
regression model by performing a month-by-month analysis from January-April,
2020 which span the period of before the pandemic till one month after the
travel restriction was in effect. In addition, we study the importance of
baseline covariates effects such as GDP, percentage of urban population, and
population size of the countries vs the network effect. The GDP and the
population of the most populated countries are in the order of $10^{12}$ and
$10^{6}$ respectively, so we scale these variables by $10^{12},10^{6}$, and
$10^{2}$, respectively, in order to stablize coefficient estimation.

With these data retrieved across different continents to model the network
effect through our porposed model in (1) we assume a directed stochastic block
model framework, where nodes correspond to countries and each block contains
countries having a large number of commercial flights traveling among them
compared to the others. Edge formation is determined by thresholding the
population-normalized count of flights arriving at the destination country.

Figure 2: Changes in air travel networks and the impact of incidence rates
over time. Monthly average number of flights per million population over the
countries and monthly estimate of $\beta$ for Jan-Apr 2020 demonstrate a
downward and an upward trend, respectively.

We let the incidence rate $y_{i}$ of the $i$th country to be the number of
cases per 1000 populations, and the covariates $\bm{x}_{i_{n\times
4}}=(\bm{1}_{n\times 1}:\bm{x}_{{i2}_{n\times 1}}:\bm{x}_{{i3}_{n\times
1}}:\bm{x}_{{i4}_{n\times 1}})$ to be a matrix whose first column has 1s as
the intercept and the second, third and fourth columns have the population
size, GDP, percentage of the urban population of the $i$th country scaled by
$10^{6}$, $10^{12}$, and $10^{2}$ respectively. For constructing the adjacency
matrix $\bm{A}$, we consider two scenarios: unweighted and weighted. For
unweighted adjacency matrices, from the directed graph dictated by the number
of flights in a particular month, we first construct a count matrix $\bm{C}$
that counts the number of flights from one country to the other. Then we
construct an unweighted adjacency matrix $\bm{A}$ with 0s and 1s where the
entry of $\bm{A}$ is 1 if it exceeds the third quartile of the elements of
$\bm{C}$, and 0 otherwise. For weighted adjacency matrices we divide the
elements of $\bm{C}$ by the total population of the destination country per
million. The main rationale behind this scaling is that one would expect more
flights traveling to a populated country compared to a less populated country.
Therefore, dividing it by the population of the country will result in the
elements of the weighted adjacency matrix being on the same scale.

With the two adjacency matrices constructed in this way, we fit model (3) and
summarize the results in Table 2. The estimates of $\beta$ increase from
February to April both for the weighted and unweighed cases suggesting that
there is an increasing association between the travel and spread of the
pandemic. To further investigate this behavior, we plot the average number of
flights per million population in Figure 2 which shows that although the
average number of flights is decreasing from January to April, the estimates
of $\beta$ are increasing. This reflects the fact that while travel bans led
to a decrease in the total number of flights in March and April, the
increasing $\beta$ in these months implies that each flight had an increased
likelihood of transmitting COVID-19 to the destination country, increasing the
correlation between the incidence rates in the two countries.

Table 2 also demonstrates that while the overall population of a country has a
negligible effect on the incidence rate leading to smaller values of
$\alpha_{2}$, the percentage of the urban population plays a crucial role
especially when the travel ban has already been in effect, for example, the
value of $\alpha_{4}$ is 62 (55 for weighted case) in April compared to 0 in
February. Moreover, the corresponding values are consistent for weighted and
unweighted networks.

Table 2 also demonstrates the utility of including baseline covariates
alongside network effect to study the effectiveness of travel bans in
mitigating the spread of COVID-19. By comparing the values of $\beta$ and
$\alpha_{4}$, one can note that for the initial months of the pandemic, the
network effect is more compared to the effect of urban population. However,
when the travel ban has already taken place during March (for most of the
countries), the effect of urban population supersedes the network effect as
its value increase to 62.33 drastically from 9.88 suggesting that urbanity is
the next important factor to consider for controlling the spread of the
pandemic after the travel bans.

We also compare our results with naive least squares in Table 3 which shows
coefficient estimates for both models. Of note, the April estimate of $\beta$
is dramatically larger in the naive least squares model, which may be a
manifestation of the increased bias we expect to see for naive least squares.
Standard errors for the coefficient estimates are uniformly smaller for least
squares as we expect, likely a representation of the overconfidence of the
least squares model that comes from ignoring network induced correlation. This
case demonstrates that using least squares would incorrectly dramatically
increase both the magnitude and degree of significance of the network effects
on incidence rates, particularly in the months post-lockdown when travel bans
were in place. In contrast, the GEE model provides a more realistic view.

Table 2: Estimates of the parameters for unweighted and (weighted) networks
under the GEE model. $\beta,\alpha_{1},\alpha_{2},\alpha_{3}$ and $\alpha_{4}$
correspond to the coefficients of the adjacency matrix (network effect),
intercept, population, GDP, and percentage of the urban population,
respectively.

Month of 2020 | $\beta$ | $\alpha_{1}$ | $\alpha_{2}$ | $\alpha_{3}$ | $\alpha_{4}$
---|---|---|---|---|---
Jan | 1.8759 (0.0038) | -0.0002 (-0.0001) | 0.0000 (0.0000) | 0.0013 (0.0013) | 0.0003 (0.0000 )
Feb | 2.1946 (0.0014) | -0.0048 (-0.0037) | 0.0000 (0.0000) | 0.0555 (0.0575) | -0.0094 (-0.0072 )
Mar | 2.7186 (0.0030) | -0.7846 (-0.7490) | -0.0005 (-0.0007) | -0.1446 (-0.0233) | 9.8844 (9.8862 )
Apr | 4.6574 (0.0328) | -4.7307 (-3.5655) | -0.0085 (-0.0125) | 0.0626 (0.4685) | 62.3300 (54.9767)

Table 3: Comparison of the naive linear regression with GEE for the weighted
networks. The numbers reported in the parenthesis correspond to our GEE method
while the others correspond to least squares.
$\beta,\alpha_{1},\alpha_{2},\alpha_{3}$ and $\alpha_{4}$ correspond to the
coefficients of the adjacency matrix, intercept, population, GDP, and
percentage of the urban population respectively.

Month of 2020 | $\beta$ | $\alpha_{1}$ | $\alpha_{2}$ | $\alpha_{3}$ | $\alpha_{4}$
---|---|---|---|---|---
Jan | 0.0424 (0.0038) | -0.0001 (-0.0001) | 0.0000 (0.0000) | 0.0013 (0.0013) | 0.0000 (0.0000 )
Feb | 0.0105 (0.0014) | -0.0037 (-0.0037) | 0.0000 (0.0000) | 0.0576 (0.0575) | -0.0073 (-0.0072 )
Mar | 0.0458 (0.0030) | -0.7060 (-0.7490) | -0.0009 (-0.0007) | 0.0017 (-0.0233) | 9.5519 (9.8862 )
Apr | 1.2655 (0.0328) | -5.5415 (-3.5655) | -0.0073 (-0.0125) | -0.2041 (0.4685) | 67.5900 (54.9767)

##  4 Discussion

We have proposed a generalized estimating equation (GEE) approach to network
regression model. Assuming independent community structure and using the
simultaneous estimation of memberships, the GEE approach allows for community
dependent covariate coefficient estimation, and thereby provides a flexible
and efficient solution to the network regression model. Moreover, this allows
us to do a hypothesis testing of the network regression parameter $\beta$
which helps us to decide the importance of including such term in our
analysis. We provided a relevant real data example of COVID-19 cases along
with baseline covariates such as GDP, population size, and percentage of urban
population of countries across different continents, and the number of
commercial flights traveling between them to study the importance of travel
bans to mitigate the spread of the COVID-19 pandemic. We have constructed the
adjacency matrix from the count of flights rendered in a network of countries
via a stochastic block model where each block contains countries having a
similar number of flights. Our proposed model has helped us to understand the
importance of the baseline covariates vs network effect. Since we have dealt
with longitudinal data, it is also instructive to note that our proposed model
offers us the flexibility of clustering both in the network space and also
over time.

One limitation of our results is the balanced community assumption made in
Proposition 3.1 which may not reflect some unbalanced networks. This
assumption can be relaxed somewhat with small departures from the balanced
design. Further, under the stochastic block model, extremely unbalanced
networks still allow for consistent estimation of $p$ and $q$, however for a
more flexible network model that allows for community-specific edge
probabilities, consistent estimation of edge probabilities requires each
community to grow asymptotically with $n$.

##  5 Appendix

###  5.1 Proof of Proposition 3.1

From the construction of the adjacency matrices, one can note that the entries
$A_{ij}$s are i.i.d. Bernoulli random variables with

$E(A_{ij})=\begin{cases}p,\;\;\ \text{if }i,j\text{ belong to the same
community}\\\ q,\;\;\ \text{otherwise}\\\ \end{cases}$

and

$Var(A_{ij})=\begin{cases}p(1-p),\;\;\ \text{if }i,j\text{ belong to the same
community}\\\ q(1-q),\;\;\ \text{otherwise}\\\ \end{cases}$

Denote by $S$ the set $\\{i,j:1\leq i<j\leq n,i,j\text{ belongs to the same
community}\\}$, and its complement by $S^{\prime}$. Therefore, for a directed
graph, one can write

$\displaystyle E\Big{(}\sum_{i,j\in S}A_{ij}\Big{)}=2K{m\choose 2}p,$
$\displaystyle Var\Big{(}\sum_{i,j\in S}A_{ij}\Big{)}=2K{m\choose
2}p(1-p)=s_{m_{p}}^{2},$ $\displaystyle E\Big{(}\sum_{i,j\in
S^{\prime}}A_{ij}\Big{)}=K(K-1)m^{2}q,\text{ and }$ $\displaystyle
Var\Big{(}\sum_{i,j\in
S^{\prime}}A_{ij}\Big{)}=K(K-1)m^{2}q(1-q)=s_{m_{q}}^{2}.$

Next, one can verify Lindeberg condition to establish the central limit
theorem:

$\displaystyle\frac{\sum_{i,j\in S}A_{ij}-2K{m\choose 2}p}{\sqrt{2K{m\choose
2}p(1-p)}}\;\;\ \text{ and }\;\;\frac{\sum_{i,j\in
S^{\prime}}A_{ij}-K(K-1)m^{2}q}{\sqrt{K(K-1)m^{2}q(1-q)}}\overset{d}{\rightarrow}N(0,1)$
(4)

as $Km^{2}p(1-p)$ and $Km^{2}q(1-q)\rightarrow\infty$.

The Lindeberg’s condition requires us to verify

$\frac{1}{s_{m_{p}}^{2}}\sum_{i,j\in
S}E\Big{(}B_{ij}^{2}I_{\\{|B_{ij}|\geq\epsilon s_{m_{p}}\\}}\Big{)}\rightarrow
0\;\;\ \text{ and }\frac{1}{s_{m_{q}}^{2}}\sum_{i,j\in
S^{\prime}}E\Big{(}B_{ij}^{2}I_{\\{|B_{ij}|\geq\epsilon
s_{m_{q}}\\}}\Big{)}\rightarrow 0,$

where $B_{ij}=A_{ij}-E(A_{ij})$. Since $|B_{ij}|\leq 1$, the above condition
is satisfied when $Km^{2}p(1-p)$ and $Km^{2}q(1-q)$ tend to $\infty$ under
which $\epsilon s_{m_{p}}$ and $\epsilon s_{m_{q}}$ are both greater than 1.

Dividing the numerator and denominator of (4) by $2K{m\choose 2}$ and
$K(K-1)m^{2}$ respectively, one obtains both

$m^{\gamma/2}\sqrt{2K{m\choose 2}}\frac{\sum_{i,j\in
S}A_{ij}/\Big{(}2K{m\choose 2}\Big{)}-p}{\sqrt{m^{\gamma}p(1-p)}}\;\;\ \text{
and }\;\;m^{\gamma/2}\sqrt{K(K-1)m^{2}}\frac{\sum_{i,j\in
S^{\prime}}A_{ij}/\Big{(}K(K-1)m^{2}\Big{)}-q}{\sqrt{m^{\gamma}q(1-q)}}$ both
converge to $N(0,1)$ distribution. Since $\hat{p}=\sum_{i,j\in
S}A_{ij}/\Big{(}K{m\choose 2}\Big{)}$, and $\hat{q}=\sum_{i,j\in
S^{\prime}}A_{ij}/\Big{(}K(K-1)m^{2}\Big{)}$, both $\hat{p}$ and $\hat{q}$ are
$m^{1+\gamma/2}$ consistent.

###  5.2 Proof of Theorem 3.2

Here we present the sketch of the proof of Theorem 3.2. One can write (3) as

$\displaystyle\sum_{k=1}^{K}U_{k}(\alpha,\bm{\beta},p,q)=\sum_{k=1}^{K}\bm{D}_{k}^{\top}\bm{V}_{k}^{-1}\bm{S}_{k}=0,$
(5)

where $\bm{S}_{k}=\bm{y}_{k}-\bm{\mu}_{k}$, and $U$ is a function of the model
parameters.

Let $\bm{b}=(\beta,\bm{\alpha})^{\top}$ denote the vector of model parameters,
and $\bm{\pi}=(p,q)^{\top}$ denote the vector of model parameters, and vector
of within and between edge probabilities respectively. Letting $\bm{b}$ fixed,
the Taylor expansion yields

$\displaystyle\frac{\sum_{k=1}^{K}U_{k}(\bm{b},\bm{\pi}^{*})}{K^{1/2}m^{1+\gamma/2}}=\frac{\sum_{k=1}^{K}U_{k}(\bm{b},\bm{\pi})}{K^{1/2}m^{1+\gamma/2}}+\frac{\sum_{k=1}^{K}\frac{\partial
U_{k}(\bm{b},\bm{\pi})}{\partial\pi}}{K^{1/2}m^{1+\gamma/2}}m^{1+\gamma/2}(\bm{\pi}^{*}-\bm{\pi})+o_{P}(1)$
(6) $\displaystyle\hskip 79.6678pt=\tilde{A}+\tilde{B}\tilde{C}+o_{P}(1),$

One can note that $\tilde{B}=o_{P}(1)$ as $\partial
U_{k}(\bm{b},\bm{\pi})/\partial\bm{\pi}$ are linear functions of
$\bm{S}_{k}$’s defined in (5) whose means are zero, and $\tilde{C}=O_{P}(1)$
thanks to Proposition 3.1. Therefore,
$\frac{\sum_{k=1}^{K}U_{k}(\bm{b},\bm{\pi}^{*})}{K^{1/2}m^{1+\gamma/2}}$ is
asymptotically equivalent to
$\frac{\sum_{k=1}^{K}U_{k}(\bm{b},\bm{\pi})}{K^{1/2}m^{1+\gamma/2}}$ whose
asymptotic distribution is multivariate normal with zero mean and covariance
is equal to $V$ as defined in Theorem 3.2. The proof is thus complete
following Liang Zeger (1986).

##  6 Acknowledgement

RG would like to thank Anupam Kundu, postdoctoral associate at Yale School of
Public Health and Thien Le, postdoctoral fellow at Harvard for many useful
discussions regarding real data processing. JP acknowledges support from
R01AI138901, and IB acknowledges support from R01MH116884.

## References

  * Aukett . (1988) aukett1988genderAukett, R., Ritchie, J. Mill, K.  1988\.  Gender differences in friendship patterns Gender differences in friendship patterns.  Sex roles19157–66.
  * Bank (2020) world2020covidBank, W.  2020\.  COVID-19 to plunge global economy into worst recession since World War II. Covid-19 to plunge global economy into worst recession since World War II.  World Bank Group Washington, DC.
  * Binkiewicz . (2017) binkiewicz2017covariateBinkiewicz, N., Vogelstein, JT. Rohe, K.  2017\.  Covariate-assisted spectral clustering Covariate-assisted spectral clustering.  Biometrika1042361–377.
  * Burton . (1998) burton1998extendingBurton, P., Gurrin, L. Sly, P.  1998\.  Extending the simple linear regression model to account for correlated responses: an introduction to generalized estimating equations and multi-level mixed modelling Extending the simple linear regression model to account for correlated responses: an introduction to generalized estimating equations and multi-level mixed modelling.  Statistics in Medicine17111261–1291.
  * Chen . (2021) chen2021analysisChen, M., Kato, K. Leng, C.  2021\.  Analysis of networks via the sparse $\beta$-model Analysis of networks via the sparse $\beta$-model.  Journal of the Royal Statistical Society: Series B (Statistical Methodology).
  * Christakis Fowler (2007) christakis2007spreadChristakis, NA. Fowler, JH.  2007\.  The spread of obesity in a large social network over 32 years The spread of obesity in a large social network over 32 years.  New England Journal of Medicine3574370–379.
  * COVID (19) covid19databaseCOVID, G.  19\.  database (2021) Retrieved from https://github. com/CSSEGISandData database (2021) retrieved from https://github. com/cssegisanddata.  COVID-19.
  * Diggle . (2002) diggle2002analysisDiggle, P., Diggle, PJ., Heagerty, P., Liang, KY., Zeger, S. .  2002\.  Analysis of longitudinal data Analysis of longitudinal data.  Oxford University Press.
  * Halekoh . (2006) geeglmHalekoh, U., Højsgaard, S. Yan, J.  2006\.  The R Package geepack for Generalized Estimating Equations The r package geepack for generalized estimating equations.  Journal of Statistical Software15/21–11.
  * P. Hoff (2021) hoff2021additiveHoff, P.  2021\.  Additive and multiplicative effects network models Additive and multiplicative effects network models.  Statistical Science36134–50.
  * PD. Hoff . (2002) hoff2002latentHoff, PD., Raftery, AE. Handcock, MS.  2002\.  Latent space approaches to social network analysis Latent space approaches to social network analysis.  Journal of the American Statistical Association974601090–1098.
  * Holland . (1983) holland1983stochasticHolland, PW., Laskey, KB. Leinhardt, S.  1983\.  Stochastic blockmodels: First steps Stochastic blockmodels: First steps.  Social Networks52109–137.
  * Holland Leinhardt (1981) holland1981exponentialHolland, PW. Leinhardt, S.  1981\.  An exponential family of probability distributions for directed graphs An exponential family of probability distributions for directed graphs.  Journal of the American Statistical Association7637333–50.
  * Igarashi . (2005) igarashi2005genderIgarashi, T., Takai, J. Yoshida, T.  2005\.  Gender differences in social network development via mobile phone text messages: A longitudinal study Gender differences in social network development via mobile phone text messages: A longitudinal study.  Journal of Social and Personal Relationships225691–713.
  * Kelejian Prucha (1998) kelejian1998generalizedKelejian, HH. Prucha, IR.  1998\.  A generalized spatial two-stage least squares procedure for estimating a spatial autoregressive model with autoregressive disturbances A generalized spatial two-stage least squares procedure for estimating a spatial autoregressive model with autoregressive disturbances.  The Journal of Real Estate Finance and Economics17199–121.
  * Le . (2022) le2022frameworkLe, TM., Raynal, L., Talbot, O., Hambridge, H., Drovandi, C., Mira, A.Onnela, JP.  2022\.  Framework for assessing and easing global COVID-19 travel restrictions Framework for assessing and easing global covid-19 travel restrictions.  Scientific Reports1211–13.
  * Liang Zeger (1986) liang1986longitudinalLiang, KY. Zeger, SL.  1986\.  Longitudinal data analysis using generalized linear models Longitudinal data analysis using generalized linear models.  Biometrika73113–22.
  * Mandel . (2021) mandel2021neuralMandel, F., Ghosh, RP. Barnett, I.  2021\.  Neural networks for clustered and longitudinal data using mixed effects models Neural networks for clustered and longitudinal data using mixed effects models.  Biometrics.
  * Mu . (2022) mu2022spectralMu, C., Mele, A., Hao, L., Cape, J., Athreya, A. Priebe, CE.  2022\.  On spectral algorithms for community detection in stochastic blockmodel graphs with vertex covariates On spectral algorithms for community detection in stochastic blockmodel graphs with vertex covariates.  IEEE Transactions on Network Science and Engineering.
  * Newman (2018) newman2018networkNewman, ME.  2018\.  Network structure from rich but noisy data Network structure from rich but noisy data.  Nature Physics146542–545.
  * Qin Rohe (2013) qin2013regularizedQin, T. Rohe, K.  2013\.  Regularized spectral clustering under the degree-corrected stochastic blockmodel Regularized spectral clustering under the degree-corrected stochastic blockmodel.  Advances in Neural Information Processing Systems26.
  * Rosvall Bergstrom (2007) rosvall2007mapsRosvall, M. Bergstrom, CT.  2007\.  Maps of information flow reveal community structure in complex networks Maps of information flow reveal community structure in complex networks.  arXiv preprint physics.soc-ph/0707.0609.
  * Schäfer . (2014) schafer2014bringingSchäfer, M., Strohmeier, M., Lenders, V., Martinovic, I. Wilhelm, M.  2014\.  Bringing up OpenSky: A large-scale ADS-B sensor network for research Bringing up opensky: A large-scale ads-b sensor network for research.  IPSN-14 Proceedings of the 13th International Symposium on Information Processing in Sensor Networks Ipsn-14 proceedings of the 13th International Symposium on Information Processing in Sensor Networks ( 83–94).
  * Shalizi Thomas (2011) shalizi2011homophilyShalizi, CR. Thomas, AC.  2011\.  Homophily and contagion are generically confounded in observational social network studies Homophily and contagion are generically confounded in observational social network studies.  Sociological Methods & Research402211–239.
  * Shrum . (1988) shrum1988friendshipShrum, W., Cheek Jr, NH. MacD, S.  1988\.  Friendship in school: Gender and racial homophily Friendship in school: Gender and racial homophily.  Sociology of Education227–239.
  * Staber (1993) staber1993friendsStaber, U.  1993\.  Friends, acquaintances, strangers: Gender differences in the structure of entrepreneurial networks Friends, acquaintances, strangers: Gender differences in the structure of entrepreneurial networks.  Journal of Small Business & Entrepreneurship11173–82.
  * Strohmeier . (2021) strohmeier2021crowdsourcedStrohmeier, M., Olive, X., Lübbe, J., Schäfer, M. Lenders, V.  2021\.  Crowdsourced air traffic data from the OpenSky Network 2019–2020 Crowdsourced air traffic data from the opensky network 2019–2020.  Earth System Science Data132357–366.
  * Tifferet (2019) tifferet2019genderTifferet, S.  2019\.  Gender differences in privacy tendencies on social network sites: A meta-analysis Gender differences in privacy tendencies on social network sites: A meta-analysis.  Computers in Human Behavior931–12.
  * VanderWeele (2011) vanderweele2011sensitivityVanderWeele, TJ.  2011\.  Sensitivity analysis for contagion effects in social networks Sensitivity analysis for contagion effects in social networks.  Sociological Methods & Research402240–255.
  * WBG (20201) capita2020availablegdpWBG.  20201\.  Available online: https://data.worldbank.org/indicator/NY.GDP.MKTP.CD Available online: https://data.worldbank.org/indicator/ny.gdp.mktp.cd.  CD (accessed on 2 September 2022).
  * WBG (20202) capita2020availabletotpopWBG.  20202\.  Available online: https://data.worldbank.org/indicator/SP.POP.TOTL Available online: https://data.worldbank.org/indicator/sp.pop.totl.  CD (accessed on 2 September 2022).
  * WBG (20203) capita2020availableurbanWBG.  20203\.  Available online: https://https://data.worldbank.org/indicator/SP.URB.TOTL.IN.ZS Available online: https://https://data.worldbank.org/indicator/sp.urb.totl.in.zs.  CD (accessed on 2 September 2022).
  * Zeger (1986) zeger1986analysisZeger, S.  1986\.  The analysis of discrete and continuous longitudinal data The analysis of discrete and continuous longitudinal data.  Biometrics42121–130.