Regression of exchangeable relational arrays

TL;DR

This paper introduces a new regression method for relational arrays that accounts for actor dependence using exchangeability, improving inference accuracy over existing methods.

Contribution

It develops a novel class of parsimonious estimators leveraging exchangeability to better handle dependence in relational array regressions.

Findings

Standard error estimators are more accurate with exchangeability assumptions.

Simulation studies show improved inference accuracy.

Application to international trade data demonstrates practical effectiveness.

Abstract

Relational arrays represent measures of association between pairs of actors, often in varied contexts or over time. Trade flows between countries, financial transactions between individuals, contact frequencies between school children in classrooms, and dynamic protein-protein interactions are all examples of relational arrays. Elements of a relational array are often modeled as a linear function of observable covariates. Uncertainty estimates for regression coefficient estimators -- and ideally the coefficient estimators themselves -- must account for dependence between elements of the array (e.g. relations involving the same actor) and existing estimators of standard errors that recognize such relational dependence rely on estimating extremely complex, heterogeneous structure across actors. This paper develops a new class of parsimonious coefficient and standard error estimators for regressions of relational arrays. We leverage an exchangeability assumption to derive standard error estimators that pool information across actors and are substantially more accurate than existing estimators in a variety of settings. This exchangeability assumption is pervasive in network and array models in the statistics literature, but not previously considered when adjusting for dependence in a regression setting with relational data. We demonstrate improvements in inference theoretically, via a simulation study, and by analysis of a data set involving international trade.

Figures

Tables

Table 1. Table S1: Twelve index patterns that form a set of equivalent bivariate distributions with the entry in Ξ Ξ \Xi corresponding to indices ( i , j , r ) 𝑖 𝑗 𝑟 (i,j,r) . The row labels correspond to those in Figure  1 and the column labels indicate index sets that share the same third dimension “slice” of Ξ Ξ \Xi with index set ( i , j , r ) 𝑖 𝑗 𝑟 (i,j,r) .

	$s=r$	$s\ne r$
(0)	$(i,j,r)$	$(i,j,s)$
(a)	$(j,i,r)$	$(j,i,s)$
(b)	$(i,l,r)$	$(i,l,s)$
(c)	$(k,j,r)$	$(k,j,s)$
(d)	$(k,i,r)$	$(k,i,s)$
(e)	$(k,l,r)$	$(k,l,s)$

Table 2. Table S2: Approximate standard deviations for exchangeable error setting.

${\sigma }_{ϵ}$	${\sigma }_a$	${\sigma }_b$	${\sigma }_{\gamma }$	${\sigma }_z$
0.866	0.957	0.677	0.677	0.677

Table 3. Table S3: Mean (across values of n 𝑛 n ) ratio of average (across covariate and error realizations) absolute bias in dyadic clustering variance estimator to average absolute bias of exchangeable variance estimator for independent and exchangeable error models.

	binary	positive	real-valued
independent	96.48	5.03	123.45
exchangeable	3.04	5.43	5.02

Table 4. Table S4: Estimated probability of a Type I error (size) and one less the probability of a Type II error (power) for the proposed permutation test.

Size	Power
0.025	0.839

Equations

$$y_{ijr}={\beta }^T{x}_{ijr}+{\xi }_{ijr},(i,j=1,\dots ,m;i\ne j;r=1,\dots ,R),$$

$${\hat{\beta }}_{GLS}=(X^T{\Omega }^{-1}X{)}^{-1}{X}^T{\Omega }^{-1}y.$$

$$\text{var}\left(\hat{\beta }\mid X\right)$$

$$\text{var}\left({\hat{\beta }}_{GLS}\mid X\right)$$

$${\hat{V}}_{DC}=(X^{T}X{)}^{-1}{X}^T{\hat{\Omega }}_{DC}X(X^{T}X{)}^{-1}.$$

$${\hat{V}}_E=(X^{T}X{)}^{-1}{X}^T{\hat{\Omega }}_{E}X(X^{T}X{)}^{-1},{\hat{\Omega }}_E=\sum _{\eta =1}^2\sum _{u=0}^d{\hat{\varphi }}_u^{(\eta )}{\mathcal{S}}_u^{(\eta )},$$

$$\displaystyle\widehat{\phi}_{b}^{(2)}=\binom{R}{2}^{-1}\frac{1}{n(n-1)(n-2)}\sum\limits_{r\neq s}\sum\limits_{i}\sum\limits_{j\neq i}e_{ijr}\Big{(}\sum\limits_{k\neq j}e_{iks}-e_{ijs}\Big{)}.$$

$$y_{ij}={\beta }_1+x_{ij}{\beta }_2+{\xi }_{ij}.$$

$$n^2\text{Bias}({\hat{V}}_{DC})+O(n^{-1/2})\le -2n^2|\text{Bias}({\hat{V}}_E)|+O(n^{-1/2})\le 0,$$

$${\hat{\varphi }}_0^{(1)}=\frac{1}{Rn(n-1)}\sum _r\sum _i\sum _{j\ne i}{e}_{ijr}^2,$$

$${\hat{\varphi }}_a^{(1)}=\frac{1}{Rn(n-1)}\sum _r\sum _i\sum _{j\ne i}{e}_{ijr}{e}_{jir},$$

$$\displaystyle\widehat{\phi}_{b}^{(1)}=\frac{1}{Rn(n-1)(n-2)}\sum_{r}\sum\limits_{i}\sum\limits_{j\neq i}e_{ijr}\Big{(}\sum\limits_{k\neq i}e_{ikr}-e_{ijr}\Big{)},$$

$$\displaystyle\widehat{\phi}_{c}^{(1)}=\frac{1}{Rn(n-1)(n-2)}\sum_{r}\sum\limits_{i}\sum\limits_{j\neq i}e_{ijr}\Big{(}\sum\limits_{k\neq j}e_{kjr}-e_{ijr}\Big{)},$$

$$\displaystyle\widehat{\phi}_{d}^{(1)}=\frac{1}{2Rn(n-1)(n-2)}\sum_{r}\sum\limits_{i}\sum\limits_{j\neq i}e_{ijr}\Big{(}\sum\limits_{k\neq i}e_{kir}+\sum\limits_{k\neq j}e_{jkr}-2e_{jir}\Big{)},$$

$${\hat{\varphi }}_0^{(2)}={\left(\genfrac{}{}{0ex}{}{R}{2}\right)}^{-1}\frac{1}{n(n-1)}\sum _{r\ne s}\sum _i\sum _{j\ne i}{e}_{ijr}{e}_{ijs},$$

$${\hat{\varphi }}_a^{(2)}={\left(\genfrac{}{}{0ex}{}{R}{2}\right)}^{-1}\frac{1}{n(n-1)}\sum _{r\ne s}\sum _i\sum _{j\ne i}{e}_{ijr}{e}_{jis},$$

$$\displaystyle\widehat{\phi}_{b}^{(2)}=\binom{R}{2}^{-1}\frac{1}{n(n-1)(n-2)}\sum_{r\neq s}\sum\limits_{i}\sum\limits_{j\neq i}e_{ijr}\Big{(}\sum\limits_{k\neq i}e_{iks}-e_{ijs}\Big{)},$$

$$\displaystyle\widehat{\phi}_{c}^{(2)}=\binom{R}{2}^{-1}\frac{1}{n(n-1)(n-2)}\sum_{r\neq s}\sum\limits_{i}\sum\limits_{j\neq i}e_{ijr}\Big{(}\sum\limits_{k\neq j}e_{kjs}-e_{ijs}\Big{)},$$

$$\displaystyle\widehat{\phi}_{d}^{(2)}=\binom{R}{2}^{-1}\frac{1}{2n(n-1)(n-2)}\sum_{r\neq s}\sum\limits_{i}\sum\limits_{j\neq i}e_{ijr}\Big{(}\sum\limits_{k\neq i}e_{kis}+\sum\limits_{k\neq j}e_{jks}-2e_{jis}\Big{)}.$$

$$\text{var}({\xi }_{ijr})={\theta }_0^{(1)},\text{cov}({\xi }_{ijr},{\xi }_{kir})={\theta }_a^{(1)},\text{cov}({\xi }_{ijr},{\xi }_{klr})=0,$$

$$\text{cov}({\xi }_{ijr},{\xi }_{ijs})={\theta }_0^{(2)},\text{cov}({\xi }_{ijr},{\xi }_{kis})={\theta }_a^{(2)},\text{cov}({\xi }_{ijr},{\xi }_{kls})=0,r\ne s.$$

$$\text{var}({\xi }_{ijr})=\text{var}\{{\xi }_{\pi (i)\pi (j)\nu (r)}\}.$$

$$\text{cov}({\xi }_{ijr},{\xi }_{klr})=\text{cov}\{{\xi }_{\pi (i)\pi (j)r},{\xi }_{\pi (k)\pi (l)r}\}.$$

$$({\xi }_{ijr},{\xi }_{ijs}),(i=1,\dots n;j=1,\dots ,n;i\ne j;r=1,\dots ,R;s=1,\dots ,R;r\ne s),$$

$$({\xi }_{ijr},{\xi }_{ilr}),$$

$$({\xi }_{ijr},{\xi }_{ils}),$$

$$\mathcal{D}=\text{span}\left\{M_{jk,lm}:(j,k)\cap (l,m)=\varnothing \right\}.$$

$$\mathcal{E}=\text{span}\left({\mathcal{S}}_i:i=1,\dots ,5\right).$$

$$({\hat{V}}_{DC}{)}_{kk}$$

$$=f(n{)}^{-2}\text{tr}\left\{(y-X\hat{\beta })(y-X\hat{\beta }{)}^T{P}_{\mathcal{D}}\left(X_k{X}_k^T\right)\right\},$$

$$0\le {\varphi }_3+{\varphi }_4\le 1,$$

$${\varphi }_3^2+{\varphi }_5^2\le {\varphi }_3,$$

$${\varphi }_4^2+{\varphi }_5^2\le {\varphi }_4.$$

$${\varphi }_3^2+{\varphi }_4^2+2{\varphi }_5^2\le 1,$$

$$n({\varphi }_3+{\varphi }_4+2{\varphi }_5)$$

$$1+{\varphi }_2-{\varphi }_3-{\varphi }_4-2{\varphi }_5$$

$$1-{\varphi }_2-{\varphi }_3-{\varphi }_4+2{\varphi }_5$$

$$n\left[{\varphi }_3+{\varphi }_4\pm {\left\{({\varphi }_3-{\varphi }_4{)}^2+4{\varphi }_5^2\right\}}^{1/2}\right]$$

$$4{\varphi }_5^2+({\varphi }_3-{\varphi }_4{)}^2$$

Full text

*Keywords: * array data; weighted network; dependent data; generalized least squares.

1 Introduction

Entries in relational arrays quantify pairwise interactions between actors that may be of multiple types or observed over time. Examples include annual flows of migrants between countries (Aleskerov et al.,, 2017) and interactions between students over the course of a semester (Han et al.,, 2016). In economics, relational arrays are used to describe monetary transfers between individuals as part of informal insurance markets  (see Bardham,, 1984; Fafchamps,, 2006; Foster and Rosenzweig,, 2001; Attanasio et al.,, 2012; Banerjee et al.,, 2013)). Other examples of data that are naturally represented as relational arrays include gene expressions (Zhang and Horvath,, 2005) and international relations (Fagiolo et al.,, 2008).

A relational array $Y=\left(y_{ijr}\right)$ where $(i,j=1,\ldots,n,i\not=j,r=1,\ldots,R)$, is composed of a series of $R$ matrices of size $(n\times n)$, each of which describes the directed pairwise relationships among $n$ actors of type $r$, e.g. time period $r$ or relation context $r$. The diagonal elements of each matrix, for example $y_{iir}$, are assumed to be undefined, as we do not consider, e.g., international relations of a country with itself. The relationship from actor $i$ to actor $j$ may differ than that from $j$ to $i$, $y_{ijr}\neq y_{jir}$, however, the methods we propose extend to the symmetric relation case (see Section S3).

The primary goal in our setting is inference for linear regressions exploring the effects of exogenous covariates on the values in the relational array, expressed as

[TABLE]

where $y_{ijr}$ is a continuous directed measure of the $r$th relation from actor $i$ to actor $j$, $x_{ijr}$ is a $(p\times 1)$ vector of covariates, and $\xi_{ijr}$ is an unobserved, scalar random error. For example, considering informal insurance markets, fafchamps2007formation examine how covariates such as geographical proximity and kinship relate to risk sharing relations after economic shocks.

A core challenge in making inference on $\beta$ arises from the innate dependencies among error relations involving the same actor. For example, dependence often exists between trade relations involving the same country or between economic transfers originating from the same individual. This dependence may arise due to variation unaccounted for in the covariates, for example, from differences in production levels between nations or from individual differences in risk aversion. Standard regression techniques may lead to poor estimates of $\beta$ and/or incorrect conclusions regarding the significance of the estimate of $\beta$. Approaches to account for error dependence in relational arrays have appeared in the statistics, biostatistics, and econometrics literatures and can be characterized into two broad classes.

The first set of approaches impose a parametric model on the errors. Specifically, they either use latent variables to model the array measurements as conditionally independent given the latent structure (see holland1983stochastic; wang1987stochastic; hoff2002latent; li2002unified; hoff2005bilinear) or model the error covariance structure directly subject to a set of simplifying assumptions  (see hoff2011separable; fosdick2014separable; hoff2015multilinear). While these methods allow for (possibly) improved estimation of $\beta$ and appropriate standard error estimators in the presence of relational dependence, the accuracy of inference on $\beta$ depends on the extent to which the true error structure is consistent with the specified parametric model. In addition, many of these models are estimated in a Bayesian paradigm using Markov chain Monte Carlo approaches, which are commonly computationally expensive to estimate.

The second set of approaches to accounting for relational dependence relies heavily on empirical estimates of the error structure based on the regression residuals, first proposed by fafchamps2007formation and based on the spatial dependence work of conley1999gmm. This framework is model agnostic, making as few assumptions as possible about the data generating process. One empirical approach estimates the regression coefficients using ordinary least squares, and then utilizes a sandwich covariance estimator – which is robust to a wide array of error structures – for the standard errors of the regression coefficients (fafchamps2007formation; aronow2015cluster). In finite samples, this estimator is hindered by the need to estimate a large number of covariance parameters with limited observations (king2014robust for a discussion in other contexts), and is the reason why wakefield2013bayesian suggests such estimators be labeled empirical rather than robust. We observe that standard errors from this empirical framework are often highly variable and are anticonservative.

In this work, we introduce an empirical estimation approach for relational arrays that incorporates an exchangeability assumption. This assumption is implicit in many of the model-based approaches discussed previously and is a hallmark of Bayesian hierarchical models within, and outside, the relational context (orbanz2015bayesian). Our key contribution is to define the covariance matrix (and a corresponding estimator) of the relational error array under exchangeability. Use of our parsimonious estimator produces superior estimates of $\beta$ and its standard errors relative to existing approaches, and our estimator is easier to compute than existing Bayesian model-based and exchangeable bootstrapping (menzel2017bootstrap; green2017bootstrapping) approaches. Reproduction code is available at https://github.com/fmarrs3/netreg_public and methods are implemented in R package netregR.

2 Inference in relational regression

2.1 Estimation of regression coefficients

We employ a least squares framework to perform inference on ${\beta}$ in the relational regression model in (1) aitkin1935least. An unbiased estimator for $\beta$ is the ordinary least squares estimator, $\widehat{\beta}=\left(X^{T}X\right)^{-1}X^{T}y$, where $X$ is an \big{(}Rn(n-1)\times p\big{)} matrix of $(p\times 1)$ covariate vectors $(x_{ijr})$ and ${y}$ is a vectorized representation of $(y_{ijr})$. The least squares estimator is the best linear unbiased estimator for $\beta$ when the covariance matrix $\Omega=\text{var}(y~{}\mid~{}X)$ is proportional to the identity matrix. Dependence is expected in relational data, e.g. between relations $(i,j,r)$ and $(i,k,r)$ which share actor $i$. If $\Omega$ were known, the best linear unbiased estimator for $\beta$ is the generalized least squares estimator of aitkin1935least,

