# Estimating the correlation in network disturbance models

A. D. Barbour111Institut für Mathematik, Universität Zürich,
Winterthurertrasse 190, CH-8057 ZÜRICH.  and Gesine Reinert222Department of
Statistics, University of Oxford, 24–29 St Giles’, OXFORD OX1 3LB, UK. GDR was
supported in part by EPSRC grant EP/R018472/1 and by the COSTNET COST Action
CA 15109.
Universität Zürich and University of Oxford

###### Abstract

The Network Disturbance Model of [3] expresses the dependency between
observations taken at the vertices of a network by modelling the correlation
between neighbouring vertices, using a single correlation parameter $\rho$. It
has been observed that estimation of $\rho$ in dense graphs, using the method
of Maximum Likelihood, leads to results that can be both biased and very
unstable. In this paper, we sketch why this is the case, showing that the
variability cannot be avoided, no matter how large the network. We also
propose a more intuitive estimator of $\rho$, which shows little bias. The
related Network Effects Model is briefly discussed.

Keywords: Network disturbance, network autocorrelation, maximum likelihood.

MRC subject classification: 91D30; 91Cxx, 62P25, 62J05

##  1 Introduction

The classical Gaussian linear regression model can be written as

$y\ =\ X\beta+\varepsilon.$ (1.1)

Here, the $n$-dimensional vector $y$ of data is described as the sum of a
structural element $X\beta$, where $X$ is a known $n\times m$ matrix of full
rank $m$, $\beta=(\beta_{1},\ldots,\beta_{m})^{\top}$ is an $m$-vector of
unknown parameters, to be estimated, and the elements
$\varepsilon_{1},\ldots,\varepsilon_{n}$ of the $n$-vector $\varepsilon$ are
realizations of independent normally distributed random variables
$\boldsymbol{\varepsilon}_{1},\ldots,\boldsymbol{\varepsilon}_{n}$, each with
zero mean and variance $\sigma^{2}$. Here and throughout, we conventionally
use bold face to denote random elements, and italic for their realizations.

If the indices $i=1,2,\ldots,n$ represent the vertices of a simple undirected
graph, whose (symmetric) adjacency matrix is denoted by $A$, it may be
reasonable to suppose that the errors $\boldsymbol{\varepsilon}_{i}$ and
$\boldsymbol{\varepsilon}_{j}$ at neighbouring vertices $i$ and $j$ are
correlated, instead of being independent. A similar idea lies at the basis of
time series analysis, where the error at a given time point is assumed to
depend on the errors at the closest previous time points. In analogy with the
moving average models of time series analysis, the Network Disturbance Model
of [3] (see also [12, (4.6)]) supposes that $\mathbf{y}\ =\
X\beta+\boldsymbol{\varepsilon}$, where the errors $\boldsymbol{\varepsilon}$
are now implicitly given by the equation

$\boldsymbol{\varepsilon}_{i}\ =\
\rho\sum_{{j=1}}^{n}W_{ij}\boldsymbol{\varepsilon}_{j}+\boldsymbol{\nu}_{i},\qquad
1\leq i\leq n,$ (1.2)

or more succinctly, by

$K_{\rho}\boldsymbol{\varepsilon}\ =\ \boldsymbol{\nu},\qquad\mbox{where}\quad
K_{\rho}\ :=\ I-\rho W.$ (1.3)

Here, the weight matrix $W$ is known, the parameter $\rho$ is unknown, and the
random variables $\boldsymbol{\nu}_{1},\ldots,\boldsymbol{\nu}_{n}$ are
independent normal random variables with mean zero and variance $\sigma^{2}$.
The weight $W_{ij}\geq 0$ can be thought of as a (relative) measure of the
influence of vertex $j$ on vertex $i$, and is assumed to take the value zero
if $i$ and $j$ are not neighbours, i.e. if $A_{ij}=0$; we also assume that
$W_{ii}=0$ for each $i$. The parameter $\rho$ can be varied, and $|\rho|$
mirrors the overall strength of the dependence between the observations at
neighbouring vertices.

The network disturbance model can be estimated by using the method of maximum
likelihood (MLE); details of the procedure are given in Section 2. In
classical settings, including repeated independent observations from a common
underlying distribution, the maximum likelihood estimator is asymptotically
unbiased, as the number of observations tends to infinity. Here, the
observations are dependent, and the standard asymptotic conclusions need no
longer be valid. Indeed, it has been widely observed in simulation studies
that, as the underlying network becomes denser, and the corresponding weight
matrix $W$ becomes less sparse, the maximum likelihood estimator
${\hat{\rho}}$ of $\rho$ becomes more variable and more negatively biased.
This has been extensively documented, largely in the context of the related
network effects models of [3], defined in (1.4) below, for instance in [10],
[16], [4], [11] and [8]. Methods for correcting the bias of the MLE, using
asymptotics appropriate in settings where consistent estimation is possible,
have since been proposed in [1], [19] and [18]. Other methods of estimation
have been developed by [7], who introduce a generalized method of moments
estimator, by [2], with a Bayesian approach, and by [17], whose techniques are
based on Hodges–Lehmann estimators. In [16], theoretical investigation shows
that, at least for networks very close to the complete graph, the MLE may lead
to misleading results: see the final example in Section 4.4 for a fuller
explanation, that makes clear why no estimator of $\rho$ can work well in such
circumstances.

In this paper, we examine the problem in greater depth, showing why the MLE
can be expected to be biased and lacking in precision, whenever the underlying
network is dense. In Section 3, we use Cramér–Rao theory to show that, if the
underlying network is sufficiently dense, then any (more or less unbiased)
estimator of $\rho$ must necessarily be unreliable. The bound is expressed as
a lower bound (3.7) on the variance of any (unbiased) estimator, and is given
explicitly as a function of the matrix $W$ of weights and of the true value of
$\rho$. The structure of the network enters by way of $W$. For instance, the
choice $W_{ij}:=A_{ij}/\deg(i)$ makes explicit the dependence of the weights
on the adjacency matrix $A$ of the network. We then, in Section 4, examine the
behaviour of the maximum likelihood estimator in enough detail to explain the
empirical observations. In particular, we are able to show that the MLE of
$\rho$ is biased when the network is sufficiently dense. In ordinary least
squares estimation, the estimator of $\beta$ has small variance if the
smallest eigenvalue of the matrix $X^{\top}X$ is large. Here, it is possible,
for rather special graphs, that the MLE of $\beta$ may still have substantial
variance, even when this condition is satisfied; see Examples 1 and 2 in
Section 4.4.

The exact formulae for the bias and variance of the MLE are in general
indigestible, but we discuss some special cases in which it is possible to
estimate their magnitudes. In Section 5, we consider a simpler and more
intuitive estimator of $\rho$, that has better properties as regards bias, but
is still lacking in precision in dense graphs, as it must be, in view of the
lower bound given in Section 3. Its properties are illustrated by some
simulations. An R function for computing the estimator, as well an empirical
measure of its precision, is available on request from the authors.

There is a related model, the Network Effects Model of [3] (see also [12,
(4.1)], which takes the form

$K_{\rho}y\ =\ X\beta+\nu,$ (1.4)

where the components of $\nu$ are realizations of independent normal random
variables with mean zero and variance $\sigma^{2}$, and $K_{\rho}$ is as
defined in (1.3). This differs from the network disturbance model because,
rewriting it as

$y\ =\ K_{\rho}^{-1}X\beta+\varepsilon,$ (1.5)

the errors in $y$ have the same dependence as before, but the model structure
explaining $y$ also changes with $\rho$. The limitations as regards estimating
$\rho$ in dense networks are analogous, and for similar reasons. We discuss
the model briefly in Section 6.

The main conclusions of the paper are sketched in Section 7, and additional
calculations can be found in the Appendix.

##  2 Details of the model

The maximum likelihood estimates of the parameters in the simple linear
regression model $\mathbf{y}\ =\ X\beta+\boldsymbol{\varepsilon}$ of (1.1),
with independent and identically distributed errors, are given by

$\hat{\beta}\ =\ \hat{\beta}(y)\ =\
(X^{\top}X)^{-1}X^{\top}y;\qquad\hat{\sigma}^{2}\ =\ \hat{\sigma}^{2}(y)\ =\
n^{-1}y^{\top}(I-H(X))y,$ (2.1)

where

$H(X)\ :=\ X(X^{\top}X)^{-1}X^{\top}$ (2.2)

is symmetric and idempotent, and represents the orthogonal projection from
$\mathbb{R}^{n}$ onto the linear span of the columns of $X$. Writing

$M_{0}\colon\ \mathbf{y}\ =\ X\beta+\boldsymbol{\varepsilon}$ (2.3)

to denote the probabilistic model from which the data are sampled, the random
quantities $\hat{\beta}(\mathbf{y})$ and $\hat{\sigma}^{2}(\mathbf{y})$ have
distributions given by

$\hat{\beta}(\mathbf{y})\ \sim\
{\mathrm{MVN}\,}_{m}(\beta,\sigma^{2}(X^{\top}X)^{-1});\qquad\hat{\sigma}^{2}(\mathbf{y})\
\sim\ n^{-1}\sigma^{2}\chi^{2}_{n-m},$ (2.4)

where ${\mathrm{MVN}\,}_{m}(\mu,\Sigma)$ denotes the $m$-dimensional
multivariate normal distribution with mean $\mu$ and covariance marix
$\Sigma$, and $\chi^{2}_{l}$ denotes the chi-squared distribution with $l$
degrees of freedom. The two estimators are also independent of one another
(see [13, Section 4.b]).

In the network disturbance model $\mathbf{y}\ =\
X\beta+\boldsymbol{\varepsilon}$, the errors $\boldsymbol{\varepsilon}$ are
given by (1.3), with dependence deriving from the weight matrix $W$. It is
often convenient to assume that the row sums of $W$ are all equal to $1$, in
which case the quantity $\sum_{j=1}^{n}W_{ij}\boldsymbol{\varepsilon}_{j}$
represents an average of the errors at the neighbours of vertex $i$; then
$\rho$ would typically be smaller than $1$ in modulus, so as to ensure that
the average influence of the neighbours is not larger than the size of a
typical error. A simple choice for $W_{ij}$ is $A_{ij}/\deg(i)$, where
$\deg(i)$ denotes the degree of the vertex $i$. If there are isolated vertices
$i$, the corresponding row sums cannot be made equal to $1$, since then
$W_{ij}=0$ for all $j\neq i$.

The implicit representation (1.3) of the errors can be written as a direct
formula for $\boldsymbol{\varepsilon}$ in terms of $\boldsymbol{\nu}$, as

$\boldsymbol{\varepsilon}\ =\ K_{\rho}^{-1}\boldsymbol{\nu},$ (2.5)

provided that the matrix $K_{\rho}$ can be inverted. The expansion

$K_{\rho}^{-1}\ =\ (I-\rho W)^{-1}\ =\ I+\sum_{l\geq 1}\rho^{l}W^{l}$ (2.6)

makes sense provided that the sum $\sum_{l\geq 1}\rho^{l}W^{l}x$ converges
absolutely for each $x\in\mathbb{R}^{n}$; that is, if ${|\rho|}$ is less than
$1/r(W)$, where $r(W)$ is the spectral radius of $W$. In particular, if $W$
has all of its row sums equal to $1$, the vector ${\bf{1}}$ with all
components equal to $1$ is a right eigenvector with eigenvalue $1$, implying
that $r(W)=1$ (see [6, Corollary 8.1.30, p.493]). Then

$\boldsymbol{\varepsilon}\ =\ K_{\rho}^{-1}\boldsymbol{\nu}\ =\
\boldsymbol{\nu}+\sum_{l\geq 1}\rho^{l}W^{l}\boldsymbol{\nu},\qquad|\rho|<1,$

represents $\boldsymbol{\varepsilon}_{i}$ as a weighted sum of the
contributions of the different elements $\boldsymbol{\nu}_{j}$, averaged along
paths in the network from $j$ to $i$, and with overall weight $\rho^{l}$ for
paths of length $l$ (possibly negative, if $\rho<0$). If $|\rho|>1/r(W)$, this
interpretation no longer holds, since the expansion (2.6) is no longer well
defined. The resulting model is then not so easy to understand, since, at
least locally, the cumulative effect on the observation at $i$ from those at
distance $d$ from $i$ increases geometrically with $d$. In time series
analysis, the analogue is that of processes that are ‘explosive’, and the
analogue of the condition $|\rho|<1/r(W)$ is used to ensure that a time series
is stationary. For this reason, we shall assume in what follows that
$|\rho|<1/r(W)$, so that the inverse of $K_{\rho}$ exists, and is explicitly
given by (2.6).

The model $\mathbf{y}\ =\ X\beta+\boldsymbol{\varepsilon}$ with errors given
by (1.2) is simple to analyze if the value of $\rho$ is known, since it can be
written in the form

$M_{\rho}\colon\quad K_{\rho}\mathbf{y}\ =\ K_{\rho}X\beta+\boldsymbol{\nu},$
(2.7)

with the components of $\boldsymbol{\nu}$ being independent normal
${\mathcal{N}}(0,\sigma^{2})$ random variables, and can thus be fitted as a
standard linear regression based on data $y^{\prime}:=K_{\rho}y$. Equations
(2.1) then give the estimates

$\displaystyle\hat{\beta}_{\rho}$ $\displaystyle=$
$\displaystyle\hat{\beta}_{\rho}(y)\ =\
((K_{\rho}X)^{\top}K_{\rho}X)^{-1}(K_{\rho}X)^{\top}K_{\rho}y;$ (2.8)
$\displaystyle{\hat{\sigma}}^{2}_{\rho}$ $\displaystyle=$
$\displaystyle{\hat{\sigma}}^{2}_{\rho}(y)\ =\
{n^{-1}\\{K_{\rho}(y-X\hat{\beta}_{\rho})\\}^{\top}\\{K_{\rho}(y-X\hat{\beta}_{\rho})\\}}\
=\ n^{-1}(K_{\rho}y)^{\top}(I-H_{\rho})K_{\rho}y,$

where, with $H$ as in (2.2),

$H_{\rho}\ :=\ H(K_{\rho}X)$ (2.9)

and so $K_{\rho}X\hat{\beta}_{\rho}=H_{\rho}K_{\rho}y$. The residuals after
fitting $M_{\rho}$ are thus

${\hat{\varepsilon}}_{\rho}(y)\ :=\ y-X{\hat{\beta}}_{\rho}(y)\ =\
K_{\rho}^{-1}(I-H_{\rho})K_{\rho}y.$ (2.10)

The distributions of $\hat{\beta}_{\rho}(\mathbf{y})$ and
$\hat{\sigma}^{2}_{\rho}(\mathbf{y})$ under the model $M_{\rho}$ are given by

${\hat{\beta}}_{\rho}(\mathbf{y})\ \sim\
{\mathrm{MVN}\,}_{m}(\beta,\sigma^{2}((K_{\rho}X)^{\top}K_{\rho}X)^{-1});\quad{\hat{\sigma}}^{2}_{\rho}(\mathbf{y})\
\sim\ n^{-1}\sigma^{2}\chi^{2}_{n-m},$ (2.11)

and the two estimators are independent of one another. In particular, both
${\hat{\beta}}$ and ${\hat{\sigma}}^{2}$ have precision much as in simple
linear regression, and (under reasonable assumptions about $X$) can be
expected to become more accurate as $n$ increases. If the true model is
$M_{\rho^{\prime}}$, for $\rho^{\prime}\neq\rho$, the expectation of
${\hat{\beta}}_{\rho}(\mathbf{y})$ is still equal to $\beta$, and, although
the variances of its components may be larger, they are typically comparable
to those on the true model, again under reasonable conditions; see Appendix
A.4 for more detail.

If $\rho$ is not known, the estimates (2.8) still maximize the log-likelihood

$\ell(\theta;y)\ :=\
-{\textstyle{\frac{n}{2}}}\\{\log(\sigma^{2})+\log(2\pi)\\}-\frac{1}{2\sigma^{2}}\\{K_{\rho}(y-X\beta)\\}^{\top}\\{K_{\rho}(y-X\beta)\\}+\log\det
K_{\rho}$ (2.12)

of the model $M_{\rho}$ of (2.7), for each given value of $\rho$, where
$\theta:=(\beta,\sigma^{2},\rho)$ denotes the full set of parameters; see [3,
Equation (1.3)]. The resulting maximum of the log-likelihood, for fixed
$\rho$, is given by

$-{\textstyle{\frac{n}{2}}}\\{\log({\hat{\sigma}}^{2}_{\rho}(y))+\log(2\pi)+1\\}+\log\det
K_{\rho};$ (2.13)

this can be deduced by observing that, with $\beta$ replaced by
$\hat{\beta}_{\rho}$, the middle term of (2.12) reduces to $-n/2$, in view of
(2.8). Indeed, the first term in (2.13) is the maximum of the likelihood for
ordinary least squares applied to the data $y^{\prime}:=K_{\rho}y$, and the
determinant enters because the probability densities $f_{\mathbf{y}^{\prime}}$
of $\mathbf{y}^{\prime}$ and $f_{\mathbf{y}}$ of $\mathbf{y}$ are related by
$f_{\mathbf{y}}(v)=\det K_{\rho}\,f_{\mathbf{y}^{\prime}}(K_{\rho}v)$. Then
$\rho$ can be chosen to (numerically) maximize the expression in (2.13), or,
equivalently, to minimize

$F(\rho;y)\ :=\
\log(\hat{\sigma}^{2}_{\rho}(y))-{\textstyle{\frac{2}{n}}}\log\det
K_{\rho},\qquad|\rho|<1/r(W).$ (2.14)

##  3 A general limit to precision

The $\rho$-component of the score function is obtained by differentiating the
log-likelihood (2.12) with respect to $\rho$. Noting that
$\frac{d}{d\rho}K_{\rho}=-W$, and that, if the eigenvalues of $W$ are denoted
by $\lambda_{1},\ldots,\lambda_{n}$, we have $\det
K_{\rho}=\prod_{i=1}^{n}(1-\rho\lambda_{i})$. This gives the $\rho$-component
of the score function as

$\displaystyle s^{(\rho)}(\theta;y)$ $\displaystyle:=$
$\displaystyle\frac{\partial}{\partial\rho}\ell(\theta;y)\ =\
\frac{1}{\sigma^{2}}(y-X\beta)^{\top}K_{\rho}^{\top}W(y-X\beta)-\sum_{l=1}^{n}\frac{\lambda_{l}}{1-\rho\lambda_{l}}$
(3.1) $\displaystyle=$
$\displaystyle\frac{1}{\sigma^{2}}\bigl{(}\\{K_{\rho}(y-X\beta)\\}^{\top}Z_{\rho}\\{K_{\rho}(y-X\beta)\\}\bigr{)}-{\mathrm{Tr}}\\{Z_{\rho}\\},$

where we define the matrix

$Z_{\rho}\ :=\ WK_{\rho}^{-1}.$ (3.2)

Now, if $C$ is any $n\times n$ matrix, and ${\bf\tilde{\boldsymbol{\nu}}}$
denotes an $n$-vector of independent standard normal random variables, then
the quadratic form
${\bf\tilde{\boldsymbol{\nu}}}^{\top}C{\bf\tilde{\boldsymbol{\nu}}}$ has mean

${\mathbb{E}}\\{{\bf\tilde{\boldsymbol{\nu}}}^{\top}C{\bf\tilde{\boldsymbol{\nu}}}\\}\
=\ {\mathrm{Tr}}\\{C\\}.$ (3.3)

Furthermore, for symmetric $C$, since
${\mbox{Cov}}\\{{\bf\tilde{\boldsymbol{\nu}}}_{i}{\bf\tilde{\boldsymbol{\nu}}}_{j},{\bf\tilde{\boldsymbol{\nu}}}_{k}{\bf\tilde{\boldsymbol{\nu}}}_{l}\\}=0$
unless either $i\neq j$ and $\\{i,j\\}=\\{k,l\\}$ or $i=j=k=l$, the variance
of ${\bf\tilde{\boldsymbol{\nu}}}^{\top}C{\bf\tilde{\boldsymbol{\nu}}}$ is
given by

$\displaystyle{\rm
Var}\\{{\bf\tilde{\boldsymbol{\nu}}}^{\top}C{\bf\tilde{\boldsymbol{\nu}}}\\}$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{n}\sum_{j\neq
i}C_{ij}(C_{ij}+C_{ji})+\sum_{i=1}^{n}C_{ii}^{2}({\mathbb{E}}\\{{\bf\tilde{\boldsymbol{\nu}}}_{1}^{4}\\}-1)$
(3.4) $\displaystyle=$ $\displaystyle 2{\mathrm{Tr}}(C^{2}),$ (3.5)

because ${\mathbb{E}}\\{{\bf\tilde{\boldsymbol{\nu}}}_{1}^{4}\\}=3$ when
${\bf\tilde{\boldsymbol{\nu}}}_{1}\sim{\mathcal{N}}(0,1)$; see [14, Chapter 2,
Theorem 4 (i) & (ii)]. If $C$ is not symmetric, the expression (3.5) can be
used with $C$ replaced by the symmetric matrix
$C^{\prime}:={\textstyle{\frac{1}{2}}}(C+C^{\top})$, since
${\bf\tilde{\boldsymbol{\nu}}}^{\top}C^{\prime}{\bf\tilde{\boldsymbol{\nu}}}={\bf\tilde{\boldsymbol{\nu}}}^{\top}C{\bf\tilde{\boldsymbol{\nu}}}$.
Hence, taking expectations on the model with parameters
$\theta_{0}:=(\beta_{0},\sigma^{2}_{0},\rho_{0})$, we have

${\mathbb{E}}_{\theta_{0}}\,s^{(\rho)}(\theta_{0};\mathbf{y})\ =\
{\mathrm{Tr}}\\{Z_{\rho_{0}}\\}-{\mathrm{Tr}}\\{Z_{\rho_{0}}\\}\ \ =\ 0,$

as it must, and

${\rm Var}_{\theta_{0}}\,s^{(\rho)}(\theta_{0};\mathbf{y})\ =\
2{\mathrm{Tr}}\bigl{\\{}\bigl{(}{\textstyle{\frac{1}{2}}}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})\bigr{)}^{2}\bigr{\\}}\
=\ {\mathrm{Tr}}\\{Z_{\rho_{0}}^{2}+Z_{\rho_{0}}Z_{\rho_{0}}^{\top}\\}.$ (3.6)

