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.
86,194 characters · 32 sections · 103 citation commands
Inference in Linear Dyadic Data Models with Network Spillovers
Dyadic data is categorized by the dependence between two sets of sampled units (dyads). For example, exports between the U.S. and Canada depend on both countries (and, plausibly, their characteristics). This contrasts to classical data in the social sciences that only depends on a single unit of observation (e.g., the GDP of the U.S., or a politician's vote in a roll-call).
The empirical relevance of dyadic data is showcased by its widespread use, which has increased over the past two decades (Graham-2020a provides an extensive review). For example, applications are found in political economy (correlation in voting behavior in Parliament across seating neighbors, Harmon_et_al-2019), international political economy and trade (export-import outcomes across countries, Anderson_and_van_Wincoop-2003), international relations (hoff04, for a salient example), among many others. In fact, dyadic data is considered to be dominant in quantitative international relations poast. In these examples, applied researchers typically model the dependence between dyadic outcomes and observable characteristics using a linear model, which they then estimate using Ordinary Least Squares (OLS). However, inference on such estimators for the linear parameters is more complex.
The main approach in recent applied work has been the use of the so-called “dyadic-robust" estimators (e.g., Cameron_et_al-2011, Aronow_et_al-2015, and Tabord-Meehan-2019, among others). Such estimators build on the widely-used assumption in dyadic data that the error terms for dyad $(i, j)$ and for dyad $(k, l)$ can only be correlated if they share a unit (see Aronow_et_al-2015 and Tabord-Meehan-2019 for a discussion; Cameron_and_Miller-2014wp for a review).
In this paper, we first argue that such an assumption does not hold in many applications using dyadic data where dyads may be indirectly connected along a network.\footnote{This is a concrete class of applied examples where the assumption fails. The possibility that cross-sectional dependence in dyadic data might be more extensive than assumed has been pointed out by Cameron_and_Miller-2014wp and cranmerdemarais.} Figure (ref) presents a simple example in the context of politicians in Congress, whose votes or decisions depend on their seating neighbors. It is completely possible that behavior across dyads $(A,B)$ and $(C,D)$ might be correlated along unobservables because they have many indirect connections (in the figure, through $A$ sitting next to $B$, who sits next to $C$).\footnote{We expand on these examples in the next section. Such spillovers could be further rationalized as individual-level unobserved heterogeneity: e.g., an unmeasured preference for voting Yes, or a preference for trading with a certain country (see Graham-2020b and references therein for details).} We show that such spillovers invalidate the assumptions for consistency of existing “dyadic-robust" variance estimators through generating interdependence, implying that they are biased for the true asymptotic variance when dyads may be correlated even when they do not share a common unit.
To deal with these issues, we develop a consistent variance estimator that explicitly accounts for such network spillovers even with dyadic data, thereby complementing existing approaches (e.g., Aronow_et_al-2015).\footnote{We provide an extensive comparison of the relative benefits of each approach in the next section. We note here, though, that neither approach subsumes the other, as they depend on different assumptions and may be more appropriate for different applications.} We prove that our proposed variance estimator is consistent for the true variance of the OLS estimator in linear models with dyadic data when the cross-sectional dependence follows an observed (exogenous) network. Our main insight is that the dependence across all dyads, including indirect spillovers, can be rewritten as correlations across a specific network over dyads. This allows us to apply the framework of Kojevnikov_et_al-2021 to such network random variables, although here it is a network over dyads, rather than individuals. Monte Carlo simulations show that our proposed estimator has good finite sample properties and outperforms other estimators for the relevant contexts.
To help practitioners, we then provide a step-by-step guideline on whether our estimator may be appropriate to their context. As we describe, this choice depends on: (i) whether spillovers from indirectly connected dyads are likely to be present, (ii) whether the researcher observes/constructs the network among dyads through which spillovers propagate, and (iii) whether those spillovers are likely to be persistent. Our variance estimator is consistent for the asymptotic variance of the OLS estimator even under (i)-(iii). And our estimator can account for decay in propagation, as Corollary (ref) and Example (ref) illustrate. On the other hand, if decay is very high, or spillovers along the network do not exist, then those same results imply that the estimator of Aronow_et_al-2015 will also be consistent.
Finally, we illustrate the extent to which neglecting network spillovers with dyadic data may bias inference results. Beyond Monte Carlo simulations, we revisit the application in Harmon_et_al-2019 of voting in the European Parliament. The authors study whether random seating arrangements (based on naming conventions) induce neighboring politicians to agree with one another in policy votes. The outcome, whether politicians $i$ and $j$ vote the same way on a policy, is dyadic in nature. However, $i$ and $j$'s votes may be positively correlated even if they are not neighbors: for instance, $i$ and $j$ may sit on either side of common neighbors $k$ and $l$, who influence them both, and this seating arrangement is observed. This chain of influences is sufficient to induce strong positive correlation across non-dyads. We show that neglecting such higher order spillovers has significant empirical consequences: their estimated variance using the estimator in Aronow_et_al-2015 is roughly $22\%$ smaller than using our consistent estimator accounting for such spillovers; while the estimate based on the Eicker-Huber-White estimator ignoring spillovers is approximately 73% smaller than our proposal, consistent with the arguments of pinto.
The use of dyadic data in Political Science has a rich history, particularly in International Relations. However, empirical challenges with such models are well known -- see poast for a historical overview. Early on, the concerns were mostly about model specification, including the error term. This includes the 2001 special issue of International Organization, mostly focusing on the use of fixed effects. More recently, pinto pointed out that ignoring dependence across dyads can lead to erroneous hypothesis testing, as computed standard errors would be too small. hoff04 and minhas19, minhas22 suggest including random coefficients and latent variables to account for dependencies across dyads. Our approach explicitly accounts for the whole network of interdependencies across dyads, which can go beyond third-order dependences (assumed in minhas19,minhas22). It does so by using asymptotic inference, rather than Bayesian (minhas19,minhas22) or randomized inference (pinto).
As a result, our paper is directly related to the literature on (asymptotic) inference in regression analysis with dyadic random variables. Aronow_et_al-2015 and Tabord-Meehan-2019 consider OLS estimation and inference in a linear dyadic regression model. Meanwhile, Graham-2020a and Graham-2020b explore a likelihood-based approach to dyadic regression models, while Graham_et_al-2022wp and Chiang_and_Tan-2022 provide results for kernel density estimation in dyadic regression models.\footnote{Similarly, our paper is related to the recent literature on inference for multiway clustering, whether OLS estimation with multiway clustering (Cameron_et_al-2011); a clustering method in high-dimensional set-ups (Chiang_et_al-2021); clustering within the time dimension (Chiang_et_al-2022wpa); clustered inference with empirical likelihood (Chiang_et_al-2022wpb); bootstrap methods in multiway clustering Davezies_et_al-2021, Menzel-2021, Mackinnon_et_al-2022wpa; clustering in the context of average treatment effects (Abadie_et_al-2022), to name but a few (see MacKinnon_et_al-2022 for review). Again, one of the common assumptions in this literature is that individual observations are divided into disjoint groups -- clusters -- and observations in different clusters are not correlated. To that end, Mackinnon_et_al-2022wpb propose measures of cluster-level influence that can be used to assess whether the underlying assumption of cluster-robust variance estimation is satisfied.} While useful to allow for correlations along time and within such groups, this separable structure may be inappropriate for environments where spillovers follow a complex form of dependence along a network. As alluded to previously, all such papers assume that dyads are uncorrelated unless they share a common unit, thereby ignoring dependence across indirectly linked dyads, while our estimator explicitly captures the higher order correlation across dyads along the dyadic network.
However, we emphasize that neither approach subsumes the other. The papers cited above leave the dependence within “clusters" (groups of dyads that share units) unrestricted, but assume independence across such clusters. This is akin to the literature with one-way clustering (e.g., Hansen_and_Lee-2019). By comparison, our approach restricts such dependence among groups of dyads that share units (i.e., dependence is assumed to follow the observed network), but allows for dependence across such “clusters" of dyads along the dyadic network.\footnote{As a result, our approach complements the existing methods similar to how inference with spatial data (e.g., Conley-1999 and Jenish_and_Prucha-2009) complements one-way clustering inference. Our approach still differs from such inference with spatial data, since the latter routinely assumes the index set to be a Euclidean metric space, whose metric relies solely on the nature of the space, and uses it to define the dependence between variables. See also Ibragimov_and_Muller-2010 and references therein. In our context, however, the index set of dyads alone does not suffice to dictate the dependence structure because indices themselves do not inform us of the network topology. Instead, we first introduce a metric on a network among dyads and our mixing condition is based on dependence as dyads grow further apart along the network.} Which approach to pursue using dyadic data depends on the researchers' applications and how they fit such assumptions.
Assume that we observe a cross section of $N\in \mathbbm{N}$ individuals located along a network -- the latter interpretable as politicians, countries, firms, or other observation units depending on the context. The dyads present in the $N$-individual network (i.e., among the $\tbinom{N}{N-1}$ possible dyads), are called active dyads, so that the dyad for two units $i$ and $j$ (e.g., politicians, countries) is denoted as some $m$. The set of active dyads is denoted $\mathcal{M}_{N}$ and $M$ denotes the cardinality of that set.
We assume that each dyad $m$ is endowed with a triplet of dyad-specific variables, forming a triangular array $\{(y_{M, m}, x_{M, m}, \varepsilon_{M, m})\}_{m\in\mathcal{M}_{N}}$ with respect to $M$, where $y_{M,m}\in \mathbbm{R}$ is a one-dimensional observable outcome, $x_{M,m}\in \mathbbm{R}^{K}$ is a $K$-dimensional vector of observable characteristics with $K\in\mathbbm{N}$, and $\varepsilon_{M,m}\in\mathbbm{R}$ is a one-dimensional random error term that is not observable to the researcher. We only consider exogenous network formation and the network is assumed to be observable. These conditions are summarized in the following assumption:
While such assumptions are standard in models of dyadic networks, they seem particularly appropriate when units or dyad pairs are linked across geographical, physical, or ex-ante social relations (e.g., family ties). This includes capturing neighboring and regional spillovers across countries, as often done in international relations, or exogenous seating arrangements in Parliament, as illustrated in the examples in the next section.
The subsequent arguments require us to distinguish between a pair of dyads who share a member (i.e., who are directly linked -- which we call, adjacent) and a pair of dyads who are directly or indirectly linked (which we call, simply, connected).
In Figure (ref), dyad (A,B) is adjacent to (B,C), and connected with, though not adjacent to, (C,D). Hence, the adjacency relationship constitutes a network structure among active dyads, and thus, a network over individuals can be transformed to one over active dyads. For example, the right-hand side panel of Figure (ref) provides a network over pairs of voting politicians (i.e., active dyads).\footnote{This corresponds to thinking about the {\sl line graph} of the original graph over individuals.} We define the geodesic distance between two connected dyads $m$ and $m'$ to be the smallest number of adjacent dyads between them. Note that adjacent dyads are a special case of connected dyads with geodesic distance equal to one.
The cross-sectional model of interest takes the form of the linear network-regression model: for any $N\in\mathbbm{N}$,
where
and $\beta$ is a $K\times 1$ vector of the regression coefficients and $X_{M}$ denotes the $M\times K$ matrix that records the observed dyad-specific characteristics, i.e., $X_{M}\coloneqq[x_{M,1},\ldots, x_{M,M}]'$.
In this paper, we assume that $\beta$ is identified, which follows from standard assumptions on strict exogeneity, lack of multicollinearity and the existence of finite second moments of $y_{M,m}$ and $x_{M,m}$. (For completeness, see Assumption (ref) and Proposition (ref) in Appendix).
We note that equation ((ref)) allows for there to be spillovers across the error terms even when dyads $m$ and $m'$ are not adjacent, as long as they are connected through indirect links. By comparison, applied researchers such as Harmon_et_al-2019 and Lustig_and_Richmond-2020 (and the estimators of Aronow_et_al-2015 and Tabord-Meehan-2019) consider a variant of the linear regression (ref) under the assumption
with $m = d(i,j)$ representing the dyad between $i$ and $j$. This specific assumption would be equivalent to setting:
Whether to allow indirect spillovers as in ((ref)) or not (as in (ref)) depends on the researchers' applications. We now present two examples where our approach may be preferable.
Throughout this paper, we focus on the Ordinary Least Squares (OLS) estimator of $\beta$, denoted by $\hat{\beta}$. Under the assumptions above, we can write
It is straightforward to verify that $\hat{\beta}$ is unbiased for $\beta$ under our identification conditions (Assumption (ref)). However, a consistency result is by no means trivial due to the dependence along the network which induces a complex form of cross-sectional dependence, hindering a na\"{i}ve application of the standard theory for independently and identically distributed ($i.i.d.$) random vectors. This is clarified in Section (ref). Before that, we provide a heuristic overview of the main point of our paper in this setting.
Inference about $\beta$ is based on a normal approximation of the distribution of $\hat{\beta}$ around $\beta$. We focus on hypothesis testing conducted using the expression:
where $\widehat{Var}(\hat{\beta})$ is a consistent estimator of the asymptotic variance of $\hat{\beta}$. Our main result in Section (ref) is providing such an appropriate estimator, which takes the form: {
} where $\kappa_{m,m'}$ is an appropriate kernel function that will formally be defined in Section (ref); $h_{m,m'}$ represents an indicator function that takes one if dyads $m$ and $m'$ are connected and zero otherwise; and $\hat{\varepsilon}_{m} \coloneqq y_{m}- x_{m}' \hat{\beta}$.
This paper derives conditions under which $\widehat{Var}(\hat{\beta})$ is consistent for the asymptotic variance of $\hat{\beta}$. Before doing so, let us compare the variance estimator ((ref)) with an often used estimator based on one-way clustering of dyad groupings.
In order to validate the test statistic ((ref)) for inference, this paper establishes the consistency and asymptotic normality of $\hat{\beta}$ and develops a consistent network-robust estimator for the asymptotic variance when the correlation structure is given by (ref).
Researchers interested in the linear specification ((ref)) may find the following guidelines useful when deciding whether to use our proposed estimator ((ref)).
Let $Y_{M,m}$ be a random vector defined as
\footnotetext{By construction, the collection of $Y_{M,m}$'s constitutes a triangular array of random vectors.} and denote $\mathcal{C}_{M} \coloneqq \{x_{M, m}\}_{m\in\mathcal{M}_{N}}$.\footnote{For the case of stochastic networks, it is defined to include information about the network topology as well as the collection of the dyad-specific attributes $\{x_{M, m}\}_{m\in\mathcal{M}_{N}}$. See Kojevnikov_et_al-2021.} From equations (ref) and (ref) we can write
Our interest lies in proving the asymptotic properties of $\hat{\beta}$ taking advantage of the expression (ref). However, the presence of $\varepsilon_{M,m}$ in $Y_{M,m}$, which is allowed to be correlated along the network over active dyads, renders our approach nonstandard. The canonical results about independently and identically distributed ($i.i.d.$) random variables are by no means applicable due to the complex form of dependence inherent to $\varepsilon_{M,m}$. Dyadic-robust estimators (e.g., Aronow_et_al-2015 and Tabord-Meehan-2019) only account for correlations between adjacent dyads. This may be appropriate in some settings, but not when there may be indirect spillovers, as in specification (ref) -- (ref). Thence, the right hand side of (ref) cannot simply be embedded in other widely-studied variants of dependent random vectors, including spatially-correlated and time-series data.
The main insight of this paper is that the spillovers across connected -- even if not adjacent -- dyads, can be rewritten as the dependence of $Y_{M,m}$'s along the network of active dyads (hereby, referred to as the “network”). This transformation allows us to embrace such complex cross-sectional dependence and appropriately rewrite the problem so that recent results on network dependent random variables Kojevnikov_et_al-2021 can be applied. Our results then allow for correlations of random variables given a common shock as explicitly dictated by distances between them on the dyadic network. Our notion of dependence restricts the covariance between the random variables for any two sets of dyads $A$ and $B$, $Y_{M,A}$ and $Y_{M,B}$, which are at a distance $s$ from one another, to be bounded and to decrease as $s$ goes to infinity. Indeed, this covariance is controlled by a sequence of (random) coefficients $\theta_{M,s}$, which capture the strength of covariation net of scaling, analogous to correlation coefficients. As the minimum distance, $s$, between $A$ and $B$ grows, the dependence $\{\theta_{M,s}\}$ between $Y_{M,A}$ and $Y_{M,B}$, tends to zero. A formal description is provided in Appendix (ref).
As will become transparent shortly, asymptotic theories for $\hat{\beta}$ rest on tradeoffs between the correlation of the network-dependent random vectors (i.e., the dependence coefficients) and the denseness of the underlying network. To measure the denseness, we first define two concepts of neighborhoods: for each $m\in\mathcal{M}_{N}$ and $s\in\mathbbm{N}\cup \{0\}$,
where $\rho_{M}(m,m')$ denotes the geodesic distance between dyads $m$ and $m'$.\footnote{Recall that we define the geodesic distance between two connected dyads $m$ and $m'$ to be the smallest number of adjacent dyads between them.} The former set collects all the $m$'s neighbors whose distance from $m$ is no more than $s$ (which we call a neighborhood), whilst the latter registers all the $m$'s neighbors whose distance from $m$ is exactly $s$ (which we call a neighborhood shell).
Next, we define two types of density measures of a network: for $k, r >0$,
where it is assumed that $\mathcal{M}_{N}(m'; -1) = \emptyset$. The former measure gauges the denseness of a network in terms of the average size of a version of the neighborhood. The latter expresses the denseness of a network as the average size of the neighborhood shell raised by $k$. Kojevnikov_et_al-2021 show that controlling the asymptotic behavior of an appropriate composite of these two measures (denoted by $c_M$ and defined in Assumption (ref)) is sufficient for the Law of Large Numbers (LLN) and Central Limit Theorem (CLT) of the network dependent random variables (Condition ND).
It is worth noting that there are two different units at play: the number of sampling units (i.e., individuals), $N$, and the number of dyads, $M$. To bridge these two units, we assume that the number of active dyads also goes to infinity as the number of sampling units is taken to infinity, eliminating the possibility of extremely sparse networks among individuals. This is empirically relevant. For example, in international trade, the entry of a new country/firm to a market will most likely increase the number of trade flows in the economy; in political economy, the more members of parliaments (MEPs) there are, the more pairs of the MEPs sitting next to each other there will be (see Section (ref)).\footnote{This assumption is similar in spirit to Assumption 2.3 of Tabord-Meehan-2019 in which the minimum degree is assumed to grow at some constant rate relative to the number of individuals. It is milder than Assumption 2.3 of Tabord-Meehan-2019 since the latter does not allow any individual to be isolated, while Assumption (ref) merely constrains the average degree. Similarly, this assumption is weaker than the assumption that the maximum degree in a network is bounded even when $N\rightarrow\infty$ (e.g., Penrose_and_Yukich-2003 and de_Paula_et_al-2018).}
We make use of the following two regularity assumptions for the proof of consistency of $\hat \beta$ for $\beta$ (Theorem (ref)) and to derive its asymptotic distribution (Theorem (ref)).\footnote{These assumptions are required for Theorem (ref), but as usual, the proof of consistency (Theorem (ref))can be derived under weaker conditions. (See Assumptions (ref) and (ref) and their associated discussion, in Appendix).}
Assumption (ref) requires that the errors are not too large once conditioned on common shocks. Together with the standard full rank assumption that guarantees identification of $\beta$\footnote{See Assumption (ref)(c) in Appendix, for completeness.}, this implies Assumption 3.1 of Kojevnikov_et_al-2021 for each $u$-th element of $Y_{M,m}^{}$, denoted by $Y_{M,m}^{u}$ with $u\in\{1,\ldots, K\}$.
Assumption (ref) is a condition that controls the tradeoff between the denseness of the underlying network and the covariability of the random vectors. If the network becomes dense, then the dependence of the associated random variables has to decay much faster. This requirement is consistent with the typical assumption in the spatial econometrics literature upon which the empirical analyses of the aforementioned examples (i.e., international trade, legislative voting, and exchange rates) are based, because it embodies the idea that spillovers decay as they propagate farther (see, e.g., Kelejian_and_Prucha-2010). For instance, Acemoglu_et_al-2015 assume that network spillovers are zero if agents are sufficiently distantly connected on a geographical network. This decay in the propagation of spillovers through the terms $\theta_{M,\cdot}$. This assumption may be violated for very dense networks with low decay of spillovers.
Theorem (ref) establishes the consistency of $\hat{\beta}$ under the scenario where the number of sampling units $N$ goes to infinity. When Assumption (ref) is dropped, the result continues to hold in terms of the number of active dyads $M$.
Let $S_{M}^{} \coloneqq \sum_{m\in\mathcal{M}_{N}} Y_{M,m}^{}$, which is present in $\hat \beta$ in equation ((ref)). Let $S_{M}^{u}$ be the $u$-th entry of $S_{M}$ for $u\in\{1,\ldots,K\}$ and denote the unconditional variance of $S_{M}^{u}$ by $\tau_{M}^{2} \coloneqq Var (S_{M}^{u} )$. Since $S_{M}^{u}$ is not a sum of independent variables, its variance cannot be simply expressed as a sum of the variances of $Y_{M,m}$. We thus need to explicitly take into account covariance between the random variables $\{Y_{M,m}^{u}\}_{m\in\mathcal{M}_{N}}$. We study the CLT for the normalized sum of $Y_{M,m}^{u}$, which is given by $\frac{S_{M}^{u}}{\tau_{M}}$.
Assumptions (ref) and (ref) are written in terms of conditional expectations, whereas we are interested in the unconditional distribution of $\frac{S_{M}}{\tau_{M}}$. Assumption (ref) below bridges the conditional and unconditional variances. The final assumption for the asymptotic normality result is a standard regularity condition guaranteeing that the asymptotic variance is well-defined. This assumption follows from both matrices in the expression being well-defined (e.g., $x_{M,k}$ having bounded support).\footnote{Further note that Theorem (ref) is proved under a weaker condition than Assumption (ref).}
Under these assumptions, the asymptotic distribution of $\hat{\beta}$ is given by:
Our objective is to consistently estimate $AVar(\hat{\beta})$ defined in Theorem (ref). As errors are mean zero, $Y_{M,m}$ is centered, i.e., $E\big[ Y_{M,m} \big]=0$ for each $m\in \mathcal{M}_{N}$.
The proposed estimator is a type of kernel estimator. Let $b_{M}$ denote the bandwidth, or the lag truncation (its choice is described in Section (ref) below) and $\omega:\overline{\mathbbm{R}}\rightarrow[-1,1]$ a kernel function such that $\omega(0)=1$, $\omega(z)=0$ whenever $|z|>1$, and $\omega(z)=\omega(-z)$ for all $z\in\overline{\mathbbm{R}}$. The feasible variance estimator of interest is {
} with $\omega_{M}(s) \coloneqq \omega\big(\frac{s}{b_{M}}\big)$ for all $s \geq 0$ and $\hat{Y}_{M,m} \coloneqq x_{M, m} \hat{\varepsilon}_{M, m}$, where $\hat{\varepsilon}_{M, m}\coloneqq y_{M, m} - x_{M, m}' \hat{\beta}$.
There are two approaches for the choice of the associated lag truncation parameter. First, the researcher may already know (or is willing to impose) the truncation, perhaps due to a theoretical/institutional motivation. For instance, Acemoglu_et_al-2015 set the lag to two in a related problem. Then, the thought exercise is that this choice will adapt as $M \to \infty$ according to the assumptions below. Alternatively, the researcher could use a data-driven choice. Assumption (ref) (c) below suggests it should depend on both the sample size and the network topology, including the average degree of the dyadic network. One such selection rule is suggested in Kojevnikov_et_al-2021 based on their proofs: $b_M = 2 \log(M)/\log(\max(average~degree, 1.05))$.
The consistency of the variance estimator requires two sets of additional assumptions. The first set is Assumption 4.1 of Kojevnikov_et_al-2021, but stated here in terms of the network over dyads.
Assumption (ref) (a) is a stronger counterpart to Assumption (ref), as it requires that a higher-order (i.e., higher than fourth order) conditional moment be well-defined. Assumption (b) posits a tradeoff between the kernel function, the denseness of a network and the dependence coefficients. Specifically, the kernel function $\omega_{M}$ is required to converge to one sufficiently fast. Kojevnikov_et_al-2021 demonstrate primitive conditions under which this requirement is fulfilled (Proposition 4.2). Assumption (c) requires that the correlation coefficients decay much faster relative to the denseness of the network. This is satisfied in the suggested choice for $b_M$ above.
Another set of conditions restricts the denseness of the network, ruling out the situation where the network becomes progressively dense: most notably, the case where every single individual unit is directly linked to every other individual.
The following theorem is the main theoretical contribution of this paper.
Theorem (ref) establishes the consistency of our proposed variance estimator accounting for network spillovers across dyads in the sense of the Frobenius norm.
It follows from Theorem (ref) that the dyadic-robust variance estimator ((ref)) is inconsistent for the true variance when the underlying network involves a non-negligible degree of far-away correlations, as suggested in the examples of the previous section.\footnote{If the network is such that there are only adjacent dyads (i.e., when equation ((ref)) holds), then the result above implies consistency of this estimator for dyadic dependence. By comparison, Lemma 1 of Aronow_et_al-2015 and Proposition 3.1 and 3.2 of Tabord-Meehan-2019 also provide consistent variance estimators for the dyadic dependence case without higher-order network spillovers. However, these results and ours do not subsume one another. Indeed, their estimators can accommodate flexible dependence within clusters of dyads that share common units, while we assume that the network of spillovers is observed even if there are only adjacent connections.} Specifically, the following corollary states that the dyadic-robust variance estimator of Aronow_et_al-2015 may not necessarily be consistent when it is na\"ively applied to the network-regression model with non-zero correlations beyond direct neighbors.
The added condition (ref) in Corollary (ref) pertains to both the network configuration of active dyads and the regression variables. It represents a setting where the spillovers from far-away neighbors are non-negligible even when $N$ is large. For instance, (ref) can hold even if there are not many neighbors, as long as the covariances between the error terms are sufficiently large. This builds on pinto that inference with dyadic data may be biased if one only partially accounts for such spillovers. On the other hand, if far-away neighbors in the network have only a negligible effect on the cross-sectional dependence, the dyadic-robust variance estimator remains a good approximation for the asymptotic variance of linear dyadic data models with network spillovers across dyads. These insights are investigated further in Section (ref) using numerical simulations.
These observations extend to settings where the spillovers decay along the network (i.e., when the correlation along unobservables decreases with the geodesic distance among dyads). Indeed, our estimator already accounts for such decay through the indirect covariances in its expression ((ref)). When such spillovers propagate and decay is not too high, then condition ((ref)) is satisfied, as such spillovers are non-neglible. On the other hand, if they decay at a very high rate (in the limit, a $100\%$ decay from adjacent to connected dyads), then our estimator will become very similar to the dyadic-robust variance estimator.
However, some researchers may be willing to tolerate some asymptotic bias to still implement the dyadic-robust estimator even with some asymptotic bias. In such cases, when should they prefer our proposed estimator? While a general answer is complex because the bias depends on both the strength of indirect spillovers and on the network configuration, we provide an intuitive criterion which is verified using: (i) a simple, yet informative, analytical example where spillovers decay exponentially with distance along the network (shown below), and (ii) extensive Monte Carlo simulations with decay and different network configurations (shown in Section (ref)). They illustrate that researchers should opt for the network-robust estimator above when they have a low tolerance for bias, when there are persistent indirect spillovers and/or when the network is not too sparse.
When the tolerated bias is small ($B \to 0$), dependence is large enough and does not decay too fast (large $\gamma$), or the network is more dense, our approach is preferable because it provides a consistent estimator even with non-negligible spillovers. Since the network is observed, $S$ and the last term are estimable and can be used for such diagnostics. Appendix (ref) provides further discussions, while the next section presents for our baseline results and discusses simulation exercises.
We compare three types of variance estimators across different specifications and network configurations. We use the Eicker-Huber-White estimator as a benchmark,\footnote{It is used in bliss98 and mansfield00, for instance.} the dyadic-robust estimator of Tabord-Meehan-2019 as a comparison accounting for the dyadic nature of the data (when inappropriately used in the presence of network spillovers), and our proposed estimator which is robust to network spillovers across dyads.
We first generate networks on which random variables are assigned. We follow Canen_et_al-2020b among others by employing two models of random graph formations. They are referred to as Specifications 1 and 2. Specification 1 uses the \citepos{Barabasi_and_Albert-1999} model of preferential attachment, with the fixed number of edges $\nu \in\{1,2,3\}$ being established by each new node.\footnote{In generating the Barab\'asi-Albert random graphs, we follow Canen_et_al-2020b by choosing the seed to be the Erd\"os-Renyi random graph with the number of nodes equal the smallest integer above $5\sqrt{N}$, where $N$ denotes the number of nodes.} Specification 2 is based on the Erd\"os-Renyi random graph Erdos_and_Renyi-1959, Erdos_and_Renyi-1960 with probability $p= \frac{\lambda}{N}$ for $N$ denoting the number of nodes and $\lambda\in\{1,2,3\}$ being a parameter that governs the probability relative to the node size. The summary statistics for the networks generated by Specification 1 and 2 are given in Appendix (ref). The maximum degree and the average degree increase monotonically as we increase the parameters in both specifications. The number of active edges (i.e., dyads) also increases with the sample size regardless of the specification. This reflects that each node tends to have more direct links as the network becomes denser. In our exercises, the number of indirectly linked dyads also increases with network denseness. However, this is due to our simulated networks being relatively sparse. In other settings, the number of indirect connections may decrease with network density.
For each of the randomly generated networks, the simulation data is generated from the following simple network-linear regression:
\footnotetext{To simplify notation, we drop the $M$ subscript, making the triangular array structure implicit.} with $m\coloneqq d(i,j)$ representing the dyad between agent $i$ and $j$. The dyad-specific regressor $x_{m}$ is defined as $x_{m}\coloneqq |z_{i} - z_{j}|$, where both $z_{i}$ and $z_{j}$ are drawn independently from $\mathcal{N}(0,1)$. The regression coefficient is fixed to $\beta=1$ across specifications.
The dyad-specific error term $\varepsilon_{m}$ is constructed to exhibit non-zero correlation with $\varepsilon_{m'}$ as long as dyads $m$ and $m'$ are connected (i.e., in the network terminology, there exists a path in the simulated network), while the strength of the correlation is assumed to decay as they are more distant. To that end, we draw $\varepsilon_{m} \coloneqq \sum_{m'} \gamma_{m,m'} \eta_{m,m'}$, where $\gamma_{m,m'}$ equals $\gamma^{s}$ if the distance between $m$ and $m'$ is $s$, and $0$ otherwise, for $\gamma\in[0,1]$\footnote{In this simulation, we focus on cases of positive spillovers, as negative spillovers can be analyzed analogously.} and $s\in\{1,\ldots,S\}$ with $S$ being the maximum geodesic distance that the spillover propagates to. Each $\eta_{m,m'}$ is drawn $i.i.d.$ from $\mathcal{N}(0,1)$. Hence, $\gamma$ controls the strength of spillover effects, representing their decay rate. If $\gamma=1$, then spillover effects are the same no matter how far the agents are apart, i.e., the spillover effects do not decay. If $\gamma = 0$, there are no spillover effects, so the dyadic-robust variance estimator should be consistent. The case of $S=2$ corresponds to a situation where up to friends of friends may matter for spillovers.
We consider three scenarios for each type of network. In the main text, we set $S=2$ and $\gamma=0.8$. The results for $S=2$ with $\gamma=0.2$ are given in Appendix (ref), and the ones for $S=1$ with $\gamma=0.8$ are in Appendix (ref). Throughout the experiments, we employ the mean-shifted (by one) rectangular kernel with the lag truncation equal to two for the sake of comparison.\footnote{These choices for the kernel and lag truncation correspond to placing unit weights on pairs of adjacent dyads and pairs of connected dyads of distance equal two, and zero otherwise. This combination is sufficient for us to accommodate all the correlations across the network since it is constructed to have $S=2$. This provides a transparent comparison to other estimators. However, it is rare for a researcher to have such information. Appendix (ref) provides further exercises where $S>2$ (so spillovers reach further dyads) and where the lag-truncation is chosen adaptively by the rule recommended in Section (ref). Those results are discussed below.}
In Table (ref) we present the coverage probability for $\beta$ and the average length of the confidence interval across simulations. To do so, we compute the $t$-statistic using the OLS estimator for $\beta$ and different variance estimators under a Normal distribution approximation.\footnote{It is well known that estimates of a variance-covariance matrix may be negative semidefinite when the sample size is very small. This occurs in four out of 5000 simulations when $N=500$. Rather than dropping such observations, we follow Cameron_et_al-2011 and augment the eigenvalues of the matrix by adding a small constant, say $0.005$, thereby obtaining a new variance estimate that is more conservative.} The finite-sample properties of the three variance estimators are further illustrated in Figure (ref) in Appendix (ref).
The results for the empirical coverage probabilities depend on two dimensions: the sample size ($N$) and the denseness of the underlying network (parametrized by $\nu$ and $\lambda$). The coverage probability for each estimator improves with the sample size. However, when spillovers are high ($\gamma = 0.8$), only our proposed network-robust variance estimator has coverage close to 95%, consistent with Theorem (ref). Meanwhile, in this set-up, both the Eicker-Huber-White and the dyadic-robust variance estimators perform poorly as the underlying network becomes denser, no matter which specification of the network is involved. For example, in Specification 1 with $\nu=3$ and the largest sample size ($N=5000$), the confidence intervals based on the Eicker-Huber-White and the dyadic-robust variance estimators do not cover the true parameter 615 and 455 times out of 5000 simulations (12.3% and 9.1%), respectively. On the other hand, the network-robust variance estimator is designed to capture higher-order correlations and, thus, its coverage remains stable across network configurations.
A similar conclusion is drawn from the average length of the confidence intervals. The confidence intervals for the Eicker-Huber-White and dyadic-robust variance estimators are typically shorter than those for our proposed estimator when $\gamma$ is large and $S=2$. This means the former undercovers the true parameter (in the presence of positive spillovers), confirming our theoretical results (Corollary (ref)). The difference is often around 10 -- 20% of the length of the network-robust estimator and it increases with even denser networks (see Table (ref) in Appendix).
However, as the magnitude of spillovers decreases (i.e. $\gamma$ tends to zero), higher-order spillovers are less pronounced, so that the biases from using the Eicker-Huber-White and dyadic-robust variance estimators disappear. This is shown in Table (ref) of Appendix (ref) for the case of $S=2$ and $\gamma = 0.2$. When $S=1$, the dyadic-robust variance estimator coincides with our proposed estimator (i.e., there are no spillovers from non-adjacent links, or spillovers fully decay immediately). This is shown in Table (ref) of Appendix (ref).
Finally, Appendix (ref) shows that the results are robust to spillovers that can reach the most distantly connected neighbors ($S=\infty$) and to choosing the lag-truncation adaptively. Rather than assuming a given value, we set $b_M = 2 \log(M)/\log(\max(average~degree, 1.05))$ as suggested above. The results are qualitatively and quantitatively very similar: when the decay is not too high (i.e., spillovers propagate), the confidence interval associated with the dyadic-robust estimator will undercover the true parameter even in moderate sample sizes.
We now turn to an empirical application to demonstrate the performance of our variance estimator with real-world data. In doing so, we revisit the work of Harmon_et_al-2019 on whether legislators who sit next to each other in Parliament tend to vote more alike on policy proposals.
They focus on the European Parliament, whose Members (MEPs) are voted in through elections in each European Union (EU) member country every five years. The Parliament convenes once or twice a month, in either Brussels or Strasbourg, to debate and vote on a series of proposals. Once elected to the European Parliament (EP), these MEPs are organized into European Political Groups (EPGs), which aggregate similar ideological members/parties across countries. As Harmon_et_al-2019 describe, these EPGs function as parties for many of the traditional party-functions in other legislatures, including coordination on policy and policy votes. Most importantly, MEPs sit within their EPG groups in the chamber. However, within each EPG group, non-party leaders traditionally sit in alphabetical order by last name. See Figure (ref) in Appendix (ref) for an example. Hence, to the extent that last names are exogenous (i.e., politicians are not selected or choose last names anticipating seating arrangements in the European Parliament), seating arrangements are used as a quasi-natural experiment.
We adopt the same dataset used in Harmon_et_al-2019, which collects the MEP-level data on votes cast in the EP. The dataset records what each MEP voted for (Yes or No), where she was seated, and a number of individual characteristics (e.g. country, age, education, gender, tenure, etc). Its sample period is the plenary sessions held in both Brussels and Strasbourg between October 2006 and November 2010. This corresponds to the sixth and seventh terms of the European Parliament.
For our empirical illustration, we restrict the sample to the policies voted in Strasbourg during the seventh term and we focus on the seating pattern between July 14th -- July 16th, 2009 (which involved 116 different proposals being voted on). The resulting sample has 2,431,261 observations, which are split into 422 politicians forming 26,099 pairs (i.e., dyads) of MEPs over 116 proposals.\footnote{There are 334 pairs of adjacent dyads and 591 pairs of connected dyads.} Further information on the construction of our sample is detailed in Appendix (ref). These restrictions keep the main set-up in Harmon_et_al-2019, while allowing us to evaluate variance estimators with smaller sample sizes.
Let $Agree_{d(i,j), t}$ be an indicator that takes one if MEP $i$ and $j$ cast the same vote on proposal $t$, and zero otherwise. Likewise, $SeatNeighbors_{d(i,j), t}$ is a binary variable that equals one if MEP $i$ and $j$ are seated next to each other when the vote for proposal $t$ is taken place, and zero otherwise. We view the seating arrangement as a network among the MEPs: i.e., two MEPs who are seated next to each other within the same political group are treated as an active dyad, which in turn constitutes a network over pairs of MEPs within each political group (see the note below Figure (ref) for details).\footnote{This accommodates the row-by-EP-by-EPG clustering implemented by the authors. This is already an extension to Harmon_et_al-2019 which assume that only seating neighbors would have correlated error terms. One could extend this further by allowing connections across parties.}
We follow Harmon_et_al-2019 in assuming that such seating arrangements are exogenously determined, i.e., the adjacency relation is exogenously formed. Their main specification is a linear model:
The authors originally conducted inference using the estimator in Aronow_et_al-2015, assuming that dyads $m = d(i,j)$ cannot be correlated with $m' = d(k,l)$ unless they share a common member (i.e., there is no correlation across errors if $i,j$ and $k,l$ do not include a common unit). We will now compare this approach to using the variance estimator introduced in Section (ref), which allows the error terms to exhibit non-zero correlations as long as they are connected on the network over dyads represented by the adjacency relation of seating arrangements in Parliament. In this analysis, we use the mean-shifted rectangular kernel with the lag truncation equal the longest path in the constructed network. This combination of the kernel density and the lag truncation enables us to accommodate all the possible correlations across connected dyads (i.e., pairs of MEPs), placing equal weight on each of them. This makes the comparisons across estimators more transparent.\footnote{In Appendix (ref), we replicate this analysis with a different choice of kernel and setting the lag-truncation parameter following the criterion suggested above/in Kojevnikov_et_al-2021. The results are very similar.}
Inspired by Harmon_et_al-2019, we consider three specifications: (I) a simple linear regression model as given in ((ref)); (II) the model ((ref)) augmented with a flexible set of other demographic variables;\footnote{Following Harmon_et_al-2019, we include indicators whether country of origins, quality of education, freshman status and gender, respectively, are the same, as well as differences in ages and tenures. See the note below Table (ref) for details.} and (III) the model ((ref)) with both a flexible set of other demographic variables and day-specific fixed effects. When fixed effects are present in their original estimation, we estimate a within-difference model via OLS.
The main results of our empirical analysis are summarized in Table (ref). Panel A displays the parameter estimates for the three different specifications. This panel shows that our point-estimates are consistent with the original estimates of Harmon_et_al-2019 (columns 6 and 7 of Table 4), as they are close to 0.006 (their original results) and stable across specifications.\footnote{Note that our dependent variable is equal to one if two MEPs vote the same and zero otherwise, while Harmon_et_al-2019 code it as one if MEPs vote differently. Hence, to compare our estimates with theirs, the signs on the estimates of $SeatNeighbors$ must be flipped.} Hence, changes to point-estimates are not due to sample selection. The positive coefficient for {\sl SeatNeighbors} suggests that the MEPs sitting together tend to vote more similarly than those sitting apart, providing evidence in favor of their original hypothesis. The coefficients on the covariates (displayed in Panel C of Appendix Table (ref)) are also quantitatively and qualitatively similar to those in their original paper. For instance, our estimates for {\sl SameCountry} are 0.056, while their estimates are around 0.051, suggesting politicians from the same country are more likely to vote similarly on policies.
Panel B shows the standard errors for the regression coefficient of {\sl SeatNeighbors} using different variance estimators. Building on Section (ref), the panel compares three different types of heteroskedasticity robust estimators: namely, the Eicker-Huber-White, the dyadic-robust estimator (used in their original work), and our proposed network-robust estimator. We note that, due to the smaller sample size, the standard errors in our exercise are larger than the original authors'. Hence, the results in this section are not directly comparable to those in Harmon_et_al-2019. Rather, we compare the different estimators within our own sample.
As foreshadowed in the Monte Carlo simulations, the Eicker-Huber-White estimates are the smallest, followed by the dyadic-robust estimates, which, in turn, are smaller than the network-robust estimates. In fact, for Specification (III), the Eicker-Huber-White estimate is roughly 73% smaller than using the estimator accounting for network spillovers across dyads, while the dyadic-robust one is 22% smaller. This fact entails two implications. First, our finding provides empirical evidence in support of the existence of indirect positive spillovers among the MEPs: even distant connections may indirectly generate correlated behavior among politicians $i$ and $j$. Second, the use of alternative estimators not accounting for such spillovers undercovers the true parameter and may generate biased hypothesis testing about the regression coefficient of {\sl SeatNeighbors}. The difference in estimates appears quantitatively meaningful in this empirical example.
With dyadic data, researchers typically assume that dyads are uncorrelated if they do not share a common unit, an assumption that is leveraged in inference. We showed that this assumption may be inappropriate in many models where interactions occur on a network: while data may be dyadic, the cross-sectional dependence may be much more complex and spillover beyond pairwise interactions. For instance, trade between countries may depend on trade between auxiliary partners and their unobservables. In political economy, whether a politician votes with a colleague may have spillovers from seating neighbors beyond one's own immediate ones. We verified this using both theoretical results and Monte Carlo simulations. We provided a new consistent variance estimator for parameters in a linear model with dyadic data with correct asymptotic coverage and good finite sample properties.
To conclude, we clarify that our goal in this exercise is neither to criticize dyadic-robust variance estimators, which are a fundamental part of the empiricist's toolkit, nor to suggest our approach should always be used. Rather, we wish to draw attention that researchers should fully specify the cross-sectional dependence in their model. If the conventional assumption of dyadic dependence correctly specifies the environment in question, or when spillovers beyond immediate neighbors might be negligible, then previous approaches suffice. However, as we have discussed above, existing applications may apply the latter method even if it is seemingly inappropriate to their setting. This includes situations where such network spillovers may be present or persistent (even with decay). In such scenarios, our estimator provides a possible solution. Those choices, though, must be guided by the application that empiricists face. Hence, building on poast, we recommend researchers to continue to fully specify their model, including full specification of their covariance structure, thereby clarifying what type of inference procedure is most appropriate for their environment.
\singlespacing \setlength\bibsep{0pt}
We thank Aimee Chin, Hugo Jales, Taisuke Otsu, Pablo Pinto, V\'{i}tor Possebom, Kevin Song, Bent E. S\o rensen, and seminar participants at the University of Houston, the 2022 Texas Econometrics Camp, the 2022 European Winter Meeting and the 2022 Asia Meeting of the Econometric Society in East and South-East Asia for valuable comments and suggestions.
Replication code and data for this article will be made available publicly online upon conditional acceptance.
Competing interests: The authors declare none.
\onehalfspacing