[TABLE]

In practice, $\Omega$ is unknown and must be estimated. Given estimator $\widehat{\Omega}$, alternating estimation of ${\Omega}$ with (2) (replacing $\Omega$ with $\widehat{\Omega}$ at each iteration) is termed feasible generalized least squares. When $\widehat{\Omega}$ is consistent, feasible generalized least squares is asymptotically efficient for $\beta$ (greene2003econometric; hansen2015econometrics).

Regardless of whether the ordinary least squares estimator or (2) is used to estimate $\beta$, uncertainty estimates are required for inference. A common approach is to approximate the distribution of the $\beta$ estimate as a multivariate normal random variable and construct confidence intervals using an estimator of its variance: in the ordinary least squares setting,

[TABLE]

and in the generalized least squares setting,

[TABLE]

where $\widetilde{\Omega}$ is the final estimate of $\Omega$ from the generalized least squares procedure. Variance estimators are often constructed by substituting an estimator for $\Omega$ in (3) and (4), and are commonly termed sandwich estimators (huber1967behavior; white1980heteroskedasticity). Thus, inference for $\beta$ requires an estimator for $\Omega$, regardless of how $\beta$ is estimated, and properties of the estimator of $\text{var}(\widehat{\beta}\mid X)$ depend strongly on the estimator of $\Omega$.

2.2 Dyadic clustering estimator

fafchamps2007formation, cameron2011robust, aronow2015cluster, and tabord2017inference propose and describe the properties of a flexible standard error estimator for relational regression which makes the sole assumption that two relations $(i,j,r)$ and $(k,l,s)$ are independent if $(i,j)$ and $(k,l)$ do not share an actor. This assumption implies that $\text{cov}(y_{ijr},y_{kls}\mid X)=\text{cov}(\xi_{ijr},\xi_{kls}\mid X)=0$ for non-overlapping relation pairs, but places no restrictions on the covariance elements for pairs of relations that share an actor. Let $\Omega_{DC}$ denote the covariance matrix of $\xi$, subject to this non-overlapping pair independence assumption. fafchamps2007formation propose estimating each nonzero entry of $\Omega_{DC}$ with a product of residuals, i.e. using $e_{ijr}e_{iks}$ to estimate $\text{cov}(\xi_{ijr},\xi_{iks})$, where $e_{ijr}=y_{ijr}-\widehat{\beta}^{T}x_{ijr}$. The estimator $\widehat{\Omega}_{DC}$ can be seen as that which takes the empirical covariance of the residuals defined by ${e}{e}^{T}$, where $e$ is a vector of the set of residuals $(e_{ijr})$, and introduces zeros to enforce the non-overlapping pair independence assumption. fafchamps2007formation propose a sandwich variance estimator for $\text{var}(\widehat{\beta}\mid X)$ in (3) based on $\widehat{\Omega}_{DC}$,

[TABLE]

We refer to $\widehat{V}_{DC}$ as the dyadic clustering estimator as it owes its derivation to the extensive literature on cluster-robust standard error estimators.

The dyadic clustering estimator in (5) has the attractive properties that it is asymptotically consistent under wide range of error dependence structures and is fast to compute. However, $\widehat{\Omega}_{DC}$ estimates ${O}(R^{2}n^{3})$ nonzero covariance elements separately based on ${O}(R^{2}n^{2})$ dependent observations, and thus, $\widehat{V}_{DC}$ is inherently quite variable. Only when there is extreme heterogeneity in the true covariance structure is the dyadic clustering method optimal and it will suffer a loss of efficiency otherwise.

3 Standard errors under exchangeability

3.1 Exchangeability in relational models

A common modeling assumption for relational and array structured errors is exchangeability. Defined by de Finetti for a univariate sequence of random variables, exchangeability was generalized to array data by hoover1979relations and aldous1981representations. The errors in a relational data model are jointly exchangeable if the probability distribution of the error array $\Xi=(\xi_{ijr})$ is invariant under any simultaneous permutation of the rows and columns, and secondary permutation of the third dimension. Mathematically, this means $\text{pr}(\Xi)=\text{pr}\{\Pi(\Xi)\},$ where $\Pi(\Xi)=\{\xi_{\pi(i)\pi(j)\nu(r)}\}$ is the error array with its indices reordered according to permutation operators $\pi$ and $\nu$. Intuitively, exchangeability in the regression context means the observed covariates are sufficiently informative such that the labels of the rows and columns in the error array are uninformative. Similarly the ordering of the third dimension of the error array is uninformative to its distribution. This assumption may be appropriate when the third dimension of the array represents different contexts of observations (such as economic trade sectors) that have no inherent ordering, or when the third dimension represents time periods, but the bulk of the temporal variation is accounted for in the covariates. Many of the conditionally independent parametric latent variable models cited in the Introduction have this joint exchangeability property (hoff2008modeling; bickel2009nonparametric).

3.2 Impact of exchangeability on covariance structure

li2002unified and hoff2005bilinear describe several particular random effects models for $R=1$. The corresponding error covariance matrices have different entries depending on the model, yet all covariance matrices have at most six unique entries. A key contribution of this paper is to formalize and extend this observation, showing that any jointly exchangeable model for relational array $\Xi$ results in an $\Omega$ of the same form, with at most six unique terms when $R=1$ and at most twelve unique terms when $R>1$.

Proposition 1**.**

If a probability model for a directed relational array $\Xi$ is jointly exchangeable and has finite second moments, then the covariance matrix of $\Xi$ contains at most twelve unique values.

The twelve (possibly) unique entries in $\Omega$ correspond to the twelve distinguishable configurations of relation pairs $(i,j,r)$ and $(k,l,s)$ with unlabeled actors. The twelve configurations can be separated into two sets of six identical configurations of relations $(i,j)$ and $(k,l)$ with unlabeled actors, as depicted in Figure 1, where each set corresponds to $r=s$ and $r\neq s$. A proof is provided in Section S4.

3.3 Covariance matrices of exchangeable relational arrays

Similar to the dyadic clustering estimator, we assume non-overlapping relation pairs are independent, such that $\text{cov}(\xi_{kls},\xi_{ijr})=0$, for any $s$ and $r$ when $(i,j,k,l)$ are distinct. This assumption sets two of the twelve parameters in $\Omega$ to zero. We introduce a new class of covariance matrices which contain ten possibly nonzero entries, $\phi_{0}^{(\eta)},\phi_{a}^{(\eta)},\phi_{b}^{(\eta)},\phi_{c}^{(\eta)},\phi_{d}^{(\eta)},(\eta=1,2)$, associated with the left panel in Figure 1. The separation of covariances by $r=s$ and $r\neq s$ implies that $\Omega$ consists of blocks of $\Omega_{1}$ and $\Omega_{2}$, each consisting of five nonzero terms for $r=s$ ($\eta=1$) and $r\neq s$ ($\eta=2$), respectively (Figure 1, right panel). We define an exchangeable covariance matrix as any covariance matrix of this form and denote it $\Omega_{E}$.

4 Exchangeable estimator definition and evaluation

4.1 Exchangeable covariance estimator

Consider relational regression models with $\Omega$ of the exchangeable form $\Omega_{E}$. The proposed exchangeable estimator of $\text{var}(\widehat{\beta}\mid X)$ is then

[TABLE]

where $\mathcal{S}^{(\eta)}_{u}$ denotes the $\{Rn(n-1)\times Rn(n-1)\}$ binary matrix with $1$’s in the entries corresponding to relation pairs of type $(u=0,a,b,c,d;\eta=1,2)$ as defined in Figure 1. We propose estimating the ten parameters in $\Omega$ by averaging the residual products that share the same index configurations, corresponding to (0)-(d) in Figure 1. For example, the estimate of ${\text{cov}}(\xi_{kls},\xi_{ijr})$, corresponding to $u=b$ and $\eta=2$, is

[TABLE]

The remaining nine estimators for $(s=0,a,\ldots,e;\eta=1,2)$, are defined analogously, and $\widehat{\Omega}_{E}$ may be interpreted as the projection of $\widehat{\Omega}_{DC}$ into the vector space over symmetric matrices of the form of $\Omega_{E}$. We provide theoretical details for the proposed exchangeable estimator, including asymtotic normality of least squares under exchangeability and consistency of $\widehat{V}_{E}$ in Sections S7- S9.

4.2 Comparison of exchangeable estimator with dyadic clustering

It is intuitive that the moment-based exchangeable estimator is consistent, and more efficient than the dyadic clustering estimator whenever the exchangeability assumption is satisfied. One might expect the highly parametrized dyadic clustering estimator to trade off high variance for reduced bias. However, we derive the result that the dyadic clustering estimator is biased downwards, and that this bias is larger than twice the bias of the exchangeable estimator. One concludes that one tradeoff for the robustness of the dyadic clustering estimator is anticonservatism. The proof of Theorem 2 is provided in Section S6. We provide additional theoretical details for the proposed exchangeable estimator (comparison of mean-square error with dyadic clustering) in Section S10.

Theorem 2**.**

Consider error vector $\xi$ and normally distributed covariate vector $x$ with mean zero, exchangeable covariance matrix, and bounded fourth moments, where

[TABLE]

Then, the dyadic clustering estimator for $\text{var}(\widehat{\beta}_{2})$ is biased downwards,

[TABLE]

noting that both $n^{2}\text{Bias}(\widehat{V}_{DC})$ and $|n^{2}\text{Bias}(\widehat{V}_{E})|$ are $O(1)$.

We conducted a simulation study to compare the bias and 95% confidence interval coverage when using the exchangeable and dyadic clustering estimators. We simulated from a model with three covariates (one each of binary, positive real, and real-valued) with exchangeable and non-exchangeable error models, for $R=1$ and $n=20,40,80,160,320$ (see Section S11 for details). Figure 2 shows the estimated mean coverage and middle 95% of coverages across various $X$ realizations. In all settings, the estimated mean coverage of the exchangeable estimator is closer to the nominal $0.95$ level than the dyadic clustering estimator. This difference is most pronounced for the binary covariate, where there is reduced signal-to-noise relative to the other covariates. The average bias of the dyadic clustering estimator is typically more than four times that of the exchangeable estimator under the exchangeable error model, confirming Theorem 2 and driving the poorer coverage performance (see Section S11.2 for details).

5 Patterns in international trade

We demonstrate the implications of using our exchangeable estimator in a study of international trade among $58$ countries over $R=T=20$ years. These data were previously analyzed and made available by westveld2011mixed. Following westveld2011mixed, ward2007persistent, and tinbergen1962shaping, we use a modified gravity mean model to represent log yearly trade between each pair of countries as linear function of seven covariates in years 1981-2000. westveld2011mixed propose a model – which we refer to as the mixed effects model – which explicitly decomposes the regression error term $\epsilon_{ijt}$ for each time period and pair of actors into time-dependent sender and receiver effects, resulting in 13 error covariance parameters which are estimated using a latent variable representation and Bayesian Markov chain Monte Carlo methodology. We propose estimating the gravity mean model using feasible generalized least squares, assuming the errors are jointly exchangeable. As noted in Section 4.1, the proposed approach estimates ten error covariance parameters. We provide a method for efficient inversion of $\Omega_{E}$ in this setting in Section S13.

We compared the exchangeable and mixed effects approaches (and ordinary least squares as a baseline) in an out-of-sample prediction study. Here we estimated the regression coefficients using the first $K$ years of trade data for $K=4,...,19$ and used the estimates to predict trade values in the following year. Figure 3 provides the coefficient of determination, $R^{2}$, for the three procedures when predicting trade flows in years 5 through 20. There is a median increase in $R^{2}$ of about 10% (30%) when using the proposed exchangeable approach relative to the mixed effects approach (and ordinary least squares, respectively). The proposed approach performs better than the other approaches for all time periods, although the gap in performance decreases as $T$ increases. These results suggest the more parsimonious exchangeable approach represents the data better than the mixed effects model, and yet, the exchangeable approach runs in a small fraction of the time of the mixed effects approach. See Section S14 for additional details.

We compared the coefficients of the ordinary least squares, mixed effects, and exchangeable approaches. The $R^{2}$ between the ordinary least squares and mixed effects coefficients is about $0.46$, while the exchangeable and mixed effects coefficients have an $R^{2}$ of $0.78$. About $40\%$ of the ordinary least squares coefficients (using dyadic clustering standard errors) were significantly different from the mixed effects coefficients, while only $1\%$ of the exchangeable coefficients were significantly different from the the mixed effects coefficients. Finally, the ordinary least squares coefficients have standard errors that are, on average, over $1.5$ times the standard errors of the exchangeable coefficients (with exchangeable standard errors). See Section S14.3 for more details. Together with the one-year-ahead prediction results, the coefficient and standard error comparisons suggest that the proposed approach can revise ordinary least squares coefficients in the direction of a higher fidelity model, giving more precise estimates of the coefficients, while requiring few modeling decisions and with limited runtime penalty. Based on the success of the exchangeable approach in the trade data analysis and simulation study, we recommend that researchers use the feasible generalized least squares estimator of the coefficient vector $\beta$ as demonstrated here (unless an unbiased estimator of the coefficient vector is specifically desired).

6 Discussion

The proposed exchangeable estimator leverages exchangeability for maximal symmetry (and thus maximal parsimony) in the covariance matrix of relational array $Y$. The exchangeability assumption may not be appropriate when the true error covariances are substantially heterogeneous. We propose using a permutation test based on the dyadic clustering estimator for testing the hypothesis of exchangeable errors. The procedure consists of generating a null distribution of $\widehat{V}_{DC}$ in (5) by randomly permuting the residual array in a manner consistent with exchangeability. If the observed estimator is extreme relative to the null distribution, this suggests the errors are non-exchangeable. Details and simulations are available in Section S15.

Supplementary material

The supplementary material contains details about the estimator and analyses, proofs, further theoretical details, and a proposed test for exchangeability.

S1 Patterns of exchangeable covariance matrices

Figure S1 shows the structure of $\Omega_{E}$ for a relational matrix with four actors $\{A,B,C,D\}$ and $R=1$.

S2 Exchangeable estimator details

The text gives an example for the estimator of $\phi_{b}^{(2)}$. We explicitly define the empirical mean estimates used in the exchangeable estimator in (6) below:

[TABLE]

Recall from the text that $\widehat{\phi}_{0}^{(1)}$ is an estimator of ${\text{var}}(\xi_{ijr})$, $\widehat{\phi}_{a}$ is an estimator of ${\text{cov}}(\xi_{ijr},\xi_{jir})$, $\widehat{\phi}_{b}^{(1)}$ is an estimator of ${\text{cov}}(\xi_{ijr},\xi_{kjr})$, $\widehat{\phi}_{c}^{(1)}$ is an estimator of ${\text{cov}}(\xi_{ijr},\xi_{ikr})$, and $\widehat{\phi}_{d}^{(1)}$ is an estimator of ${\text{cov}}(\xi_{ijr},\xi_{jkr})$ (and equivalently an estimator of ${\text{cov}}(\xi_{ijr},\xi_{kir})$). The estimator $\widehat{\phi}_{i}^{(2)}$, for $i=0,a,b,c,d$, is the analogous estimator to $\widehat{\phi}_{i}^{(1)}$, only when $r\neq s$.

S3 Undirected arrays

This section specializes the results presented in the manuscript to undirected relational data. Consider the case when $R=1$ and suppose the relational data contains the relations among $n$ actors. The covariance of the errors $\Omega$ contains six possibly unique elements

[TABLE]

As in the directed case, we assume covariances corresponding to relations which share no common actor are zero. We again propose estimating the remaining nonzero covariances using the corresponding means of residuals.

S4 Proof of Proposition 1

Proof.