This, together with [5, Theorem 6.2] and [9, Theorem 6.6 and (6.25)], implies
the following limit on the precision of estimators of $\rho$.

Cramér–Rao lower bound. The variance of any unbiased estimator of $\rho$ is no
smaller than

$\gamma_{\rho_{0}}^{2}\ :=\
\frac{1}{{\mathrm{Tr}}\\{Z_{\rho_{0}}^{2}+Z_{\rho_{0}}Z_{\rho_{0}}^{\top}\\}}\,,$
(3.7)

where $Z_{\rho}:=WK_{\rho}^{-1}$ and $K_{\rho}$ is as in (2.5).

In practice, the bound implies a precision that cannot be significantly
improved upon for any estimator of $\rho$, without introducing substantial
bias; even the best ‘sensible’ estimator of $\rho$ has to have random
variability on the scale of $\gamma_{\rho_{0}}$.

Note that, for any $n\times n$ matrix $M$,

$|{\mathrm{Tr}}\\{M^{2}\\}|\ =\
\Bigl{|}\sum_{i=1}^{n}\sum_{j=1}^{n}M_{ij}M_{ji}\Bigr{|}\ \leq\
\sum_{i=1}^{n}\sum_{j=1}^{n}M_{ij}^{2}\ =\ {\mathrm{Tr}}\\{MM^{\top}\\},$
(3.8)

by Cauchy–Schwarz, so that the two traces appearing in $\gamma_{\rho_{0}}^{2}$
are not in general the same, and that

${\mathrm{Tr}}\\{M^{2}+MM^{\top}\\}\ =\
\sum_{i=1}^{n}\sum_{j=1}^{n}\\{M_{ij}M_{ji}+M_{ij}^{2}\\}\ =\
{\textstyle{\frac{1}{2}}}\sum_{i=1}^{n}\sum_{j=1}^{n}(M_{ij}+M_{ji})^{2}\
\geq\ 0.$ (3.9)

#### Consequences of (3.7): the case $\rho_{0}=0$.

As a simple example of what is to be expected, suppose that $\rho_{0}=0$, and
that the underlying graph has no isolated vertices, so that $\deg(i)\geq 1$
for all $i$. Let $W:=D^{-1}A$, where $D$ is the diagonal matrix with elements
$D_{ii}:=\deg(i)$, $1\leq i\leq n$. Then $Z_{\rho_{0}}=Z_{0}=W$, and we can be
specific about the traces in (3.7):

$\displaystyle{\mathrm{Tr}}\\{W^{2}\\}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{A_{ij}}{\deg(i)\deg(j)}$
(3.10) $\displaystyle\leq$ $\displaystyle{\mathrm{Tr}}\\{WW^{\top}\\}\ =\
\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{A_{ij}}{\deg(i)^{2}}\ =\
\sum_{i=1}^{n}\frac{1}{\deg(i)}\ \leq\ n.$

It follows from the equality in (3.10) that ${\mathrm{Tr}}\\{WW^{\top}\\}\geq
n/{\overline{d}}$, where ${\overline{d}}:=n^{-1}\sum_{i=1}^{n}\deg(i)$ denotes
the average degree of a vertex, and hence that
$\gamma_{0}^{2}\leq{\overline{d}}/n$. However, if the vertex degrees are
reasonably homogeneous, in the sense all the vertices have degrees at least
$c{\overline{d}}$, for some $0<c<1$, then it follows from (3.10) that
${\mathrm{Tr}}\\{WW^{\top}\\}\leq n/c{\overline{d}}$, and thus, from (3.7),
that $\gamma_{0}^{2}\geq c{\overline{d}}/n$. Hence, in these simple
circumstances, the variability of any sensible estimator of $\rho$ is on a
scale of at least $({\overline{d}}/n)^{1/2}$.

These calculations have particular implications for dense networks. If
$W=D^{-1}A$, and if all vertices have degrees exceeding $cn$, for some fixed
$c>0$, it follows that $\gamma_{0}\geq\sqrt{c}$, and it is not possible to
obtain an accurate estimate of $\rho$ by any method, no matter how large the
size $n$ of the network, if in fact $\rho_{0}=0$. For more general choices of
$W$, ${\mathrm{Tr}}\\{WW^{\top}\\}=\sum_{i,j=1}^{n}W_{ij}^{2}$, and if all the
elements of $W$ are smaller than $C/n$, for some $C<\infty$, it follows that
${\mathrm{Tr}}\\{WW^{\top}\\}\leq C^{2}$, and hence, from (3.7), that
$\gamma_{0}^{2}>1/C^{2}$. Hence, here too, any estimator of $\rho$ has
significant variability on the scale of $1/C$, whatever the size $n$ of the
network, if in fact $\rho_{0}=0$.

#### Consequences of (3.7): general $\rho_{0}$.

If $\rho_{0}\neq 0$, analogous results are true. Supposing that
$\sum_{j=1}^{n}W_{ij}=1$ for each $i$, it follows that, for all $l\geq 1$,
$\sum_{j=1}^{n}(W^{l})_{ij}=1$ for all $i$ also. Hence
$(W^{l})_{ij}=\sum_{k=1}^{n}(W^{l-1})_{ik}W_{kj}$ is an average of the values
$W_{1j},\ldots,W_{nj}$, Thus, if all the elements of $W$ are smaller than
$C/n$, it follows that $(W^{l})_{ij}\leq C/n$ also. Hence it is immediate
that, for $l,l^{\prime}\geq 1$,

${\mathrm{Tr}}\\{W^{l}\\}\ \leq\
C\quad\mbox{and}\quad{\mathrm{Tr}}\\{W^{l}(W^{\top})^{l^{\prime}}\\}\ \leq\
C.$ (3.11)

Thus, for $\rho\geq 0$, we have

$\displaystyle{\mathrm{Tr}}\\{Z_{\rho}^{2}\\}$ $\displaystyle\leq$
$\displaystyle{\mathrm{Tr}}\\{Z_{\rho}Z_{\rho}^{\top}\\}\ =\
{\mathrm{Tr}}\\{W(I-\rho W)^{-1}(I-\rho W^{\top})^{-1}W^{\top}\\}$ (3.12)
$\displaystyle\leq$ $\displaystyle C\sum_{l\geq 0}(l+1)\rho^{l}\ =\
\frac{C}{(1-\rho)^{2}}\,,$

and so

$\gamma_{\rho_{0}}\ \geq\ C^{-1/2}(1-\rho).$

Once again, there is substantial variability in the estimate of $\rho$ unless
$C:=n\max_{i,j}W_{ij}$ is large. If $\rho<0$, the bounds in (3.11) can be used
to justify similar conclusions, by bounding the difference
${\mathrm{Tr}}\\{Z_{\rho}Z_{\rho}^{\top}\\}-{\mathrm{Tr}}\\{WW^{\top}\\}$,
provided that $|\rho|$ is not too large.

In particular, in uniformly dense networks, if $W=D^{-1}A$, there is
substantial variability in any estimator of $\rho$ if all vertex degrees
exceed $n/C$, and $C$ is not large.

##  4 Maximum likelihood estimation

The discussion in the previous section shows that there is an upper limit on
the precision of any (unbiased) estimator of $\rho$, and that, in uniformly
dense networks, any such estimator must have substantial variability. In this
section, we go into more detail about the properties of the MLE.

###  4.1 General theory

An equation for the MLE

To minimize the quantity $F(\rho;y)$ of (2.14), it is usual to look for values
of $\rho$ such that

$U(\rho;y)\ :=\ \frac{d}{d\rho}F(\rho;y)\ =\ 0.$ (4.1)

Since $F(\rho;y)$ is explicitly given as a function of $\rho$, its derivative
can be calculated, giving

$\displaystyle U(\rho;y)$ $\displaystyle=$
$\displaystyle\frac{1}{{\hat{\sigma}}^{2}_{\rho}(y)}\,\frac{d{\hat{\sigma}}^{2}_{\rho}(y)}{d\rho}-\frac{2}{n}\,\frac{d}{d\rho}\,\log\det
K_{\rho}$ (4.2) $\displaystyle=$
$\displaystyle-\frac{1}{n{\hat{\sigma}}^{2}_{\rho}(y)}\,(K_{\rho}y)^{\top}(I-H_{\rho})(Z_{\rho}^{\top}+Z_{\rho})(I-H_{\rho})K_{\rho}y+\frac{2}{n}\sum_{l\geq
1}\rho^{l}{\mathrm{Tr}}\\{W^{l+1}\\};\phantom{XX}$

the details are given in Appendix A.1. The maximum likelihood estimate
${\hat{\rho}}_{\mathrm{ML}}$ of $\rho$ is now typically obtained by solving
$U(\rho;y)=0$; Example 3 in Section 4.4 shows that this may not always work.

Note that the data enter through a ratio of quadratic forms in $y$. Under
$M_{\rho_{0}}$, as observed in (2.11),
${\hat{\sigma}}^{2}_{\rho_{0}}(\mathbf{y})\ \sim\
n^{-1}\sigma^{2}\chi^{2}_{n-m}$, having mean $\sigma^{2}(1-m/n)$ and variance
less than $2n^{-1}\sigma^{4}$, so that, to a reasonable approximation,
${\hat{\sigma}}^{2}_{\rho_{0}}(y)$ can be replaced by $\sigma^{2}$ in (4.2) if
$M_{\rho_{0}}$ is true, provided that $n$ is large. As a result, $U(\rho;y)$
depends on the data $y$ largely through a single quadratic form in $y$.

Approximating the distribution of $U(\rho;\mathbf{y})$ under $M_{\rho_{0}}$

In order to gauge the typical values of $U(\rho;y)$ to be obtained under the
true model $M_{\rho_{0}}$, we express $U(\rho;\mathbf{y})$ in terms of the
(not directly observable) errors $\boldsymbol{\nu}$, rather than in terms of
$\mathbf{y}$. Under $M_{\rho_{0}}$,
$\mathbf{y}=X\beta+K_{\rho_{0}}^{-1}\boldsymbol{\nu}$, and it thus follows
from (2.14) and (4.2) that

$\displaystyle U(\rho;\mathbf{y})$ $\displaystyle=$
$\displaystyle\frac{1}{{\hat{\sigma}}^{2}_{\rho}(\mathbf{y})}\,\frac{d{\hat{\sigma}}^{2}_{\rho}(\mathbf{y})}{d\rho}-\frac{2}{n}\,\frac{d}{d\rho}\,\log\det
K_{\rho}$ (4.3) $\displaystyle=$
$\displaystyle-\frac{1}{n{\hat{\sigma}}^{2}_{\rho}(K_{\rho_{0}}^{-1}\boldsymbol{\nu})}\,\boldsymbol{\nu}^{\top}Q(\rho;\rho_{0})\boldsymbol{\nu}+\frac{2}{n}\sum_{l\geq
1}\rho^{l}{\mathrm{Tr}}\\{W^{l+1}\\},$

where

$Q(\rho;\rho_{0})\ :=\
(K_{\rho}K_{\rho_{0}}^{-1})^{\top}(I-H_{\rho})(Z_{\rho}^{\top}+Z_{\rho})(I-H_{\rho})K_{\rho}K_{\rho_{0}}^{-1}.$
(4.4)

Note that $\beta$ does not appear in (4.3) because
$(I-H_{\rho})^{T}=I-H_{\rho}$ and $(I-H_{\rho})K_{\rho}X=0$, implying in
particular, using (2.8), that
${\hat{\sigma}}^{2}_{\rho}(\mathbf{y}-X\beta)={\hat{\sigma}}^{2}_{\rho}(\mathbf{y})$.
Hence, and because
${\hat{\sigma}}^{2}_{\rho}(K_{\rho_{0}}^{-1}\boldsymbol{\nu})$ is close to
$\sigma^{2}$ for large $n$, the random part of $U(\rho;\mathbf{y})$ is given,
to a close approximation, by
$-\boldsymbol{\nu}^{\top}Q(\rho;\rho_{0})\boldsymbol{\nu}/(n\sigma^{2})$.

In the remaining argument, we make frequent use of elementary properties of
the trace, that ${\mathrm{Tr}}(AB)={\mathrm{Tr}}(BA)$ and that
${\mathrm{Tr}}(A^{\top})={\mathrm{Tr}}(A)$, as well as noting that
$z^{\top}Mz=z^{\top}M^{\top}z$ for any compatible $z$ and $M$, and that
$H_{\rho}$ and $I-H_{\rho}$ are idempotent. The expectation of
$\boldsymbol{\nu}^{\top}Q(\rho_{0},\rho_{0})\boldsymbol{\nu}$, using (3.3), is
given by

$\displaystyle{\mathbb{E}}\\{\boldsymbol{\nu}^{\top}Q(\rho_{0},\rho_{0})\boldsymbol{\nu}\\}$
$\displaystyle=$
$\displaystyle\sigma^{2}{\mathrm{Tr}}\\{Q(\rho_{0},\rho_{0})\\}\ =\
{2\sigma^{2}{\mathrm{Tr}}\\{(I-H_{\rho_{0}})Z_{\rho_{0}}(I-H_{\rho_{0}})\\}}$
(4.5) $\displaystyle=$ $\displaystyle
2\sigma^{2}{\mathrm{Tr}}\\{(I-H_{\rho_{0}})Z_{\rho_{0}}\\}\ =\
2\sigma^{2}\bigl{(}{\mathrm{Tr}}\\{Z_{\rho_{0}}\\}-{\mathrm{Tr}}\\{H_{\rho_{0}}Z_{\rho_{0}}\\}\bigr{)},$

because $K_{\rho_{0}}K_{\rho_{0}}^{-1}=I$, and using
${\mathrm{Tr}}(AB)={\mathrm{Tr}}(BA)$ and
$(I-H_{\rho_{0}})^{2}=I-H_{\rho_{0}}$. Then, using ${\mathrm{Tr}}\\{W\\}=0$,

${\mathrm{Tr}}\\{Z_{\rho_{0}}\\}\ =\ {\mathrm{Tr}}\\{W(I-\rho_{0}W)^{-1}\\}\
=\ \sum_{l\geq 1}\rho_{0}^{l}{\mathrm{Tr}}\\{W^{l+1}\\}.$ (4.6)

Thus we have

$-{\mathbb{E}}\bigl{\\{}\sigma^{-2}\boldsymbol{\nu}^{\top}Q(\rho_{0};\rho_{0})\boldsymbol{\nu}\bigr{\\}}+2\sum_{l\geq
1}\rho^{l}{\mathrm{Tr}}\\{W^{l+1}\\}\ =\
2{\mathrm{Tr}}\\{H_{\rho_{0}}Z_{\rho_{0}}\\}.$

Hence, from (4.3) and (4.5), neglecting the small variability in
${\hat{\sigma}}^{2}_{\rho}(\mathbf{y})$, the quantity
$nU(\rho_{0};\mathbf{y})$ under the model $M_{\rho_{0}}$ satisfies

${\mathbb{E}}_{\rho_{0}}\\{nU(\rho_{0};\mathbf{y})\\}\ \approx\ \mu_{0}\ :=\
\mu(\rho_{0})\ :=\ 2{\mathrm{Tr}}\\{H_{\rho_{0}}Z_{\rho_{0}}\\}.$ (4.7)

The principal source of variability in $nU(\rho_{0};\mathbf{y})$ arises from
$\sigma^{-2}\boldsymbol{\nu}^{\top}Q(\rho_{0};\rho_{0})\boldsymbol{\nu}$,
whose variance under $M_{\rho_{0}}$, from (3.5), is given by
$2{\mathrm{Tr}}\\{Q(\rho_{0};\rho_{0})^{2}\\}$. This yields

$\displaystyle{\rm Var}_{\rho_{0}}(nU(\rho_{0};\mathbf{y}))$
$\displaystyle\approx$ $\displaystyle
2{\mathrm{Tr}}\\{Q(\rho_{0};\rho_{0})^{2}\\}$ (4.8) $\displaystyle=$
$\displaystyle{2{\mathrm{Tr}}\\{(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{\rho_{0}})(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{\rho_{0}})\\}}$
$\displaystyle=:$ $\displaystyle\tau_{0}^{2}\ :=\ \tau_{0}^{2}(\rho_{0}).$

Combining this with (4.7), we reach the following

Conclusion 1: the distribution of $nU(\rho_{0};\mathbf{y})$ can be written in
the form

$nU(\rho_{0};\mathbf{y})\ =\ \mu_{0}+\tau_{0}L(\mathbf{y}),$ (4.9)

where $L(\mathbf{y})$ has mean approximately $0$ and variance approximately
$1$.

Consequences for the distribution of the MLE

Using the information above, we can now make approximate deductions about the
MLE ${\hat{\rho}}:={\hat{\rho}}_{\mathrm{ML}}$, obtained by setting
$U(\rho;y)$ equal to zero. We begin by noting that $U({\hat{\rho}},y)=0$,
whereas, as in (4.9), $nU(\rho_{0};y)=\mu_{0}+\tau_{0}L(y)$, where
$L(\mathbf{y})$ approximately has mean zero and variance $1$. Writing
$U^{\prime}(\rho;y):=\frac{d}{d\rho}U(\rho;y)$, this implies that

