Dyadic dependence

Let $\left[n\right]\equiv \{1,2,\dots ,n\}$ be a set of individuals. Suppose I have data $\left\{\right(y_{ij},x_{ij}):i,j\in [n]\text{with}i<j\}$ on pairs in $\left[n\right]$ generated by the process $$y_{ij}=x_{ij}\beta +{\epsilon }_{ij},$$ where $x_{ij}$ is a row vector of pair $\{i,j\}$'s characteristics, $\beta$ is a vector of coefficients to be estimated, and ${\epsilon }_{ij}$ is a random error term with zero mean and zero correlation with the $x_{ij}$. For example, $\left[n\right]$ could be the nodes in a network, $x_{ij}$ the dimensions along which nodes $i$ and $j$ interact, and $y_{ij}$ the outcome of such interaction.

We can rewrite the data-generating process (DGP) in matrix form as $$y=X\beta +\epsilon ,$$ where $y$ is the vector of outcomes, $X$ is the design matrix, and $\epsilon$ is the vector of errors. Here $X$ has $$N\equiv \frac{n(n-1)}{2}$$ rows, each corresponding to a(n unordered) pair of individuals in $\left[n\right]$. Since the $x_{ij}$ and ${\epsilon }_{ij}$ are uncorrelated, the ordinary least squares estimator $$\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}y$$ of $\beta$ is unbiased. However, $\hat{\beta }$ may not be efficient because the errors ${\epsilon }_{ij}$ may be correlated. For example, if $${\epsilon }_{ij}=u_i+u_j+v_{ij}$$ with $u_i$, $u_j$, and $v_{ij}$ independent then $$\mathrm{Cov}({\epsilon }_{ij},{\epsilon }_{jk})=\mathrm{Var}\left(u_j\right).$$ Intuitively, the pairs $\{i,j\}$ and $\{j,k\}$ are linked through individual $j$, and so any errors specific to that individual affect the errors for both pairs. Consequently, the homoskedastic estimator $${\widehat{\mathrm{Var}}}_{\text{Hom.}}\left(\hat{\beta }\right)={\hat{\sigma }}^2(X^{T}X{)}^{-1}$$ with $${\hat{\sigma }}^2=\frac{1}{N}\sum _{ij}{\hat{\epsilon }}_{ij}^2$$ and $${\hat{\epsilon }}_{ij}=y_{ij}-x_{ij}\hat{\beta }$$ will typically under-estimate the variance in $\hat{\beta }$ by failing to account for linked pairs having dependent errors.

So, how can we account for such dependence? Consider the “sandwich” form $$\mathrm{Var}\left(\hat{\beta }\right)=BMB$$ of the (co)variance matrix for $\hat{\beta }$, where $B=(X^{T}X{)}^{-1}$ is the “bread” matrix and $M=X^{T}VX$ is the “meat” matrix with $V=\mathrm{Var}\left(\epsilon \right)$ the error (co)variance matrix. We need to estimate $M$ because we don’t observe the ${\epsilon }_{ij}$. Indexing pairs by $p$, the homoskedastic estimator defined above uses $$\begin{array}{cc}{\hat{M}}_{\text{Hom.}}& ={\hat{\sigma }}^2{X}^{T}X\\ & ={\hat{\sigma }}^2\sum _{p=1}^N{x}_p^T{x}_p,\end{array}$$ which assumes all errors have equal variance. In contrast, White (1980) suggests using $$\begin{array}{cc}{\hat{M}}_{\text{White}}& =X^{T}diag\left({\hat{\epsilon }}_p^2\right)X\\ & =\sum _{p=1}^N{\hat{\epsilon }}_p^2{x}_p^T{x}_p,\end{array}$$ which allows for unequal error variances (heteroskedasticity). But neither ${\hat{M}}_{\text{Hom.}}$ nor ${\hat{M}}_{\text{White}}$ allow for dyadic dependence among the errors. To that end, Aronow et al. (2017) suggest augmenting White’s estimator via $$\begin{array}{cc}{\hat{M}}_{\text{Aronow}}& ={\hat{M}}_{\text{White}}+\sum _{p=1}^N\sum _{q\in D\left(p\right)}{\hat{\epsilon }}_p{\hat{\epsilon }}_q{x}_p^T{x}_q,\end{array}$$ where $D\left(p\right)$ is the set of pairs $q\ne p$ linked to $p$ by a shared individual. We can express ${\hat{M}}_{\text{Aronow}}$ in matrix form as $${\hat{M}}_{\text{Aronow}}=X^T(D\odot \hat{\epsilon }{\hat{\epsilon }}^T)X,$$ where $D=\left(d_{pq}\right)$ is the dyadic dependence matrix with $$d_{pq}=\{\begin{array}{cc}1& \text{if pairs}p\text{and}q\text{are linked}\\ 0& \text{otherwise},\end{array}$$ and where $\odot$ denotes element-wise multiplication. Aronow et al. show that, under mild conditions, $B{\hat{M}}_{\text{Aronow}}B$ is a consistent estimator for $\mathrm{Var}\left(\hat{\beta }\right)$ when the data exhibit dyadic dependence.¹

To see Aronow et al.‘s estimator in action, suppose the DGP is given by the system $$\begin{array}{cc}{y}_{ij}& =\beta x_{ij}+{\epsilon }_{ij}\\ x_{ij}& =z_i+z_j\\ {\epsilon }_{ij}& =u_i+u_j+v_{ij},\end{array}$$ where $z_i$, $z_j$, $u_i$, $u_j$ and $v_{ij}$ are iid standard normal, and $\beta =1$ is the (scalar) coefficient to be estimated. Both the $x_{ij}$ and the ${\epsilon }_{ij}$ exhibit dyadic dependence, so we expect the homoskedastic and White estimators to under-estimate the true variance in $\hat{\beta }$. Indeed, the box plots below show that Aronow et al.‘s estimator is less biased than the homoskedastic and White estimators, and gets more accurate as the number of individuals $n$ grows.

Aronow et al.‘s estimator can also be applied to generalized linear models. For example, suppose $$y_{ij}=\{\begin{array}{cc}1& \text{if nodes}i\text{and}j\text{are adjacent}\\ 0& \text{otherwise}\end{array}$$ is an indicator for the event in which nodes $i$ and $j$ are adjacent in a network. We can model the link formation process as $$Pr(y_{ij}=1)={\Lambda }^{-1}(x_{ij}\beta +{\epsilon }_{ij}),$$ where $\Lambda \left(x\right)\equiv \log \left(x/\right(1-x\left)\right)$ is the logit link function. The logistic regression estimate $\hat{\beta }$ of $\beta$ reveals how the observable characteristics $x_{ij}$ of nodes $i$ and $j$ determine their probability of being adjacent. We can estimate the variance of $\hat{\beta }$ consistently by letting ${\hat{P}}_{ij}={\Lambda }^{-1}\left(x_{ij}\hat{\beta }\right)$ be the predicted probability for pair $\{i,j\}$, replacing the bread matrix $B=(X^{T}X{)}^{-1}$ with $$\hat{B}={\left(X^{T}diag\left({\hat{P}}_{ij}(1-{\hat{P}}_{ij})\right)X\right)}^{-1},$$ and computing $\hat{B}{\hat{M}}_{\text{Aronow}}\hat{B}$. My co-authors and I use this approach in “Research Funding and Collaboration:” we estimate how grant proposal outcomes determine the probability with which pairs of researchers co-author, and we compare ${\hat{\sigma }}^2\hat{B}$ and $\hat{B}{\hat{M}}_{\text{Aronow}}\hat{B}$ to show that our inferences are robust to dyadic dependence.