Consider a probability model for a directed relational array $\Xi$ that satisfies the joint exchangeability and second moment criteria defined in the proposition. Now consider arbitrary permutations $\pi(\cdot)$ of the actor set $(1,...,n)$ and $\nu(\cdot)$ of the actor set $(1,...,R)$. Exchangeability implies the that the probability distribution of $\xi_{ijr}$ is the same as $\xi_{\pi(i)\pi(j)\nu(r)}$. Thus, every entry in $\Xi$ is marginally identically distributed, and

[TABLE]

This is one of the entries in $\Omega$, the covariance matrix of $\Xi$.

We now move onto the covariances in $\Omega$; it remains to show that there are 11 unique entries beyond the variance. This result follows from two facts. The first fact is that exchangeability implies that the bivariate distribution of the pair $(\xi_{ijr},\xi_{klr})$ must be the same as bivariate distribution of $\{\xi_{\pi(i)\pi(j)r},\xi_{\pi(k)\pi(l)r}\}$, and thus,

[TABLE]

The second fact is that permutations preserve equality and inequality, meaning that if $i=j$, then $\pi(i)=\pi(j)$, and if $i\neq j$, then $\pi(i)\neq\pi(j)$, for any permutation $\pi$. We will enumerate the 11 patterns of indices that have the same bivariate distributions, as implied by the equality-inequality preservation of permutations, being mindful of the fact that the first two indices in $\xi_{ijr}$ refer to the same index set $(1,\ldots,n)$.

First take the case when $i=k$ and $j=l$. When $r=s$, we have the variance previously discussed. When $r\neq s$, we obtain the second set of equivalent bivariate distributions (after the variance):

[TABLE]

where it is understood that $i\neq j$ due to the network setting. Moving on, consider the case when $i=k$ and $j\neq l$, noting that, by assumption, $i\neq j$ and $k\neq l$. Then, we obtain two sets with equivalent bivariate distributions, pertaining to $r=s$ and $r\neq s$:

[TABLE]

Thus far, we have enumerated four sets of equivalent bivariate distributions under exchangeability. The remaining eight may be determined by cycling through combinations of indices in $(i,j)$ equal to, and not equal to, indices in $(k,l)$, for both $r=s$ and $r\neq s$. Table S1 summarizes the twelve index patterns that result; the first and third row in the table have already been discussed. The second row, first column, states that pairs of entries in $\Xi$ with the index pattern $(\xi_{ijr},\xi_{jir})$ form a set with the same bivariate distribution.

Finally, one might at first think that we have missed the bivariate distribution of $(\xi_{ijr},\xi_{jlr})$. However, we note that the bivariate distribution of random variables $(A,B)$ is equal to that of $(B,A)$, and thus, the bivariate distribution of $(\xi_{ijr},\xi_{jlr})$ is captured in the first column of row $(d)$ of Table S1. Since we have enumerated all entries in $\Omega$, and there are twelve possibly unique sets with equivalent covariances (or variances), the result is shown for $R>1$. When $R=1$, it is sufficient to consider only the first column of Table S1, which establishes that there are at most six unique entries in $\Omega$ for $R=1$. ∎

S5 Eigenvalues of exchangeable covariance matrix

S5.1 Positive definite $\widehat{\Omega}_{E}$

Since the entries in the exchangeable covariance matrix estimator $\widehat{\Omega}_{E}$ are empirical averages, it is possible the estimate $\widehat{\Omega}_{E}$ is not positive definite. Here we briefly investigate the constraints on the parameters that guarantee the resulting covariance matrix is positive definite when $R=1$. Note that for computing the sandwich estimator variance of $\widehat{\beta}$ and making inference on $\widehat{\beta}$, positive definiteness of $\widehat{\Omega}_{E}$ is not necessary. However, if generalized least squares is employed to estimate $\beta$, the inverse of the covariance matrix estimator is required, and hence it is desirable that $\widehat{\Omega}_{E}$ be positive definite.

S5.2 Undirected relational data

We first focus on the undirected case, where the exchangeable covariance matrix contains two distinct nonzero entries: a variance $\sigma^{2}$ and a parameter $\phi$ in the off-diagonal representing the correlation between any pairs of relations that share an actor. Below we consider the correlation matrix, rather than the covariance matrix, which contains only a single nonzero correlation value. We denote this value by $a=\theta_{a}^{(1)}/\theta_{0}^{(1)}$.

Based on a thorough empirical investigation, we conjecture that the exchangeable correlation matrix corresponding to an undirected set of relations among $n$ actors, which has nonzero value $a$ in select off-diagonal entries, has exactly three eigenvalues as given below.

[TABLE]

The correlation matrix is positive definite if and only if all eigenvalues are positive. Thus, if $a\in\left(1/({2(n-2)}),1/2\right)$, the correlation matrix is positive definite. Notice that the upper bound on $a$ does not vary with $n$. Using the relation between $a$ and $(\theta_{0}^{(1)},\theta_{a}^{(1)})$, this constraint can be re-expressed as a constraint on the covariance parameters.

S5.3 Directed relational data

Since matrices of the form of $\Omega_{E}$ form a five-dimensional linear subspace over square symmetric matrices, it is intuitive that the eigenspace of such matrices is limited. We find empirically that square symmetric matrices of the form of $\Omega_{E}$ have five unique eigenvalues (when $R=1$). As in the undirected case, we again focus on the exchangeable correlation matrix which contains four nonzero off-diagonal elements $\{a,b,c,d\}$ corresponding, respectively, in placement to $\{\phi_{a}^{(1)},\phi_{b}^{(1)},\phi_{c}^{(1)},\phi_{d}^{(1)}\}$ in the exchange covariance matrix $\Omega_{E}$. Note $a=\phi_{a}^{(1)}/\phi_{0}^{(1)}$, $b=\phi_{b}^{(1)}/\phi_{0}^{(1)}$, and so on. Based on empirical studies, we conjecture the eigenvalues for the exchangeable correlation matrix associated with a directed set of relations among $n$ actors has exactly five eigenvalues as given below.

[TABLE]

where $\alpha=(c^{2}+b^{2})(n^{2}-2n+1)+4d^{2}(n^{2}-6n+9)+2bc(1-n^{2}+2n)$ and $\beta=ad(8n-24)+(b+c)d(12-4n)+4a\{a-(b+c)\}$. As in the undirected case, these constraints can be re-expressed as constraints on the original five covariance parameters.

S6 Proof of Theorem 2: bias of the exchangeable estimator

S6.1 Notation and proof preliminaries

In this section we prove that the dyadic clustering estimator is biased downwards, and that this bias is greater in absolute value than the bias of the exchangeable estimator, in the case when $R=1$.

First, we make some definitions. Recall that the dyadic clustering estimator is based on the outer product of residuals, with entries corresponding to pair of relations that do not share an actor set to zero. This operation may be viewed as a projection onto a vector space of square, symmetric matrices of dimension $n(n-1)$ with zeros in the appropriate entries. The dimension of this vector space is $O(n^{3})$. We define the vector space below:

Definition S3**.**

Vector space $\mathcal{D}$: Let $M_{jk,lm}$ be the square matrix of dimension $n(n-1)$ with a ‘1’ in the entry corresponding to relations $jk$ and $lm$, and zeros elsewhere. Then,

[TABLE]

We define the orthogonal projection of matrix $A$ onto $\mathcal{D}$ as $P_{\mathcal{D}}\left(A\right)$.

A similar definition to Definition S3 may be made concerning the space of covariance matrices corresponding to jointly exchangeable random variables, as defined below.

Definition S4**.**

Vector space $\mathcal{E}$: Let $\mathcal{S}_{i}$ be an indicator matrix of entries in square, symmetric matrices of dimension $n(n-1)$ corresponding to pairs of relations that share an actor in the $i\text{th}$ manner, for $i=1,\ldots,5$. Then,

[TABLE]

We define the orthogonal projection of matrix $A$ onto $\mathcal{E}$ as $P_{\mathcal{E}}\left(A\right)$.

It will be useful to write the dyadic clustering and exchangeable estimators as functions of the projections in Definitions S3 and S4, supposing that $X$ is orthogonal, that is, $X^{T}X=f(n)I$, where $f(n)$ is a function of the number of actors $n$. The dyadic clustering estimator of $\text{var}(\widehat{\beta}_{k})$ is then

[TABLE]

where $X_{k}$ is column $k$ of covariate matrix $X$. The second equality in (S1) results from the application of a property of inner products, which we apply to the inner product over matrices, where, for any projection operator $P(.)$, we have that $\text{tr}\left\{AP(B)\right\}=\text{tr}\left\{P(A)B\right\}=\text{tr}\left\{P(A)P(B)\right\}$. An analogous definition to (S1) holds for the exchangeable estimator of $\text{var}(\widehat{\beta}_{k})$ when replacing the projections onto $\mathcal{D}$ in (S1) with projections onto $\mathcal{E}$.

Finally, to prove Theorem 2, we need some relationships between the error covariance parameters $(\phi_{3}^{(1)},\phi_{4}^{(1)},\phi_{5}^{(1)})$ when $\Omega$ is positive definite. We summarize these relationships in the inequalities in the following proposition. We drop the superscript ‘$(1)$’ to lighten notation.

Proposition S5**.**

Let $\Omega\in\mathcal{E}$ be an exchangeable correlation matrix, with $\Omega={I}+\sum_{i=2}^{5}\phi_{i}\mathcal{S}_{i}$. Then,

[TABLE]

Proof.

The inequalities result from the requirement that eigenvalues of $\Omega$ are positive for all $n$. Taking $n$ arbitrarily large, the eigenvalues of $\Omega$ are

[TABLE]

The second and third eigenvalues give that $\phi_{3}+\phi_{4}\leq 1$. The fourth and fifth eigenvalues imply $\phi_{3}+\phi_{4}\geq 0$.

Using the fourth and fifth eigenvalue and rearranging,

[TABLE]

Then, we have that

[TABLE]

which establishes the second result. The third result follows from an analogous argument to (S3) for $\phi_{4}^{2}+\phi_{5}^{2}$. The final result follows from summing the previous two, such that

[TABLE]

∎

S6.2 Proof of Theorem 2

In this proof, we start by establishing the expression for the true variance of $\widehat{\beta}_{2}$. Then, we write the bias of the dyadic clustering and exchangeable estimators as a function of the variance parameters of $x$ and $\xi$. Then, we write the bias of the dyadic clustering estimator as a function of the bias of the exchangeable estimator, and then bound the absolute value of the bias of the exchangeable estimator, which gives the result. Throughout, we assume that the bias of the exchangeable and dyadic clustering estimators, scaled by $n^{2}$, are $O(1)$. We verify this assumption at the end of the proof. Without loss of generality, we prove for the scaled problem where $E(xx^{T})=\Psi=\sum_{i=1}^{5}\psi_{i}\mathcal{S}_{i}$ with $\psi_{1}=1$.

We begin by expressing the true variance of $\widehat{\beta}_{2}$. We require that $n^{1/2}(\widehat{\beta}_{2}-\beta_{2})$ converges in distribution to a normal random variable with finite variance, which may be expected from Proposition 3.2 in tabord2017inference. By assumption, $\widehat{\beta}_{2}=(x^{T}x)^{-1}x^{T}y$. The true variance of $\widehat{\beta}_{2}$ is then

[TABLE]

where we use that $x^{T}x=n^{2}+O(n^{-1})$ by assumption. Using (S1), it will be useful to rewrite the true variance

[TABLE]

where we use the fact that $\Omega$ is a member of both spaces $\mathcal{D}$ and $\mathcal{E}$ by assumption.

We now write the expectation of the dyadic clustering estimator, taking only the leading terms. First, the expectation

[TABLE]

where $P_{X}$ is the projection onto the column space of the covariate matrix $X$. By assumption, we have that

[TABLE]

where ${1}$ is a vector of 1’s of appropriate length (in this case, length $n(n-1)$). Using the expression of the dyadic clustering estimator in (S1), its expectation is

[TABLE]

The second term in (S4) is zero in probability, as the row sums of $\Omega$ are all the same, and are $O(n)$. Thus,

[TABLE]

Then, the second term is

[TABLE]

since $n^{-2}1^{T}x=O(n-1)$ by assumption, which implies that the second term vanishes relative to the rest.

The third term in (S4) may be written

[TABLE]

where $c_{n}=O(n)$ is the row sum of $\Omega$, which again is the same for all rows. Further, using the Cauchy-Scharz inequality,

[TABLE]

Substituting, the third term is

[TABLE]

We now bound the order of $||P_{\mathcal{D}}\left(xx^{T}\right)||_{F}$. By assumption, $x^{T}x=n(n-1)+O(n^{-1})$, and thus

[TABLE]

Using this bound, the third term is

[TABLE]

and the third term vanishes relative to the other terms in (S4). Further, the sixth term in (S4) is of the same order,

[TABLE]

and may also be neglected.

Using (S5), (S6), and (S7), the expression for the bias of the dyadic clustering estimator based on (S4) becomes

[TABLE]

Since, by the properties of projections, $||P_{\mathcal{E}}\left(xx^{T}\right)||_{F}\leq||P_{\mathcal{D}}\left(xx^{T}\right)||_{F}$ since $\mathcal{E}\subset\mathcal{D}$, we may write the bias of the exchangeable estimator by analogy:

[TABLE]

We now explicitly determine the value of the exchangeable bias in (S9). By assumption,

[TABLE]

recalling that $E(xx^{T})=\Psi=\sum_{i=1}^{5}\psi_{i}\mathcal{S}_{i}$. Then, the projection of $xx^{T}$ onto the exchangeable space is

[TABLE]

and the scaled true variance is

[TABLE]

Using (S10) and (S11), the first term in the exchangeable bias in (S9) is

[TABLE]

The second term in the exchangeable bias in (S9) is

[TABLE]

Now we must determine $n^{-4}x^{T}\mathcal{S}_{i}\mathcal{S}_{j}x$ for each $i$ and $j$. Recalling that $\mathcal{S}_{i}$ is a matrix that indicates pairs of relations $jk$ and $lm$ that share an actor in the $i^{th}$ manner,

[TABLE]

Then, the sums (S12) may be determined by examining the relationships between the actors in relations $tu$ and $ab$. First, for any $i$ and $j$, $tu$ must share at least one actor with $ab$, as the expectation $E(x_{tu}x_{ab})$ is zero otherwise. This result implies that every sum in (S12) is at most $O(n^{4})$. As an example, we consider $i=3$ and $j=3$. In this case, $tu$ is of the form $ru$ and $ab$ is of the form $rb$. Whenever $b\neq u$, the expectation $E(x_{ru}x_{rb})=\psi_{3}$. When $b=u$, the expectation $E(x_{ru}x_{rb})=\psi_{1}$, although there are an order of magnitude fewer of these terms. Finally, similar to the expectation argument, the variance of the sums in (S12) consists of at most $O(n^{6})$ finite terms, and thus tends to zero. Then, the sum tends to the limit of its expectation. Summarizing, we have that

[TABLE]

Similar counting operations give the following result, for $i\leq j$,

[TABLE]

Substituting, we have the following expression for the bias of the exchangeable estimator

[TABLE]

We now return to the bias of the dyadic clustering estimator in (S8). The first term in the bias of the dyadic clustering estimator is

[TABLE]

As with the exchangeable estimator, the sum in (S14) converges to the limit of its expectation, as the variance tends to zero. Then, using Iserlis’ theorem,

[TABLE]

using the assumption that $\psi_{1}=1$. Then, the first term in the bias of the dyadic clustering estimator in (S8) is

[TABLE]

We now turn to the second term in the bias of the dyadic clustering estimator in (S8).

[TABLE]

where ‘$\circ$’ represents element-wise multiplication. As with the exchangeable estimator, the sum in (S15) has variance that tends to zero, and converges to its expectation. Writing in sum form and using Isserlis’ theorem,

[TABLE]

Substituting into (S15), the second term in the bias of the dyadic clustering estimator is

[TABLE]