$-n^{-1}(\mu_{0}+\tau_{0}L(y))\ =\ \\{U({\hat{\rho}};y)-U(\rho_{0};y)\\}\ =\
({\hat{\rho}}-\rho_{0})U^{\prime}(\rho_{0};y)+o(|{\hat{\rho}}-\rho_{0}|),$

where the final error term is relatively small if
$(|\mu_{0}|+\tau_{0})/n|U^{\prime}(\rho_{0};y)|$ is small enough, leading to
the approximate expression for the error ${\hat{\rho}}-\rho_{0}$ in estimating
$\rho_{0}$:

${\hat{\rho}}_{\mathrm{ML}}-\rho_{0}\ \approx\
-\frac{\mu_{0}+\tau_{0}L(y)}{nU^{\prime}(\rho_{0};y)}\,.$ (4.10)

This is essentially a one–step Newton–Raphson approximation to the value of
$\rho$ that solves $U(\rho;y)=0$, starting from the guess $\rho_{0}$. Thus it
is necessary also to examine the magnitude of $nU^{\prime}(\rho_{0};y)$. Since
$U(\rho;y)$ is explicitly given, there are exact formulae for its derivative,
given in Appendix A.2; however, they are not easy to understand. Instead, it
is possible to approximate the expectation of
$nU^{\prime}(\rho_{0};\mathbf{y})$ on the model $M_{\rho_{0}}$, giving

$n{\mathbb{E}}_{\rho_{0}}U^{\prime}(\rho_{0};\mathbf{y})\ \approx\
\Delta_{0},$ (4.11)

where

$\displaystyle\Delta_{0}$ $\displaystyle=$
$\displaystyle{\textstyle{\frac{1}{2}}}\tau_{0}^{2}+2{\mathrm{Tr}}\\{-Z_{\rho_{0}}H_{\rho_{0}}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{\rho_{0}})+H_{\rho_{0}}Z_{\rho_{0}}^{2}\\};$
(4.12)

the details are again to be found in Appendix A.2. Thus, from (4.10), we have
the following

Conclusion 2: the difference between ${\hat{\rho}}_{\mathrm{ML}}$ and
$\rho_{0}$ can be expected to be of magnitude roughly

$\frac{|\mu_{0}|{\pm}\tau_{0}}{\Delta_{0}}\,,$ (4.13)

with the term $|\mu_{0}|/\Delta_{0}$ corresponding to systematic bias, and
$\psi_{0}:=\tau_{0}/\Delta_{0}$ reflecting the scale of the randomness in the
estimator.

Note that the measure $\psi_{0}$ should be treated only as a rough guide to
the actual variability of ${\hat{\rho}}_{\mathrm{ML}}$ under $M_{\rho_{0}}$,
since (4.10) is only a linear approximation to the solution of a non-linear
equation, and because there may be added variability arising from the
variation of $nU^{\prime}(\rho_{0};\mathbf{y})$ about its mean.

###  4.2 Estimating $\rho$ when there is no structural element

To get an idea of the relative magnitudes of the terms $\mu_{0}$, $\tau_{0}$
and $\Delta_{0}$ in (4.13), it is helpful to look at some special cases. For
the first, suppose that there is no structural element at all — not even a
mean to be fitted — which corresponds to taking $X=0$; in this case, the
formulae (4.7), (4.8) and (4.12) should be interpreted as having $H_{\rho}=0$,
the projection onto the space $\\{0\\}$. Then $\mu_{0}=0$, so that any bias is
small,

$\tau_{0}^{2}\ =\
4{\mathrm{Tr}}\bigl{\\{}Z_{\rho_{0}}^{2}+Z_{\rho_{0}}Z_{\rho_{0}}^{\top}\bigr{\\}},\quad\mbox{where}\quad
Z_{\rho}\ =\ W(I-\rho W)^{-1},$

and

$\Delta_{0}\ =\ \tau_{0}^{2}/2,$

giving

$\psi_{0}\ =\ \frac{\tau_{0}}{\Delta_{0}}\ =\
\frac{1}{\sqrt{{\mathrm{Tr}}\\{Z_{\rho_{0}}^{2}+Z_{\rho_{0}}Z_{\rho_{0}}^{\top}\\}}}.$
(4.14)

This is the same as $\gamma_{\rho_{0}}$ in (3.7). Thus, from (4.13), the
typical magnitude of $|{\hat{\rho}}_{\mathrm{ML}}-\rho_{0}|$ is of order
$\gamma_{\rho_{0}}$. This is the best order of error that could theoretically
have been hoped for. Nonetheless, as observed above, it may not lead to
accurate estimation of $\rho$, for instance if $W=D^{-1}A$ and the underlying
network is sufficiently dense, even though the bias is small.

###  4.3 Estimating $\rho$ when structure is present

To accommodate the effect of having structure in the model, suppose first that
only an overall mean is to be fitted, so that $X$ is the $n$-vector
${\bf{1}}$, having all elements equal to $1$. Suppose also that there are no
isolated vertices, and that $W$ is chosen to have $\sum_{j=1}^{n}W_{ij}=1$ for
all $j$, so that $W{\bf{1}}={\bf{1}}$ and the spectral radius $r(W)=1$. Then
it is immediately verified that
$H_{\rho}=H_{*}:=n^{-1}{\bf{1}}{\bf{1}}^{\top}$ for all $\rho$, because, for
instance, $K_{\rho}X=(I-\rho W){\bf{1}}=(1-\rho){\bf{1}}$. Hence
$WH_{\rho}=H_{*}$ and $Z_{\rho}H_{\rho}=(1-\rho)^{-1}H_{*}$. As a result, by
arguments that are made precise in Appendix A.3, (4.8) and (4.12) yield

$\displaystyle\tau_{0}^{2}$ $\displaystyle=$ $\displaystyle
4{\mathrm{Tr}}\\{(Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2})(I-H_{*})\\}{+2\\{(1-\rho_{0})^{-2}-n^{-1}|Z_{\rho_{0}}^{\top}{\bf{1}}|^{2}\\}}$
(4.15) $\displaystyle\leq$
$\displaystyle{4{\mathrm{Tr}}\\{(Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2})(I-H_{*})\\}};$
$\displaystyle\Delta_{0}$ $\displaystyle=$
$\displaystyle{\textstyle{\frac{1}{2}}}\tau_{0}^{2}+2(1-\rho_{0})^{-2};$

here, and in what follows, $|x|$ for a vector $x$ denotes the Euclidean norm.
Hence, if $\rho_{0}$ is fixed and less than $1$, $\psi_{0}$ is of order
$1/(\tau_{0}+1)$, and this, from (4.13), represents the scale of the random
variation in the estimate of $\rho_{0}$. Now $1/(\tau_{0}+1)$ is small if
$\tau_{0}^{2}$ is large compared to $1$, but can be significantly larger than
the lower bound $\gamma_{\rho_{0}}$ in (3.7), if $\tau_{0}^{2}$ is much
smaller than
$1/\gamma_{\rho_{0}}^{2}={\mathrm{Tr}}\\{Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2}\\}$.
As in the case where there is no structure, if $W=D^{-1}A$ and the graph is
sufficiently dense, then $\gamma_{\rho_{0}}$ need not be small, in which case
the MLE has substantial variability.

In contrast to the situation without structure, the principal element
$\mu_{0}$ of the bias in the estimation of $\rho$ need not be zero. If
$X={\bf{1}}$, it follows from (4.7) that

$\mu_{0}\ =\ 2{\mathrm{Tr}}\\{Z_{\rho_{0}}H_{0}\\}\ =\ 2(1-\rho_{0})^{-1}\ >\
0.$

In view of Conclusion 2, the bias becomes important when $\Delta_{0}$ is not
large, which, from (4.15), is when $\tau_{0}^{2}$ is not large; as discussed
above, this occurs for example if $W=D^{-1}A$ and the graph is uniformly
dense. Note also that, from (4.10) and (4.13), the value of the difference
${\hat{\rho}}_{\mathrm{ML}}-\rho_{0}$ resulting from $\mu_{0}$ is roughly
$-2/\\{\Delta_{0}(1-\rho_{0})\\}$, which typically implies a negative bias, as
demonstrated in the results of the simulations referred to in the
introduction.

If there is more structure in the regression, but still assuming that the
first column of $X$ is ${\bf{1}}$ and that $W{\bf{1}}={\bf{1}}$, the detailed
formulae become more complicated, but the basic message remains the same; the
details are given in Appendix A.3. The matrix $H_{\rho}$ can be represented in
the form

$H_{\rho}\ =\ H_{*}+\sum_{l=2}^{m}x_{\rho}^{(l)}(x_{\rho}^{(l)})^{\top},$
(4.16)

where $H_{*}:=n^{-1}{\bf{1}}{\bf{1}}^{\top}$ is as before, and
$x_{\rho}^{(2)},\ldots,x_{\rho}^{(m)}$ are orthonormal vectors orthogonal to
${\bf{1}}$ that, with ${\bf{1}}$, span the space generated by the columns of
$K_{\rho}X$. The expression for $\tau_{0}^{2}$ derived in (A.3), in addition
to the quantity

$v(\rho_{0})\ :=\
{\mathrm{Tr}}\\{(Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2})(I-H_{\rho_{0}})\\},$
(4.17)

also involves the traces of the matrices $Z_{\rho_{0}}^{2}H_{\rho_{0}}$,
$Z_{\rho_{0}}^{\top}Z_{\rho_{0}}H_{\rho_{0}}$,
$(Z_{\rho_{0}}H_{\rho_{0}})^{2}$ and
$Z_{\rho_{0}}^{\top}H_{\rho_{0}}Z_{\rho_{0}}H_{\rho_{0}}$. Similarly, from
(4.12), the quantity $\Delta_{0}$ differs from
${\textstyle{\frac{1}{2}}}\tau_{0}^{2}$ by a quantity that is also expressed
in terms of these traces. It is shown in Appendix A.3 that if, for instance,
it is assumed that each of the vectors $x_{\rho}^{(l)}$, $2\leq l\leq m$, has
components bounded in modulus by $cn^{-1/2}$, for some fixed $c\geq 1$ — as is
the case for $x^{(1)}=n^{-1/2}{\bf{1}}$ — then each of the traces listed above
is bounded in modulus by $m^{2}c^{4}/(1-|\rho|)^{2}$. This gives rise to the
following

Conclusion 3: If, in the setting of this section, $m$ and $c$ are also fixed
for all $n$, and $|\rho_{0}|$ is bounded away from $1$, then the formulae for
$\tau_{0}^{2}$ and $\Delta_{0}$ given in (4.8) and (4.12) differ from
$4v(\rho_{0})$ and $2v(\rho_{0})$ by quantities that are uniformly bounded in
$n$ as $n\to\infty$, where $v(\rho_{0})$ is as defined in (4.17). The
variability in the MLE of $\rho$ is then on the scale
$1/\max\\{1,\sqrt{v(\rho_{0})}\\}$.

Again, this scale of variation may be much larger than $\gamma_{0}$ if
$v(\rho_{0})$ is much smaller than
${\mathrm{Tr}}\\{Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2}\\}$.

However, the principal difference is in the quantity $\mu_{0}$ that is
responsible for the bias in estimating $\rho$, which, from (4.7), (4.10),
(4.11) and (A.16), and because $(I-\rho
W)^{-1}{\bf{1}}=(1-\rho)^{-1}{\bf{1}}$, can be expressed in the following

Bias formula: assuming that $W{\bf{1}}={\bf{1}}$, and that
$x_{\rho}^{(2)},\ldots,x_{\rho}^{(m)}$ are the remaining vectors in the
expression (4.16) for $H_{\rho}$, the bias in the MLE is approximately given
by $-\mu_{0}/\Delta_{0}$, where

$\mu_{0}{\ =\ 2{\mathrm{Tr}}\\{Z_{\rho_{0}}H_{0}\\}}\ =\
2(1-\rho_{0})^{-1}+\sum_{l=2}^{m}(x_{\rho_{0}}^{(l)})^{\top}Z_{\rho_{0}}x_{\rho_{0}}^{(l)}.$

In the simulations in [10], the remaining columns of $X$ were realized from
independent standard normally distributed random variables, and therefore had
no particular relationship with $A$ or $W$, or with each other. In such
circumstances, they make no appreciable contribution to the bias. To see this,
consider, for simplicity, what happens if $\rho_{0}=0$. Then, for $2\leq l\leq
m$, the orthonormal vectors $x_{0}^{(l)}$ can be thought of roughly as being
columns of $X$, divided by $|x^{(l)}|\approx\sqrt{n}$ so as to have norm $1$;
‘roughly’ involves neglecting the small adjustments that are needed, so as to
be made orthogonal to ${\bf{1}}$ and to each other. As a result, using (3.3)
and (3.5) (with
$C:={\textstyle{\frac{1}{2}}}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})$), it follows
that typical values of
$(x_{0}^{(l)})^{\top}Z_{0}x_{0}^{(l)}=(x_{0}^{(l)})^{\top}Wx_{0}^{(l)}$ are
roughly centred around $n^{-1}{\mathrm{Tr}}\,W=0$, with standard deviation
$n^{-1}\bigl{(}{\mathrm{Tr}}\\{{\textstyle{\frac{1}{2}}}(W+W^{\top})^{2}\\}\bigr{)}^{1/2}$.
Hence they are individually negligible when compared to
$\tau_{0}\asymp\\{{\mathrm{Tr}}\\{WW^{\top}\\}\\}^{1/2}$, and thus have no
appreciable effect on the bias. The same is true for their sum, if $m$ is not
large. A similar, but more complicated, argument leads to the same conclusion
for $\rho_{0}\neq 0$. This phenomenon is confirmed in the simulations of [10].
However, if the underlying graph encodes some community structure, with more
edges between members of the same community than outside, and if the
regression is designed to investigate community related effects, then the
resulting contributions to $\mu_{0}$ can be as important as that from fitting
an overall mean, and can have a corresponding effect on the bias. Model $M2$
of Section 5.2 is of this form.

###  4.4 Ill conditioned graphs

If the weight matrix associated with the underlying network is ill
conditioned, there are consequences for estimation, both of $\rho$ and of
$\beta$. Taking $\beta$ first, the distribution of $\hat{\beta}_{\rho_{0}}$
from (2.8) on the model $M_{\rho_{0}}$ has covariance matrix
$\sigma^{2}(X^{\top}S_{\rho_{0}}X)^{-1}$, as in (2.11), where
$S_{\rho}:=K_{\rho}^{\top}K_{\rho}$. It is shown in Appendix A.4 that the
variance of the linear combination $a^{\top}\beta$, for any given
$m$-dimensional unit vector $a$, is bounded above by

$\frac{\sigma^{2}}{\lambda_{\mathrm{min}}(X^{\top}S_{\rho_{0}}X)}\ \leq\
\frac{\sigma^{2}}{\lambda_{\mathrm{min}}(X^{\top}X)\lambda_{\mathrm{min}}(S_{\rho_{0}})},$

where $\lambda_{\mathrm{min}}(M)$ and $\lambda_{\mathrm{max}}(M)$ denote the
smallest and largest eigenvalues of a symmetric matrix $M$. Thus, for fixed
$\sigma^{2}$, provided that $\lambda_{\mathrm{min}}(S_{\rho_{0}})$ is not
close to zero, linear combinations of $\beta$ can be accurately estimated if
$\lambda_{\mathrm{min}}(X^{\top}X)$ is large, which is the usual condition for
ordinary least squares. And, even if $\lambda_{\mathrm{min}}(S_{\rho_{0}})$ is
small, it may still be the case that
$\lambda_{\mathrm{min}}(X^{\top}S_{\rho_{0}}X)$ is comparable to
$\lambda_{\mathrm{min}}(X^{\top}X)$, if the space spanned by the columns of
$X$ does not contain vectors that are substantially shrunk by
$S_{\rho_{0}}^{1/2}$.

If $\rho_{0}$ is not accurately estimated, it will typically be the case that
$\hat{\beta}_{\rho}$ is used to estimate $\beta$, for some value
$\rho\neq\rho_{0}$. In such cases, it is shown in Appendix A.4 that linear
combinations of $\beta$ can be estimated using $\hat{\beta}_{\rho}$ to an
accuracy comparable to that obtained using $\hat{\beta}_{\rho_{0}}$, provided
that $\lambda_{\mathrm{max}}(S(\rho,\rho_{0}))$ is not large, where
$S(\rho,\rho_{0}):=S_{\rho}^{1/2}S_{\rho_{0}}^{-1}S_{\rho}^{1/2}$. However, if
$\lambda_{\mathrm{max}}(S(\rho,\rho^{\prime}))$ can be large for values of
$\rho$ and $\rho^{\prime}$ that are plausible, estimation of $\beta$ may be
rather less accurate. In particular, it is shown in Appendix A.4 that a large
value of $\lambda_{\mathrm{max}}(W^{\top}W)$ indicates that problems with the
estimation of $\beta$ may also arise; see Example 1 in Section 4.4.

For bias and variability in the estimation of $\rho$, it has already been
shown that problems can easily arise in dense networks, without $W$ having any
unusual structure. However, these problems can be accentuated by particular
networks and design matrices. Again assuming that $W{\bf{1}}={\bf{1}}$, this
is because the trace
${v(\rho_{0})=}{\mathrm{Tr}}\\{(Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2})(I-H_{\rho_{0}})\\}$
in (4.17) can be very much smaller than
${\mathrm{Tr}}\\{Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2}\\}=1/\gamma_{\rho_{0}}^{2}$.
If $\tau_{0}\geq 1$, so that $\psi_{0}:=\tau_{0}/\Delta_{0}\asymp
1/(\tau_{0}+1)$, this leads to a value of $\psi_{0}$ that is much larger than
the theoretical lower bound $\gamma_{\rho_{0}}$. Taking $\rho_{0}=0$ for
simplicity, the comparison is between
${\mathrm{Tr}}\\{(WW^{\top}+W^{2})(I-H_{0})\\}$ and
${\mathrm{Tr}}\\{WW^{\top}\\}$.

Example 1. As a first example, take that of the $(n-1)$-star, for which
$A_{1j}=A_{j1}=1$, $2\leq j\leq n$, and $A_{ij}=0$ otherwise, and take
$W:=D^{-1}A$. Then the matrix $W^{2}$ has $(W^{2})_{11}=1$, and all other
elements less than or equal to $1/(n-1)$. If $X$ consists only of the single
column ${\bf{1}}$, then $H_{0}:=H_{*}:=n^{-1}{\bf{1}}{\bf{1}}^{\top}$,
$WH_{0}=H_{*}$, and the matrix $W^{2}H_{0}=H_{*}$ has all its elements equal
to $1/n$. Thus the traces of $W^{2}$ and $W^{2}H_{0}$ are bounded for all $n$;
direct calculation gives ${\mathrm{Tr}}\\{W^{2}\\}=2$ and
${\mathrm{Tr}}\\{W^{2}H_{0}\\}=1$. The matrix $WW^{\top}$ consists of an
$(n-1)\times(n-1)$ block of $1$’s, for $2\leq i,j\leq n$, together with
$1/(n-1)$ at position $(1,1)$, and the remaining elements zero; the matrix
$WW^{\top}H$ has a corresponding $(n-1)\times(n-1)$ block, now with values
$1-1/n$, all elements in the first row have the value $1/n(n-1)$, and the
remaining elements in the first column have value $1-1/n$. Thus the diagonal
elements of $WW^{\top}(I-H_{0})$ all have values equal to $1/n$, and so its
trace is exactly $1$, whereas ${\mathrm{Tr}}\\{WW^{\top}\\}=n-1+1/(n-1)\geq
n-1$. Hence the lower bound $\gamma_{0}$ at $\rho=0$ is less than $1/(n-1)$,
but, with the simple structure corresponding to estimating the overall mean in
the model, using the formula (4.15) from Section 4.3,