Then, combining (S14) and (S16), we obtain the following expression for the bias of the dyadic clustering estimator,

[TABLE]

wherein we recognize the bias of the exchangeable estimator, and

[TABLE]

Since we have $\Omega$ positive definite and

[TABLE]

Thus, we have shown that the bias of the dyadic clustering estimator is less than twice the bias of the exchangeable estimator. To establish that this same property holds for the absolute values of the biases, it remains to show that $n^{2}\text{Bias}\{(\widehat{V}_{E})_{22}\}\leq x^{T}\Omega x+O(n^{-1/2})$, that is, that the right hand side of (S17) is negative and bounded away from zero. Combining terms in (S13), the bias of the exchangeable estimator is

[TABLE]

Now, using Proposition S5, we have that

[TABLE]

and thus,

[TABLE]

To establish a lower bound on the bias of the exchangeable estimator, using (S2) and (S3) in Proposition S5,

[TABLE]

and a similar argument holds for $-\psi_{4}(1+\psi_{5}^{2}/\psi_{4}^{2})$. Then, recalling from Proposition S5 that

[TABLE]

and combining with (S20), we obtain the following bound on the key expressions in (S18),

[TABLE]

To obtain the final inequality in (S21), we use the fact that $0\leq\phi_{3}+\phi_{4}\leq 1$ to find that

[TABLE]

and that

[TABLE]

A similar argument to that in (S21) can be obtained starting with $\psi_{3}^{2}+\psi_{4}^{2}+2\psi_{5}^{2}-2\psi_{4}(1+\psi_{5}^{2}/\psi_{4}^{2})$. Then, combining (S21) with the complete expression for the bias in (S18) and the upper bound in (S19),

[TABLE]

Substituting the bound on the absolute value of the exchangeable bias (S22) into (S17), we have that

[TABLE]

which establishes the result. To complete the proof, we note that by the explicit expression for $n^{2}\text{Bias}\{(\widehat{V}_{E})_{22}\}$ in (S13), the exchangeable bias scaled by $n^{2}$ is $O(1)$. Similarly, $n^{2}\text{Bias}\{(\widehat{V}_{DC})_{22}\}$ is $O(1)$ by (S17). ∎

S6.3 A corollary to Theorem 2

In the following corollary we expand the conditions of the bias theorem slightly, such that the error array is not strictly in the exchangeable, similar to Corollary S12.

Corollary S6**.**

Theorem 2 holds under the relaxation exchangeability of $\xi$, where the first two entries in the covariance matrix of $\xi$ are allowed to be heterogenous, that is, $E(\xi_{jk}^{2})=\sigma^{2}_{jk}$ and $E(\xi_{jk}\xi_{kj})=\theta_{jk}$, where $|\theta_{jk}|<\sigma^{2}_{jk}<C<\infty$. Similarly, the conditions hold when $E(x_{jk}x_{jk}^{T})=A_{jk}$ and $E(x_{jk}x_{kj}^{T})=D_{jk}$, where the entries in matrices $A_{jk}$ and $D_{jk}$ are bounded.

The result follows directly from noting that the proof of Theorem 2 depends only on the values of $E(x_{jk}x_{lm})$ and $E(\xi_{jk}\xi_{lm})$ for $(jk,lm)\in\Theta_{i}$ for $i\in\{3,4,5\}$, at least up to a vanishing constant. Thus, so long as the expectations $E(x_{jk}x_{lm})$ and $E(\xi_{jk}\xi_{lm})$ for $(jk,lm)\in\Theta_{i}$ for $i\in\{1,2\}$ are bounded, Theorem 2 holds. Corollary S6 states that even under heteroskedasticity in the covariance matrix of the error vector $\xi$, the exchangeable estimator still outperforms the dyadic clustering estimator in bias. Only when the entries in $\Omega$ corresponding to configurations (b), (c), and (d) are heterogeneous does the theory fail to hold.

S7 Theoretical properties of the exchangeable estimator

In this section, we provide additional theoretical results for the proposed exchangeable estimator. Namely, we show that the proposed exchangeable estimator is consistent and provides an improvement in mean-square error over the dyadic clustering estimator, when the assumed model is correct. First, we provide conditions under which the theory holds, and prove asymptotic normality of ordinary least squares (which is necessary for the results that follow). For brevity, we focus on the ordinary least squares coefficient estimator and the case of a singe time period or array observation, where $R=1$.

S7.1 Conditions for theory

We define the conditions under which we establish the theory, starting with a formal definition of the class of exchangeable covariance matrices.

Definition S7**.**

An exchangeable covariance matrix is defined as $\Omega_{E}=E[\xi\xi^{T}]$ arising from mean-zero random vector $\xi={\rm vec}(\Xi)$, where $\Xi$ is a jointly exchangeable random matrix with entries $\xi_{ij}$ independent $\xi_{kl}$ whenever $\{i,j\}\cap\{k,l\}\not=\varnothing$.

For the theoretical assessment of the exchangeable estimator, we take $X$ random, but still evaluate $\widehat{\beta}$, $\widehat{V}_{DC}$, and $\widehat{V}_{E}$ conditional on $X$. We assume the rows of the matrix $X$ are jointly exchangeable, meaning that a reordering of the rows $\{x_{jk}\}_{j,k=1}^{n}$ to $\{x_{\pi(j)\pi(k)}\}_{j,k=1}^{n}$ leaves the distribution of matrix $X$ invariant for any permutation $\pi(.)$. As with the dependence in the errors, we assume that two rows of $X$ that correspond to relations which do not share an actor are independent, that is row $x_{ij}^{T}$ is independent row $x_{kl}^{T}$ whenever $\{i,j\}\cap\{k,l\}=\varnothing$. This dependence in the rows of matrix $X$ (along with some assumptions on the finiteness of its moments) implies the following:

[TABLE]

where $\Theta_{i}$ is the set of pairs of relations $(jk,lm)$ that share a member in the $i$th manner and ‘[math]’ refers to self-relation (i.e. variance).

The theoretical setting is as follows:

Condition S8**.**

Define the following data generating process:

1. (A1)

The true data generating model is ${y}=X\beta+{\xi}$, where the errors ${\xi}$ are mean-zero with exchangeable covariance matrix as in Definition S7. 2. (A2)

At least one of $\{\phi_{b},\phi_{c},\phi_{d}\}$ is nonzero.

In addition, consider the following regularity conditions:

1. (B1)

The covariate matrix $X$ has rows that are jointly exchangeable with at least one of $\{M_{i}\}_{i\in\{b,c,d\}}$ in (S24) nonzero, and where row $x_{ij}^{T}$ is independent row $x_{kl}^{T}$ whenever $\{i,j\}\cap\{k,l\}\not=\varnothing$. 2. (B2)

The fourth moments of the errors and the eight moments of the covariates are bounded: $E(|\xi_{jk}|^{4})<C<\infty$ and $\max_{l\in\{1,2,\ldots,p\}}E(|x_{jk}^{(l)}|^{8})<C^{\prime}<\infty$ where $x_{jk}=[x_{jk}^{(1)},\ldots,x_{jk}^{(p)}]^{T}$. 3. (B3)

The errors ${\Xi}$ and covariates $X$ are independent. 4. (B4)

$X$* is full rank.*

S7.2 Additional theoretical results

The following theorem establishes asymptotic normality of $\widehat{\beta}$, and supports the normal approximation to $\widehat{\beta}$ to produce confidence intervals.

Theorem S9**.**

Under Conditions S8,

[TABLE]

where ‘$\rightarrow_{d}$’ denotes convergence in distribution.

Theorem S9 is similar to Proposition 3.2 in tabord2017inference, although with unbounded errors. It may be seen as an extension of central limit theorems for sums over infinitely exchangeable arrays (fortini2012central) to regression scenarios, and belongs to a broader literature on distributions of random processes that satisfy a symmetry property (kallenberg2006probabilistic; austern2018limit).

We establish consistency of the proposed estimator, and its improvement in performance over the dyadic clustering estimator, in the following three theorems.

Theorem S10**.**

Under Conditions S8, the exchangeable covariance estimator is consistent in the sense that

[TABLE]

where ‘$\rightarrow_{p}$’ denotes convergence in probability.

Theorem S11**.**

Under Conditions S8, the mean-square error of the exchangeable estimator is less than that of the dyadic clustering estimator with probability approaching 1, that is

[TABLE]

S7.3 Corollary to theory

The proofs of Theorem S9 through Theorem S11 depend only on the uniformity of the covariance entries in $\Omega$ corresponding to $\phi_{b}$, $\phi_{c}$, and $\phi_{d}$. Then, we can relax the assumption of exchangeability of the error vector $\xi$ to allow some heterogeneity, namely heteroskedasticity, under which the theory still holds.

Corollary S12**.**

Theorem S9 through Theorem S11 hold under the relaxation of (A2) and (B1), where the first two entries in the covariance matrix of $\xi$ are allowed to be heterogenous, that is, $E(\xi_{jk}^{2})=\sigma^{2}_{jk}$ and $E(\xi_{jk}\xi_{kj})=\theta_{jk}$, where $|\theta_{jk}|<\sigma^{2}_{jk}<C<\infty$. Similarly, the conditions hold when $E(x_{jk}x_{jk}^{T})=A_{jk}$ and $E(x_{jk}x_{kj}^{T})=D_{jk}$, where the entries in matrices $A_{jk}$ and $D_{jk}$ are bounded.

S8 Proof of asymptotic normality of ordinary least squares

S8.1 Notation and outline of proof

For this proof, and throughout the remainder of the supplementary material, we adopt slightly different notation to simplify the representation of the exchangeable covariance estimator. Recall that the exchangeable covariance estimator for the OLS estimating equations is

[TABLE]

where $\widehat{\Omega}_{E}$ is the exchangeable estimate of the error covariance matrix, consisting of five averages of residual products. Here we express $\widehat{\Omega}_{E}$ as

[TABLE]

This amounts to mapping $\sigma^{2}\mapsto\phi_{1}$, $\phi_{a}\mapsto\phi_{2}$, …, $\phi_{d}\mapsto\phi_{5}$, and re-indexing the $\mathcal{S}$ and $M$ matrices accordingly. Additionally, when we consider sequences of jointly exchangeable random variables $\{W_{ij}\}_{i,j=1}^{n}$, it is understood that the sequence arises from a relational array such that entries with $i=j$ are undefined. Thus, sums over the sequence are of $n(n-1)$ terms and we define $\sum_{ij}W_{ij}=\sum_{i\neq j}W_{ij}$.

We work in the asymptotic regime where actors are added incrementally to the relational data set, i.e. $n$ is continually increasing. To establish asymptotic normality of $\widehat{\beta}$, we wish to show

[TABLE]

where $\{M_{i}\}_{i\in\{1,3,4,5\}}$ are as in (S24) and ‘$\rightarrow_{d}$’ denotes element-wise convergence in distribution.

The motivation for the proof argument follows from the expression

[TABLE]

We state that $\left(\frac{\sum_{jk}x_{jk}x_{jk}}{n(n-1)}\right)^{-1}$ converges in probability to $M_{1}^{-1}$, and then show asymptotic normality of the second multiplicative term in (S26).

In analyzing $\{x_{ij}\xi_{ij}\}_{i,j=1}^{n}$, by condition (B1), the joint exchangeability and independence of non-overlapping pairs of the sequence $\{\xi_{ij}\}_{i,j=1}^{n}$ extends to the component sequences in the vectors $\{x_{ij}\xi_{ij}\}_{i,j=1}^{n}$. Thus, to prove asymptotic normality of $\widehat{\beta}$, we first prove a theorem stating that the average of a mean-zero sequence of jointly exchangeable random variables is asymptotically normal. Specifically, for $\{W_{ij}\}_{i,j=1}^{n}$ mean zero and jointly exchangeable, we show

[TABLE]

for some normalizing constant $\sigma$ and fixed sequence $k_{n}\rightarrow 0$ as $n\rightarrow\infty$.

To prove (S27), we rely on a result from bolthausen1982central, as well as a supporting lemma which we present here. Below we outline the significance of these results in the proof.

• •

Lemma S13 (bolthausen1982central): Provides a sufficient condition for asymptotic normality of a sequence of measures based on the standard normal characteristic function.

• •

Lemma S14: Provides a bound for a variance that surfaces in the proof of asymptotic normality in (S27).

From (S27), we immediately have the marginal asymptotic normality of the sample mean of each of the vector components in the sequence $\{x_{ij}\xi_{ij}\}_{i,j=1}^{n}$. To establish joint asymptotic normality, we employ the Cramér-Wold device (cramer1936some), where asymptotic normality of $\{{v}^{T}x_{ij}\xi_{ij}\}_{i,j=1}^{n}$, for all ${v}\in\mathbb{R}^{p}$ with $||{v}||=1$, establishes joint normality. To achieve the asymptotic normality of this inner product, we simply recognize that this inner product is itself the mean of an exchangeable sequence of random variables. Joint asymptotic normality of the mean of the sequence of vectors $\{x_{ij}\xi_{ij}\}_{i,j=1}^{n}$ establishes joint asymptotic normality of $\widehat{\beta}$ via (S26).

S8.2 Lemmas and theorem in support of Theorem S9

The following is Lemma 2 in bolthausen1982central and provides a sufficient condition for asymptotic normality. We abuse notation slightly, letting $i$ be the imaginary unit where appropriate.

Lemma S13** (bolthausen1982central).**

Let $\nu_{n}$ be a sequence of probability distributions over $\mathbb{R}$ which satisfies

$\ \sup_{n}\int x^{2}d\nu_{n}(x)<\infty$, and 2. 2.

for all $\lambda\in\mathbb{R}$, $\ \lim_{n}\int(i\lambda-x){\rm e}^{i\lambda x}d\nu_{n}(x)=0.$

Then, $\nu_{n}\rightarrow_{d}{\rm N}(0,1).$

To provide intuition for Lemma S13, the integral in condition (2) is identically zero when $\nu_{n}$ is the standard normal distribution.

The next lemma provides a sufficient condition on the dependence structure in $\{W_{ij}\}_{i,j=1}^{n}$ necessary for the proof of asymptotic normality in (S27). Reecall that terms in $\{W_{ij}\}_{i,j=1}^{n}$ with $i=j$ are undefined.

Lemma S14**.**

Let $\{W_{ij}\}_{i,j=1}^{n}$ be a sequence of jointly exchangeable random variables as in Definition S7 with $||W_{ij}||_{4}<L<\infty$, where $||W_{ij}||_{p}:=E\left(|W_{ij}|^{p}\right)^{1/p}$ for $p>0$. Then,

[TABLE]

for some $C<\infty$, where $\Theta_{ij}$ is the set of ordered pairs $(k,l)$ that share an index with $(i,j)$.

Proof.

By definition we write

[TABLE]

Each covariance of (S28) is bounded by $L^{4}$. To bound the variance, we will show the number of nonzero entries in the sum is ${O}(n^{5})$. For ${\rm cov}(W_{ij}W_{kl},W_{rs}W_{tu})\neq 0$, there must be overlap between the index sets $\{i,j,k,l\}$ and $\{r,s,t,u\}$. Further, the sum in (S28) is taken over index sets that themselves contain overlap, i.e. $\{i,j\}\cap\{k,l\}\neq\varnothing$ and $\{r,s\}\cap\{t,u\}\neq\varnothing$. For example, the index sets $\{i,j,i,l\}$ and $\{i,s,i,u\}$ have nonzero covariance in (S28). Since there are 5 unique indices in the union of the sets $\{i,j,i,l\}$ and $\{i,s,i,u\}$, there are ${O}(n^{5})$ such index set pairs of this form in total. There are 96 pairs of index sets that result in nonzero covariance ${\rm cov}(W_{ij}W_{kl},W_{rs}W_{tu})$. For example, another such pair of index sets is $\{i,j,i,l\}$ and $\{i,j,i,j\}$. Each of these 96 pairs of index sets is ${O}(n^{5})$. Thus, the sum of covariances in (S28) is over ${O}(n^{5})$ bounded elements. ∎

It is worth noting that we repeat the counting argument in the proof of Lemma S14 in many of the following proofs, including those in later sections. Now that we have Lemma S13 and S14, we prove that a general sequence of mean-zero exchangeable random variables is asymptotically normal.

Theorem S15**.**

Let $\{W_{ij}\}_{i,j=1}^{n}$ be a mean-zero sequence of jointly exchangeable random variables with at least one of $\{\phi_{3},\phi_{4},\phi_{5}\}$ nonzero. If $||W_{ij}||_{4}<L<\infty$, then

[TABLE]

Proof.

We first show that $\phi_{3}+\phi_{4}+2\phi_{5}$ is the correct limiting variance. Writing the variance of the expression on the left hand side of (S29) explicitly and recalling that entries such that $i=j$ are undefined, we see

[TABLE]

by the properties of joint exchangeability of $\{W_{ij}\}_{i,j=1}^{n}$ as described in Section 3.3. This variance is finite and nonzero by assumption. To prove (S29), it is sufficient to show

[TABLE]

Define the limiting variance as $\sigma_{n}^{2}=n^{3}(\phi_{3}+\phi_{4}+2\phi_{5})$ and the sum $S_{n}=\sum_{ij}W_{ij}$.

To establish (S30), we employ Lemma S13, where $\nu_{n}$ is the probability measure corresponding to $\bar{S}_{n}$ for all $n$. The first condition of Lemma S13 is satisfied since

[TABLE]

for $C<\infty$ and all $n$. Thus, to prove (S30), it is sufficient to show the second condition of Lemma S13: for all $\lambda\in\mathbb{R}$, as $n\rightarrow\infty$,

[TABLE]

We decompose the term in the expectation as in guyon1995random and lumley2003asymptotics:

[TABLE]

[TABLE]

where $\bar{S}_{ij,n}=S_{ij,n}/\sigma_{n}$. To satisfy (S31) it remains to be shown that $\lim_{n\rightarrow\infty}E(A_{m})=0$ for each $m\in\{1,2,3\}$.

Beginning with ${A_{1}}$, $|{\rm e}^{i\lambda\bar{S}_{n}}|\leq 1$. Using this fact and Lemma S14,

[TABLE]

for all real-valued $\lambda$. $E(|A_{1}|)^{2}$ limiting to zero implies $E(|A_{1}|)$ limits to zero, and hence $E(A_{1})$ limits to zero.

Now, for ${A_{2}}$, by Taylor expansion of ${\rm e}^{-i\lambda\bar{S}_{ij,n}}$, we can write

[TABLE]

for some $0<c<\infty$ and all $n,\lambda$. Using this bound and the fact that $|\Theta_{ij}|=4n-6$, we evaluate $E(|A_{2}|)$ directly below:

[TABLE]

for all real $\lambda$. As $E(|A_{2}|)$ limits to zero, so does $E(A_{2})$.

Finally, for ${A_{3}}$, the expression ${S}_{ij,n}$ sums all terms in the sequence $\{W_{ij}\}_{i,j=1}^{n}$ that depend upon $W_{ij}$, including $W_{ij}$ itself. Thus, $W_{ij}$ and $\bar{S}_{n}-\bar{S}_{ij,n}$ are independent. It follows immediately that

[TABLE]

since $E\left(W_{ij}\right)=0$ for all ordered pairs $(i,j)$.

Hence, $\lim_{n\rightarrow\infty}E(A_{m})=0$ for each $m\in\{1,2,3\}$ and we have the convergence in (S31), implying $\bar{S}_{n}\rightarrow_{d}{\rm N}(0,1)$ by Lemma S13, which gives the desired result in (S29). ∎

S8.3 Proof of Theorem S9

We begin by writing

[TABLE]

again emphasizing that entries in the sum with $j=k$ are undefined and omitted. Addressing the first multiplicative term in (S32), we recall that the inverse map is continuous. Then, by (S24) and the continuous mapping theorem, we have

[TABLE]

We now analyze the second multiplicative term in (S32). Showing asymptotic normality of this term is sufficient to show asymptotic normality of the expression on the left hand side of (S32). Recall $x_{jk}^{T}=[x_{jk}^{(1)},x_{jk}^{(2)},\ldots,x_{jk}^{(p)}]$. We wish to show that the sum of vectors

[TABLE]

for some limiting variance $\Sigma$. By the Cramér-Wold device (cramer1936some), ${U}_{n}$ is asymptotically normal with asymptotic variance $\Sigma$ if and only if ${v}^{T}{U}_{n}$ is asymptotically normal with asymptotic variance ${v}^{T}\Sigma{v}$ for every vector ${v}\in\mathbb{R}^{p}$ such that $||{v}||=1$. Clearly,

[TABLE]

where we define $\widetilde{x}_{jk}={v}^{T}x_{jk}$. We wish to apply Theorem S15 to the sequence $\{\widetilde{x}_{jk}\xi_{jk}\}_{j,k=1}^{n}$. First, the condition of finite moments in (B2) of Theorem S9 and $||{v}||=1$ implies that $||\widetilde{x}_{jk}\xi_{jk}||_{4}<L$ for some finite $L<\infty$. Secondly, by the independence of $X$ and $\Xi$ in (B3) of Theorem S9, the sequence $\{\widetilde{x}_{jk}\xi_{jk}\}_{j,k=1}^{n}$ is a mean-zero exchangeable sequence of scalar random variables. Taking the variance directly,

[TABLE]

Then, we apply Theorem S15 with $\sigma_{n}^{2}=\text{var}(\sum_{jk}\widetilde{x}_{jk}\xi_{jk})$ in (S35), which gives that

[TABLE]

Thus, by the Cramér-Wold device, we get the desired joint asymptotic normality

[TABLE]

Combining the convergence in probability in (S33) and the asymptotic normality of (S36), we obtain the desired result. ∎

S9 Proof of consistency of the exchangeable estimator

S9.1 Notation and outline of proof

For the proof of the consistency of the exchangeable estimator $\widehat{V}_{E}$, we adopt the same change in notation as in Section S8, defined in (S25). We deviate slightly in that we denote $\Theta_{i}$ to denote dyadic pairs $(j,k)$ and $(l,m)$ that share a member in the $i$th manner. For example, for $i=3$ we must have $j=l$ and $m\neq k$. We use the same assumptions as in Theorem S9.

This proof is outlined as follows. We initially prove that the exchangeable estimator $\widehat{V}_{E}$ is consistent if the exchangeable parameter estimates $\{\widehat{\phi}_{i}:i=1,\ldots,5\}$ are consistent for the true parameters. We then prove consistency of $\{\widehat{\phi}_{i}\}$ in two steps: (a) we show parameter estimates $\{\widetilde{\phi}_{i}\}$ based on the unobserved true errors $\Xi$ are consistent and then (b) we show that the parameter estimates $\{\widehat{\phi}_{i}\}$ are asymptotically equivalent to $\{\widetilde{\phi}_{i}\}$. We require the consistency of $\widehat{\beta}$ result (implied by Theorem S9) for this last step.

S9.2 Proof of Theorem S10: consistency of the exchangeable estimator

First, from Theorem S9, the order of convergence of $\widehat{\beta}$ is ${n}^{1/2}$. Thus, we choose the rate $n$ as our asymptotic regime for consistency of $\widehat{V}_{E}$. We wish to show that

[TABLE]

First we show that to prove consistency of $\widehat{V}_{E}$, it is sufficient to prove the consistency of the parameter estimates $\{\widehat{\phi}_{i}\}$ for the true parameters. We begin by writing the difference of variances $n\widehat{V}_{E}-n\text{var}(\widehat{\beta})$ in (S37), as

[TABLE]

where $c_{i,n}=|\Theta_{i}|/n(n-1)^{2}$ and $h_{i,n}(X)$ contains the remaining terms which are functions of $X$. By the counting argument used to show Lemma S14, each $|\Theta_{i}|$ is at most ${O}(n^{3})$, so each $c_{i,n}\rightarrow d_{i}$ for some finite constant $d_{i}$. Namely, $d_{i}=0$ for $i\in\{1,2\}$, $d_{i}=1$ when $i\in\{3,4\}$, and $d_{i}=2$ for $i=5$. To obtain the result in (S37), it is sufficient then to show $\widehat{\phi}_{i}-\phi_{i}\rightarrow_{p}0$ and $h_{i,n}(X)$ converges in probability to some constant for all $i$. The latter comes easily, that is, by assumption and Slutsky’s theorem,

[TABLE]

The continuous mapping theorem allows us to take the probability limit of ${X^{T}X}/({n(n-1))}$ before inversion, as the inversion map is continuous.

We now consider consistency of the parameter estimates $\widehat{\phi}_{i}$. First, define error averages $\{\widetilde{\phi}_{i}:i=1,\ldots,5\}$ analogous to the parameter estimates, such that for each $i$,

[TABLE]

We will show $\widetilde{\phi}_{i}-\phi_{i}$ converges in probability to zero, and then do the same for $\widehat{\phi}_{i}-\widetilde{\phi}_{i}$. This is sufficient for showing $\widehat{\phi}_{i}-{\phi}_{i}\rightarrow_{p}0$ as $\widehat{\phi}_{i}-{\phi}_{i}=(\widehat{\phi}_{i}-\widetilde{\phi}_{i})+(\widetilde{\phi}_{i}-\phi_{i})$.

To show convergence in probability of $\widetilde{\phi}_{i}-\phi_{i}$ to zero, we use the argument that the bias and variance both tend to zero. By assumption (A1), $E[\xi_{jk}\xi_{\ell m}]=\phi_{i}$ for every relation pair $(jk,\ell m)\in\Theta_{i}$. Thus, $E[\widetilde{\phi}_{i}-\phi_{i}]=0$ for all $n$ and $i\in\{1,\ldots,5\}$. We now turn to the variance:

[TABLE]

We again make a counting argument similar to that in Lemma (S14). By condition (B2), each of the $|\Theta_{i}|^{2}$ covariances in the sum above are bounded. The covariance between $\xi_{jk}\xi_{\ell m}$ and $\xi_{rs}\xi_{tu}$ is nonzero only if there is overlap between their two index sets. This reduces the number of nonzero covariances from the maximum possible $|\Theta_{i}|^{2}$ by a factor of at least $n$. Again, consider the case of $i=3$ where $|\Theta_{3}|={O}(n^{3})$. Each pair of relations in $\Theta_{3}$ must be of the form $(jk,jm)$, and thus the second set of indices must be of the form $(js,ju)$, for example, for the covariance to be nonzero. The set of indices $\{j,k,j,m,j,s,j,u\}$ is of order ${O}(n^{5})=|\Theta_{3}|^{2}n^{-1}$. There are other forms of relation pairs in the second sum that give rise to nonzero covariance, such as $(ks,ku)$ and so on. However, there are nine such forms, each of which is ${O}(n^{5})$. Thus, the number of nonzero covariances is ${O}(n^{5})$, and hence, we have

[TABLE]

This same argument holds for all $i$, and thus, we have the desired consistency: $\widetilde{\phi}_{i}-\phi_{i}\rightarrow_{p}0$ for $i=1,\ldots,5$.

We now show that $\widehat{\phi}_{i}-\widetilde{\phi}_{i}$ converges in probability to zero. We first write the expression in terms of the estimated coefficients $\widehat{\beta}$:

[TABLE]

By Theorem S9, $\widehat{\beta}-\beta$ converges to zero in probability. By Slutsky’s theorem, if the terms in (S40) involving elements of $X$ and $\xi$ converge in probability to any constant, then $\widehat{\phi}_{i}-\widetilde{\phi}_{i}$ converges in probability to zero. By (B1) and (S24) we have the convergence in probability of the term involving $x_{jk}x_{\ell m}^{T}$. Furthermore, by condition (B3), we have that $E(\xi_{jk}x_{\ell m})=E(\xi_{\ell m}x_{jk})=0$. It remains to be shown that the variance of the error-covariate averages tend to zero. Consider the variance of the first error-covariate averages:

[TABLE]

In writing (S41), we use condition (B3) and simplify by conditioning on $X$ and using the law of total variance. By the same counting arguments used to establish (S39), there are $|\Theta_{i}|^{2}{O}(n^{-1})$ nonzero bounded covariances in (S41). Thus, we have

[TABLE]

Since the expectation and variance both tend to zero, we have

[TABLE]

The same argument applies to the second error-covariate term in (S40). Thus, we have shown that consistency of $\widehat{\beta}$ implies

[TABLE]

∎

S10 Proof of Theorem S11: mean-square error of the exchangeable and dyadic clustering estimators

In this section, we prove that the mean-square error of the estimator of $\text{var}(\widehat{\beta})$, conditional on $X$, is lower when using the exchangeable estimator than that when using the dyadic clustering estimator with high probability in $X$, assuming that the error structure is exchangeable. Before proving the theorem, we provide a lemma that states that the mean-square error of the each estimator is asymptotically equivalent to the mean-square error of each estimator based on the true errors, which vastly simplifies the proof of the mean-square error theorem. Even so, we must consider higher order moments of $\xi$ than the covariances ${\rm cov}\left(\xi_{jk},\xi_{lm}\right)$. So, we also provide a lemma in which we define the covariance of any pair of product of error relations ${\rm cov}\left(\xi_{jk}\xi_{lm},\xi_{rs}\xi_{tu}\right)$ and define the limiting values of the covariance of the error averages, $n{\rm cov}(\widetilde{\phi}_{v},\widetilde{\phi}_{w})$, for every pair $(v,w)\in\{1,\ldots,5\}\times\{1,\ldots,5\}$.

In this Section, we use the notation ${O}(n^{a})$ and ${\Theta}(n^{a})$, for some $a\in\mathbb{R}$, to denote the convergence a sequence of numbers to a constant (possibly zero) and a nonzero constant, respectively, as $n$ grows to infinity. In other words, $X_{n}={O}(n^{a})$ means that the sequence $n^{-a}X_{n}$ converges to a constant that may be zero. The notation $X_{n}={\Theta}(n^{a})$ means that the sequence $n^{-a}X_{n}$ converges to a nonzero constant. Lastly, it follows that $X_{n}={O}(n^{a-\epsilon})$ means that the sequence $n^{-a}X_{n}$ converges to zero.

We use similar notation for convergence of sequences of random variables. The notation $X_{n}={O}_{p}(n^{a})$ for some $a\in\mathbb{R}$ means that the sequence $n^{-a}X_{n}$ converges in distribution to a random variable (possibly a constant). The notation $X_{n}=o_{p}(n^{a})$ for some $a\in\mathbb{R}$ means that the sequence $n^{-a}X_{n}$ converges in probability to zero. Finally, we define $X_{n}={\Theta}_{p}(n^{a})$ to mean that $n^{-a}X_{n}$ converges in distribution to a random variable with distribution that is not a point mass at zero, and thus possibly a nonzero constant (as will always be the case in this section).

S10.1 Lemmas in support of Theorem S11

The first lemma describes the covariances of parameter estimates based on the errors, which arise in the proof of the mean-square error theorem. Of interest are the covariances ${\rm cov}(\widetilde{\phi}_{v},\widetilde{\phi}_{w})$ for $(v,w)\in\{3,4,5\}\times\{3,4,5\}$, as there are ${\Theta}(n)$ times as many of these covariances in $\widehat{V}_{E}$ as those covariances where at least one of $v$ or $w$ is in $\{1,2\}$. However, we provide limiting values of all covariances for completeness. The proof of this lemma follows from recognizing that $\widetilde{\phi}_{v}$ is a sample average and from defining all possible covariances that make up ${\rm cov}(\widetilde{\phi}_{v},\widetilde{\phi}_{w})$ and their multiplicities.

Lemma S16**.**