$\tau_{0}^{2}\ =\ 4{\mathrm{Tr}}\\{(WW^{\top}+W^{2})(I-H_{*})\\}\ =\ 4(1+2-1)\
=\ 8$

is comparable to $1$, implying that the MLE does not give accurate estimation
of $\rho$.

For estimating $\beta$, note that, from (2.8) and (2.11), the variance of
$\hat{\beta}_{\rho}$ on the model $M_{\rho}$ is given by
$\sigma^{2}({\bf{1}}^{\top}K_{\rho}^{\top}K_{\rho}{\bf{1}})^{-1}=\sigma^{2}/\\{n(1-\rho)^{2}\\}$,
because $K_{\rho}{\bf{1}}=(1-\rho){\bf{1}}$. Hence, if $\rho$ is known,
$\beta$ can be estimated accurately by $\hat{\beta}_{\rho}$ if $n$ is large
and $\rho$ is not too close to $1$. However, in this example, if $e^{(1)}$
denotes the vector $(1,0,\ldots,0)^{\top}$, then
$\lambda_{\mathrm{max}}(W^{\top}W)\geq(e^{(1)})^{\top}W^{\top}We^{(1)}=n-1$ is
large, suggesting that using the wrong value of $\rho$ in estimating $\beta$
may cause problems. Applying (A.22) from Appendix A.4 below, and taking
$\rho_{0}=0$ to simplify the calculations, since then $S_{\rho_{0}}=I$, we
obtain the explicit expression

$\displaystyle{\rm Var}_{0}\hat{\beta}_{\rho}$ $\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{n^{2}(1-\rho)^{2}}{\bf{1}}^{\top}K_{\rho}K_{\rho}^{\top}{\bf{1}}\
=\
\frac{\sigma^{2}}{n^{2}(1-\rho)^{2}}\\{n(1-2\rho)+\rho^{2}{\bf{1}}^{\top}WW^{\top}{\bf{1}}\\}$
$\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{n^{2}(1-\rho)^{2}}\Bigl{\\{}n(1-2\rho)+\rho^{2}\Bigl{(}(n-1)^{2}+\frac{1}{n-1}\Bigr{)}\Bigr{\\}}.$

Hence the variance of $\hat{\beta}_{\rho}$ grows like
$\sigma^{2}\\{\rho/(1-\rho)\\}^{2}$ away from $\rho_{0}=0$, and is thus
substantial when the wrong value of $\rho$ is used. As a result, since $\rho$
cannot be reliably estimated, nor can $\beta$.

Conclusion 4: This example shows that having a large number of vertices of low
degree — in this case, $n-1$ vertices of degree $1$ — is not of itself enough
to imply accurate estimation of $\rho$ or $\beta$, even though it ensures that
${\mathrm{Tr}}\\{WW^{\top}\\}$ is large, and hence that $\gamma_{0}$ is small.

Adding further columns to $X$ results in a projection matrix
$H=H_{*}+{\widetilde{H}}$, with $H_{*}{\widetilde{H}}=0$,
${\widetilde{H}}^{\top}={\widetilde{H}}$ and
${\widetilde{H}}^{2}={\widetilde{H}}$, replacing the trace
${\mathrm{Tr}}\\{WW^{\top}(I-H_{*})\\}$ by
${\mathrm{Tr}}\\{WW^{\top}(I-H)\\}$; the latter is smaller, because

$\displaystyle{\mathrm{Tr}}\\{WW^{\top}(I-H_{*})\\}-{\mathrm{Tr}}\\{WW^{\top}(I-H)\\}$
(4.18) $\displaystyle\ =\
{\mathrm{Tr}}\\{WW^{\top}(I-H_{*})\\}-{\mathrm{Tr}}\\{(I-H_{*}-{\widetilde{H}})WW^{\top}(I-H_{*}-{\widetilde{H}})\\}$
$\displaystyle\ =\ {\mathrm{Tr}}\\{WW^{\top}(I-H_{*})\\}$
$\displaystyle\qquad\qquad\mbox{}-{\mathrm{Tr}}\\{(I-H_{*})WW^{\top}(I-H_{*})-2{\widetilde{H}}WW^{\top}(I-H_{*})+{\widetilde{H}}WW^{\top}{\widetilde{H}}\\}$
$\displaystyle{\ =\
{\mathrm{Tr}}\\{2(I-H_{*}){\widetilde{H}}WW^{\top}-{\widetilde{H}}WW^{\top}\\}\
=\ {\mathrm{Tr}}\\{{\widetilde{H}}WW^{\top}\\}}$ $\displaystyle\ =\
{\mathrm{Tr}}\\{({\widetilde{H}}W)({\widetilde{H}}W)^{\top}\\}\ \geq\ 0.$

The value of ${\mathrm{Tr}}\\{W^{2}H\\}$ is not materially changed, because,
apart from $W_{11}=1$, the elements of $W^{2}$ are uniformly of order
$O(1/n)$. Hence, even with the more realistic model structure, the value of
$v(\rho_{0})$ is still comparable to $1$ when $\rho_{0}=0$, and, in view of
Conclusion 3, the variability of the MLE of $\rho$ is substantial.

Example 2. The same sort of discussion can be carried out, if the network
above is generalized to a small number of stars, each of which has high
degree, the remaining vertices being leaves, and the centres of the stars
being connected among themselves. The elements of $W^{2}$, apart from the
diagonal elements corresponding to the centres of the stars, are then
uniformly small, as in (3.10), typically resulting in only moderate values of
${\mathrm{Tr}}\\{W^{2}\\}$ and ${\mathrm{Tr}}\\{W^{2}H\\}$. The block
structure of $WW^{\top}$ in the previous example is replaced by a collection
of blocks of $1$’s, one for the leaves of each star. As a result, by arguments
similar to those in the example above, the trace
${\mathrm{Tr}}\\{WW^{\top}(I-H)\\}$ becomes small if the design matrix $X$
contains columns identifying the blocks; the elements of the column
identifying a block should take the value $1$ at positions within the
corresponding block, and zero within the other blocks, the remaining elements
being chosen to ensure that the sum of the columns is ${\bf{1}}$, so that the
overall mean belongs to the design. In consequence, the value of $v(0)$ is
comparable to $1$, and, in view of Conclusion 3, the variability of the MLE of
$\rho$ is substantial in this case as well.

Example 3. A further, more extreme example, is the one considered by [16], in
which $A={\bf{1}}{\bf{1}}^{\top}-I$ is the adjacency matrix of the complete
graph on $n\geq 2$ vertices. Then, for $W=D^{-1}A$, we have

$W\ =\ W^{\top}\ =\ \frac{1}{n-1}\\{{\bf{1}}{\bf{1}}^{\top}-I\\}\ =\
H_{*}-\frac{1}{n-1}\\{I-H_{*}\\},$ (4.19)

where, as before, $H_{*}=n^{-1}{\bf{1}}{\bf{1}}^{\top}$, implying, since
$H_{*}$ and $I-H_{*}$ are idempotent and $H_{*}(I-H_{*})=0$, that

$W^{r}\ =\
H_{*}+\Bigl{(}\frac{-1}{n-1}\Bigr{)}^{r}\\{I-H_{*}\\}\quad\mbox{and}\quad
Z_{\rho}\ =\ \frac{H_{*}}{1-\rho}-\frac{I-H_{*}}{n-1+\rho}{\ =\
Z_{\rho}^{\top}}.$ (4.20)

Hence, in particular,

${{\mathrm{Tr}}\\{Z_{\rho}Z_{\rho}^{\top}+Z_{\rho}^{2}\\}=2\Bigl{\\{}\frac{1}{(1-\rho)^{2}}+\frac{n-1}{(n-1+\rho)^{2}}\Bigr{\\}},}$

and so, in view of (3.7), estimation of $\rho$ cannot be accurate if $\rho$ is
not close to $1$. However, if the model consists only of an intercept term, so
that $X={\bf{1}}$, and recalling that ${\bf{1}}^{\top}W={\bf{1}}^{\top}$, the
estimate of $\beta$ from (2.8) is given by

$\hat{\beta}_{\rho}(y)\ =\ ({\bf{1}}^{\top}(I-\rho
W)^{2}{\bf{1}})^{-1}{\bf{1}}^{\top}(I-\rho W)^{2}y\ =\
n^{-1}\sum_{j=1}^{n}y_{j},$

which is the same for all values of $\rho$, and its variance, obtained using
the general formula in (A.22) below, reduces to $1/\\{n(1-\rho_{0})^{2}\\}$ on
$M_{\rho_{0}}$, and becomes small as $n$ increases. Thus, although $\rho$ may
not be accurately estimable, there is no problem in estimating $\beta$.

Returning to the estimation of $\rho$, still with $X={\bf{1}}$, note also
that, using (4.20),

${\mathrm{Tr}}\\{(I-H_{0})WW^{\top}\\}\ =\ {\mathrm{Tr}}\\{(I-H_{*})W^{2}\\}\
=\ \frac{1}{(n-1)^{2}}{\mathrm{Tr}}\\{I-H_{*}\\}\ =\ \frac{1}{n-1}.$

Thus ${\mathrm{Tr}}\\{(I-H_{0})WW^{\top}\\}$ is an order of magnitude smaller
than ${\mathrm{Tr}}\\{WW^{\top}\\}$. This in turn, from (4.15), implies that
the variability in the MLE of $\rho$ is on the scale of $\tau_{0}\ll 1$, in
apparent contradiction to the lower bound in (3.7); this can only be the case
if the bias in the MLE is extremely large. We now investigate what happens in
more detail.

First, allowing a more general design matrix $X$ with $m$ columns, observe
that

$WX\ =\ (n-1)^{-1}\\{{\bf{1}}{\bf{1}}^{\top}X-X\\},$

so that the $i$-th column of $WX$ is a linear combination of ${\bf{1}}$ and
the $i$-th column of $X$. If $X$ contains the column ${\bf{1}}$, this implies
that the columns of $K_{\rho}X=X-\rho WX$ span the same space as the columns
of $X$, and hence that the orthogonal projection $H_{\rho}$ is the same for
all $\rho$: $H_{\rho}=H_{0}=H(X)$, and $H_{0}{\bf{1}}={\bf{1}}$. As a result,
it follows that

$H_{0}H_{*}\ =\ H_{*}H_{0}\ =\ H_{*}\quad\mbox{and}\quad(I-H_{0})H_{*}\ =\
H_{*}(I-H_{0})\ =\ 0,$ (4.21)

so that, from (4.8), (4.20) and (4.21),

$\tau_{0}^{2}\ =\
8{\mathrm{Tr}}\\{Z_{\rho_{0}}(I-H_{0})Z_{\rho_{0}}(I-H_{0})\\}\ =\
\frac{8}{(n-1+\rho_{0})^{2}}{\mathrm{Tr}}\\{I-H_{0}\\}\ =\
\frac{8(n-m)}{(n-1+\rho_{0})^{2}}\ \approx\ \frac{8}{n}\,.$

Similar calculations show that, using (4.20) and (4.21) in (4.12), there is
great simplification, yielding

$\displaystyle\Delta_{0}$ $\displaystyle=$
$\displaystyle{\textstyle{\frac{1}{2}}}\tau_{0}^{2}+2{\mathrm{Tr}}\\{-Z_{\rho_{0}}H_{0}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{0})+H_{0}Z_{\rho_{0}}^{2}\\}$
$\displaystyle=$
$\displaystyle{\textstyle{\frac{1}{2}}}\tau_{0}^{2}+2\Bigl{\\{}\frac{1}{(1-\rho_{0})^{2}}+\frac{m-1}{(n-1+\rho_{0})^{2}}\Bigr{\\}}{\
=\ \frac{2}{(1-\rho_{0})^{2}}+O(n^{-1})},$

since ${\mathrm{Tr}}\\{H_{0}-H_{*}\\}=m-1$; and then

$\displaystyle\mu_{0}$ $\displaystyle=$ $\displaystyle
2{\mathrm{Tr}}\\{H_{0}Z_{\rho_{0}}\\}\ =\
\frac{2}{1-\rho_{0}}{-\frac{m-1}{n-1+\rho}}.$

Hence it follows that $\mu_{0}/\Delta_{0}$ is not small, and that the bias in
the MLE is thus substantial.

As it happens, using (4.19) and (4.21),

$(I-H_{\rho})K_{\rho}y\ =\ (I-H_{0})(I-\rho W)y\ =\
\Bigl{\\{}1+\frac{\rho}{n-1}\Bigr{\\}}(I-H_{0})y,$

and, from (2.8), and because $(I-H_{\rho})=(I-H_{\rho})^{2}=(I-H_{0})^{2}$,
this implies that

${\hat{\sigma}}^{2}_{\rho}(y)\ =\
n^{-1}(K_{\rho}y)^{\top}(I-H_{\rho})K_{\rho}y\ =\
\Bigl{\\{}1+\frac{\rho}{n-1}\Bigr{\\}}^{2}n^{-1}y^{\top}(I-H_{0})y.$

Thus ${\hat{\sigma}}^{2}_{\rho}(y)$ varies very little as $\rho$ ranges over
the interval $(-1,1)$. What is worse, $\log({\hat{\sigma}}^{2}_{\rho}(y))$ is
the sum of a function of $\rho$ alone and a function of the data alone, as is
therefore the function $F(\rho;y)$ in (2.14) also. Thus, in this setting, the
maximum likelihood estimator, choosing $\rho\in(-1,1)$ to minimize
$F(\rho;y)$, represents a choice on which the data have no influence, so that
its variability is actually zero, and the same fixed choice of $\rho$ would
always be made.

Conclusion 5: For the complete graph, the likelihood gives no information at
all about the value of $\rho$.

For this choice of $W$, $K_{\rho}$ has eigenvalues $1-\rho$, together with
$1+\rho/(n-1)$ repeated $n-1$ times, and so the function

$\displaystyle F(\rho;y)$ $\displaystyle=$
$\displaystyle\log\\{n^{-1}y^{\top}(I-H_{0})y\\}+2\log\\{1+\rho/(n-1)\\}$
$\displaystyle\qquad\mbox{}-{\textstyle{\frac{2}{n}}}\,\bigl{\\{}\log(1-\rho)+(n-1)\log\\{1+\rho/(n-1)\\}\bigr{\\}}$
$\displaystyle=$
$\displaystyle\log\\{n^{-1}y^{\top}(I-H_{0})y\\}+{\textstyle{\frac{2}{n}}}\,\bigl{\\{}\log\\{1+\rho/(n-1)\\}-\log(1-\rho)\bigr{\\}}$

defined in (2.14) has strictly positive derivative throughout the (open)
interval $(-\infty,1)$. Hence setting its derivative equal to zero to find a
minimum would also fail here.

##  5 Estimation using a quadratic form

The conclusions of Section 4 indicate that, in dense graphs, the MLE of the
correlation parameter $\rho$ in the network disturbance model may quite
generally exhibit significant bias and variability. The ill–conditioned
networks of Section 4.4 indicate that these effects can be dramatic. Because
of the Cramér–Rao bound given in (3.7), it is not possible to avoid the
variability, whatever estimator of $\rho$ is used. However, estimators with
smaller bias can be envisaged; see [1], [19] and [18] for progress in the
context of correcting the bias in the MLE for network effects models, when the
variability is small. In this section, we move away from the MLE, and instead
propose a simple and intuitive estimator of $\rho$, which suffers much less
from bias, even in circumstances in which all estimators have substantial
variability.

###  5.1 A simple estimator of $\rho$

An important element in $n\frac{dF(\rho;y)}{d\rho}$, the quantity that is set
equal to zero for the MLE, is the quadratic form

$(K_{\rho}y)^{\top}(I-H_{\rho})(Z_{\rho}+Z_{\rho}^{\top})(I-H_{\rho})K_{\rho}y$

in $y$, as can be seen in (4.2). Here, we replace this quadratic form with a
simpler one. Since, under the model $M_{\rho_{0}}$ from (2.7), the innovations
$\boldsymbol{\nu}_{i}=(K_{\rho_{0}}\boldsymbol{\varepsilon})_{i}$, $1\leq
i\leq n$, are independent, even at neighbouring vertices, the expectation of
$(K_{\rho_{0}}\boldsymbol{\varepsilon})_{i}(K_{\rho_{0}}\boldsymbol{\varepsilon})_{j}$
is zero for all $i\neq j$. On the other hand, for $\rho\neq\rho_{0}$, the
expectation of
$(K_{\rho}\boldsymbol{\varepsilon})_{i}(K_{\rho}\boldsymbol{\varepsilon})_{j}$
under $M_{\rho_{0}}$ is typically not zero, if $i$ and $j$ are neighbours. For
example, if $\rho_{0}=0$, then, for $i\neq j$,

${\sigma^{-2}}{\mathbb{E}}_{0}\\{(K_{\rho}\boldsymbol{\varepsilon})_{i}(K_{\rho}\boldsymbol{\varepsilon})_{j}\\}\
=\ \\{(I-\rho W)(I-\rho W)^{\top}\\}_{ij}\ =\
-\rho(W_{ij}+W_{ji})+\rho^{2}(WW^{\top})_{ij},$ (5.1)

which is not zero whenever $\rho\neq 0$ is small enough, if one of $W_{ij}$
and $W_{ji}$ is non-zero. This suggests building sums

$T^{C}_{\rho}(y)\ :=\
(K_{\rho}{\hat{\varepsilon}}_{\rho}(y))^{\top}CK_{\rho}{\hat{\varepsilon}}_{\rho}(y)$
(5.2)

of products of the form
$(K_{\rho}{\hat{\varepsilon}}_{\rho}(y))_{i}(K_{\rho}{\hat{\varepsilon}}_{\rho}(y))_{j}$,
where ${\hat{\varepsilon}}_{\rho}(y):=K_{\rho}^{-1}(I-H_{\rho})K_{\rho}y$ is
the vector of residuals after fitting $M_{\rho}$, as derived in (2.10); thus

$T^{C}_{\rho}(y)\ =\
y^{\top}K_{\rho}^{\top}(I-H_{\rho})C(I-H_{\rho})K_{\rho}y.$ (5.3)

The matrices $C$, having non-negative coefficients that do not depend on the
data or on $\rho$, and with all diagonal elements equal to zero, are to be
suitably chosen. Any such sum $T^{C}_{\rho}(y)$, for fixed $C$ and $y$, should
take values close to zero when $\rho$ is close to $\rho_{0}$, and this can
potentially be exploited to gain information about the value of $\rho$. As
indicated in (5.1), the elements of $C$ should be chosen so as to give larger
weight to pairs $(i,j)$ where $W_{ij}+W_{ji}$ is large, and a natural choice
is thus to take $C={\textstyle{\frac{1}{2}}}(W+W^{\top})$, or, equivalently,
to take $C=W$. Another intuitive choice, reflecting the structure of the
underlying graph, rather than the model weights, is to take $C=A$, roughly
analogous to the lag-one autocovariance in time series analysis.

Recalling that, on $M_{\rho_{0}}$,
$\mathbf{y}=X\beta+K_{\rho_{0}}^{-1}\boldsymbol{\nu}$, it follows that, on
$M_{\rho_{0}}$,

${\hat{\varepsilon}}_{\rho}(\mathbf{y})\ =\
K_{\rho}^{-1}(I-H_{\rho})K_{\rho}K_{\rho_{0}}^{-1}\boldsymbol{\nu},$

since $(I-H_{\rho})K_{\rho}X=0$, and hence that, on $M_{\rho_{0}}$,

$T^{C}_{\rho}(\mathbf{y})\ =\
\boldsymbol{\nu}^{\top}Q_{C}(\rho,\rho_{0})\boldsymbol{\nu},$ (5.4)

where, much as in (4.4),