If ${\xi}$ is a mean zero random vector with positive definite covariance matrix in the exchangeable class, $\Omega=\sum_{i=1}^{5}\phi_{i}\mathcal{S}_{i}$, and $E[\xi_{jk}^{4}]<L<\infty$, then the covariance $n{\rm cov}\left(\widetilde{\phi}_{v},\widetilde{\phi}_{w}\right)$ for $(v,w)\in\{1,\ldots,5\}\times\{1,\ldots,5\}$ converges to

[TABLE]

and $C\left(v,w\right)_{i}$, $D\left(v,w\right)_{i}$, and $F\left(v,w\right)_{i}$ are unknown finite constants equal to ${\rm cov}\left(\xi_{jk}\xi_{lm},\xi_{rs}\xi_{tu}\right)$ for various configurations of the sets $\{j,k,l,m\}$ and $\{r,s,t,u\}$.

Proof.

By definition,

[TABLE]

The sum is over $|\Theta_{v}||\Theta_{w}|$ terms. Whenever $\{j,k,l,m\}\cap\{r,s,t,u\}=\varnothing$, the covariance is zero. This removes a power of $n$ from the sum in (S43), such that the sum is over ${O}(|\Theta_{v}||\Theta_{w}|n^{-1})$ possibly nonzero covariances. The scaled sum in (S43) converges – provided that the number of values that ${\rm cov}\left(\xi_{jk}\xi_{lm},\xi_{rs}\xi_{tu}\right)$ can take is finite – as each covariance is finite by assumption and the sequence of covariances is homogeneous as $n$ grows by exchangeability. In the remainder of the proof, we enumerate and define the covariances ${\rm cov}\left(\xi_{jk}\xi_{lm},\xi_{rs}\xi_{tu}\right)$ in (S43) for particular pairs $(v,w)\in\{1,\ldots,5\}\times\{1,\ldots,5\}$, showing that the number of values that ${\rm cov}\left(\xi_{jk}\xi_{lm},\xi_{rs}\xi_{tu}\right)$ can take is finite. This is sufficient to establish convergence.

We begin by analyzing the case of interest, that is when both $v$ and $w$ are members of $\{3,4,5\}$. As an example, we focus on $v=3$ and $w=4$, where the first product of error relations corresponds to the same-sender covariance (b) in Figure 1 and the second corresponds to the same-receiver covariance (c) in Figure 1. In this case, both $|\Theta_{3}|=|\Theta_{4}|={\Theta}(n^{3})$. When $v=3$ and $w=4$, the covariance in (S43) becomes

[TABLE]

where ‘$=_{a}$’ denotes equality in the limit as $n$ grows to infinity.

Only pairs of relation products that share a single actor will remain in the limit, as there are an order of $n$ fewer covariances resulting from pairs of relation products that share two actors. One such pair of relation products that share a single actor correspond to the case when $s=j$, i.e. ${\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{jr}\xi_{tr}\right)$, of which there are ${\Theta}(n^{5})$ in the sum in (S44). There are ${\Theta}(n^{4})$ covariances corresponding to the case when $s=j$ and $r=k$, i.e. ${\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{jk}\xi_{tk}\right)$. The values of all covariances in (S44) are finite by assumption and not equal in general. However, by exchangeability, covariances resulting from pairs of relations that share an actor in the same way are equal. For example, the covariance corresponding to $s=j$ is the same regardless of the node labeling, that is ${\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{jr}\xi_{tr}\right)={\rm cov}\left(\xi_{ab}\xi_{ac},\xi_{ad}\xi_{ed}\right)$ for any set $\{a,b,c,d,e\}\subset\{1,\ldots,n\}$ with $|\{a,b,c,d,e\}|=5$.

There are nine ways that we may have $|\{j,k,l\}\cap\{r,s,t\}|=1$ in (S44), i.e. there are nine ways that exactly one of $\{j,k,l\}$ equals exactly one of $\{r,s,t\}$. However, these reduce into four unique covariance values for each pair $(v,w)\in\{3,4,5\}\times\{3,4,5\}$. As an example, when $t=j$ the covariance is the same as that when $s=j$, that is ${\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{jr}\xi_{tr}\right)={\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{sr}\xi_{jr}\right)$. Now we define these four covariance values and their multiplicities out of the nine possible ways that exactly one of $\{j,k,l\}$ equals exactly one of $\{r,s,t\}$:

• •

When $r=j$, we define the covariance $C(3,4)_{1}={\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{sj}\xi_{tj}\right)$, of which there is one out of nine possible;

• •

When $s=j$, the covariance is the same as when $t=j$ (multiplicity two), and we define this covariance $C(3,4)_{2}={\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{jr}\xi_{tr}\right)$;

• •

When $r=k$, the covariance is the same as when $r=l$ (multiplicity two), and we define this covariance $C(3,4)_{3}={\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{sk}\xi_{tk}\right)$;

• •

We define the covariance when $s=k$ to be $C(3,4)_{4}={\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{kr}\xi_{tr}\right)$, of which there are four, the remaining terms of which correspond to $t=k$, $s=l$, and $t=l$.

Now, noting that there are $n^{5}+{\Theta}(n^{4})$ covariances in the sum (S44) corresponding to each of the nine possible ways that exactly one of $\{j,k,l\}$ equals exactly one of $\{r,s,t\}$, we see that

[TABLE]

where ‘$\rightarrow$’ denotes convergence in the limit as $n$ goes to infinity. Under appropriate definition of $C(v,w)_{i}$ for $i\in\{1,2,3,4\}$, the same argument applies when both $v$ and $w$ are one of $\{3,4\}$. When $w=5$ (relation products of the form $\{\xi_{jk}\xi_{kl}\}$ and $\{\xi_{jk}\xi_{lj}\}$) and $v=3$, however, we then must consider covariances ${\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{rs}\xi_{tr}\right)$ and ${\rm cov}\left(\xi_{jk}\xi_{jl},\xi_{sr}\xi_{rt}\right)$ from $w=5$, which doubles the coefficients in (S45). This accounts for $\beta_{w}=2$ when $w=5$ and $\beta_{w}=1$ otherwise in (S42). The same argument applies when $v=5$.

We now analyze both $v$ and $w$ in $\{1,2\}$, corresponding to variance and the reciprocal covariance (a) in Figure 1. In this case, both $|\Theta_{v}|=|\Theta_{w}|=n(n-1)$. Taking $v=1$ and $w=1$ as an example, in the limit, the covariance in (S43) is

[TABLE]

Again, we only consider covariances corresponding to pairs of relation products that share a single actor as only these covariances survive in the limit. There are four possible ways that $\{j,k\}$ shares exactly one actor with $\{r,s\}$. We define the three unique covariances and their multiplicities corresponding to the four ways that $\{j,k\}$ shares exactly one actor with $\{r,s\}$ as follows:

• •

When $r=j$, we define the covariance $F(1,1)_{1}={\rm cov}\left(\xi_{jk}^{2},\xi_{js}^{2}\right)$, of which there is one out of the four possibilities;

• •

When $s=k$, we define the covariance $F(1,1)_{2}={\rm cov}\left(\xi_{jk}^{2},\xi_{rk}^{2}\right)$, of which there is one;

• •

When $s=j$, we define the covariance $F(1,2)_{3}={\rm cov}\left(\xi_{jk}^{2},\xi_{rj}^{2}\right)$, which is the same as when $r=k$, accounting for the remaining two possibilities.

Now, the fact that there are $n^{3}+{\Theta}(n^{2})$ covariances in the sum (S46) corresponding to each of the four possible ways that $\{j,k\}$ shares exactly one actor with $\{r,s\}$ gives that

[TABLE]

Of course, the same argument applies to any $v$ and $w$ both in $\{1,2\}$. In the case where $v=1$ and $w=2$, by symmetry, $F(1,2)_{1}=F(1,2)_{2}$. Similarly, for $v=w=2$, all $F(2,2)_{1}=F(2,2)_{2}=F(2,2)_{3}$.

Similar counting arguments to those in the previous paragraphs apply when one of $v,w$ is in $\{3,4,5\}$ and the other is in $\{1,2\}$. As an example, consider $v=1$ and $w=3$. Once again, only pairs of relations that share a single actor will remain in the limit. Then, in the limit, the covariance in (S43) becomes

[TABLE]

Now, there are six ways in which the first pair of relations share an actor with the second pair, i.e. all sets with exactly one actor from $\{j,k\}$ equal to exactly one other from $\{r,s,t\}$. We define the covariances corresponding to the six possibilities below:

• •

When $r=j$, we define the covariance $D(1,3)_{1}={\rm cov}\left(\xi_{jk}^{2},\xi_{js}\xi_{jt}\right)$, of which there is one out of the six possibilities;

• •

When $r=k$, we define the covariance $D(1,3)_{2}={\rm cov}\left(\xi_{jk}^{2},\xi_{ks}\xi_{kt}\right)$, of which there is one;

• •

The overlaps where $s=j$ and $t=j$ result in the same covariance (multiplicity two), which we define $D(1,3)_{3}={\rm cov}\left(\xi_{jk}^{2},\xi_{rj}\xi_{rt}\right)$;

• •

The overlaps where $s=k$ and $t=k$ result in the same covariance (multiplicity two), which we define $D(1,3)_{4}={\rm cov}\left(\xi_{jk}^{2},\xi_{rk}\xi_{rt}\right)$.

Then, noting that there are $n^{4}+{\Theta}(n^{3})$ covariances in the sum (S47) corresponding to each of the six possible ways that exactly one actor from $\{j,k\}$ is equal to exactly one other from $\{r,s,t\}$, we have that $n{\rm cov}\left(\widetilde{\phi}_{2},\widetilde{\phi}_{3}\right)$ converges to

[TABLE]

When $D(v,w)_{k}$ for $k\in\{1,2,3,4\}$ is appropriately defined, the same argument applies for all settings where one of $v,w$ is in $\{3,4\}$ and the other is in $\{1,2\}$. When $w=5$ (relation products of the form $\{\xi_{jk}\xi_{kl}\}$ and $\{\xi_{jk}\xi_{lj}\}$), however, we then must consider covariances in (S47) ${\rm cov}(\xi_{jk}^{2},\xi_{rs}\xi_{tr})$ and ${\rm cov}(\xi_{jk}^{2},\xi_{sr}\xi_{rs})$, which doubles the coefficients in (S48). This accounts for $\beta_{w}=2$ when $w=5$ and $\beta_{w}=1$ otherwise in (S42). When $v=2$, for example, we have the simplification that $D(2,3)_{1}=D(2,3)_{2}$ and $D(2,3)_{3}=D(2,3)_{4}$. ∎

The expressions for the estimators based on the errors are simpler to analyze than those based on the residuals. For example, when comparing the mean-square errors of the exchangeable and dyadic clustering estimators, it is desirable to analyze $MSE_{\xi}(\widetilde{V}_{E}\mid X)$ instead of $MSE_{\xi}(\widehat{V}_{E}\mid X)$. The following lemma allows us to do just this. This lemma states that mean-square errors of the estimators based on the errors are asymptotically equivalent to the mean-square error of those based on the residuals. The proof consists of first evaluating the mean-square error conditional on $X$. We then show that $n^{3}MSE_{\xi}(\widetilde{V}_{E}\mid X)$ converges in $X$-probability to a nonzero constant in general and that $n^{3}MSE_{\xi}(\widetilde{V}_{E}\mid X)-n^{3}MSE_{\xi}(\widehat{V}_{E}\mid X)$ converges in $X$-probability to zero, implying that the difference between $MSE_{\xi}(\widetilde{V}_{E}\mid X)$ and $MSE_{\xi}(\widehat{V}_{E}\mid X)$ is asymptotically negligible. We repeat the procedure for $MSE_{\xi}(\widetilde{V}_{DC}\mid X)$ and $MSE_{\xi}(\widehat{V}_{DC}\mid X)$.

Lemma S17**.**

Assuming $E\left(|x_{jk}^{(l)}|^{8}\right)<L^{8}<\infty$ for all $(l=1,\ldots,p)$ and under the assumptions of Theorem S9, the mean-square error for both the exchangeable and dyadic clustering estimators based on the residuals is asymptotically equivalent to the mean-square error of each respective estimator based on the errors. That is,

[TABLE]

and analogously for dyadic clustering.

Proof.

We will focus on the exchangeable estimator first, and then the dyadic clustering estimator. Throughout, we drop the conditioning on $X$ in the mean-square error as it is understood, for example $MSE_{\xi}(\widehat{V}_{E})=MSE_{\xi}(\widehat{V}_{E}\mid X)$.

We begin with the exchangeable estimator. By definition, the mean-square error of the exchangeable estimator is

[TABLE]

where $V^{*}=\text{var}(\widehat{\beta})$, the true variance of $\widehat{\beta}$. By the Cauchy-Schwarz inequality,

[TABLE]

If we show that $n^{3}MSE_{\xi}\left(\widetilde{V}_{E}\right)$ converges in $X$-probability to a constant, i.e. $MSE_{\xi}\left(\widetilde{V}_{E}\right)={O}_{p}(n^{-3})$, and that $E\left\{\left(\widehat{V}_{E}-\widetilde{V}_{E}\right)^{2}\ \mid\ X\right\}={O}_{p}(n^{-4})$, then (S51) implies that the third additive term of (S50) is ${O}_{p}(n^{-7/2})$. This is sufficient to establish (S49). We begin with showing $n^{3}MSE_{\xi}\left(\widetilde{V}_{E}\right)={O}_{p}(1)$. By definition, the scaled mean-square error is n^{3}MSE_{\xi}\left(\widetilde{V}_{E}\right)=n^{3}E\left\{{\rm tr}\left(\widetilde{V}_{E}^{2}\right)\ \big{|}\ X\right\}, which is equal to

[TABLE]

By Lemma S16, $n{\rm cov}(\widetilde{\phi}_{v},\widetilde{\phi}_{w})$ converges to a finite constant for every $(v,w)\in\{1,\ldots,5\}\times\{1,\ldots,5\}$. The convergence in probability of each multiplicative term in (S52) containing $X$ is defined by assumption (B1); only those with both $v$ and $w$ in $\{3,4,5\}$ survive in the limit as these have $|\Theta_{v}|={\Theta}(n^{3})$ whereas $|\Theta_{v}|={\Theta}(n^{2})$ for $v\in\{1,2\}$. Thus, we have that

[TABLE]

which is finite., and where ‘$\rightarrow_{pr_{X}}$’ denotes convergence in $X$-probability.

It remains to show that $E\left\{\left(\widehat{V}_{E}-\widetilde{V}_{E}\right)^{2}\ \mid\ X\right\}={O}_{p}(n^{-4})$. To establish this fact, it is sufficient to show that $\widehat{V}_{E}-\widetilde{V}_{E}={O}_{p}(n^{-2})$, and then, by the continuous mapping theorem, $\left(\widehat{V}_{E}-\widetilde{V}_{E}\right)^{2}={O}_{p}(n^{-4})$, which implies the desired result. Writing directly,

[TABLE]

By assumption (B1), the multiplicative terms involving $X$ converge in probability to constants. To establish $\widehat{V}_{E}-\widetilde{V}_{E}={O}_{p}(n^{-2})$, it is sufficient show that $\widehat{\phi}_{v}-\widetilde{\phi}_{v}={O}_{p}(n^{-1})$ for all $v\in\{1,\ldots,5\}$. Writing this expression directly, the difference $\widehat{\phi}_{v}-\widetilde{\phi}_{v}$ is

[TABLE]

By Theorem S9, $\widehat{\beta}-\beta={O}_{p}(n^{-1/2})$. Also, by assumption (B1), the sum involving $X$ in the second term converges in probability to a constant; thus, the second additive term in (S55) is ${O}_{p}(n^{-1})$. Turning to the first additive term, we notice its expectation is zero, that is $E(x_{jk}\xi_{lm})=0$ for all relations $jk$ and $lm$. The variance is

[TABLE]