$Q_{C}(\rho,\rho_{0})\ :=\
(K_{\rho}K_{\rho_{0}}^{-1})^{\top}(I-H_{\rho})C(I-H_{\rho})K_{\rho}K_{\rho_{0}}^{-1}.$
(5.5)

From (5.4) and (5.5), we have

${\mathbb{E}}_{\rho_{0}}T^{C}_{\rho_{0}}(\mathbf{y})\ =\
\sigma^{2}{\mathrm{Tr}}\\{(I-H_{\rho_{0}})C\\}\ =\
-\sigma^{2}{\mathrm{Tr}}\\{H_{\rho_{0}}C\\}\ =:\
\sigma^{2}{\hat{\mu}}^{C}_{\rho_{0}},$ (5.6)

since $C$ has zero diagonal, and so ${\mathrm{Tr}}\\{C\\}=0$. Equation (5.6)
suggests the following ‘method of moments’ procedure for estimating $\rho$.

Elementary estimate of $\rho$: find $\rho={\hat{\rho}}_{C}$ so that

$T^{C}_{\rho}(y)+{\hat{\sigma}}_{\rho}^{2}(y){\mathrm{Tr}}\\{H_{\rho}C\\}\ =\
0.$ (5.7)

Equation (5.7) is used to correct the bias that would result, in view of
(5.6), if, in estimating $\rho$, $T^{C}_{\rho}(y)$ were itself set equal to
zero.

Procedure:

Input:

a network with weight matrix $W$, a linear model with design matrix $X$ and
observed data $y$, and a non-negative coefficient matrix $C$ with zero
diagonal.

Function:

for any value of $\rho$:

  1. 1.

Calculate the trace ${\mathrm{Tr}}\\{H_{\rho}C\\}$, using (2.9) and (2.2);

  2. 2.

Calculate the parameter estimates ${\hat{\beta}}_{\rho}(y)$ and
${\hat{\sigma}}^{2}_{\rho}(y)$ using (2.8);

  3. 3.

Calculate ${\hat{\varepsilon}}_{\rho}(y)$, using ${\hat{\beta}}_{\rho}(y)$ and
(2.10);

  4. 4.

Calculate the value of $T^{C}_{\rho}(y)$ from (5.2).

Solution:

use a root finding algorithm to find the zero of
$T^{C}_{\rho}(y)+{\hat{\sigma}}_{\rho}^{2}(y){\mathrm{Tr}}\\{H_{\rho}C\\}$,
which is the estimate ${\hat{\rho}}_{C}$ of the correlation $\rho$ using
$T^{C}$.

In this paper, we do not attempt to give detailed conditions under which the
procedure using (5.7) gives good estimates ${\hat{\rho}}_{C}$ of $\rho$. It is
clear from the example of the complete graph in Section 4.4 that this cannot
always be so; indeed, for the complete graph, (5.7) has no solution. Instead,
we carried out some simulations, illustrating that ${\hat{\rho}}_{C}$ can be
effective, within the limitations of what is possible. These are discussed in
Section 5.2 below. In the remainder of this section, we give a rough guide to
its precision from a theoretical standpoint, as well as proposing a
permutation procedure for assessing its precision based on the data to hand.

The scale of variation of ${\hat{\rho}}_{C}$
To determine the scale of variation about $\rho_{0}$ of the resulting
estimator ${\hat{\rho}}_{C}$, write

$U^{C}_{\rho}(y)\ :=\
T^{C}_{\rho}(y)+{\hat{\sigma}}_{\rho}^{2}(y){\mathrm{Tr}}\\{H_{\rho}C\\}.$
(5.8)

As in the discussion of the MLE around (4.10) and (4.11) in Section 4, a
plausible measure is given by
${\hat{\psi}}^{C}_{\rho_{0}}:={\hat{\tau}}^{C}_{\rho_{0}}/|{\widehat{\Delta}}^{C}_{\rho_{0}}|$,
to be estimated by
${\hat{\tau}}^{C}_{{\hat{\rho}}_{C}}/|{\widehat{\Delta}}^{C}_{{\hat{\rho}}_{C}}|$,
where

$({\hat{\tau}}^{C}_{\rho_{0}})^{2}\ :=\ \sigma^{-4}{\rm
Var}_{\rho_{0}}(U^{C}_{\rho_{0}}(\mathbf{y}))\quad\mbox{and}\quad{\widehat{\Delta}}^{C}_{\rho_{0}}\
:=\
{-}\sigma^{-2}{\mathbb{E}}_{\rho_{0}}\Bigl{\\{}\frac{dU^{C}_{\rho}(\mathbf{y})}{d\rho}\Bigr{|}_{\rho=\rho_{0}}\Bigr{\\}}.$

The quantity $({\hat{\tau}}^{C}_{\rho_{0}})^{2}$ is dominated by the trace of
a matrix,

${\mathrm{Tr}}\\{(I-H_{\rho_{0}})(C+C^{\top})(I-H_{\rho_{0}})(C+C^{\top})\\},$
(5.9)

which is reminiscient of the definition of $\tau^{2}_{0}$ in (4.8), and the
quantity ${\widehat{\Delta}}^{C}_{\rho_{0}}$ has leading term

${\mathrm{Tr}}\bigl{\\{}(I-H_{\rho_{0}}){[Z_{\rho_{0}}^{\top}(I-H_{\rho_{0}})-Z_{\rho_{0}}H_{\rho_{0}}]}(C+C^{\top})\bigr{\\}}.$
(5.10)

The derivations, and more detailed expressions, are given in Section A.5. As
before, the scale of variation

${\hat{\psi}}^{C}_{\rho_{0}}\ :=\
{\hat{\tau}}^{C}_{\rho_{0}}/|{\widehat{\Delta}}^{C}_{\rho_{0}}|,$ (5.11)

determined from (5.9) and (5.10), need not correspond precisely to the
standard deviation of ${\hat{\rho}}_{C}$ on $M_{\rho_{0}}$, because
$U^{C}_{\rho}(y)$ is not linear as a function of $\rho$, and because the
derivative of $U^{C}_{\rho}(y)$ may fluctuate appreciably about its
expectation.

To get an idea of the magnitude of ${\hat{\psi}}^{C}_{\rho_{0}}$, let $C=W$,
and suppose that $W{\bf{1}}={\bf{1}}$ and that $X={\bf{1}}$, so that only the
overall mean is being fitted. Then, as in Section 4.3,
$K_{\rho}X=(1-\rho){\bf{1}}$, implying that, for all $\rho$,
$H_{\rho}=H_{*}:=n^{-1}{\bf{1}}{\bf{1}}^{\top}=WH_{\rho}$ and that
$Z_{\rho}H_{\rho}=(1-\rho)^{-1}H_{*}$. As a result, we easily deduce that

$({\hat{\tau}}^{W}_{\rho_{0}})^{2}\ =\
{2}({\mathrm{Tr}}\\{W^{2}+WW^{\top}-H_{*}WW^{\top}\\}-1)+O(n^{-1}),$

and that

${\widehat{\Delta}}^{C}_{\rho_{0}}\ =\
-{\mathrm{Tr}}\\{WZ_{\rho_{0}}+WZ_{\rho_{0}}^{\top}-H_{*}WZ_{\rho_{0}}^{\top}\\}+(1-\rho)^{-1}+n^{-1}{\mathrm{Tr}}\\{(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{*})\\},$

where we have used the full expressions in (A.25) and (A.30) instead of (5.9)
and (5.10) to obtain the smaller order terms. In particular, if $\rho_{0}=0$,
so that $Z_{\rho_{0}}=W$,
${{\textstyle{\frac{1}{2}}}}({\hat{\tau}}^{W}_{0})^{2}=-{\widehat{\Delta}}^{C}_{0}+O(n^{-1})$,
and the magnitude of ${\hat{\psi}}^{C}_{0}$ is roughly
$2/{\hat{\tau}}^{W}_{0}$; this remains the case for $\rho_{0}$ not too far
from zero. As in the discussion in Section 3, this implies that if, for
instance, $W_{ij}\leq c/n$ for all $i,j$, then the scale of variation is at
least $1/\sqrt{c}$, and is appreciable if $c$ is not large. For $W:=D^{-1}A$,
this corresponds to having a dense underlying graph, as before.

The justification for using (5.7) to estimate $\rho$ relies only on the first
two joint moments of the vector $\boldsymbol{\nu}$ of errors, and does not
assume that the components are independent or normally distributed, though
they should be uncorrelated. However, the formula (A.25), used in the
theoretical discussion of the scale of random variation
${\hat{\psi}}^{C}_{\rho_{0}}$ of the estimator ${\hat{\rho}}_{C}$, is already
based on stronger assumptions. For instance, if the errors $\boldsymbol{\nu}$
are independent and have zero mean, but if
${\mathbb{E}}\\{\sigma^{-4}\boldsymbol{\nu}_{1}^{4}\\}\neq 3$, the expression
(3.4) for the variance of a quadratic form needs to be used in place of (3.5).

A permutation procedure for assessing the accuracy of ${\hat{\rho}}_{C}$
Instead of assessing the precision of ${\hat{\rho}}_{C}$ using
${\hat{\psi}}^{C}_{\rho_{0}}$, some form of simulation procedure can be used.
One possibility is the following permutation algorithm. Start by computing the
esimated $\boldsymbol{\nu}$-residuals

${\hat{\nu}}(y)\ :=\
K_{{\hat{\rho}}_{C}}{\hat{\varepsilon}}_{{\hat{\rho}}_{C}}(y)\ =\
(I-H_{{\hat{\rho}}_{C}})K_{{\hat{\rho}}_{C}}y.$ (5.12)

Since the components of $\boldsymbol{\nu}$ are independent and identically
distributed, this suggests that synthetic data could be derived by replacing
the $\boldsymbol{\nu}$-residuals $({\hat{\nu}}_{i}(y),\,1\leq i\leq n)$ by
permuted residuals

${\hat{\nu}}^{\pi}_{i}(y)\ :=\ {\hat{\nu}}_{\pi(i)}(y),\qquad 1\leq i\leq n,$
(5.13)

where $\pi$ is a randomly chosen permutation of $\\{1,2,\ldots,n\\}$. These
can be used to determine residuals

${\hat{\varepsilon}}^{\pi}(y)\ :=\ K_{{\hat{\rho}}_{C}}^{-1}{\hat{\nu}}^{\pi}$
(5.14)

for the original model, and hence to generate synthetic data
$y^{\pi}:=X{\hat{\beta}}+{\hat{\varepsilon}}^{\pi}(y)$. These can in turn be
analyzed in the same way as the original data, to provide an estimate
${\hat{\rho}}_{C}^{\pi}$ of ${\hat{\rho}}_{C}(y)$. Repeating this procedure
for a large number $M$ of independently chosen random permutations $\pi_{j}$,
$1\leq j\leq M$, yields a selection of estimates ${\hat{\rho}}_{C}^{\pi_{j}}$
of ${\hat{\rho}}_{C}(y)$, whose empirical variability should be much like that
of the variability of ${\hat{\rho}}_{C}$ about the model value $\rho_{0}$,
when $\rho_{0}={\hat{\rho}}_{C}(y)$.

###  5.2 Simulations

In this section, we illustrate the effectiveness of ${\hat{\rho}}_{C}$ using
simulated data, showing that it substantially reduces the bias exhibited by
the MLE, while its precision remains constrained by the lower bound (3.7). In
Table 1, data were simulated at the vertices of a Bernoulli graph $G(n,p)$,
with $n=100$ and with three choices of $p$. The data were simulated from the
model $M_{\rho}$ with six different values of $\rho$. The parameters
$\beta=(1,0.5,0.4,0.3)$ were held fixed, and the design matrix $X$ was chosen
to have ${\bf{1}}$ as its first column, with its remaining columns consisting
of independent realizations of standard normal random variables. For each
combination of $p$ and $\rho$, the elements of $X$ and the network $G$ were
sampled anew; then 100 replicates were generated, starting from realizations
of independent standard normal random variables $\boldsymbol{\nu}$, and using
(1.2) to determine the values of $\varepsilon$, taking
$W_{ij}:=D_{i}^{-1}A_{ij}$ if $D_{i}\geq 1$, and $W_{ij}=0$ if $i$ is an
isolated vertex. Each cell in Table 1 contains the mean of the corresponding
100 values of ${\hat{\rho}}_{W}$, calculated using (5.7) with $C=W$, together
with the standard error of a single estimate. It is clear that the standard
errors are substantial. The means are close to the truth, within the limits of
variability; for $100$ replicates, the standard error of the mean is one tenth
of that of a single estimate. The estimates of $\beta$ were uniformly good,
and are not presented. In a Bernoulli graph $G(n,p)$, the degrees of the
vertices are reasonably homogeneous, being concentrated around their mean
$(n-1)p$ when $p\gg 1/n$. Thus, from (3.10), calculating as if all the degrees
were exactly $(n-1)p$, the lower bound (3.7) on the variance of an estimator
when $\rho=0$ gives $\gamma_{\rho}$ approximately equal to $\sqrt{p/2}$. The
actual values of $\gamma_{\rho}$ in the simulations depend on the value of
$\rho$ and on the particular realization of the Bernoulli random graphs that
were sampled. The value of $\rho$ made little difference, but the randomness
in the degrees was still enough, for $n=100$, to reduce the value of
$\gamma_{\rho}$ by about 10%; typical values are quoted in the final line of
Table 1(a). They indicate that the observed variablility of the quadratic form
estimator is not far from the best that could be achieved in this setting.

$p$
---
|  | 0.0975 | 0.19 | 0.36
| -0.2 | -0.21 (0.20) | -0.24 (0.30) | -0.24 (0.44)
| -0.1 | -0.14 (0.23) | -0.10 (0.35) | -0.11 (0.46)
$\rho$ | 0.0 | 0.01 (0.23) | -0.03 (0.37) | 0.02 (0.47)
| 0.1 | 0.08 (0.24) | 0.09 (0.31) | 0.08 (0.46)
| 0.2 | 0.21 (0.23) | 0.16 (0.34) | 0.11 (0.50)
| 0.3 | 0.28 (0.21) | 0.23 (0.32) | 0.27 (0.47)
CRLB  | 0.0 | (0.20) | (0.29) | (0.39)
(a) Means of 100 simulated estimates of $\rho$ (with the standard error of a
single estimate) calculated using (5.7) with $C=W$. The Cramér–Rao Lower Bound
(CRLB) line gives the values of the lower bound $\gamma_{\rho_{0}}$ derived
using (3.7).

$p$
---
|  | 0.0975 | 0.19 | 0.36
| -0.2 | -0.26 (0.21) | -0.28 (0.28) | -0.33 (0.39)
| -0.1 | -0.13 (0.20) | -0.22 (0.34) | -0.28 (0.37)
$\rho$ | 0.0 | -0.03 (0.21) | -0.10 (0.30) | -0.19 (0.42)
| 0.1 | 0.02 (0.21) | -0.05 (0.30) | -0.11 (0.41)
| 0.2 | 0.13 (0.20) | 0.14 (0.31) | -0.05 (0.42)
| 0.3 | 0.19 (0.20) | 0.16 (0.30) | 0.02 (0.33)
CRLB  | 0.0 | (0.20) | (0.29) | (0.39)
(b) Means of 100 simulated MLE estimates of $\rho$. The standard errors of a
single observation are consistent with the values given in the CRLB line.

Table 1: Comparison between ${\hat{\rho}}_{C}$ and
${\hat{\rho}}_{\mathrm{ML}}$ for a network disturbance model based on
$G(100,p)$ with $W=D^{-1}A$.

For comparison, we give the results in Table 1(b) of using the MLE in the same
setting. The standard deviations were uniformly close to the values given in
the CRLB line in the table above, and were thus a little smaller than those in
the experiments using (5.7) with $C=W=D^{-1}A$. However, the MLE exhibits a
systematic bias, which becomes progressively worse as the density of the graph
increases. Even with edge density $p=0.0975$, the smallest value investigated,
the means of the estimates of $\rho$ were consistently too small, by amounts
comparable to $0.05$, whereas the whole range of values of $\rho$ being
investigated was only $0.5$.

In the second experiment, the graph $G(n,p)$ was replaced by an Erdős–Rényi
mixture model, consisting of two classes with $50$ vertices in each; the
within class probabilities were taken to be $2p(1-p)$ and the between class
probabilities to be $p$, for $p=0.05,0.1$ and $0.2$. The value of $\rho$ was
kept constant at $0.1$, and the parameters $\beta$ were as in the first
experiment. The results of estimation using ${\hat{\rho}}_{C}$ are given in
Table 2(a). In the first row, model $M1$, the design matrix was also
constructed as before, with the first column taken to be ${\bf{1}}$, and with
the remaining columns chosen to have independent standard normal entries. For
the second row, model $M2$, the second column of the design matrix had
elements $+1$ for the vertices in the first of the classes, and $-1$ for
vertices in the second, corresponding to an inter-block contrast, which
significantly changes the bias correction ${\mathrm{Tr}}\\{H_{\rho}C\\}$. In
each cell, the results of 1,000 replicates are given. At this accuracy, a
small negative bias in the estimates of $\rho$ using ${\hat{\rho}}_{C}$ can be
discerned, although it is negligible when compared to the standard errors of
the estimates. The approximate ideal standard errors derived from (3.7) for
$\rho=0.1$ are given in the final row. In Table 2(b), the results of using
maximum likelihood are given for comparison; here, the bias is clearly
visible.

$p$ | 0.05 | 0.1 | 0.2
---|---|---|---
$M1$ | 0.102 (0.20) | 0.104 (0.28) | 0.111 (0.43)
$M2$ | 0.094 (0.21) | 0.082 (0.29) | 0.098 (0.42)
CRLB | (0.18) | (0.26) | (0.35)
(a) Means of 1,000 simulated estimates of $\rho$ (with the standard error of a
single estimate) using (5.7) with $C=W$. The CRLB line gives the values of the
lower bound $\gamma_{\rho_{0}}$ derived using (3.7).

$p$ | 0.05 | 0.1 | 0.2
---|---|---|---
$M1$ | 0.06 (0.18) | 0.01 (0.26) | -0.09 (0.37)
$M2$ | 0.05 (0.19) | -0.02 (0.27) | -0.14 (0.37)
CRLB | (0.18) | (0.26) | (0.35)
(b) Means of 1,000 simulated estimates of $\rho$ (with the standard error of a
single estimate) using the MLE. The standard errors of a single observation
are consistent with the values given in the CRLB line.

Table 2: Comparison between ${\hat{\rho}}_{C}$ and
${\hat{\rho}}_{\mathrm{ML}}$ for two network disturbance models based on an
Erdős–Rényi mixture graph, each with $W:=D^{-1}A$ and with $\rho=0.1$.

Table 3 gives the results of assessing the variability of ${\hat{\rho}}_{W}$
using the permutation algorithm based on (5.12)–(5.14). The experiments were
carried out in the context of the network disturbance model $M2$ based on an
Erdős–Rényi mixture graph. In order to make a sensible comparison with the
variability of the estimates ${\hat{\rho}}_{W}$ derived from repeatedly
simulating data from the model, it was necessary to apply the permutation
algorithm to a realization of the model for which the original estimate
${\hat{\rho}}_{W}$ was reasonably close to the true value $\rho_{0}=0.1$.
Because of the substantial variability in ${\hat{\rho}}_{W}$, this was not the
case in most realizations. In each experiment, the simulated variability of
${\hat{\rho}}_{W}^{\pi}$ agreed well with the variability of
${\hat{\rho}}_{W}$ under the model, and the value of
${\hat{\psi}}^{W}_{\rho_{0}}$ was also quite similar.

$p$ | ${\hat{\rho}}_{W}$ | $E({\hat{\rho}}_{W}^{\pi})$ | $SE({\hat{\rho}}_{W}^{\pi})$ | ${\hat{\psi}}^{W}_{\rho_{0}}$ | $SE({\hat{\rho}}_{W})$
---|---|---|---|---|---
0.05 | 0.094 | 0.092 | 0.19 | 0.17 | 0.21
0.1 | 0.107 | 0.112 | 0.29 | 0.26 | 0.29
0.2 | 0.091 | 0.073 | 0.43 | 0.42 | 0.42
|  |  |  |  |
Table 3: Means $E({\hat{\rho}}_{W}^{\pi})$ and standard errors
$SE({\hat{\rho}}_{W}^{\pi})$ of 1,000 simulated estimates of
${\hat{\rho}}_{W}$, using the permutation procedure, for the model $M2$ with
$\rho=0.1$, together with the scale of variation
${\hat{\psi}}^{W}_{\rho_{0}}$, calculated using (5.11), (A.25) and (A.30), and
the corresponding standard errors $SE({\hat{\rho}}_{W})$ of ${\hat{\rho}}_{W}$
from Table 2.

In Figure 1, a histogram of $1,000$ values of ${\hat{\rho}}_{W}$, simulated
from the model $M2$, but now with $\rho=0.3$, is shown next to a histogram of
$1,000$ values of ${\hat{\rho}}_{W}^{\pi}$, generated using the permutation
method described above from a realization of the model $M2$ in which
${\hat{\rho}}_{W}=0.302$. The two histograms are very similar, and neither is
substantially skewed, even though the value $0.3$ chosen for $\rho$ is close
enough to the upper limit of $1$ for some compression in the upper tail to be
expected.

(a) Histogram of ${\hat{\rho}}_{W}$

(b) Histogram of ${\hat{\rho}}_{W}^{\pi}$

Figure 1: Histogram of $1,000$ values of ${\hat{\rho}}_{W}$ for the model $M2$
with $\rho=0.3$, together with a histogram of $1,000$ values of the
permutation estimates ${\hat{\rho}}_{W}^{\pi}$ derived from a realization of
the same model in which ${\hat{\rho}}_{W}=0.302$

##  6 Network effects models

The network effects model (1.4) can be written in the form

$y\ =\ K_{\rho}^{-1}X\beta+\varepsilon,$ (6.1)

where the errors $\varepsilon$ in $y$ have the same dependence as in the
network disturbance model, but the model structure explaining $y$ also changes
with $\rho$. The corresponding log-likelihood, as in [3, Equation (1.7)], is
given by

${l(\theta;y)\ :=\
}-{\textstyle{\frac{n}{2}}}\\{\log(\sigma^{2})+\log(2\pi)\\}-\frac{1}{2\sigma^{2}}(K_{\rho}y-X\beta)^{\top}(K_{\rho}y-X\beta)+\log\det
K_{\rho},$ (6.2)

where $\theta:=(\beta,\sigma^{2},\rho)$. Since the structural element also
contains information about $\rho$, this is reflected in the estimation
procedures. In particular, the estimate of $\beta$ using maximum likelihood,
for a given value of $\rho$, is

${\hat{\beta}}_{\rho}(y)\ :=\ (X^{\top}X)^{-1}XK_{\rho}y,$

as can be deduced from ordinary least squares (2.1) applied to the data
$y^{\prime}:=K_{\rho}y$. As a result, if the true model is
$M_{\rho^{\prime}}$, for $\rho^{\prime}\neq\rho$, we have

${\mathbb{E}}_{\rho^{\prime}}{\hat{\beta}}_{\rho}(\mathbf{y})\ =\
(X^{\top}X)^{-1}XK_{\rho}K_{\rho^{\prime}}^{-1}X\beta.$ (6.3)

This expectation is in general different from $\beta$, in which case the MLE
of $\beta$ is likely to be biased if $\rho$ is incorrectly estimated. Hence,
in order to get good estimates of the structural component $\beta$, it is
important to have an accurate estimate of $\rho$.

Now, differentiating (6.2) with respect to $\rho$, the $\rho$-component of the
score function analogous to (3.1) takes the form

$\displaystyle{\tilde{s}}^{(\rho)}(\theta;y)$ $\displaystyle=$
$\displaystyle\sigma^{-2}(K_{\rho}y-X\beta)^{\top}Z_{\rho}K_{\rho}y-{\mathrm{Tr}}\\{Z_{\rho}\\},$
(6.4)

where $Z_{\rho}:=WK_{\rho}^{-1}$. Its variance at $\theta_{0}$ under the model
with $\theta=\theta_{0}$ is now given by

$V_{1}\ :=\ {\rm
Var}_{\theta_{0}}\,{\tilde{s}}^{(\rho)}(\theta_{0};\mathbf{y})\ =\
{\mathrm{Tr}}\\{Z_{\rho_{0}}^{2}+Z_{\rho_{0}}Z_{\rho_{0}}^{\top}\\}+\sigma_{0}^{-2}|Z_{\rho_{0}}X\beta_{0}|^{2},$
(6.5)

containing an extra term, as compared to (3.6), that involves the parameter
$\beta_{0}/\sigma_{0}$. This implies that the lower bound $1/V_{1}$ for the
variance of an unbiased estimator is smaller than in the analogous network
disturbance model. This does not, of itself, mean that estimation of $\rho$
can be made more precise, since there is no guarantee that the lower bound can
be attained; instead, a closer analysis is needed.

As an example, suppose that $X$ consists solely of the column ${\bf{1}}$, with
parameter $\beta$, and that the weight matrix $W$ satisfies
$W{\bf{1}}={\bf{1}}$. Then the additional term is just
$n(\beta_{0}/(1-\rho_{0})\sigma_{0})^{2}$, a quantity that grows with $n$. On
the other hand, the network effects model (6.1) translates into

$y\ =\ (1-\rho)^{-1}\beta{\bf{1}}+\varepsilon\ =\
{\tilde{\beta}}{\bf{1}}+\varepsilon,$

which can also be viewed as a network disturbance model, with a new structural
parameter ${\tilde{\beta}}=\beta/(1-\rho)$ in place of $\beta$. Thus
${\tilde{\beta}}$ can be accurately estimated, if the network is large enough,
as can be deduced from (2.11), but, as before, any estimator of $\rho$
exhibits substantial variability if the network is dense. This explains the
simulation results observed in [10], in which the maximum likelihood
estimators of both $\rho$ and the coefficient of the constant exhibited
substantial variability, and were strongly negatively correlated. In the same
experiments, the remaining $m-1$ columns $X^{\prime}$ of $X$ had independent,
normally distributed entries with mean zero and variance $1$. In this case,
with $\beta_{1}/(1-\rho)$ replaced by a parameter ${\tilde{\beta}}$ and with
$\beta^{\prime}:=(\beta_{2},\ldots,\beta_{m})^{T}$, the extra term in the
variance of the score function reduces to

$\sigma^{-2}{\mathbb{E}}\\{|Z_{\rho}X^{\prime}\beta^{\prime}|^{2}\\}\ =\
(|\beta^{\prime}|/\sigma)^{2}{\mathrm{Tr}}\\{Z_{\rho}Z_{\rho}^{\top}\\}$

(where the expectation is with respect to the random choice of $X^{\prime}$),
and, as in (3.12), this does not become large with $n$ for uniformly dense
graphs. Hence there are still substantial limitations on the precision of
estimation of $\rho$, for Cramér–Rao reasons, and, in consequence, on the
estimation of $\beta_{1}={\tilde{\beta}}(1-\rho)$ also.

The same phenomenon makes itself felt if the model structure $X$ is more
complicated, but the columns of $WX$ all belong to the linear span of the
columns of $X$; $WX=X\Gamma$, for an $m\times m$ matrix $\Gamma$. In this
case, the model structure $K_{\rho}^{-1}X\beta$ can be expressed as

$\displaystyle K_{\rho}^{-1}X\beta$ $\displaystyle=$ $\displaystyle(I-\rho
W)^{-1}X\beta\ =\ X\beta+\sum_{l\geq 1}\rho^{l}W^{l}X\beta$ (6.6)
$\displaystyle=$ $\displaystyle X\beta+\sum_{l\geq 1}\rho^{l}X\Gamma^{l}\beta\
=\ X(I-\rho\Gamma)^{-1}\beta\ =:\ X\gamma,$

where the new parameter $\gamma:=(I-\rho\Gamma)^{-1}\beta$. Once again, this
yields a network disturbance model, in which the parameter $\gamma$ can be
accurately estimated if $n$ is large, whereas, in dense networks, the
estimator of $\rho$ may exhibit substantial variability. For the original
parameter $\beta$, this would express itself in variable, but strongly
dependent, estimators of the components of $\beta=(I-\rho\Gamma)\gamma$, which
are linked through the common value of the estimator of $\rho$.

In general, $WX$ can be written in the form

$WX\ =\ X\Gamma+E,$ (6.7)

where $E$ is an $n\times m$ matrix, whose columns are orthogonal to those of
$X$, which measures how far $WX$ lies from the linear span of the columns of
$X$. In practice, except for the model in which only the overall mean is to be
estimated, and $W{\bf{1}}={\bf{1}}$ by construction, it is usual to expect
that $E$ is not precisely the zero matrix. It is shown in Appendix A.6 that,
if the quantity

$\lambda_{\mathrm{max}}\\{E^{\top}E\\}+{\mathrm{Tr}}\\{W^{2}+WW^{\top}\\}\ =\
\lambda_{\mathrm{max}}\\{(WX)^{\top}(I-H)WX\\}+{\mathrm{Tr}}\\{W^{2}+WW^{\top}\\}$

is not large, then neither $\rho$ nor $\beta$ can be accurately estimated.

Conclusion 6: In the network effects model, if the columns of $WX$ (almost)
belong to the space generated by the columns of $X$, then estimating $\rho$
and the parameters $\beta$ in uniformly dense networks may both be
problematic.

##  7 Conclusion

Maximum likelihood estimation of the correlation parameter $\rho$ in the
network disturbance and network effects models suffers from substantial bias
and variability when the underlying graph is dense, as has been documented,
for instance, in [10], [16], [4], [11] and [8]. In the context of the network
disturbance model (2.7), we have shown that, in dense networks, there is a
lower bound to the variability of any reasonable estimator of $\rho$, as a
consequence of Cramér–Rao theory, and that increasing the size of the network
does not reduce the variability, if the density remains constant. We also
demonstrate that, under similar circumstances, the maximum likelihood
estimator is significantly biased. For certain ill conditioned networks, these
effects are particularly marked. In Section 5.1, we propose a simple estimator
of $\rho$ that has little bias, though its variability in dense networks is
still substantial, because all estimators are subject to the limitations
imposed by the Cramér–Rao lower bound. We also suggest a permutation algorithm
for assessing the variability of the estimator, on the basis of data
consisting of a single network. The effectiveness of the new estimator is
illustrated by simulation in a few simple settings in Section 5.2. There is a
discussion of network effects models in Section 6. Here, in circumstances in
which $\rho$ cannot be accurately estimated, the estimates of the structural
parameters $\beta$ are also affected. An explicit quantity is given that
indicates whether such problems are likely to occur. It is expressed in terms
of the design matrix $X$ and the weight matrix $W$.

## Acknowledgement

ADB thanks the mathematics departments of the University of Melbourne and
Monash University, for their kind hospitality while part of the work was
undertaken. GR was supported in part by EP/T018445/1 and EP/R018472/1. Both
authors express their warm appreciation of the many helpful suggestions by the
referees; these have greatly improved the presentation.

## References

  * [1] Y. Bao (2013) Finite sample bias of the QMLE in spatial autoregressive models. Econometric Theory 29, 68–89.
  * [2] D. Dittrich, R. T. A. J. Leenders & J. Mulder (2017) Bayesian estimation of the network autocorrelation model. Social Networks 48, 213–236.
  * [3] P. Doreian (1989) Network autocorrelation models: problems and prospects. In: Spatial Statistics: Past, Present, Future, D. A. Griffin (Ed.), Michigan Document Services, Ann Arbor: 369–389.
  * [4] S. Farber, A. Páez & E. Volz (2010) Topology, dependency tests and estimation bias in network autoregressive models. In: Progress in Spatial Analysis, A. Páez, J. Gallo, R. N. Bullung & S. Dall’erba (Eds), Springer, Berlin: 29–57.
  * [5] W. K, Härdle & L. Simar (2015) Applied Multivariate Statistical Analysis, 4th Edn. Springer, Heidelberg.
  * [6] R. A. Horn & C. R. Johnson (2012) Matrix analysis. Cambridge University Press.
  * [7] H. H. Kelejian & I. R. Prucha (1999) A generalized moments estimator for the autoregressive parameter in a spatial model. International Economic Review 40, 509–533.
  * [8] M. La Rocca, G. C. Porzio, M. P. Vitale & P. Doreian (2018) Finite sample behavior of MLE in network autocorrelation models. In: Classification, (Big) Data Analysis and Statistical Learning, F. Mola, C. Conversano & M. Vichi (Eds), Springer, Cham: 43–50.
  * [9] E. L. Lehmann & G. Casella (2006) Theory of point estimation, 2nd Edn. Springer, New York.
  * [10] M. S. Mizruchi & E. J. Neuman (2008) The effect of density on the level of bias in the network autocorrelation model. Social Networks 30, 190–200.
  * [11] E. J. Neuman & M. S. Mizruchi (2010) Structure and bias in the network autocorrelation model. Social Networks 32, 290–300.
  * [12] K. Ord (1975) Estimation methods for models of spatial interaction. Journal of the American Statistical Association 70, 120–126.
  * [13] C. R. Rao (1973) Linear statistical inference and its applications, 2nd Edn. Wiley, New York.
  * [14] S. R. Searle & M. H. J. Gruber (1971) Linear models. Wiley Online Library, vol 10.
  * [15] E. Seneta (2006) Non-negative matrices and Markov chains, Springer Science & Business Media.
  * [16] T. E. Smith (2009) Estimation bias in spatial models with strongly connected weight matrices. Geographical Analysis 41, 307–332.
  * [17] C. Strumann (2019) Hodges–Lehmann estimation of static panel models with spatially correlated disturbances. Computational Econometrics 53, 141–168.
  * [18] Z. Yang (2015) A general method for third-order bias and variance corrections on a nonlinear estimator. Journal of Econometrics 186, 178–200.
  * [19] D. Yu, P. Bai & C. Ding (2015) Adjusted quasi-maximum likelihood estimator for mixed regressive, spatial autoregressive model and its small sample bias. Computational Statistics and Data Analysis 87, 116–135.

## Appendix

###  A.1 Rewriting the likelihood equation for $\rho$

To minimize the quantity $F(\rho;y)$ of (2.14), it is usual to look for values
of $\rho$ such that $\frac{d}{d\rho}F(\rho;y)=0$. First, as for (3.1),

$\frac{d}{d\rho}\,\log\det K_{\rho}\ =\ -{\mathrm{Tr}}\\{Z_{\rho}\\}.$ (A.1)

Then the derivative of ${\hat{\sigma}}_{\rho}^{2}$ in (2.8) is given by

$n\,\frac{d{\hat{\sigma}}^{2}_{\rho}(y)}{d\rho}\ =\
y^{\top}\Bigl{\\{}-W^{\top}(I-H_{\rho})K_{\rho}-K_{\rho}^{\top}\Bigl{(}\frac{dH_{\rho}}{d\rho}\Bigr{)}K_{\rho}-K_{\rho}^{\top}(I-H_{\rho})W\Bigr{\\}}y.$
(A.2)

To compute $\frac{d}{d\rho}H_{\rho}$, note that, for any family of invertible
matrices $(Y_{x},\,x\in\mathbb{R})$ with differentiable components,
$0=\frac{d}{dx}\\{Y_{x}^{-1}Y_{x}\\}$, implying that

$\frac{d}{dx}\\{Y_{x}^{-1}\\}\ =\
-Y_{x}^{-1}\frac{d}{dx}\\{Y_{x}\\}Y_{x}^{-1}.$ (A.3)

Hence, writing $Y_{\rho}:=((K_{\rho}X)^{\top}K_{\rho}X)$, and using
$\frac{dK_{\rho}}{d\rho}=-W$, we have

$\displaystyle\frac{dY_{\rho}^{-1}}{d\rho}$ $\displaystyle=$ $\displaystyle
Y_{\rho}^{-1}\\{X^{\top}W^{\top}K_{\rho}X+X^{\top}K_{\rho}^{\top}WX\\}Y_{\rho}^{-1}.$

Thus, noting for the last line that $W=Z_{\rho}K_{\rho}$ and that
$K_{\rho}XY_{\rho}^{-1}(K_{\rho}X)^{\top}=H_{\rho}$, we deduce that

$\displaystyle\frac{dH_{\rho}}{d\rho}$ $\displaystyle=$
$\displaystyle\frac{d}{d\rho}\Bigl{\\{}K_{\rho}X\,Y_{\rho}^{-1}(K_{\rho}X)^{\top}\Bigr{\\}}$
(A.4) $\displaystyle=$ $\displaystyle-
WXY_{\rho}^{-1}(K_{\rho}X)^{\top}-K_{\rho}XY_{\rho}^{-1}(WX)^{\top}$
$\displaystyle\qquad\mbox{}+K_{\rho}XY_{\rho}^{-1}\\{X^{\top}W^{\top}K_{\rho}X+X^{\top}K_{\rho}^{\top}WX\\}Y_{\rho}^{-1}(K_{\rho}X)^{\top}$
$\displaystyle=$
$\displaystyle-(I-H_{\rho})Z_{\rho}H_{\rho}-H_{\rho}Z_{\rho}^{\top}(I-H_{\rho}).$

Note also that, because $H_{\rho}$ is idempotent, which implies that
$(I-H_{\rho})H_{\rho}=0$, we have

${\mathrm{Tr}}((Z_{\rho}H_{\rho})^{\top}(I-H_{\rho}))\ =\
{\mathrm{Tr}}((I-H_{\rho})Z_{\rho}H_{\rho})\ =\
{\mathrm{Tr}}(Z_{\rho}H_{\rho}(I-H_{\rho}))\ =\ 0.$ (A.5)

Combining (A.2) and (A.4), it follows that

$n\,\frac{d{\hat{\sigma}}^{2}_{\rho}(y)}{d\rho}\ =\
-(K_{\rho}y)^{\top}(I-H_{\rho})(Z_{\rho}^{\top}+Z_{\rho})(I-H_{\rho})K_{\rho}y,$
(A.6)

and thus, writing $U(\rho;y)\ :=\ \frac{d}{d\rho}F(\rho;y)$ and using (A.1),

$\displaystyle U(\rho;y)$ $\displaystyle=$
$\displaystyle\frac{1}{{\hat{\sigma}}^{2}_{\rho}(y)}\,\frac{d{\hat{\sigma}}^{2}_{\rho}(y)}{d\rho}-\frac{2}{n}\,\frac{d}{d\rho}\,\log\det
K_{\rho}$ (A.7) $\displaystyle=$
$\displaystyle-\frac{1}{n{\hat{\sigma}}^{2}_{\rho}(y)}\,(K_{\rho}y)^{\top}(I-H_{\rho})(Z_{\rho}^{\top}+Z_{\rho})(I-H_{\rho})K_{\rho}y+\frac{2}{n}\sum_{l\geq
1}\rho^{l}{\mathrm{Tr}}\\{W^{l+1}\\},\phantom{XX}$ (A.8)

as given in (4.2).

###  A.2 The magnitude of $nU^{\prime}(\rho;y)$

From (A.7) and (A.8), and since, for any $n\times n$ matrix $M$ and $n$-vector
$y$, $y^{\top}My=y^{\top}M^{\top}y$, it follows that