where we use the fact that $E\left(x_{jk}x_{rs}^{T}\xi_{lm}\xi_{tu}\right)$ is only nonzero when relation $lm$ shares an actor with relation $tu$ since $E\left(\xi_{lm}\xi_{tu}\right)=0$ whenever $\{j,k\}\cap\{l,m\}=\varnothing$ and $X$ is independent ${\xi}$ by assumption (B3). This fact removes at least a factor of $n$ from the sum. Thus, we have that $\sum_{(jk,lm)\in\Theta_{v}}{x_{jk}\xi_{lm}+x_{lm}\xi_{jk}}/{|\Theta_{v}|}={O}_{p}(n^{-1/2})$, which gives that $\widehat{\phi}_{v}-\widetilde{\phi}_{v}={O}_{p}(n^{-1})$ and

[TABLE]

which establishes (S49) for the exchangeable estimator.

The same argument following (S50) applies to the dyadic clustering estimator. To establish (S49) for the dyadic clustering estimator, it is thus sufficient to show that $n^{3}MSE_{\xi}\left(\widetilde{V}_{DC}\right)$ converges in $X$-probability to a constant and that $E\left\{\left(\widehat{V}_{DC}-\widetilde{V}_{DC}\right)^{2}\ \mid\ X\right\}={O}_{p}(n^{-4})$. We begin with the former. By definition, $n^{3}MSE_{\xi}\left(\widetilde{V}_{DC}\right)$ is

[TABLE]

where $\Theta_{0}$ is the set of relation pairs that share an actor in any manner. Then, substituting the asymptotic values for ${\rm cov}\left(\xi_{jk}\xi_{lm},\xi_{rs}\xi_{tu}\right)$ from Lemma S16 and separating the sum by the five ways that two relations may share an actor,

[TABLE]

where ‘$=_{a}$’ denotes equality in the limit and $T(v,w)_{i}$ is the set of relations $(jk,lm,rs,tu)$ such that $(jk,lm)\in\Theta_{v}$ and $(rs,tu)\in\Theta_{w}$ and such that the pairs of relations $(jk,lm)$ and $(rs,tu)$ share a single actor as appropriate for $C(v,w)_{i}$ in Lemma S16. In (S57), we substitute the limit of $\left({X^{T}X}/{n^{2}}\right)^{-2}$ from assumption (B1). Also in (S57), only terms with $v$ and $w$ both in $\{3,4,5\}$ survive in the limit as the set $|T(v,w)_{i}|={\Theta}(n^{5})$ (as detailed in Lemma S16), while the order is less for either $v$ or $w$ in $\{1,2\}$, so these terms vanish in the limit. Evaluating the vector products, the expression on the right hand side of (S57) equal to

[TABLE]

where $\left(m_{1}^{-2}\right)_{ab}$ is the $(a,b)$ entry in $M_{1}^{-2},$ e.g., and $x_{jk}^{(a)}$ is the entry in $X$ pertaining to column $a$ and relation $jk$. Further, the variance $\text{var}\left({n^{-5}}\sum_{T(v,w)_{i}}x_{lm}^{(a)}x_{rs}^{(b)}x_{tu}^{(c)}x_{jk}^{(d)}\right)$ is

[TABLE]

where $U(v,w)_{i}=T(v,w)_{i}$ and $(lm,rs,tu,jk)$ indexes the first sum and $(ef,gh,np,yz)$ indexes the second sum. The convergence is the result of the independence portion of assumption in (B1) and the bounded moment assumption on $X$. The variance in (S59) converges to zero for every set of covariates $\{a,b,c,d\}$, every relation type $v$ and $w$ both in $\{3,4,5\}$, and every covariance type $i\in\{1,2,3,4\}$. Thus, provided that the expectation of $n^{-5}\sum_{T(v,w)_{i}}x_{lm}^{(a)}x_{rs}^{(b)}x_{tu}^{(c)}x_{jk}^{(d)}$ converges to a constant, this expression converges in probability to that same constant. This expectation is

[TABLE]

where $\left(m_{v}\right)_{ij}$ is the $(i,j)$ entry in $M_{v}$ and we use the symmetry of $M_{v}$ for all $v\in\{3,4,5\}$. Unlike ${\xi}$, for a given $i\in\{1,2,3,4\}$, the covariances in (S60) may not be the same for two relation sets in $T(v,w)_{i}$. However, by assumption (B1) and taking $i=1$, $v=3$, and $w=4$ for example, we still have that

[TABLE]

for $|\{j,k,l,s,t,e,f,g,h,p\}|=10$. That is, covariances that share an actor in the same way are still equal. So, for fixed $i\in\{1,2,3,4\}$ and pair of $v$ and $w$ both in $\{3,4,5\}$, we may collect the $\alpha_{i}$ possible covariances and average them to attain the convergent value. We thus define the limit

[TABLE]

where $W(v,w)_{i}$ is the set of $\alpha_{i}$ ways that $(jk,lm,rs,tu)$ correspond to $T(v,w)_{i}$. For example, when $i=4$, $v=3$, and $w=4$, $W(v,w)_{i}$ contains four index sets corresponding to the four multiplicities of $C(v,w)_{4}$ as defined in Lemma (3). The convergence of (S62) results from (S61). As the average over $W(v,w)_{i}$ is over a finite number of terms since each $\alpha_{i}$ is bounded, there is no possibility of divergence. The covariances in (S62) are finite by assumption (B2) and $\beta_{v}\beta_{w}$ arises from the asymptotic limit of $n^{-5}|T(v,w)_{i}|$. Taking (S62) together with (S60), we have the convergence of the expectation of $n^{-5}\sum_{T(v,w)_{i}}x_{lm}^{(a)}x_{rs}^{(b)}x_{tu}^{(c)}x_{jk}^{(d)}$. Along with (S59), convergence of the expectation of $n^{-5}\sum_{T(v,w)_{i}}x_{lm}^{(a)}x_{rs}^{(b)}x_{tu}^{(c)}x_{jk}^{(d)}$ establishes the convergence of $n^{-5}\sum_{T(v,w)_{i}}x_{lm}^{(a)}x_{rs}^{(b)}x_{tu}^{(c)}x_{jk}^{(d)}$ to the same limit.

Now, for a particular $v$ and $w$ both in $\{3,4,5\}$ and $i\in\{1,2,3,4\}$, we collect the set of $\{C_{X}^{(d,a,b,c)}(v,w)_{i}\}$ for every $a,b,c,$ and $d$ in $\{1,2,\ldots,p\}$ from (S63) into a $p^{2}\times p^{2}$ matrix defined $D_{X}(v,w)_{i}$. Substituting this definition into (S58) while noting each $\{M_{v}\}_{v=1}^{5}$ is symmetric, the convergent value for the dyadic clustering estimator is

[TABLE]

Since the convergent value in (S64) is a finite constant, it remains to show that

[TABLE]

As with the exchangeable estimator, it is sufficient to show that $\widehat{V}_{DC}-\widetilde{V}_{DC}={O}_{p}(n^{-2})$. Using the residual definition $e_{jk}=\xi_{jk}+x_{jk}^{T}(\beta-\widehat{\beta})$, the expression for $\widehat{V}_{DC}-\widetilde{V}_{DC}$ is

[TABLE]

where we substitute the convergence of $\left({X^{T}X}/{n^{2}}\right)^{-1}$ to $M_{1}^{-1}$ and have used the exchangeability property to get the factor of two on the second additive term in the center of (S65). Analyzing the first additive term in the center of (S65),

[TABLE]

recalling the notation that $\Theta_{jk}$ is the set of all relations that share an actor with relation $jk$ and that $|\Theta_{jk}|={\Theta}(n)$. We attain the convergence rate by noting that the $X$-terms in (S66) converge in probability to constants by assumption (B1) and $\widehat{\beta}-\beta={O}_{p}(n^{-1/2})$ by Theorem S9. The convergences in (S66) are for $p\times p$ matrices; these convergences are element-wise.

We now analyze the convergence rate of the second additive term in the center of (S65),

[TABLE]

Again, the convergence of the first multiplicative term is a result of assumption (B1) and $\widehat{\beta}-\beta={O}_{p}(n^{-1/2})$ by Theorem S9. The mean $n^{-2}\sum_{lm}x_{lm}\xi_{lm}$ is expectation zero and ${O}_{p}(n^{-1/2})$ by previous arguments, for example in (S34). Thus, we have $\widehat{V}_{DC}-\widetilde{V}_{DC}={O}_{p}(n^{-2})$, and the dyadic clustering estimator satisfies the relation in (S49). ∎

S10.2 Proof of Theorem S11

We now establish that the mean-square error of the exchangeable estimator is less than that of the dyadic clustering estimator with high probability. To do so, we show that the value to which the difference in mean-square errors converges is nonnegative. Throughout the proof, we drop the conditioning on $X$ in the mean-square error as it is understood, for example $MSE_{\xi}(\widehat{V}_{E})=MSE_{\xi}(\widehat{V}_{E}\mid X)$.

The asymptotic difference in mean-square errorss is as follows, where we substitute the expressions for the estimators based on the errors in (S53) and (S64), as justified by Lemma S17:

[TABLE]

It remains to show that this is a nonnegative constant. To do so, we show that the matrix in the quadratic form in (S67) is the limit of a variance matrix, and thus positive semi-definite. We will show that the scaled variance

[TABLE]

converges to the desired matrix. First, the sum in (S68) is ${\Theta}(n^{5})$ as the relations $jk$ and $lm$ must share at least one actor with the relations $rs$ and $tu$ for the covariance to be nonzero. Then, by the arguments in Lemma S16, only pairs of relations that share a single actor survive in the limit. Finally, by assumption (B3), $X$ is independent ${\xi}$ and the variance in (S68) is asymptotically equivalent to

[TABLE]

where only terms with both $v$ and $w$ in $\{3,4,5\}$ survive in the limit. Then, by assumptions (A1) and (B1) and applying Lemma S16, the variance converges to

[TABLE]

where we substitute the definition of $D_{X}(v,w)_{i}$ following (S63). Thus, the matrix in (S69) is positive semi-definite. Now, (S67) becomes

[TABLE]

∎

S11 Simulation study

S11.1 Details of simulation study

We simulated data from a linear regression model with an exchangeable error model, independent and identically distributed errors, and a non-exchangeable error model. For each error setting, we employed the following three-covariate regression model:

[TABLE]

where $1_{(.)}$ denotes the indicator function. In this model, $\beta_{1}$ is an intercept; $\beta_{2}$ is a coefficient on a binary indicator of whether individuals $i$ and $j$ both belong to a pre-specified class $C\subset(1,\ldots,n)$; $\beta_{3}$ is a coefficient on the absolute difference of a continuous, actor-specific covariate $x_{3i}$; and $\beta_{4}$ is that for a pair-specific continuous covariate $x_{4ij}$. For all simulations, we fixed true coefficients $\beta=(1,1,1,1)^{T}$. We drew each $x_{2i}$ from a Bernoulli$(1/2)$ distribution independently. In the rare event that $x_{2i}=x_{2j}$ for all $(i,j)$ pairs, one realization $x_{2k}$ was randomly flipped to a $1$ or [math]. All $x_{3i}$ and $x_{4ij}$ were drawn independently from a standard normal distribution.

Each error setting was specified to have the same total variance:

[TABLE]

This variance was chosen so that the variance of the error would be similar to that of the regression mean model $\beta^{T}x_{ij}$. In the independent and identifically distributed errors setting, $\xi_{ij}\sim_{iid}{\rm N}(0,3)$ for all $(i,j)$. To generate the non-exchangeable errors, a mean-zero random effect was added to the upper left quadrant of $\text{var}(\xi)$ under independent and identically distributed errors. The errors for the non-exchangeable error setting may be written

[TABLE]

where ‘${1}_{(j\leq\left\lfloor{n/2}\right\rfloor)}$’ is an indicator of index $i$ less than or equal to the floor of $n/2$, for example. Since noise was added to actor relations in the same position in $\Xi$ in each simulation run, the distribution of the relations was not exchangeable: the distribution of the errors would be different for a reordering of the rows and columns of the error array.

Finally, the distribution of the exchangeable (bilinear mixed effects model) error setting is defined in hoff2005bilinear and detailed in full below. The bilinear mixed effects model is a generalization of the “additive common shocks” error structure used in simulation studies to justify the dyadic clustering estimator (cameron2014robust; aronow2015cluster; tabord2017inference).

[TABLE]

where $a_{i}$, $b_{j}$, $z_{i}$, $z_{j}$, and $\epsilon_{ij}$ are independent. We selected the dimension of the latent space to be $d=2$, the correlation between sender and receiver effects as $\rho_{ab}=1/2$, and the sender variance to be twice that of receiver variance: $\sigma_{a}^{2}=2\sigma_{b}^{2}$. We further specified $\sigma_{z}=\sigma_{\gamma}=\sigma_{b}$. Finally, we selected $\sigma^{2}_{\epsilon}=3/4$. With the aforementioned choices, the restriction $\sum_{ij}\text{var}(\xi_{ij})=3n(n-1)$ generated a quadratic equation in $\sigma_{b}$. The standard deviations that resulted from solving this quadratic equation are shown in Table S2.

We generated 500 random realizations of covariates for relational data sets with numbers of actors $n=20,40,80,160,320$. For each covariate realization, 1,000 random error realizations were generated for each of the three error settings: independent and identically distributed, exchangeable, and non-exchangeable. Using the model (S70), a simulated data set was created from each covariate realization and error realization pair. The regression model was fit to each data set using ordinary least squares, and standard errors were estimated using the exchangeable and dyadic clustering sandwich variance estimators. Confidence interval coverage was estimated for each covariate realization by counting the fraction of confidence intervals that contained the true coefficient. Bias and variance of the standard error estimators was estimated for each covariate realization relative to the known true standard errors of the ordinary least squares estimator, conditional on the covariate realization.

S11.2 Bias in simulation study

In addition to the estimated coverages in Figure 2, we plot the estimated bias of the standard error estimators in Figure S2. We note that the bias of the dyadic clustering estimator is generally substantially larger in absolute value, especially in the independent and exchangeable error cases, and most standard error estimators are biased downwards.

To verify Theorem 2, we compared the biases of the dyadic clustering and exchangeable estimators of the variances of the regression coefficients in the simulation studies with exchangeable and independent errors. We computed the mean bias for each variance estimator across all simulations for each error structure, data set size, coefficient, and estimator. Then, we took the ratio of mean absolute biases of the dyadic clustering variance estimator to that of the exchangeable estimator. The ratios, averaged for all $n$, are given in Table S3. We note that, in all cases, the mean bias of the dyadic clustering estimator was negative, and larger in absolute value than that of the exchangeable estimator. The fact that every ratio is greater than 2 confirms Theorem 2.

S12 Dyadic clustering covariance matrix invertibility

Ideally, for a covariance matrix estimate $\widehat{\Omega}$ to be of utmost utility, it must be invertible. For example, if we wish to employ feasible generalized least squares, the estimate of the the covariance matrix must be nonsingular. However, in many cases the dyadic clustering estimator is singular and hence cannot be inverted. In cases when the dyadic clustering estimator $\widehat{\Omega}_{DC}$ is singular, it can still be used to produce the dyadic clustering estimator $\widehat{V}_{DC}$ in (5).

Theorem S18**.**

The dyadic clustering estimate of the error variance, $\widehat{\Omega}_{DC}$, is singular for directed data with $R=1$.

Proof.

The dyadic clustering estimator can be written as the Hadamard product between the outer product of the residuals and a matrix of indicators of whether the dyad indices share a member,

[TABLE]

where ${1}{\{(i,j)\cap(k,l)\neq\varnothing\}}$ is an $n(n-1)\times n(n-1)$ matrix of indicators of indices $(i,j)$ sharing an index with $(k,l)$. The rank of the outer product of the residuals is one: ${\rm rank}(ee^{T})=1$. The rank of the indicator matrix is at most $n(n-1)/2$, since the indices $(i,j)$ share a member with an arbitrary pair $(k,l)$ if and only if the indices $(j,i)$ do as well. Thus, the column of ${1}{\{(i,j)\cap(k,l)\neq\varnothing\}}$ corresponding to $(i,j)$ is the same as that corresponding to $(j,i)$.