$\displaystyle U^{\prime}(\rho;y)$ $\displaystyle=$
$\displaystyle\frac{1}{{\hat{\sigma}}^{2}_{\rho}(y)}\,\frac{d^{2}{\hat{\sigma}}^{2}_{\rho}(y)}{d\rho^{2}}-\Bigl{\\{}\frac{1}{{\hat{\sigma}}^{2}_{\rho}(y)}\,\frac{d{\hat{\sigma}}^{2}_{\rho}(y)}{d\rho}\Bigr{\\}}^{2}+\frac{2}{n}\sum_{l\geq
1}l\rho^{l-1}{\mathrm{Tr}}\\{W^{l+1}\\}$ (A.9) $\displaystyle=$
$\displaystyle-\frac{2}{n{\hat{\sigma}}^{2}_{\rho}(y)}\frac{d}{d\rho}\bigl{\\{}y^{\top}K_{\rho}^{\top}\,\\{(I-H_{\rho})Z_{\rho}(I-H_{\rho})\\}\,K_{\rho}y\bigr{\\}}$
$\displaystyle\qquad\mbox{}-\Bigl{\\{}U(\rho;y)-\frac{2}{n}\sum_{l\geq
1}\rho^{l}{\mathrm{Tr}}\\{W^{l+1}\\}\Bigr{\\}}^{2}+\frac{2}{n}\sum_{l\geq
1}l\rho^{l-1}{\mathrm{Tr}}\\{W^{l+1}\\},\phantom{XX}$

invoking (A.6) for the first term. Using
$\frac{d}{d\rho}K_{\rho}=-W{=-Z_{\rho}K_{\rho}}$ and
$\frac{d}{d\rho}Z_{\rho}=W\frac{d}{d\rho}\\{K_{\rho}^{-1}\\}=Z_{\rho}^{2}$,
from (A.3), and again using $y^{\top}My=y^{\top}M^{\top}y$, we can express the
derivative of the quadratic form in (A.9) in terms of the quantities
$H_{\rho}$, $K_{\rho}$, $Z_{\rho}$ and $\frac{d}{d\rho}H_{\rho}$ as

$\displaystyle\frac{d}{d\rho}\\{y^{\top}K_{\rho}^{\top}\,\\{(I-H_{\rho})Z_{\rho}(I-H_{\rho})\\}\,K_{\rho}y\\}$
(A.10) $\displaystyle=$ $\displaystyle
y^{\top}K_{\rho}^{\top}\Bigl{\\{}-Z_{\rho}^{\top}(I-H_{\rho})(Z_{\rho}+Z_{\rho}^{\top})(I-H_{\rho})-\frac{dH_{\rho}}{d\rho}(Z_{\rho}+Z_{\rho}^{\top})(I-H_{\rho})$
$\displaystyle\qquad\mbox{}+(I-H_{\rho})Z_{\rho}^{2}(I-H_{\rho})\Bigr{\\}}K_{\rho}y.$

The middle term in (A.9) is expressed in terms of the same quantities
$H_{\rho}$, $K_{\rho}$, $Z_{\rho}$ and $\frac{d}{d\rho}H_{\rho}$, together
with $W$, using (A.2). The final term in (A.9), much as for (4.6), gives
$2n^{-1}{\mathrm{Tr}}\\{Z_{\rho}^{2}\\}$, and, from (A.4),

$\displaystyle
y^{\top}K_{\rho}^{\top}\frac{dH_{\rho}}{d\rho}(Z_{\rho}+Z_{\rho}^{\top})(I-H_{\rho})K_{\rho}y$
(A.11) $\displaystyle\ =\
-y^{\top}K_{\rho}^{\top}\\{(I-H_{\rho})Z_{\rho}H_{\rho}+H_{\rho}Z_{\rho}^{\top}(I-H_{\rho})\\}(Z_{\rho}+Z_{\rho}^{\top})(I-H_{\rho})K_{\rho}y.$

The resulting expressions for $U^{\prime}(\rho;y)$ are not easy to understand.

It is somewhat easier to consider, as a typical value, the expectation
$n{\mathbb{E}}_{\rho_{0}}U^{\prime}(\rho_{0};\mathbf{y})$ on the model
$M_{\rho_{0}}$. Approximating ${\hat{\sigma}}^{2}_{\rho}(\mathbf{y})$ by
$\sigma^{2}$, the expectation of $n$ times the first and third terms in (A.9)
is approximated by

${\mathbb{E}}_{\rho_{0}}\Bigl{\\{}-\frac{2}{\sigma^{2}}\frac{d}{d\rho}\\{\mathbf{y}^{\top}K_{\rho}^{\top}\,\\{(I-H_{\rho})Z_{\rho}(I-H_{\rho})\\}\,K_{\rho}\mathbf{y}\\}{\Bigr{|}_{\rho=\rho_{0}}}\Bigr{\\}}+2{\mathrm{Tr}}\\{Z_{\rho_{0}}^{2}\\}.$

This expression, using (3.3) together with (A.10) and (A.11) and properties of
the trace, and also observing that $(I-H_{\rho})K_{\rho}X=0$, gives

$\displaystyle
2{\mathrm{Tr}}\\{Z_{\rho_{0}}^{\top}(I-H_{{\rho_{0}}})(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{{\rho_{0}}})$
(A.12)
$\displaystyle\qquad\mbox{}-\\{(I-H_{\rho_{0}})Z_{\rho_{0}}H_{\rho_{0}}+H_{\rho_{0}}Z_{\rho_{0}}^{\top}(I-H_{\rho_{0}})\\}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{{\rho_{0}}})\phantom{XXXX}$
$\displaystyle\qquad\qquad\mbox{}-(I-H_{{\rho_{0}}})Z_{\rho_{0}}^{2}(I-H_{{\rho_{0}}})+Z_{\rho_{0}}^{2}\\}$
$\displaystyle\ =\
{\mathrm{Tr}}\\{(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{{\rho_{0}}})(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{{\rho_{0}}})\\}$
$\displaystyle\qquad\mbox{}+2{\mathrm{Tr}}\\{-Z_{\rho_{0}}H_{\rho_{0}}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{\rho_{0}})+H_{\rho_{0}}Z_{\rho_{0}}^{2}\\}\
=:\ \Delta_{0},$

with the first trace equal to ${\textstyle{\frac{1}{2}}}\tau_{0}^{2}$, as
given in (4.8). The middle term in (A.9), multiplied by $n$, gives

$n^{-1}{\rm
Var}_{\rho_{0}}(nU(\rho;\mathbf{y}))+n^{-1}(n{\mathbb{E}}_{\rho_{0}}U(\rho;\mathbf{y})-2{\mathrm{Tr}}\\{Z_{0}\\})^{2}\
\approx\
n^{-1}(\tau_{0}^{2}+4{\mathrm{Tr}}\\{(I-H_{\rho_{0}})Z_{\rho_{0}}\\}^{2}),$

using (4.7) and (4.8). In the circumstances discussed in Section 4.3,
$\Delta_{0}$ and $\tau_{0}^{2}$ are seen to be typically comparable to
${\mathrm{Tr}}\\{Z_{\rho_{0}}Z_{\rho_{0}}^{\top}\\}$, so that this last
contribution is of relative order $O(n^{-1})$ as compared to $\Delta_{0}$.
Hence, in such circumstances, the expectation
$n{\mathbb{E}}_{\rho_{0}}U^{\prime}(\rho_{0};\mathbf{y})$ on the model
$M_{\rho_{0}}$ is approximately given by $\Delta_{0}$.

###  A.3 Details of the MLE when structure is present

Suppose that $X$ consists solely of the $n$-vector ${\bf{1}}$, having all
elements equal to $1$, and that there are no isolated vertices. Suppose also
that $W$ is chosen to have $\sum_{j=1}^{n}W_{ij}=1$ for all $j$, so that
$W{\bf{1}}={\bf{1}}$. Then $H_{\rho}=H_{*}:=n^{-1}{\bf{1}}{\bf{1}}^{\top}$ for
all $\rho$, because $K_{\rho}X=(I-\rho W){\bf{1}}=(1-\rho){\bf{1}}$, and hence
$WH_{\rho}=H_{*}$ and $Z_{\rho}H_{\rho}=(1-\rho)^{-1}H_{*}$. The matrix
product in (4.8) can be multiplied out to give

$\displaystyle\tau_{0}^{2}$ $\displaystyle=$ $\displaystyle
2{\mathrm{Tr}}\\{(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{\rho_{0}})(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{\rho_{0}})\\}$
$\displaystyle=$ $\displaystyle
2{\mathrm{Tr}}\\{2(Z_{\rho_{0}}^{2}+Z_{\rho_{0}}Z_{\rho_{0}}^{\top})-2H_{\rho_{0}}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})^{2}+[H_{\rho_{0}}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})]^{2}\\}$
$\displaystyle=$ $\displaystyle
4{\mathrm{Tr}}\\{(Z_{\rho_{0}}^{2}+Z_{\rho_{0}}Z_{\rho_{0}}^{\top})(I-H_{\rho_{0}})\\}-2{\mathrm{Tr}}\\{H_{\rho_{0}}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})^{2}-[H_{\rho_{0}}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})]^{2}\\}.$

With the particular choices of $X$ and $W$, and repeatedly using the
properties of the trace, this gives

$\displaystyle\tau_{0}^{2}$ $\displaystyle=$ $\displaystyle
4{\mathrm{Tr}}\\{(Z_{\rho_{0}}^{2}+Z_{\rho_{0}}Z_{\rho_{0}}^{\top})(I-H_{*})\\}+{2}\Bigl{\\{}\frac{1}{(1-\rho_{0})^{2}}-\frac{1}{n}|Z_{\rho_{0}}^{\top}{\bf{1}}|^{2}\Bigr{\\}},$

where, for instance, the final term is just
$-{\mathrm{Tr}}\\{H_{*}Z_{\rho_{0}}Z_{\rho_{0}}^{\top}\\}$, whereas, for
example,

${\mathrm{Tr}}\\{H_{*}Z_{\rho_{0}}^{\top}Z_{\rho_{0}}\\}\ =\
\frac{1}{1-\rho_{0}}{\mathrm{Tr}}\\{H_{*}Z_{\rho_{0}}\\}\ =\
\frac{1}{1-\rho_{0}}{\mathrm{Tr}}\\{Z_{\rho_{0}}H_{*}\\}\ =\
\frac{1}{(1-\rho_{0})^{2}}\,.$

In much the same way, and in particular because
$Z_{\rho}H_{*}=(1-\rho)^{-1}H_{*}$ and $H_{*}(I-H_{*})=0$, we obtain

$2{\mathrm{Tr}}\\{-Z_{\rho_{0}}H_{*}(Z_{\rho_{0}}+Z_{\rho_{0}}^{\top})(I-H_{*})+H_{*}Z_{\rho_{0}}^{2}\\}\
=\ 2{\mathrm{Tr}}\\{Z_{\rho_{0}}^{2}H_{*}\\}\ =\
\frac{2}{(1-\rho_{0})^{2}}\,,$

giving the formula for $\Delta_{0}$ in (4.15).

Note that

$\frac{1}{n}|Z_{\rho_{0}}^{\top}{\bf{1}}|^{2}\ \geq\
\Bigl{\\{}\frac{1}{n}{\bf{1}}^{\top}Z_{\rho_{0}}^{\top}{\bf{1}}\Bigr{\\}}^{2}\
=\ \frac{1}{(1-\rho_{0})^{2}}\,,$

by Cauchy–Schwarz, giving the upper bound for $\tau_{0}^{2}$ in (4.15). Note
also that, again using $Z_{\rho}H_{\rho}=(1-\rho)^{-1}H_{*}$, and because
${\mathrm{Tr}}\\{H_{*}\\}=1$ and $H_{*}^{2}=H_{*}$, it follows that

$\displaystyle{\textstyle{\frac{1}{4}}}\tau_{0}^{2}$ $\displaystyle=$
$\displaystyle{\mathrm{Tr}}\\{Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2}\\}-{\mathrm{Tr}}\\{H_{*}Z_{\rho_{0}}(H_{*}Z_{\rho_{0}})^{\top}\\}-{\mathrm{Tr}}\\{Z_{\rho_{0}}^{2}H_{*}\\}$
(A.14) $\displaystyle\leq$
$\displaystyle{\mathrm{Tr}}\\{Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2}\\}-\frac{1}{(1-\rho_{0})^{2}}\,,$

so that $\tau_{0}^{2}$ is indeed smaller than
$4{\mathrm{Tr}}\\{Z_{\rho_{0}}Z_{\rho_{0}}^{\top}+Z_{\rho_{0}}^{2}\\}$.

Suppose now that there is more structure in the regression, so that $X$
consists of columns in addition to the first column ${\bf{1}}$. If we still
assume that $W{\bf{1}}={\bf{1}}$, the basic message remains the same. Because
the matrix $H_{\rho}$ is idempotent, its spectral decomposition can be written
as

$H_{\rho}\ =\ H_{*}+\sum_{l=2}^{m}x_{\rho}^{(l)}(x_{\rho}^{(l)})^{\top},$
(A.15)

where $H_{*}:=n^{-1}{\bf{1}}{\bf{1}}^{\top}$ is as before, and
$x_{\rho}^{(2)},\ldots,x_{\rho}^{(m)}$ are orthonormal vectors orthogonal to
${\bf{1}}$ that, with ${\bf{1}}$, span the space generated by the columns of
$K_{\rho}X$. If, for instance, it is assumed that each of the vectors
$x_{\rho}^{(l)}$, $2\leq l\leq m$, has components bounded in modulus by
$cn^{-1/2}$, for some fixed $c\geq 1$ — as is the case for
$x^{(1)}=n^{-1/2}{\bf{1}}$ — then the traces of the matrices
$Z_{\rho_{0}}^{2}H_{\rho_{0}}$, $Z_{\rho_{0}}^{\top}Z_{\rho_{0}}H_{\rho_{0}}$,
$(Z_{\rho_{0}}H_{\rho_{0}})^{2}$ and
$Z_{\rho_{0}}^{\top}H_{\rho_{0}}Z_{\rho_{0}}H_{\rho_{0}}$ can all be shown to
be bounded by $m^{2}c^{4}/(1-|\rho|)^{2}$, which is not large if $m$ and $c$
are moderate and $|\rho|$ is not too close to $1$. The argument runs as
follows.

In view of (4.16), and writing $x^{(1)}:=n^{-1/2}{\bf{1}}$, the trace
${\mathrm{Tr}}\\{MH_{\rho}\\}$, for any $n\times n$ matrix $M$, can be written
in the form

${\mathrm{Tr}}\\{MH_{\rho}\\}\ =\
\sum_{l=1}^{m}{\mathrm{Tr}}\\{Mx_{\rho}^{(l)}(x_{\rho}^{(l)})^{\top}\\}\ =\
\sum_{l=1}^{m}(x_{\rho}^{(l)})^{\top}Mx_{\rho}^{(l)}.$ (A.16)

Thus, using (2.6), it follows that, for any $1\leq l,l^{\prime}\leq m$,

$\displaystyle|Z_{\rho}x_{\rho}^{(l)}|$ $\displaystyle=$
$\displaystyle\Bigl{|}\sum_{r\geq 0}W^{r+1}\rho^{r}x_{\rho}^{(l)}\Bigr{|}$
(A.17) $\displaystyle\leq$ $\displaystyle\sum_{r\geq
0}|\rho|^{r}cn^{-1/2}|W^{r+1}{\bf{1}}|\ =\
\frac{c}{n^{1/2}(1-|\rho|)}|{\bf{1}}|\ =\ \frac{c}{1-|\rho|}\,;$

that

$\displaystyle|(x_{\rho}^{(l)})^{\top}Z_{\rho}x_{\rho}^{(l^{\prime})}|$
$\displaystyle=$ $\displaystyle\Bigl{|}\sum_{r\geq
0}(x_{\rho}^{(l)})^{\top}W^{r+1}\rho^{r}x_{\rho}^{(l^{\prime})}\Bigr{|}$
(A.18) $\displaystyle\leq$ $\displaystyle{\Bigl{|}\sum_{r\geq
0}|\rho|^{r}(cn^{-1/2}{\bf{1}}^{\top})W^{r+1}(cn^{-1/2}{\bf{1}})\Bigr{|}}\ =\
c^{2}n^{-1}{{\bf{1}}^{\top}{\bf{1}}}\sum_{r\geq 0}|\rho|^{r}$ $\displaystyle=$
$\displaystyle\frac{c^{2}}{1-|\rho|}\,;$

and that

$\displaystyle|(x_{\rho}^{(l)})^{\top}Z_{\rho}^{2}x_{\rho}^{(l^{\prime})}|$
$\displaystyle=$ $\displaystyle|(x_{\rho}^{(l)})^{\top}W^{2}(1-\rho
W)^{-2}x_{\rho}^{(l^{\prime})}|\ =\ \Bigl{|}\sum_{r\geq
0}(x_{\rho}^{(l)})^{\top}(r+1)W^{r+2}\rho^{r}x_{\rho}^{(l^{\prime})}\Bigr{|}$
(A.19) $\displaystyle\leq$ $\displaystyle
c^{2}n^{-1}{\bf{1}}^{\top}{\bf{1}}\sum_{r\geq 0}(r+1)|\rho|^{r}\ =\
\frac{c^{2}}{(1-|\rho|)^{2}}\,.$

Then, for instance, from (A.16) and (A.19), we have

$\displaystyle{\mathrm{Tr}}\\{Z_{\rho_{0}}^{2}H_{\rho_{0}}\\}$
$\displaystyle=$
$\displaystyle\sum_{l=1}^{m}(x_{\rho}^{(l)})^{\top}Z_{\rho_{0}}^{2}x_{\rho}^{(l)}\
\leq\ \frac{mc^{2}}{(1-|\rho|)^{2}}\,;$

from (A.16) and (A.17), we have

$\displaystyle{\mathrm{Tr}}\\{Z_{\rho_{0}}^{\top}Z_{\rho_{0}}H_{\rho_{0}}\\}$
$\displaystyle=$
$\displaystyle\sum_{l=1}^{m}(x_{\rho}^{(l)})^{\top}Z_{\rho_{0}}^{\top}Z_{\rho_{0}}x_{\rho}^{(l)}\
=\ \sum_{l=1}^{m}|Z_{\rho_{0}}x_{\rho}^{(l)}|^{2}\ \leq\
\frac{mc^{2}}{(1-|\rho|)^{2}}\,;$

and, from (A.16) and (A.18), we have

$\displaystyle{\mathrm{Tr}}\\{Z_{\rho_{0}}^{\top}H_{\rho_{0}}Z_{\rho_{0}}H_{\rho_{0}}\\}$
$\displaystyle=$
$\displaystyle\sum_{l=1}^{m}(x_{\rho}^{(l)})^{\top}Z_{\rho_{0}}^{\top}H_{\rho_{0}}Z_{\rho_{0}}x_{\rho}^{(l)}\
=\
\sum_{l=1}^{m}{\mathrm{Tr}}\\{Z_{\rho_{0}}x_{\rho}^{(l)}(x_{\rho}^{(l)})^{\top}Z_{\rho_{0}}^{\top}H_{\rho_{0}}\\}$
$\displaystyle=$
$\displaystyle\sum_{l=1}^{m}\sum_{l^{\prime}=1}^{m}(x_{\rho}^{(l^{\prime})})^{\top}Z_{\rho_{0}}x_{\rho}^{(l)}(x_{\rho}^{(l)})^{\top}Z_{\rho_{0}}^{\top}x_{\rho}^{(l^{\prime})}\
\leq\ \frac{m^{2}c^{4}}{(1-|\rho|)^{2}}\,.$

By similar calculations, each of the traces listed above is bounded in modulus
by the quantity $m^{2}c^{4}/(1-|\rho|)^{2}$.

###  A.4 Estimation of $\beta$

From (2.8) and (2.11), the distribution of $\hat{\beta}_{\rho}$ on the model
$M_{\rho}$ has covariance matrix $\sigma^{2}(X^{\top}S_{\rho}X)^{-1}$, where
$S_{\rho}:=K_{\rho}^{\top}K_{\rho}$. The variance of a linear combination
$a^{\top}\hat{\beta}_{\rho}$, for a given $m$-dimensional unit vector $a$, is
thus

$\sigma^{2}a^{\top}(X^{\top}S_{\rho}X)^{-1}a\ \leq\
\sigma^{2}\lambda_{\mathrm{max}}\\{(X^{\top}S_{\rho}X)^{-1}\\}\ =\
\sigma^{2}/\lambda_{\mathrm{min}}(X^{\top}S_{\rho}X),$ (A.20)

where $\lambda_{\mathrm{min}}(M)\leq\lambda_{\mathrm{max}}(M)$ denote the
smallest and largest eigenvalues of $M$, both positive if $M$ is positive
definite symmetric. Letting ${\tilde{u}}_{X}$ denote the unit $n$-vector
$Xu/|Xu|$, and using Rayleigh–Ritz twice, we have

$\displaystyle\lambda_{\mathrm{min}}(X^{\top}S_{\rho}X)$ $\displaystyle=$
$\displaystyle\min_{u\in\mathbb{R}^{m}\colon|u|=1}u^{\top}X^{\top}S_{\rho}Xu\
=\
\min_{u\in\mathbb{R}^{m}\colon|u|=1}|Xu|^{2}{\tilde{u}}_{X}^{\top}S_{\rho}{\tilde{u}}_{X}$
(A.21) $\displaystyle\geq$
$\displaystyle\lambda_{\mathrm{min}}(S_{\rho})\min_{u\in\mathbb{R}^{m}\colon|u|=1}|Xu|^{2}\
=\ \lambda_{\mathrm{min}}(X^{\top}X)\lambda_{\mathrm{min}}(S_{\rho}).$

Hence, for fixed $\sigma^{2}$, provided that
$\lambda_{\mathrm{min}}(S_{\rho})$ is not close to zero, (A.20) and (A.21)
imply that accurate estimation of all linear combinations of $\beta$ depends
on having a large value of $\lambda_{\mathrm{min}}(X^{\top}X)$, which is the
usual condition for ordinary least squares. In particular, if the network size
$n$ becomes large, but $\lambda_{\mathrm{min}}(S_{\rho})$ remains uniformly
bounded away from zero, the estimate of $\beta$ is consistent if
$\lambda_{\mathrm{min}}(X^{\top}X)\to\infty$. Example 1 in Section 4.4 shows
that the latter condition alone is not enough to guarantee consistent
estimation of $\beta$.

If the true value of $\rho$ is $\rho_{0}$, the distribution of
$\hat{\beta}_{\rho}(\mathbf{y})$, as determined by substituting
$\mathbf{y}=X\beta+K_{\rho_{0}}^{-1}\boldsymbol{\nu}$ for $y$ into (2.8), is
multivariate normal, with mean $\beta$ and covariance matrix
$\sigma^{2}\Sigma(\rho,\rho_{0})$, where

$\displaystyle\Sigma(\rho,\rho_{0})$ $\displaystyle:=$
$\displaystyle(X^{\top}S_{\rho}X)^{-1}(K_{\rho}X)^{\top}K_{\rho}K_{\rho_{0}}^{-1}(K_{\rho_{0}}^{-1})^{\top}K_{\rho}^{\top}K_{\rho}X(X^{\top}S_{\rho}X)^{-1}$
(A.22) $\displaystyle=$
$\displaystyle(X^{\top}S_{\rho}X)^{-1}X^{\top}S_{\rho}(S_{\rho_{0}})^{-1}S_{\rho}X(X^{\top}S_{\rho}X)^{-1}.$

Writing $\widetilde{X}_{\rho}:=S_{\rho}^{1/2}X(X^{\top}S_{\rho}X)^{-1}$, where
$M^{1/2}$ denotes the non-negative square root of a non-negative definite
matrix $M$, and setting
$S(\rho,\rho_{0}):=S_{\rho}^{1/2}(S_{\rho_{0}})^{-1}S_{\rho}^{1/2}$, (A.22)
gives

$\Sigma(\rho,\rho_{0})\ =\
\widetilde{X}_{\rho}^{\top}S(\rho,\rho_{0})\widetilde{X}_{\rho}.$ (A.23)

Using Rayleigh–Ritz, for any unit $m$-vector $a$, we thus have

$\displaystyle a^{\top}\Sigma(\rho,\rho_{0})a\ =\
(\widetilde{X}_{\rho}a)^{\top}S(\rho,\rho_{0})\widetilde{X}_{\rho}a$
$\displaystyle\leq$
$\displaystyle|\widetilde{X}_{\rho}a|^{2}\lambda_{\mathrm{max}}(S(\rho,\rho_{0}))$
$\displaystyle=$
$\displaystyle\\{a^{\top}(X^{\top}S_{\rho}X)^{-1}a\\}\lambda_{\mathrm{max}}(S(\rho,\rho_{0})).$

Thus, if the true model is $M_{\rho_{0}}$, estimating a linear combination
$a^{\top}\beta$ using $\hat{\beta}_{\rho}$ instead of $\hat{\beta}_{\rho_{0}}$
gives an estimator $a^{\top}\hat{\beta}_{\rho}$ whose variance is at most
$\lambda_{\mathrm{max}}(S(\rho,\rho_{0}))$ times the value, given in (A.20),
that it would have had, were $M_{\rho}$ the true model. Hence, if estimation
on the model $M_{\rho}$, using $\hat{\beta}_{\rho}$, is accurate for any
$\rho$, then estimating using $\hat{\beta}_{\rho}$ on the model $M_{\rho_{0}}$
is also accurate, provided that $\lambda_{\mathrm{max}}(S(\rho,\rho_{0}))$ is
not large.

Note that, if $n(W):=\sqrt{\lambda_{\mathrm{max}}(W^{\top}W)}$ denotes the
operator norm of the matrix $W$, then, by Rayleigh–Ritz, because $|Wy|\leq
n(W)|y|$ for all $y$, the largest eigenvalue of
$S(\rho,0)=I-\rho(W+W^{\top})+\rho^{2}W^{\top}W$ is at least
$(|\rho|n(W)-1)^{2}$, and this is large, for moderate non-zero values of
$\rho$, if $n(W)$ is large. Hence, if $\lambda_{\mathrm{max}}(W^{\top}W)$ is
large, this indicates that there may possibly be problems with estimating
$\beta$.

###  A.5 The scale of precision of estimation using $T^{C}_{\rho}$

As discussed in Section 5, a plausible measure of the scale of precision of
estimation using $T^{C}_{\rho}$ is given by
${\hat{\psi}}^{C}_{\rho_{0}}:={\hat{\tau}}^{C}_{\rho_{0}}/|{\widehat{\Delta}}^{C}_{\rho_{0}}|$,
to be estimated by
${\hat{\tau}}^{C}_{{\hat{\rho}}_{C}}/|{\widehat{\Delta}}^{C}_{{\hat{\rho}}_{C}}|$,
where

$({\hat{\tau}}^{C}_{\rho_{0}})^{2}\ :=\ \sigma^{-4}{\rm
Var}_{\rho_{0}}(U^{C}_{\rho_{0}}(\mathbf{y}))\quad\mbox{and}\quad{\widehat{\Delta}}^{C}_{\rho_{0}}\
:=\
{-}\sigma^{-2}{\mathbb{E}}_{\rho_{0}}\Bigl{\\{}\frac{dU^{C}_{\rho}(\mathbf{y})}{d\rho}\Bigr{|}_{\rho=\rho_{0}}\Bigr{\\}},$

and $U^{C}$ is as defined in (5.8).

To find $({\hat{\tau}}^{C}_{\rho_{0}})^{2}$, note that, from (5.4) and (2.8),
since $\mathbf{y}=X\beta+K_{\rho_{0}}^{-1}\boldsymbol{\nu}$ under
$M_{\rho_{0}}$, we have

$\displaystyle U^{C}_{\rho_{0}}(\mathbf{y})$ $\displaystyle=$
$\displaystyle\boldsymbol{\nu}^{\top}(Q_{C}(\rho_{0},\rho_{0})+n^{-1}(I-H_{\rho_{0}}){\mathrm{Tr}}\\{H_{\rho_{0}}C\\})\boldsymbol{\nu}\
=\ \boldsymbol{\nu}^{\top}{\widehat{Q}}(C,\rho_{0})\boldsymbol{\nu},$

where $Q_{C}(\rho,\rho_{0})$ is defined in (5.5), and

${\widehat{Q}}(C,\rho_{0})\ :=\
(I-H_{\rho_{0}})C(I-H_{\rho_{0}})+n^{-1}(I-H_{\rho_{0}}){\mathrm{Tr}}\\{H_{\rho_{0}}C\\}.$
(A.24)

Hence, replacing $C$ by ${\textstyle{\frac{1}{2}}}(C+C^{\top})$, which does
not change $T^{C}$, and using (3.5), it follows that

$\displaystyle({\hat{\tau}}^{C}_{\rho_{0}})^{2}$ $\displaystyle=$
$\displaystyle\sigma^{-4}{\rm Var}_{\rho_{0}}(U^{C}_{\rho_{0}}(\mathbf{y}))\
=\
{\textstyle{\frac{1}{2}}}{\mathrm{Tr}}\bigl{\\{}\\{{\widehat{Q}}(C+C^{\top},\rho_{0})\\}^{2}\bigr{\\}}.$
(A.25)

In view of the factor $n^{-1}$ in the second element of
${\widehat{Q}}(C,\rho_{0})$, the dominant contribution is that from

${\textstyle{\frac{1}{2}}}{\mathrm{Tr}}\bigl{\\{}\\{Q_{C}(\rho_{0},\rho_{0})\\}^{2}\bigr{\\}}\
=\ {\mathrm{Tr}}\\{(I-H_{\rho_{0}})C(I-H_{\rho_{0}})(C+C^{\top})\\},$

once more of the form ${\mathrm{Tr}}\\{Z^{2}+ZZ^{\top}\\}$, now with
$Z:=(I-H_{\rho_{0}})C(I-H_{\rho_{0}})$.

To find ${\widehat{\Delta}}^{C}_{\rho_{0}}$, observe that, since
$\frac{d}{d\rho}K_{\rho}=-W=-Z_{\rho}K_{\rho}$ and $\frac{d}{d\rho}H_{\rho}$
is as in (A.4), we have

$\displaystyle\frac{dT^{C}_{\rho}(\mathbf{y})}{d\rho}$ $\displaystyle=$
$\displaystyle\frac{d}{d\rho}\\{(K_{\rho}\mathbf{y})^{\top}(I-H_{\rho})C(I-H_{\rho})K_{\rho}\mathbf{y}\\}$
(A.26) $\displaystyle=$
$\displaystyle-(K_{\rho}\mathbf{y})^{\top}\Bigl{(}Z_{\rho}^{\top}(I-H_{\rho})+\frac{dH_{\rho}}{d\rho}\Bigr{)}(C+C^{\top})(I-H_{\rho})K_{\rho}\mathbf{y}\phantom{XXXXXXXXXX}$
$\displaystyle=$
$\displaystyle-({X\beta+}K_{\rho_{0}}^{-1}\boldsymbol{\nu})^{\top}K_{\rho}^{\top}\bigl{\\{}Z_{\rho}^{\top}(I-H_{\rho})-(I-H_{\rho})Z_{\rho}H_{\rho}$
$\displaystyle\qquad\qquad\qquad\qquad\mbox{}-H_{\rho}Z_{\rho}^{\top}(I-H_{\rho})\bigr{\\}}(C+C^{\top})(I-H_{\rho})K_{\rho}K_{\rho_{0}}^{-1}\boldsymbol{\nu},$

where $X\beta$ is not present at the second appearance of $\boldsymbol{\nu}$
because $(I-H_{\rho})K_{\rho}X=0$. From (A.26) and (3.3), and because
$H_{\rho}(I-H_{\rho})=0$, it thus follows that

$-{\mathbb{E}}_{\rho_{0}}\Bigl{\\{}\frac{dT^{C}_{\rho}(\mathbf{y})}{d\rho}\Big{|}_{\rho=\rho_{0}}\Bigr{\\}}\
=\
\sigma^{2}{\mathrm{Tr}}\bigl{\\{}(I-H_{\rho_{0}}){[Z_{\rho_{0}}^{\top}(I-H_{\rho_{0}})-Z_{\rho_{0}}H_{\rho_{0}}]}(C+C^{\top})\bigr{\\}}\
=:\ \sigma^{2}{\widehat{\Delta}}^{C}_{\rho_{0},1},$ (A.27)

say, with $X\beta$ disappearing because ${\mathbb{E}}\boldsymbol{\nu}=0$.
Then, from (A.4),

$-\frac{d}{d\rho}{\mathrm{Tr}}\\{H_{\rho}C\\}\ =\
-{\mathrm{Tr}}\Bigl{\\{}C\,\frac{dH_{\rho}}{d\rho}\Bigr{\\}}\ =\
{\mathrm{Tr}}\bigl{\\{}C\\{(I-H_{\rho})Z_{\rho}H_{\rho}+H_{\rho}Z_{\rho}^{\top}(I-H_{\rho})\\}\bigr{\\}}\
=:\ {\widehat{\Delta}}^{C}_{\rho_{0},2},$ (A.28)

say, and, from (A.6), and on $M_{\rho_{0}}$,

$-\frac{d{\hat{\sigma}}^{2}_{\rho_{0}}(\mathbf{y})}{d\rho}\ =\
n^{-1}\boldsymbol{\nu}^{\top}(I-H_{\rho_{0}})(Z_{\rho_{0}}^{\top}+Z_{\rho_{0}})(I-H_{\rho_{0}})\boldsymbol{\nu},$

so that

$-\sigma^{-2}{\mathbb{E}}_{\rho_{0}}\Bigl{\\{}\frac{d{\hat{\sigma}}^{2}_{\rho_{0}}(\mathbf{y})}{d\rho}\Bigr{\\}}{\mathrm{Tr}}\\{H_{\rho_{0}}C\\}\
=\
2n^{-1}{\mathrm{Tr}}\\{(I-H_{\rho_{0}})Z_{\rho_{0}}(I-H_{\rho_{0}})\\}{\mathrm{Tr}}\\{H_{\rho_{0}}C\\}\
=:\ {\widehat{\Delta}}^{C}_{\rho_{0},3},$ (A.29)

say. Then

${\widehat{\Delta}}^{C}_{\rho_{0}}\ :=\
\sum_{l=1}^{3}{\widehat{\Delta}}^{C}_{\rho_{0},l},$ (A.30)

where ${\widehat{\Delta}}^{C}_{\rho_{0},l}$, $1\leq l\leq 3$, are given by
(A.27), (A.28) and (A.29).

###  A.6 Details of network effects estimation

In the network effects model, except when only the overall mean is to be
estimated, and $W{\bf{1}}={\bf{1}}$ by construction, it is usual to expect
that

$WX\ =\ HWX+(I-H)WX\ =\ X\Gamma+E,$ (A.31)

for a non-zero $n\times m$ matrix $E$. Here, $H:=X(X^{\top}X)^{-1}X^{\top}$ is
the projection onto the space spanned by the columns of $X$; $HWX$ is the
corresponding projection of the columns of $WX$, and can thus be written in
the form $X\Gamma$; and $E:=(I-H)WX$ has columns orthogonal to the columns of
$X$. Then, as in the calculation leading to (6.6), we obtain

$\displaystyle K_{\rho}^{-1}X\beta$ $\displaystyle=$ $\displaystyle
X(I-\rho\Gamma)^{-1}\beta+\rho K_{\rho}^{-1}E\beta$ $\displaystyle=$
$\displaystyle X\gamma+\rho K_{\rho}^{-1}E(I-\rho\Gamma)\gamma\ =:\
X_{\rho}\gamma,$

say, parametrized in terms of $\gamma$. The log-likelihood, in terms of
$\gamma$, becomes

$\ell(\gamma,\sigma^{2},\rho;y)\ =\
-{\textstyle{\frac{n}{2}}}\\{\log(\sigma^{2})+\log(2\pi)\\}-\frac{1}{2\sigma^{2}}\\{K_{\rho}(y-X_{\rho}\gamma)\\}^{\top}\\{K_{\rho}(y-X_{\rho}\gamma)\\}+\log\det
K_{\rho},$

and differentiating with respect to $\rho$ gives the analogue of the score
function in (6.4) to be

$\displaystyle{\tilde{s}}^{(\rho)}(\gamma,\sigma^{2},\rho;y)$ $\displaystyle=$
$\displaystyle\frac{1}{\sigma^{2}}\bigl{(}\\{K_{\rho}(y-X_{\rho}\gamma)\\}^{\top}Z_{\rho}\\{K_{\rho}(y-X_{\rho}\gamma)\\}\bigr{)}$
(A.32)
$\displaystyle\qquad\mbox{}+\frac{1}{\sigma^{2}}\\{K_{\rho}(y-X_{\rho}\gamma)\\}^{\top}G_{\rho}\gamma-{\mathrm{Tr}}\\{Z_{\rho}\\},$

where

$\displaystyle G_{\rho}$ $\displaystyle:=$ $\displaystyle
K_{\rho}\frac{d}{d\rho}\\{\rho K_{\rho}^{-1}E(I-\rho\Gamma)\\}$ (A.33)
$\displaystyle=$ $\displaystyle E(I-\rho\Gamma)-\rho
K_{\rho}^{-1}WE(I-\rho\Gamma)-\rho E\Gamma.$

The variance of ${\tilde{s}}^{(\rho)}(\gamma_{0},\sigma_{0}^{2},\rho_{0};y)$
on the model with parameters $(\gamma_{0},\sigma_{0}^{2},\rho_{0})$ is now

$V_{2}\ :=\
{\mathrm{Tr}}\\{Z_{\rho_{0}}^{2}+Z_{\rho_{0}}Z_{\rho_{0}}^{\top}\\}+\frac{1}{\sigma_{0}^{2}}|G_{\rho_{0}}\gamma|^{2},$
(A.34)

differing from (6.5) in the second term. The corresponding Cramér–Rao lower
bound for the variance of any unbiased estimator of $\rho$ is thus at least as
big as $1/V_{2}$.

Just as for the network disturbance model, the first term in (A.34) may not be
large in dense networks. To get an idea of the magnitude of the second term,
consider the case where $\rho_{0}=0$, so that $G_{0}=E$ and
$|G_{0}\gamma|^{2}=\gamma^{\top}E^{\top}E\gamma$. Then

$\frac{1}{\sigma_{0}^{2}}|G_{0}\gamma|^{2}\ \leq\
\frac{|\gamma|^{2}}{\sigma_{0}^{2}}\lambda_{\mathrm{max}}(E^{\top}E),$

where $\lambda_{\mathrm{max}}(E^{\top}E)$ denotes the largest eigenvalue of
the positive definite symmetric $m\times m$ matrix $E^{\top}E$. For instance,
if the elements of $E$, that represent to what extent the columns of $WX$ do
not belong to the linear span of the columns of $X$, are uniformly of
magnitude at most $c/\sqrt{n}$, then the elements of $E^{\top}E$ are at most
of magnitude $c^{2}$, and thus $\lambda_{\mathrm{max}}(E^{\top}E)$ is not
large if $c$ is not large. In such circumstances, any reasonable estimator of
$\rho$ has to have substantial variability. In practice, the model
decomposition (A.31) can be applied to yield the value of

$\lambda_{\mathrm{max}}\\{E^{\top}E\\}\ =\
\lambda_{\mathrm{max}}\\{(WX)^{\top}(I-H)WX\\},$

and, if neither this nor ${\mathrm{Tr}}\\{W^{2}+WW^{\top}\\}$ is large, then
neither $\rho$ nor $\beta$ can be accurately estimated.