For any two square matrices of equal size $A$ and $B$, ${\rm rank}(A\circ B)\leq{\rm rank}(A){\rm rank}(B)$. Thus,

[TABLE]

$\widehat{\Omega}_{DC}$ is therefore not full rank. ∎

Remark S19**.**

Theorem S18 does not hold for undirected data when $R=1$. If the data are undirected, then the bound does not guarantee singularity of $\widehat{\Omega}_{DC}$ since the dimension of $\widehat{\Omega}_{DC}$ is exactly $n(n-1)/2$. In practice, we find that $\widehat{\Omega}_{DC}$ may be full rank in this case.

Remark S20**.**

The result of Theorem S18 holds for both directed and undirected data when $R>1$. In this case, the column in the indicator matrix ${1}{\{(i,j)\cap(k,l)\neq\varnothing\}}$ corresponding to the indices $(i,j,s)$ is the same as that column corresponding to $(i,j,r)$ for all values of $r=1,\ldots,R$. Thus, again $\widehat{\Omega}_{DC}$ is not full rank.

S13 Efficient inversion of the exchangeable covariance matrix

To perform the generalized least squares estimation procedure as described in Section 5, we must invert the exchangeable covariance matrix ${\Omega}_{E}$ as defined in Figure S1. For now, we work in the case where $R=1$. Since ${\Omega}_{E}$ is a real symmetric matrix, its inverse is real and symmetric as well. However, we can say more about the patterns in the inverse ${\Omega}_{E}^{-1}$. Recall that ${\Omega}_{E}$ has at most six unique terms; call these parameters $\phi$. We find that the inverse ${\Omega}_{E}^{-1}$ has at most six unique terms as well. If we define the parameters in ${\Omega}_{E}^{-1}$ as ${p}$, we can write

[TABLE]

where $I$ is the $n(n-1)\times n(n-1)$ identity. Lastly, we make the conjecture that the parameter pattern in ${\Omega}_{E}^{-1}$ is exactly the same as that in ${\Omega}_{E}$; we find this conjecture to be true in practice. One caveat is that the locations in which we assume zeros in ${\Omega}_{E}$ are not zero in ${\Omega}_{E}^{-1}$ in general.

We can find the inverse parameters ${p}$ from ${\phi}$ without inverting the entire matrix ${\Omega}_{E}$ by instead solving the following linear system

[TABLE]

where $C({\phi},n)$ is a set of six linear equations based on the parameters ${\phi}$ and the number of actors $n$ and is depicted in Figure S3. Thus, we replace the need to invert the $n(n-1)\times n(n-1)$ matrix $\Omega_{E}$ by the inversion of the $6\times 6$ matrix $C({\phi},n)$. Using this procedure, the computational cost associated with finding the inverse of $\Omega_{E}$ is independent of the number of actors $n$.

Now consider the case of array data with $R>1$. Inversion of the exchangeable covariance $\Omega_{E}$ in Figure 1 requires consideration of the patterns in the block matrices. Recalling that $\Omega_{E}$ is parametrized by twelve terms when $R>1$, we denote the first six parameters as ${\phi}^{(1)}$ and the second six ${\phi}^{(2)}$, corresponding to $\Omega_{1}$ and $\Omega_{2}$ respectively. Again the inverse $\Omega^{-1}$ has the exact same block matrix pattern as $\Omega$. Thus, the inverse may be parametrized by ${p}^{(1)}$ and ${p}^{(2)}$, each with length six, defined by the following linear equations.

[TABLE]

This is twelve linear equations in ${p}^{(1)}$ and ${p}^{(2)}$. In this formulation we reduce a $Rn(n-1)\times Rn(n-1)$ inversion to a $12\times 12$ inversion for calculation of $\Omega_{E}^{-1}$. Again, there is no dependence of the complexity of the inversion on the array dimensions $n$ and $R$.

S14 Trade data analysis

S14.1 Trade model details

We analyzed international trade among 58 countries over 20 years. The data were analyzed and made available (DOI:10.1214/10-AOAS403) by westveld2011mixed, which used a modified gravity mean model to represent the logarithm of yearly trade between each pair of countries as linear function of seven covariates in years 1981-2000. The gravity model, proposed by tinbergen1962shaping, posits that the total trade between countries is proportional to overall economic activity of the countries weighted by the inverse of the distance between them (raised to a power). Following ward2007persistent we added an indicator for whether the nations’ militaries cooperated in the given year and a measure of democracy, i.e. polity, which ranges from [math] (highly authoritarian) to $20$ (highly democratic). The complete model has the form:

[TABLE]

where $\text{log}\mbox{Trade}_{ijt}$ is the (log) volume of trade sent from country $i$ to country $j$ in year $t$; $\text{log}\text{GDP}_{it}$ and $\text{log}\text{GDP}_{jt}$ are the (log) Gross Domestic Product of nations $i$ and $j$, respectively; $\text{log}\text{D}_{ijt}$ is the (log) geographic distance between nations; $\text{CC}_{ijt}$ is the measure of cooperation in conflict (coded as +1 if nations were on the same side of a dispute and -1 if they were on opposing sides); and $\text{Pol}_{it}$ and $\text{Pol}_{jt}$ are the polity measures for $i$ and $j$, respectively, where polity ranges from [math] (highly authoritarian) to $20$ (highly democratic).

S14.2 Estimation

Let $\beta$ denote the length $8T$ vector of regression coefficients defined in (S74), where $t=1,\ldots,T$ and $T=20$ is the final time period. We estimated $\beta$ using feasible generalized least squares as in (2), where $y$ is a vector consisting of $(\text{log}\mbox{Trade}_{ijt})$ where $(i=1,\ldots,n;j=1,\ldots,n;i\neq j;t=1,\ldots,T)$, $X$ is a matrix of corresponding covariates defined by the model in (S74), and we assume $\Omega$ is jointly exchangeable. We initialized $\beta$ using ordinary least squares. Then, we estimated $\Omega$ using the residuals and the joint exchangeability assumption, as in (6) and (7). We then used (2) with $\Omega=\widehat{\Omega}$ to re-estimate $\beta$. We iterated estimates of $\Omega$ and $\beta$ until convergence. We defined convergence as an absolute change in weighted residual inner product being less than a particular tolerance $\epsilon>0$. Explicitly, convergence is defined to occur when

[TABLE]

where $\widehat{\beta}^{\gamma}$ and $\widehat{\Omega}^{\gamma}$ are the estimators at iteration $\gamma$ of the regression coefficients and error covariance matrix, respectively. For the trade data analysis, we take $\epsilon=10^{-6}$. After convergence, we set $\Omega=\widetilde{\Omega}$ in (4) to obtain standard errors for the estimated regression coefficients, which is the standard approach in feasible generalized least squares.

S14.3 In-sample evaluation

We compared the coefficients and confidence intervals resulting from the proposed feasible generalized least squares approach to the estimated coefficients and credible intervals from the mixed effects model of westveld2011mixed, and to coefficients estimated by ordinary least squares and confidence intervals using the dyadic clustering estimator $\widehat{V}_{DC}$. The estimated coefficients and corresponding 95% confidence intervals and posterior credible intervals are shown in Figure S4.

S14.4 Prediction study

To compare the ability of the proposed exchangeable approach and the model from westveld2011mixed – which we refer to as the hierarchical, longitudinal mixed effects model – to represent the trade data, we examined the out-of-sample predictive performance of the estimators. We detail the computation of predicted values here. First, we establish some notation for this section. Recall that the trade data set has covariate measures in $X$ that vary by year. Thus, we rewrite the linear model in (S74) as

[TABLE]

where $y_{t}$ represents the $n(n-1)$ vector of relations among the $n=58$ countries in year $t$, $X_{t}$ is a matrix eight covariates corresponding to year $t$, and $\Xi_{t}$ is a vector of errors for year $t$.

We estimated both models on the first 4 through 19 years of data, and then used these estimates to predict the trade in the following year. To generate predictions from each model, we computed the conditional expectation $E({y}_{T}\mid\{y_{t}\}_{t=1}^{T-1})$ based on the assumption that $y_{T}$ and $(y_{r})$, for $r=1,\ldots,T-1$, are jointly normal. The exchangeable approach and mixed effects model correspond to different models of the variance-covariance matrices of $\Xi_{t}$ and the covariances matrices between vectors $\Xi_{t}$ and $\Xi_{t+h}$. As a baseline, we included ordinary least squares, assuming independence of each year, i.e. $\Xi_{t}$ independent of $\Xi_{t+h}$.

We first discuss the ordinary least squares estimation and prediction procedure used for this study, as this approach is the simplest. We estimated the coefficients in the prediction time period, $\beta_{T}$, with the coefficients from the previous time period and assumed independent and identically distributed entries in all $\Xi_{t}$. Thus, the ordinary least squares estimator of the trade in year $T$ is

[TABLE]

For the exchangeable procedure, we again set $\widehat{\beta}_{T}=\widehat{\beta}_{T-1}$. Based on the model underlying the exchangeable procedure (Figure 6(a)), the variance $\text{var}(y_{t})=\Omega_{1}$ for all $t=1,\ldots,T$ and the covariance $\text{cov}(y_{t},y_{t+h})=\Omega_{2}$ for all $h$. Further, it can be shown that the precision corresponding to the concatenated vector $z_{T-1}^{T}=(y_{1}^{T},y_{2}^{T},\ldots,y_{T-1}^{T})$ is of the same pattern as the variance $\text{var}(z_{T-1})$, which has $\Omega_{1}$ along the diagonal blocks and $\Omega_{2}$ in the off-diagonal blocks (as in Figure 6(a)). We define the diagonal blocks of $\text{var}(z_{T-1})^{-1}$ as $\Psi_{1}$ and the off-diagonal blocks $\Psi_{2}$. Then, when the relations $\{y_{t}\}_{t=1}^{T}$ are jointly normally distributed, we define the prediction from the exchangeable procedure is

[TABLE]

Finally, to detail predictions from the mixed effects model, we first review some terms used in westveld2011mixed. The proposed model is

[TABLE]

where $s_{i,t}$ and $r_{j,t}$ are sender and receiver random effects, respectively, and $g_{ij,t}$ is a reciprocal random effect. These effects evolve according to autoregressive order one processes, such that

[TABLE]

where $\Phi_{sr}$ is a $2\times 2$ autoregressive matrix and $\Phi_{g}$ is a symmetric autoregressive matrix. The error terms $\varepsilon_{i,t}$ and ${e}_{ij,t}$ are mean-zero independent bivariate Gaussian random variables. Please see westveld2011mixed for more details of the model. With the notation and model defined, and again setting $\widehat{\beta}_{T}=\widehat{\beta}_{T-1}$, algebra reveals that the prediction from the mixed effects model is

[TABLE]

To estimate the mixed effects model, westveld2011mixed generated a Markov chain of length 55,000, of which the first 10,000 were discarded and every 20th iteration saved, giving 2,250 samples to approximate the posterior distributions of $\beta_{T-1}$, $\Phi_{sr}$, $\Phi_{g}$, and the random effects at $t=T-1$ in the above equations. The mean of these 2,250 samples from the joint posterior distribution, for each pair $i\neq j$, was used to construct $E(y_{T}\mid\{y_{t}\}_{t=1}^{T-1})_{MEM}$.

We evaluate performance in the predicted year using the mean square prediction error, defined

[TABLE]

where the expectation is replaced by the expectation of the appropriate prediction estimator (ordinary least squares, exchangeable, or mixed effects).

S15 A test for exchangeability

S15.1 Motivation and details

In this section we briefly motivate and describe a permutation test for assumption of exchangeability of the errors. For simplicity, we focus on the case where ordinary least squares is used to estimate regression coefficients $\beta$, and where $R=1$.

If the errors are truly non-exchangeable in a way that impacts inference, then the dyadic clustering estimator should be extreme relative to the distribution of the dyadic clustering estimator under the assumption of exchangeability. The idea of the proposed test is to generate the null (exchangeable) distribution using permutations that mimic the exchangeability assumption.

If the exchangeable assumption is correct, then by definition, a simultaneous permutation of the rows and columns of the error matrix leaves the distribution unchanged. That is,

[TABLE]

where $pr(\xi_{ij})$ denotes the joint distribution of the set $(\xi_{ij})$ over all relations such that $i\neq j$. Although the errors are unobservable, the residuals are observable and can be used to approximate the errors. Thus, the residuals should be approximately exchangeable, such that

[TABLE]

where $e_{ij}=y_{ij}-\widehat{\beta}^{T}x_{ij}$ is the residual for relation $(i,j)$. We rewrite the dyadic clustering estimator of the variance $\text{var}(\widehat{\beta}\mid X)$ as

[TABLE]

where $\Theta_{0}$ represents the set of pairs of relations $(jk,lm)$ that share an actor and where $\widetilde{x}_{jk}=(X^{T}X)^{-1}x_{jk}$. Applying the joint exchangeability assumption to the residuals, one may generate an approximate sample of the dyadic clustering estimator under the null assumption of exchangeability

[TABLE]

We use an $L_{2}$ norm to evaluate the discrepancy between the observed dyadic clustering estimator based on the observed residual matrix in (S75) and the distribution of dyadic clustering estimators under the exchangeable null in (S76). The squared distance can be computed using

[TABLE]

where $\bar{V}_{DC}^{(\pi)}$ is the element-wise sample mean dyadic clustering estimator of the exchangeable permutations of $\widetilde{V}_{DC}^{(\pi)}$ in (S76). A distribution of these squared distances under the exchangeable null may be obtained by replacing the dyadic clustering estimator based on the observed residuals with the permuted samples $\widetilde{V}_{DC}^{(\pi)}$,

[TABLE]

If $\bar{d}^{(\pi)}$ in (S77) is extreme relative to the distribution of ${d}^{(\pi)}$ in (S78) under the exchangeable null, then this constitutes evidence that the errors are non-exchangeable. We obtain an empirical $p$-value by computing the fraction of the permuted, squared distances ${d}^{(\pi)}$ that are greater than $\bar{d}^{(\pi)}$.

S15.2 Evaluating the test

We simulated from exchangeable and non-exchangeable relational regression models to assess the ability of our proposed hypothesis testing procedure to detect non-exchangeability. We evaluated the power and size of the test by simulating data from null and alternative models and evaluating the test on each simulated dataset.

We begin by specifying the null exchangeable model. In addition to an intercept, we include a single relational covariate, the coefficient of which is the target of inference. The model is specified as

[TABLE]

Such a model produces true variance parameters $\phi^{T}=[1,1/2,1/4,1/4,-1/4,0]$. Further, $\text{var}(x_{ij})=1$, such that the variance of the errors and variance of the covariate $x_{ij}$ are equivalent. We take $\beta_{0}=\beta_{1}=1$.

The non-exchangeable model we consider can be expressed as follows:

[TABLE]

Under this model, both $X$ and the errors are non-exchangeable and the variance of the errors errors depend on values of $X$. The numerical prefactors were chosen such that $\text{var}(x_{ij})\approx 1$ and $\text{var}(\xi_{ij})\approx 1$ as in the exchangeable model in (S79).

To evaluate the power and size of the test, we simulated 1,000 datasets from each of the null and alternative models in (S79) and (S80), respectively. For each simulated dataset, we evaluated the empirical two-sided $p$-value using 1,000 permutations of the form of (S76), and rejected the null hypothesis at the $\alpha=0.05$ level. The test results in Table S4 show approximate size, i.e. the probability of rejecting the null when the null is true, of $0.025$. Under the considered alternative distribution, Table S4 gives the estimate approximate power of about 0.839. Thus, for the given alternative model at moderate relational data size ($n=50$), the proposed permutation test has a good chance of detecting the non-exchangeability in the alternative model, while failing to reject a true null with probability less than $0.05$.