The SAR Model for Very Large Datasets: A Reduced Rank Approach

Abstract

The SAR model is widely used in spatial econometrics to model Gaussian processes on a discrete spatial lattice, but for large datasets, fitting it becomes computationally prohibitive, and hence, its usefulness can be limited. A computationally-efficient spatial model is the spatial random effects (SRE) model, and in this article, we calibrate it to the SAR model of interest using a generalisation of the Moran operator that allows for heteroskedasticity and an asymmetric SAR spatial dependence matrix. In general, spatial data have a measurement-error component, which we model, and we use restricted maximum likelihood to estimate the SRE model covariance parameters; its required computational time is only the order of the size of the dataset. Our implementation is demonstrated using mean usual weekly income data from the 2011 Australian Census.

1. Introduction

The spatial autoregressive (SAR) model, originally proposed by Whittle [1] (see also, [2]), is widely used in spatial econometrics to model Gaussian processes on a discrete spatial lattice (which may be irregular). It has an autoregressive structure that models the spatial dependence of the attributes through a precision matrix that is typically a function of the proximity between elements of the lattice. This spatial dependence matrix is generally sparse. An important application of the SAR model is when the lattice is made up of n spatially-contiguous areas, such that for two areas that are “far apart”, the corresponding entry in the spatial dependence matrix is zero.

Despite its popularity, the SAR model has some limitations for the analysis of large datasets, because maximum likelihood estimation of its parameters is not scalable in the size of the spatial dataset. This limits the size of the problems for which SAR models can be effectively used. Algorithmic approaches to reduce the computational burden include sparse-matrix techniques [3,4] and approximate methods, such as Taylor series and Chebyshev approximations [4,5].

Rather than modelling spatial dependence through a precision matrix that is a function of the spatial proximity of areas, it could be modelled directly through a covariance matrix; a reduced-rank covariance model could then be used to lower the order of computations (e.g., [6]) to $O\left(n\right)$. Reduced-rank approaches include latent-factor-analysis models [7], predictive-process models for Gaussian data [8,9] and non-Gaussian data [10], spatial random effects (SRE) models [11,12] and spatial generalised linear mixed models (SGLMMs) [13,14,15]. The advantage of the SRE approach is a substantial improvement in computational efficiency for a broad class of flexible, spatial, non-stationary covariance functions [12,16]. The spatial dependence is captured through the choice of basis functions and through an $r\times r$ covariance matrix of the coefficients of the basis functions. By reducing the rank of the spatial covariance function from n to r through the use of r spatial basis functions, the SRE model has a computational cost of $O\left(n{r}^2\right)$, which is $O\left(n\right)$ for $r\ll n$ [12]. Our challenge in this paper is to calibrate and fit a computationally-efficient SRE model so that its parameters are directly interpretable in terms of the familiar SAR model.

An important aspect of model definition is specifying appropriate basis functions. Several classes of basis functions have previously been used for SRE models, including empirical orthogonal functions (EOFs) (summarised in, e.g., [6,17]), Fourier basis functions [18], smoothing splines [6], radial basis functions (e.g., [19,20], pp. 186–187), and multi-resolutional basis functions (i.e., basis functions obtained from functions defined at multiple resolutions), such as bisquare spatial basis functions [12,21] and wavelets [16,22], that capture the spatial dependence at multiple scales. Selecting an appropriate class of basis functions is generally problem dependent, although the use of the conditional Akaike information criterion and the generalised-degrees-of-freedom criterion have been considered [23]. The spatial basis functions selected for the model need not be orthogonal, but they are recommended to be multi-resolutional and fast to evaluate [12]. More recently, a spectral decomposition of the matrix in the numerator of Moran’s I statistic, called the Moran operator [14], has been used to identify a set of orthogonal eigenvectors associated with dominant patterns of positive spatial dependence (e.g., [24,25,26]). These eigenvectors have been used to obtain spatial basis functions for SGLMMs [27], and substantial computational efficiency has been achieved for both SGLMMs [14] and SRE models [28] by including only the first $r\ll n$ basis functions in the model.

Moran’s I can be motivated from a simple SAR model with homoskedastic errors and a symmetric spatial dependence matrix given by nearest-neighbour adjacencies. In this article, we consider a heteroskedastic SAR model with an asymmetric spatial dependence matrix and define a generalised Moran operator (GMO) whose eigenvectors define our spatial basis functions. In particular, we propose to model spatial data for small areas using an SRE model with these spatial basis functions calibrated from the proposed SAR model. The SAR-calibrated SRE model, for which computational efficiencies are readily available, is well suited to very large datasets, and its parameters are interpretable in terms of the original SAR parameters.

The plan for the rest of the article is as follows. Section 2 outlines the main features of the SAR and SRE models. A motivation is given in Section 3 for spatial basis functions obtained from a generalised Moran operator. Dimension reduction using these basis functions is developed in Section 4. Parameter estimation under the SAR-calibrated SRE model is presented in Section 5. An example of their use for Australian Census data is presented in Section 6, and the paper concludes with a discussion and conclusions in Section 7.

2. SAR Model and SRE Model Specifications

Consider a set of n geographic areas $\{{\mathit{A}}_{\mathit{i}}:\mathit{i}=1,...,\mathit{n}\}$ in a spatial domain, and identify the location of a well-defined feature in each area, denoted ${\mathbf{s}}_1,...,{\mathbf{s}}_{\mathit{n}}$, respectively. Define an irregular, finite spatial lattice ${\mathit{D}}_{\mathit{s}}\equiv \{{\mathbf{s}}_1,...,{\mathbf{s}}_{\mathit{n}}\}$ over the region of interest, ${\bigcup }_{\mathit{i}=1}^{\mathit{n}}{A}_{\mathit{i}}$, and consider the attributes $\{\mathit{Y}\left({\mathbf{s}}_{\mathit{i}}\right):\mathit{i}=1,...,\mathit{n}\}$ associated with the $\mathit{n}$ areas.

We are interested in making inference on a spatial process $Y(·)$, which we henceforth assume is Gaussian. However, $Y(·)$ is a latent process: rather than observing $\{\mathit{Y}\left({\mathbf{s}}_{\mathit{i}}\right):\mathit{i}=1,...,\mathit{n}\}$, we observe data $\{\mathit{Z}\left({\mathbf{s}}_{\mathit{i}}\right):\mathit{i}=1,...,\mathit{n}\}$, which offer an imperfect view of the latent process due to limitations in the measurement process (including averaging, the use of proxy measures rather than direct measurement and observation error).

In this section, we specify hierarchical models that allow us to make inference on the latent process $\mathbf{Y}\equiv {(\mathit{Y}\left({\mathbf{s}}_1\right),...,\mathit{Y}\left({\mathbf{s}}_{\mathit{n}}\right))}^{\prime }$ using the observed data $\mathbf{Z}\equiv {(\mathit{Z}\left({\mathbf{s}}_1\right),...,\mathit{Z}\left({\mathbf{s}}_{\mathit{n}}\right))}^{\prime }$. We first define a SAR model for $\mathbf{Y}$ and then a measurement-error model for $\mathbf{Z}$ given $\mathbf{Y}$. Subsequently, we consider an SRE model for $\mathbf{Y}$ and discuss the implications of equating the SAR and SRE models’ respective means, variances and covariances.

2.1. SAR Model

A SAR model for $\mathbf{Y}$ over the lattice $D_s$, which includes a spatially-varying mean, is defined as:

$$\begin{array}{cc}\hfill Y\left({\mathbf{s}}_i\right)& =\text{μ}\left({\mathbf{s}}_i\right)+\sum _{j=1}^n{b}_{ij}(Y\left({\mathbf{s}}_j\right)-\text{μ}\left({\mathbf{s}}_j\right))+\text{ν}\left({\mathbf{s}}_i\right),i=\text{1},...,n,\hfill \end{array}$$

(1)

where $\mathbf{\mu }\equiv {(\text{μ}\left({\mathbf{s}}_{\text{1}}\right),...,\text{μ}\left({\mathbf{s}}_n\right))}^{\prime }$ is an n-dimensional vector representing the mean of $\mathbf{Y}$ and $b_{ij}$ is a measure of the spatial dependence between the locations ${\mathbf{s}}_i$ and ${\mathbf{s}}_j$ with the convention that $b_{ii}=0$; $i=\text{1},...,n$. We define the SAR’s $n\times n$ spatial dependence matrix as $\mathbf{B}\equiv \{b_{ij}:i,j=1,...,n\}$. In general, $b_{ij}\ne b_{ji}$; that is, $\mathbf{B}$ is an asymmetric square matrix. The elements of the error vector $\mathbf{\text{ν}}\equiv {(\text{ν}\left({\mathbf{s}}_1\right),...,\text{ν}\left({\mathbf{s}}_n\right))}^{\prime }$ are independent Gaussian random variables with mean zero and variance ${\text{σ}}_{\text{ν}i}^2$, $i=\text{1},...,n$. Their joint distribution can be written as $\mathbf{\text{ν}}\sim \text{Gau}(\mathbf{0},{\mathbf{\Sigma }}_{\text{ν}})$, where ${\mathbf{\Sigma }}_{\text{ν}}\equiv \text{diag}({\text{σ}}_{\text{ν}1}^2,...,{\text{σ}}_{\text{ν}n}^2)$. For many applications, ${\mathbf{\Sigma }}_{\text{ν}}$ is parametrised as ${\text{σ}}_{\text{ν}i}^2={\text{σ}}_{\text{ν}}^2{V}_{\text{ν}}\left({\mathbf{s}}_i\right)$, where ${\text{σ}}_{\text{ν}}^2>0$, and $V_{\text{ν}}(·)>0$ is a known function of location. We write ${\mathbf{\Sigma }}_{\text{ν}}\equiv {\text{σ}}_{\text{ν}}^2{\mathbf{V}}_{\text{ν}}$, where the diagonal matrix ${\mathbf{V}}_{\text{ν}}\equiv \text{diag}(V_{\text{ν}}\left({\mathbf{s}}_1\right),...,V_{\text{ν}}\left({\mathbf{s}}_n\right))$ may be the identity or may model heteroskedasticity (e.g., the inverse of the population size in each area; see Section 6).

Here, we model the mean as a linear function of p known covariates; that is $\mathbf{\mu }\equiv \mathbf{X}\mathbf{\beta }$, where $\mathbf{X}\equiv {(\mathbf{X}\left({\mathbf{s}}_{\text{1}}\right),...,\mathbf{X}\left({\mathbf{s}}_n\right))}^{\prime }$, for $\mathbf{X}\left({\mathbf{s}}_i\right)\equiv {(X_{\text{1}}\left({\mathbf{s}}_i\right),...,X_p\left({\mathbf{s}}_i\right))}^{\prime }$, and $\mathbf{\beta }\equiv {({\beta }_{\text{1}},...,{\beta }_p)}^{\prime }$ is the associated p-dimensional vector of regression coefficients. In this case, (1) can be written equivalently as,

$$\begin{array}{cc}\hfill \mathbf{Y}& \equiv \mathbf{X}\mathbf{\beta }+{(\mathbf{I}-\mathbf{B})}^{-1}\mathbf{\nu }.\hfill \end{array}$$

(2)

The spatial dependence matrix $\mathbf{B}$ has zeros down the diagonal; generally it is asymmetric, and it is often defined in terms of the proximity of the lattice elements. For (2) to be well defined, we consider only those $\mathbf{B}$ for which $(\mathbf{I}-\mathbf{B})$ is invertible. Generally, $\mathbf{B}$ (and hence, $(\mathbf{I}-\mathbf{B})$) is sparse. In many applications, $\mathbf{B}$ is parametrised as $\text{ρ}\mathbf{H}$, where $\mathbf{H}\equiv \{h_{ij}:i,j=1,...,n\}$ is a known, in general asymmetric, $n\times n$ spatial-weights matrix with zeros down the diagonal, and ρ is a spatially-dependent parameter. The case of $\text{ρ}=0$ corresponds to “no spatial dependence,” since $\mathbf{B}=\text{ρ}\mathbf{H}=\mathbf{0}$.

Henceforth, we parametrise our model as $\mathbf{B}\equiv \text{ρ}\mathbf{H}$ and ${\mathbf{\Sigma }}_{\text{ν}}\equiv {\text{σ}}_{\text{ν}}^2{\mathbf{V}}_{\text{ν}}$, where $\mathbf{H}$ and ${\mathbf{V}}_{\text{ν}}$ are known. Consequently, there are two covariance parameters, ρ and ${\text{σ}}_{\text{ν}}^2$. From (2), the mean vector and the covariance matrix of $\mathbf{Y}$ are, respectively, given by,

$$\begin{array}{cc}\hfill {\text{E}}_{SAR}\left(\mathbf{Y}\right)& \equiv \mathbf{X}\mathbf{\beta }\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\text{cov}}_{SAR}\left(\mathbf{Y}\right)& \equiv {\text{σ}}_{\text{ν}}^2{(\mathbf{I}-\text{ρ}\mathbf{H})}^{-1}{\mathbf{V}}_{\text{ν}}{(\mathbf{I}-\text{ρ}{\mathbf{H}}^{\prime })}^{-1}.\hfill \end{array}$$

(3)

A statistical model for the observations, $\mathbf{Z}$, is defined in terms of the latent process, $\mathbf{Y}$, by assuming a conditional Gaussian distribution for the observed data. The resulting data model can be written in terms of additive measurement error, as follows:

$$Z\left({\mathbf{s}}_i\right)=Y\left({\mathbf{s}}_i\right)+\epsilon \left({\mathbf{s}}_i\right),i=1,...,n.$$

(4)

The term $\epsilon (·)$ in (4) represents an uncorrelated Gaussian process with mean zero and individual variances given by $\text{var}\left(\epsilon \left({\mathbf{s}}_i\right)\right)={\text{σ}}_{\epsilon }^2{V}_{\epsilon }\left({\mathbf{s}}_i\right)$, for $i=1,...,n$, where ${\text{σ}}_{\epsilon }^2>0$, and $V_{\epsilon }(·)>0$ is a known function of location; $\epsilon (·)$ is also independent of $Y(·)$. Hence,

$$\mathbf{Z}|\mathbf{Y}\sim \text{ind.Gau}(\mathbf{Y},{\text{σ}}_{\epsilon }^2{\mathbf{V}}_{\epsilon }),$$

where “Gau” is the Gaussian (or normal) distribution, and ${\mathbf{V}}_{\epsilon }\equiv \text{diag}(V_{\epsilon }\left({\mathbf{s}}_1\right),...,V_{\epsilon }\left({\mathbf{s}}_n\right))$ may be the identity matrix or it may represent heterogeneous conditions associated with the measurement process. Together, (2) and (4) yield the following model for $\mathbf{Z}$,

$$\mathbf{Z}\equiv \mathbf{X}\mathbf{\beta }+{(\mathbf{I}-\text{ρ}\mathbf{H})}^{-1}\mathbf{\nu }+\mathbf{\epsilon },$$

(5)

where $\mathbf{\epsilon }\equiv {(\epsilon \left({\mathbf{s}}_1\right),...,\epsilon \left({\mathbf{s}}_n\right))}^{\prime }\sim \text{Gau}(\mathbf{0},{\text{σ}}_{\epsilon }^2{\mathbf{V}}_{\epsilon })$, independent of $\mathbf{\nu }\sim \text{Gau}(\mathbf{0},{\text{σ}}_{\text{ν}}^2{\mathbf{V}}_{\text{ν}})$. Traditional spatial econometrics has fitted a SAR model directly to the spatial data, $\mathbf{Z}$, as follows (e.g., [2]): $\mathbf{Z}\equiv \mathbf{X}\mathbf{\beta }+\text{ρ}\mathbf{H}\mathbf{Z}+\mathbf{\nu }$. We believe that (5) presents the three components of variability, namely trend, SAR spatial dependence and measurement error, in an unambiguous and interpretable manner.

Motivated by Besag et al. [29], we have assumed a SAR model for the spatial errors, $\mathbf{Y}-\mathbf{X}\mathbf{\beta }$, to avoid confounding $\text{ρ}\mathbf{H}$ with β. Furthermore, to recognise that the measurement is never perfect, we have added measurement error to the latent process $\mathbf{Y}$. Details regarding these models and their estimation can be found in, for example, Anselin [2] Cliff and Ord [30], Cressie [31], Le Sage [32], Anselin [33], as well as LeSage and Pace [34].

2.2. SRE Model

The SRE model is a reduced-rank spatial statistical model that can be written as follows:

$$Y\left({\mathbf{s}}_i\right)=\mathbf{X}{\left({\mathbf{s}}_i\right)}^{\prime }\mathbf{\beta }+\mathbf{S}{\left({\mathbf{s}}_i\right)}^{\prime }\mathbf{\eta }+\xi \left({\mathbf{s}}_i\right),i=\text{1},...,n,$$

(6)

where $\mathbf{S}(·)\equiv {(S_{\text{1}}(·),...,S_r(·))}^{\prime }$ is a known vector of $r\ll n$ spatial basis functions and $\mathbf{\eta }\equiv {({\eta }_{\text{1}},...,{\eta }_r)}^{\prime }$ is its corresponding r-dimensional vector of random effects. The random effects are assumed to have a multivariate Gaussian distribution with mean $\mathbf{0}$ and $r\times r$ covariance matrix $\mathbf{K}$, so that $\mathbf{\eta }\sim \text{Gau}(\mathbf{0},\mathbf{K})$. The last term in (6), $\xi (·)$, represents fine-scale variation, independent of $\mathbf{Y}$, which is assumed here to follow,

$$\xi \left({\mathbf{s}}_i\right)\sim \text{ind.}\text{Gau}(0,{\text{σ}}_{\xi }^2{V}_{\xi }\left({\mathbf{s}}_i\right)),i=\text{1},...,n,$$

where ${\text{σ}}_{\xi }^2>0$, and $V_{\xi }(·)>0$ is a known function of location. Hence, $\mathbf{\xi }\sim \text{Gau}(\mathbf{0},{\text{σ}}_{\xi }^2{\mathbf{V}}_{\xi })$, where $\mathbf{\xi }\equiv {(\xi \left({\mathbf{s}}_{\text{1}}\right),...,\xi \left({\mathbf{s}}_n\right))}^{\prime }$, and the diagonal matrix ${\mathbf{V}}_{\xi }\equiv \text{diag}(V_{\xi }\left({\mathbf{s}}_1\right),...,V_{\xi }\left({\mathbf{s}}_n\right))$ may be the identity matrix or may represent heterogeneity, for example through ${\left(\text{population}\text{size}\right)}^{-1}$ in each area; see Section 6. Using vector notation, (6) can also be written as,

$$\mathbf{Y}\equiv \mathbf{X}\mathbf{\beta }+\mathbf{S}\mathbf{\eta }+\mathbf{\xi },$$

(7)

where $\mathbf{S}\equiv {(\mathbf{S}\left({\mathbf{s}}_{\text{1}}\right),...,\mathbf{S}\left({\mathbf{s}}_n\right))}^{\prime }$. The mean vector and the covariance matrix of $\mathbf{Y}$ for the SRE model are, respectively, given by,

$$\begin{array}{cc}\hfill {\text{E}}_{SRE}\left(\mathbf{Y}\right)& \equiv \mathbf{X}\mathbf{\beta }\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\text{cov}}_{SRE}\left(\mathbf{Y}\right)& \equiv {\mathbf{SKS}}^{\prime }+{\text{σ}}_{\xi }^2{\mathbf{V}}_{\xi }.\hfill \end{array}$$

(8)

As was done for the SAR model, we define $\mathbf{Z}$ in terms of the latent process $\mathbf{Y}$ through the data model, (4). Together, (7) and (4) yield the following model for $\mathbf{Z}$:

$$\mathbf{Z}\equiv \mathbf{X}\mathbf{\beta }+\mathbf{S}\mathbf{\eta }+\mathbf{\xi }+\mathbf{\epsilon }.$$

(9)

We calibrate the specification of $\mathbf{S}$, $\text{cov}\left(\mathbf{\eta }\right)$ and $\text{cov}\left(\mathbf{\xi }\right)$ in the SRE model to the spatial dependence in the SAR model, by initially matching their first two moments. Now, the mean of the SRE model (7) is the same as the mean of the SAR model (2). Then, after setting the measurement-error terms in (5) and (9) to be the same, we set ${\text{cov}}_{SRE}\left(\mathbf{Y}\right)\equiv {\text{cov}}_{SAR}\left(\mathbf{Y}\right)$ (e.g., [35]). That is, we set:

$${\mathbf{SKS}}^{\prime }+{\text{σ}}_{\xi }^2{\mathbf{V}}_{\xi }={\text{σ}}_{\text{ν}}^2{(\mathbf{I}-\rho \mathbf{H})}^{-1}{\mathbf{V}}_{\text{ν}}{(\mathbf{I}-\rho {\mathbf{H}}^{\prime })}^{-1}.$$

(10)

Notice that in the SRE model (left-hand side), the contribution to spatial dependence is from the reduced-rank matrix ${\mathbf{SKS}}^{\prime }$, whilst in the SAR model (right-hand side), the contribution to spatial dependence comes from the full-rank matrix, $\mathbf{I}-\rho \mathbf{H}$. When ${\mathbf{SKS}}^{\prime }=\mathbf{0}=\text{ρ}\mathbf{H}$, we respectively obtain an SRE model and an SAR model with no spatial dependence. In this case, the SRE model and the SAR model are equivalent if:

$${\text{σ}}_{\xi }^2{\mathbf{V}}_{\xi }={\text{σ}}_{\text{ν}}^2{\mathbf{V}}_{\text{ν}}.$$

(11)

This says that heteroskedasticity in the SRE’s fine-scale variation component ξ should reflect the heteroskedasticity in the SAR model’s ν, when there is no spatial dependence.

When spatial dependence is present in the SAR model, the spatial dependence component of the SRE model should, according to (10), be equal to:

$${\text{σ}}_{\text{ν}}^2{(\mathbf{I}-\text{ρ}\mathbf{H})}^{-1}{\mathbf{V}}_{\text{ν}}{(\mathbf{I}-\text{ρ}{\mathbf{H}}^{\prime })}^{-1}-{\text{σ}}_{\xi }^2{\mathbf{V}}_{\xi }.$$

(12)

Critically, (12) should also be positive-definite with rank $r\ll n$. Ensuring that these restrictions hold leads to our new methodology, which is presented in Section 4.

3. Motivation for the Spatial Basis Functions

We use this section to motivate the construction of our SAR-calibrated spatial basis functions for use in the SRE model. The motivation is clearest for the case where ${\text{σ}}_{\epsilon }^2=0$, and hence, in this section, we temporarily assume $\mathbf{Z}=\mathbf{Y}$. The SAR model parameters are $(\mathbf{\beta },\text{ρ},{\text{σ}}_{\text{ν}}^2)$. Minus twice the log-likelihood function for the SAR model is:

$$\begin{array}{cc}\hfill -2\mathcal{L}& (\mathbf{\beta },\text{ρ},{\text{σ}}_{\text{ν}}^2;\mathbf{Y})\hfill \\ & =nlog\left(2\pi \right)-log\left(\left|{\text{σ}}_{\text{ν}}^{-2}(\mathbf{I}-\text{ρ}{\mathbf{H}}^{\prime }){\mathbf{V}}_{\text{ν}}^{-1}(\mathbf{I}-\text{ρ}\mathbf{H})\right|\right)+{\text{σ}}_{\text{ν}}^{-2}{(\mathbf{Y}-\mathbf{X}\mathbf{\beta })}^{\prime }(\mathbf{I}-\text{ρ}{\mathbf{H}}^{\prime }){\mathbf{V}}_{\text{ν}}^{-1}(\mathbf{I}-\text{ρ}\mathbf{H})(\mathbf{Y}-\mathbf{X}\mathbf{\beta })\hfill \\ & =nlog\left(2\pi \right)-log\left(\left|{\text{σ}}_{\text{ν}}^{-2}({\mathbf{V}}_{\text{ν}}^{-1/2}-\text{ρ}{\mathbf{H}}^{\prime }{\mathbf{V}}_{\text{ν}}^{-1/2})({\mathbf{V}}_{\text{ν}}^{-1/2}-{\mathbf{V}}_{\text{ν}}^{-1/2}\text{ρ}\mathbf{H})\right|\right)+\hfill \\ \mathrm{(13)}& & {\displaystyle {\text{σ}}_{\text{ν}}^{-2}{(\mathbf{Y}-\mathbf{X}\mathbf{\beta })}^{\prime }({\mathbf{V}}_{\text{ν}}^{-1/2}-\text{ρ}{\mathbf{H}}^{\prime }{\mathbf{V}}_{\text{ν}}^{-1/2})({\mathbf{V}}_{\text{ν}}^{-1/2}-{\mathbf{V}}_{\text{ν}}^{-1/2}\text{ρ}\mathbf{H})(\mathbf{Y}-\mathbf{X}\beta ).}\hfill \end{array}$$

Consider the transformation $\tilde{\mathbf{Y}}={\mathbf{V}}_{\text{ν}}^{-1/2}\mathbf{Y}$. The spatial process $\tilde{\mathbf{Y}}$ is distributed as $\tilde{\mathbf{Y}}\sim \text{Gau}(\tilde{\mathbf{X}}\mathbf{\beta },{\text{σ}}_{\text{ν}}^2{(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}})}^{-1}{(\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })}^{-1})$, where $\tilde{\mathbf{X}}={\mathbf{V}}_{\text{ν}}^{-1/2}\mathbf{X}$, and $\tilde{\mathbf{H}}={\mathbf{V}}_{\text{ν}}^{-1/2}\mathbf{H}{\mathbf{V}}_{\text{ν}}^{1/2}$ (e.g., see [36], for details). Hence (13) can be written as,

$$\begin{array}{cc}\hfill -2\mathcal{L}& (\beta ,\text{ρ},{\text{σ}}_{\text{ν}}^2;\mathbf{Y})\hfill \\ & =nlog\left(2\pi \right)-log\left(\left|{\text{σ}}_{\text{ν}}^{-2}({\mathbf{V}}_{\text{ν}}^{-1/2}-{\mathbf{V}}_{\text{ν}}^{-1/2}\text{ρ}{\tilde{\mathbf{H}}}^{\prime })({\mathbf{V}}_{\text{ν}}^{-1/2}-\text{ρ}\tilde{\mathbf{H}}{\mathbf{V}}_{\text{ν}}^{-1/2})\right|\right)+\hfill \\ & {\displaystyle {\text{σ}}_{\text{ν}}^{-2}{({\mathbf{V}}_{\text{ν}}^{1/2}\tilde{\mathbf{Y}}-{\mathbf{V}}_{\text{ν}}^{1/2}\tilde{\mathbf{X}}\beta )}^{\prime }({\mathbf{V}}_{\text{ν}}^{-1/2}-{\mathbf{V}}_{\text{ν}}^{-1/2}\text{ρ}{\tilde{\mathbf{H}}}^{\prime })({\mathbf{V}}_{\text{ν}}^{-1/2}-\text{ρ}\tilde{\mathbf{H}}{\mathbf{V}}_{\text{ν}}^{-1/2})({\mathbf{V}}_{\text{ν}}^{1/2}\tilde{\mathbf{Y}}-{\mathbf{V}}_{\text{ν}}^{1/2}\tilde{\mathbf{X}}\beta )}\hfill \\ & =nlog\left(2\pi \right)-2log\left|{\mathbf{V}}_{\text{ν}}^{-1/2}\right|-log\left(\left|{\text{σ}}_{\text{ν}}^{-2}(\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}})\right|\right)+{\text{σ}}_{\text{ν}}^{-2}{(\tilde{\mathbf{Y}}-\tilde{\mathbf{X}}\beta )}^{\prime }(\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}})(\tilde{\mathbf{Y}}-\tilde{\mathbf{X}}\beta )\hfill \\ \mathrm{(14)}& & =-2\mathcal{L}(\beta ,\text{ρ},{\text{σ}}_{\text{ν}}^2;\tilde{\mathbf{Y}})+log\left|{\mathbf{V}}_{\text{ν}}\right|.\hfill \end{array}$$

The SAR model maximum likelihood estimates, $(\widehat{\mathbf{\beta }},\widehat{\text{ρ}},{\widehat{\text{σ}}}_{\text{ν}}^2)$, are those values that minimise $-2\mathcal{L}(\mathbf{\beta },\text{ρ},{\text{σ}}_{\text{ν}}^2;\mathbf{Y})$ or, equivalently, they are those values that minimise $-2\mathcal{L}(\mathbf{\beta },\text{ρ},{\text{σ}}_{\text{ν}}^2;\tilde{\mathbf{Y}})$.

Maximum likelihood estimates for the covariance parameters ρ and ${\text{σ}}_{\text{ν}}^2$ are well known to exhibit negative bias (e.g., [37]). To reduce the bias, estimation commonly proceeds using restricted (or residual) maximum likelihood (REML) (e.g., see [38], for a summary of its use for spatial data). REML estimation is carried out by performing maximum likelihood estimation on $\mathbf{M}\tilde{\mathbf{Y}}$, where $\mathbf{M}$ is chosen, such that $\text{E}\left(\mathbf{M}\tilde{\mathbf{Y}}\right)=\mathbf{0}$, rank$\left(\mathbf{M}\right)=n-p$, and p is the rank of the matrix of covariates, $\tilde{\mathbf{X}}$. Suppose we choose $\mathbf{M}={\tilde{\mathbf{P}}}^{\perp }\equiv \mathbf{I}-\tilde{\mathbf{X}}{\left({\tilde{\mathbf{X}}}^{\prime }\tilde{\mathbf{X}}\right)}^{-1}{\tilde{\mathbf{X}}}^{\prime }$. Note that ${\tilde{\mathbf{P}}}^{\perp }$ is symmetric (i.e., ${\left({\tilde{\mathbf{P}}}^{\perp }\right)}^{\prime }={\tilde{\mathbf{P}}}^{\perp }$) and idempotent (i.e., ${\left({\tilde{\mathbf{P}}}^{\perp }\right)}^2={\tilde{\mathbf{P}}}^{\perp }$). Effectively, we are substituting an ordinary-least-squares estimator of β (based on $\tilde{\mathbf{Y}}$ and $\tilde{\mathbf{X}}$) into (14), and then, we minimise this empirical version of the likelihood to obtain REML estimates. The “data” are now ${\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{Y}}$, with $E\left({\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{Y}}\right)=\mathbf{0}$ and $\text{cov}\left({\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{Y}}\right)={\tilde{\mathbf{P}}}^{\perp }\text{cov}\left(\tilde{\mathbf{Y}}\right){\tilde{\mathbf{P}}}^{\perp }$. Hence, we minimise the restricted likelihood,

$$\begin{array}{cc}\hfill -2\mathcal{L}(\mathbf{\beta },\text{ρ},{\text{σ}}_{\text{ν}}^2;{\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{Y}})& =nlog\left(2\pi \right)-log\left(\left|{\text{σ}}_{\text{ν}}^{-2}{\tilde{\mathbf{P}}}^{\perp }(\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}}){\tilde{\mathbf{P}}}^{\perp }\right|\right)+\hfill \\ & {\displaystyle {\text{σ}}_{\text{ν}}^{-2}{({\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{Y}}-{\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{X}}\mathbf{\beta })}^{\prime }{\tilde{\mathbf{P}}}^{\perp }(\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}}){\tilde{\mathbf{P}}}^{\perp }({\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{Y}}-{\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{X}}\mathbf{\beta })}\hfill \\ \mathrm{(15)}& \hfill & =nlog\left(2\pi \right)-log\left(\left|{\text{σ}}_{\text{ν}}^{-2}{\tilde{\mathbf{P}}}^{\perp }(\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}}){\tilde{\mathbf{P}}}^{\perp }\right|\right)+{\text{σ}}_{\text{ν}}^{-2}{\tilde{\mathbf{Y}}}^{\prime }{\tilde{\mathbf{P}}}^{\perp }(\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}}){\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{Y}},\hfill \end{array}$$

resulting in the REML estimates. From (15), these estimates depend critically on the non-negative-definite matrix,

$${\text{σ}}_{\text{ν}}^{-2}{\tilde{\mathbf{P}}}^{\perp }(\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}}){\tilde{\mathbf{P}}}^{\perp }.$$

(16)

Expanding (16) gives the equivalent expression,

$${\text{σ}}_{\text{ν}}^{-2}\left({\tilde{\mathbf{P}}}^{\perp }{\tilde{\mathbf{P}}}^{\perp }-2\text{ρ}\frac{{\tilde{\mathbf{P}}}^{\perp }({\tilde{\mathbf{H}}}^{\prime }+\tilde{\mathbf{H}}){\tilde{\mathbf{P}}}^{\perp }}{2}+{\text{ρ}}^2{\tilde{\mathbf{P}}}^{\perp }{\tilde{\mathbf{H}}}^{\prime }\tilde{\mathbf{H}}{\tilde{\mathbf{P}}}^{\perp }\right).$$

(17)

In (17), the dominant term in ρ is proportional to:

$${\tilde{\mathbf{P}}}^{\perp }\left\{\frac{{\tilde{\mathbf{H}}}^{\prime }+\tilde{\mathbf{H}}}{2}\right\}{\tilde{\mathbf{P}}}^{\perp }.$$

(18)

Notice that Moran’s I statistic for data $\tilde{\mathbf{Y}}$ and covariates $\tilde{\mathbf{X}}$ is,

$$\mathbf{I}\left(\mathbf{A}\right)=\frac{n}{{\mathbf{1}}^{\prime }\mathbf{A}\mathbf{1}}\frac{{\tilde{\mathbf{Y}}}^{\prime }{\tilde{\mathbf{P}}}^{\perp }\mathbf{A}{\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{Y}}}{{\tilde{\mathbf{Y}}}^{\prime }{\tilde{\mathbf{P}}}^{\perp }\tilde{\mathbf{Y}}},$$

(19)

where $\mathbf{A}$ is the adjacency matrix whose entries are one (when two entries are “adjacent”) and zero (otherwise); see [24]. When an asymmetric matrix $\mathbf{B}$ is used in place of $\mathbf{A}$, Tiefelsdorf ([39], p. 29) has suggested that $\mathbf{B}$ be “symmetrised” by replacing $\mathbf{A}$ in (19) with the symmetric matrix $({\mathbf{B}}^{\prime }+\mathbf{B})/2$.

We have seen that (18) plays an important role in the restricted likelihood (15), and it takes the form of the numerator of Tiefelsdorf’s modification to the Moran statistic for an asymmetric spatial dependence matrix. Hence, we use it to define a generalised Moran operator. Recall that ${\mathbf{V}}_{\text{ν}}$ captures heterogeneity and $\mathbf{H}$ captures asymmetric spatial dependence in the SAR model. We can write ${\tilde{\mathbf{P}}}^{\perp }\equiv \mathbf{I}-{\mathbf{V}}_{\text{ν}}^{-1/2}\mathbf{X}{\left({\mathbf{X}}^{\prime }{\mathbf{V}}_{\text{ν}}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }{\mathbf{V}}_{\text{ν}}^{-1/2}$ and, since $\tilde{\mathbf{H}}={\mathbf{V}}_{\text{ν}}^{-1/2}\mathbf{H}{\mathbf{V}}_{\text{ν}}^{1/2}$,

$$\begin{array}{cc}\hfill {\tilde{\mathbf{P}}}^{\perp }({\tilde{\mathbf{H}}}^{\prime }+\tilde{\mathbf{H}}){\tilde{\mathbf{P}}}^{\perp }& ={\tilde{\mathbf{P}}}^{\perp }{\mathbf{V}}_{\text{ν}}^{1/2}({\mathbf{H}}^{\prime }{\mathbf{V}}_{\text{ν}}^{-1}+{\mathbf{V}}_{\text{ν}}^{-1}\mathbf{H}){\mathbf{V}}_{\text{ν}}^{1/2}{\tilde{\mathbf{P}}}^{\perp }\hfill \\ & =\mathbf{Q}\left({\mathbf{V}}_{\text{ν}}\right)({\mathbf{H}}^{\prime }{\mathbf{V}}_{\text{ν}}^{-1}+{\mathbf{V}}_{\text{ν}}^{-1}\mathbf{H})\mathbf{Q}{\left({\mathbf{V}}_{\text{ν}}\right)}^{\prime },\hfill \end{array}$$

where:

$$\mathbf{Q}\left({\mathbf{V}}_{\text{ν}}\right)\equiv {\mathbf{V}}_{\text{ν}}^{1/2}-{\mathbf{V}}_{\text{ν}}^{-1/2}\mathbf{X}{\left({\mathbf{X}}^{\prime }{\mathbf{V}}_{\text{ν}}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }.$$

(20)

Hence, we define the generalized Moran operator (GMO),

$$\mathbf{G}(\mathbf{X},\mathbf{H},{\mathbf{V}}_{\text{ν}})\equiv \mathbf{Q}\left({\mathbf{V}}_{\text{ν}}\right)\left\{\frac{{\mathbf{H}}^{\prime }{\mathbf{V}}_{\text{ν}}^{-1}+{\mathbf{V}}_{\text{ν}}^{-1}\mathbf{H}}{2}\right\}\mathbf{Q}{\left({\mathbf{V}}_{\text{ν}}\right)}^{\prime },$$

(21)

where $\mathbf{Q}\left({\mathbf{V}}_{\text{ν}}\right)$ is defined by (20).

Notice that for ${\mathbf{V}}_{\text{ν}}=\mathbf{I}$ and $\mathbf{H}=\mathbf{A}$ (the symmetric adjacency matrix), we obtain $\mathbf{Q}\left(\mathbf{I}\right)={\mathbf{P}}^{\perp }\equiv \mathbf{I}-\mathbf{X}{\left({\mathbf{X}}^{\prime }\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }$, and hence, (21) is equal to ${\mathbf{P}}^{\perp }\mathbf{A}{\mathbf{P}}^{\perp }$. This is the Moran operator given by Griffith [24] and Hughes and Haran [14]. Spectral decomposition of the Moran operator has been shown to yield natural spatial basis functions for spatial modelling (e.g., see [24,27,40,41] and, more recently, [14,28]), and we propose the same for the GMO given by (21). For spatial generalised linear models where heteroskedasticity appears through a mean-variance relationship, Reich et al. [27] suggest modifying the Moran operator using weights obtained from iteratively reweighted least squares. This bears a formal resemblance to ${\tilde{\mathbf{P}}}^{\perp }$ and lends support to our proposed generalisation. Notice that the GMO allows for both heteroskedasticity (through ${\mathbf{V}}_{\text{ν}}$) and asymmetry (through $\mathbf{H}$) simultaneously.

4. Specification of the SAR-Calibrated SRE Model

In this section, we define an SAR-calibrated SRE model for the latent process $\mathbf{Y}$, using the SRE model that is defined in Section 2.2. Our approach is to use $r\ll n$ eigenvectors from the spectral decomposition of $\mathbf{G}(\mathbf{X},\mathbf{H},{\mathbf{V}}_{\text{ν}})$ given by (21), to define the r basis functions for the SRE model, (8). Now, the spectral decomposition of (21) yields,

$$\begin{array}{c}\hfill \mathbf{G}(\mathbf{X},\mathbf{H},{\mathbf{V}}_{\text{ν}})={\tilde{\mathbf{P}}}^{\perp }\left\{\frac{{\tilde{\mathbf{H}}}^{\prime }+\tilde{\mathbf{H}}}{2}\right\}{\tilde{\mathbf{P}}}^{\perp }\equiv \mathbf{U}\mathbf{\Omega }{\mathbf{U}}^{\prime },\end{array}$$

(22)

where the eigenvectors $\mathbf{U}$ of $\mathbf{G}(\mathbf{X},\mathbf{H},{\mathbf{V}}_{\text{ν}})$ satisfy ${\mathbf{U}}^{\prime }\mathbf{U}=\mathbf{U}{\mathbf{U}}^{\prime }=\mathbf{I}$; and $\mathbf{\Omega }\equiv \text{diag}({\text{ω}}_1,...,{\text{ω}}_n)$ is an $n\times n$ diagonal matrix with elements in the diagonal matrix ordered so that $|{\text{ω}}_1|\ge |{\text{ω}}_2|\ge ...\ge |{\text{ω}}_n|\ge 0$, for ${\text{ω}}_1,...,{\text{ω}}_n$, the eigenvalues of $\mathbf{G}(\mathbf{X},\mathbf{H},{\mathbf{V}}_{\text{ν}})$. Define the $n\times r$ matrix ${\mathbf{U}}_1\equiv ({\mathbf{u}}_1,...,{\mathbf{u}}_r)$, which corresponds to the r largest absolute eigenvalues of $\mathbf{G}$, namely the diagonal elements in ${\mathbf{\Omega }}_1\equiv \text{diag}({\text{ω}}_1,...,{\text{ω}}_r)$. Similarly, ${\mathbf{U}}_2$ and ${\mathbf{\Omega }}_2$ correspond to the $n-r$ smallest absolute eigenvectors. Then,

$$\begin{array}{cc}\hfill \mathbf{U}\mathbf{\Omega }{\mathbf{U}}^{\prime }=({\mathbf{U}}_1,{\mathbf{U}}_2)\left(\begin{array}{cc}{\mathbf{\Omega }}_1& \mathbf{0}\\ \mathbf{0}& {\mathbf{\Omega }}_2\end{array}\right)\left(\begin{array}{c}{\mathbf{U}}_1^{\prime }\\ {\mathbf{U}}_2^{\prime }\end{array}\right)& ={\mathbf{U}}_1{\mathbf{\Omega }}_1{\mathbf{U}}_1^{\prime }+{\mathbf{U}}_2{\mathbf{\Omega }}_2{\mathbf{U}}_2^{\prime }.\hfill \end{array}$$

(23)

Notice that eigenvectors corresponding to large positive and negative eigenvalues may be chosen; Hughes and Haran [14] only choose eigenvectors that correspond to large positive eigenvalues, which limits their approach to capturing only positive spatial dependence. In the application given in Section 6, several large absolute eigenvectors that corresponded to negative dependence were chosen.

From (3), using a second-order Taylor-series expansion for ρ around zero, the SAR model covariance matrix for $\tilde{\mathbf{Y}}$ is:

$${\text{σ}}_{\text{ν}}^2{(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}})}^{-1}{(\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })}^{-1}\simeq {\text{σ}}_{\text{ν}}^2\left\{\left[\mathbf{I}+\text{ρ}({\tilde{\mathbf{H}}}^{\prime }+\tilde{\mathbf{H}})+{\text{ρ}}^2\left({\left({\tilde{\mathbf{H}}}^{\prime }\right)}^2+{\tilde{\mathbf{H}}}^{\prime }\tilde{\mathbf{H}}+{\tilde{\mathbf{H}}}^2\right)\right]\right\}.$$

(24)

Hence, a candidate for the $r\times r$ matrix $\mathbf{K}$ in the SRE model (8) is obtained by pre- and post-multiplying (24) by ${\tilde{\mathbf{U}}}_1$, where ${\tilde{\mathbf{U}}}_1\equiv {\mathbf{V}}_{\text{ν}}^{-1/2}{\mathbf{U}}_1$:

$$\begin{array}{c}\hfill {\text{σ}}_{\text{ν}}^2\left\{{\tilde{\mathbf{U}}}_1^{\prime }\left[\mathbf{I}+\text{ρ}({\tilde{\mathbf{H}}}^{\prime }+\tilde{\mathbf{H}})+{\text{ρ}}^2\left({\left({\tilde{\mathbf{H}}}^{\prime }\right)}^2+{\tilde{\mathbf{H}}}^{\prime }\tilde{\mathbf{H}}+{\tilde{\mathbf{H}}}^2\right)\right]{\tilde{\mathbf{U}}}_1\right\}.\end{array}$$

(25)

The candidate $\mathbf{K}$ given by (25) is fast to evaluate for given values of ρ and ${\text{σ}}_{\text{ν}}^2$. However, while the $r\times r$ matrix in (25) is symmetric, it is not necessarily positive-definite. Higham [42] shows how to obtain the nearest symmetric positive approximant to any real square matrix $\mathbf{B}$. Call this positive-definite matrix ${\mathcal{A}}^{+}\left(\mathbf{B}\right)$, and hence, define,

$$\begin{array}{c}\hfill {\mathbf{R}}^{+}\left(\text{ρ}\right)\equiv {\mathcal{A}}^{+}\left({\tilde{\mathbf{U}}}_1^{\prime }\left[\mathbf{I}+\text{ρ}({\tilde{\mathbf{H}}}^{\prime }+\tilde{\mathbf{H}})+{\text{ρ}}^2\left({\left({\tilde{\mathbf{H}}}^{\prime }\right)}^2+{\tilde{\mathbf{H}}}^{\prime }\tilde{\mathbf{H}}+{\tilde{\mathbf{H}}}^2\right)\right]{\tilde{\mathbf{U}}}_1\right).\end{array}$$

(26)

Then, since ${\text{σ}}_{\text{ν}}^2>0$, the $r\times r$ matrix ($r\ll n$),

$$\begin{array}{c}\hfill {\mathbf{K}}^{+}(\text{ρ},{\text{σ}}_{\text{ν}}^2)\equiv {\text{σ}}_{\text{ν}}^2{\mathbf{R}}^{+}\left(\text{ρ}\right),\end{array}$$

(27)

is positive-definite. In conclusion, we shall use the $r\times r$ positive-definite matrix ${\mathbf{K}}^{+}$ given by (27) to define our SAR-calibrated SRE model. What remains in the calibration is specification of ${\text{σ}}_{\xi }^2{\mathbf{V}}_{\xi }$ in (8). Since ${\mathbf{V}}_{\xi }$ is assumed known, we only need to specify ${\text{σ}}_{\xi }^2$, which will generally be a function of ρ and ${\text{σ}}_{\text{ν}}^2$. We use the closeness of the inverses of the SAR and SRE covariance matrices to achieve this, which is concomitant with the eigenvectors in ${\mathbf{U}}_1$ coming from large absolute eigenvalues.

Recall that the SAR model’s inverse covariance matrix is easy to evaluate as is the SRE model’s inverse covariance matrix using the Sherman–Morrison–Woodbury identity (e.g., [43]): for generic matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{U}$ and $\mathbf{V}$ of appropriate order,

$${(\mathbf{A}+\mathbf{U}\mathbf{B}\mathbf{V})}^{-1}={\mathbf{A}}^{-1}-{\mathbf{A}}^{-1}\mathbf{U}{({\mathbf{B}}^{-1}+\mathbf{V}{\mathbf{A}}^{-1}\mathbf{U})}^{-1}\mathbf{V}{\mathbf{A}}^{-1},$$

where the inverses on the right-hand side are fast to evaluate. Consequently, for an SRE model for $\tilde{\mathbf{Y}}$ with the form of (8) and with $\tilde{\mathbf{S}}={\tilde{\mathbf{U}}}_1$, diagonal ${\tilde{\mathbf{V}}}_{\xi }={\mathbf{V}}_{\text{ν}}^{-1/2}{\mathbf{V}}_{\xi }{\mathbf{V}}_{\text{ν}}^{-1/2}$ and positive-definite $\mathbf{K}\equiv {\mathbf{K}}^{+}(\text{ρ},{\text{σ}}_{\text{ν}}^2)={\text{σ}}_{\text{ν}}^2{\mathbf{R}}^{+}\left(\text{ρ}\right)$, the inverse covariance matrix is,

$$\begin{array}{cc}\hfill {(\tilde{\mathbf{S}}\mathbf{K}{\tilde{\mathbf{S}}}^{\prime }+{\text{σ}}_{\xi }^2{\tilde{\mathbf{V}}}_{\xi })}^{-1}& ={({\tilde{\mathbf{U}}}_1{\text{σ}}_{\text{ν}}^2{\mathbf{R}}^{+}\left(\text{ρ}\right){\tilde{\mathbf{U}}}_1^{\prime }+{\text{σ}}_{\xi }^2{\tilde{\mathbf{V}}}_{\xi })}^{-1}\hfill \\ & ={\text{σ}}_{\text{ν}}^{-2}\{\zeta {\tilde{\mathbf{V}}}_{\xi }^{-1}-{\zeta }^2{\tilde{\mathbf{V}}}_{\xi }^{-1}{\tilde{\mathbf{U}}}_1{\left({\mathbf{R}}^{+}{\left(\text{ρ}\right)}^{-1}+\zeta {\tilde{\mathbf{U}}}_1^{\prime }{\tilde{\mathbf{V}}}_{\xi }^{-1}{\tilde{\mathbf{U}}}_1\right)}^{-1}{\tilde{\mathbf{U}}}_1^{\prime }{\tilde{\mathbf{V}}}_{\xi }^{-1}\}\hfill \\ \mathrm{(28)}& & \equiv {\text{σ}}_{\text{ν}}^{-2}\mathbf{J}(\text{ρ},\zeta ),\hfill \end{array}$$

where $\zeta \equiv {\text{σ}}_{\text{ν}}^2/{\text{σ}}_{\xi }^2$. Consider the Frobenius norm between the inverse of the SAR covariance matrix of $\tilde{\mathbf{Y}}$ and (28); this is proportional to:

$$\begin{array}{c}\hfill \parallel (\mathbf{I}-\text{ρ}{\tilde{\mathbf{H}}}^{\prime })(\mathbf{I}-\text{ρ}\tilde{\mathbf{H}})-\mathbf{J}(\text{ρ},\text{ζ}){\parallel }_F^2,\end{array}$$

(29)

as a function of ζ. The final step in the calibration is to minimise (29) with respect to $\text{ζ}>0$, for all ρ in its parameter space; call the minimised value $\text{ζ}\left(\text{ρ}\right)$.

Hence, the SAR-calibrated SRE covariance matrix of $\mathbf{Y}$ is,

$$\begin{array}{cc}\hfill \Sigma (\text{ρ},{\text{σ}}_{\text{ν}}^2)& \equiv {\mathbf{V}}_{\text{ν}}^{1/2}\text{cov}\left(\tilde{\mathbf{Y}}\right){\mathbf{V}}_{\text{ν}}^{1/2}\hfill \\ \mathrm{(30)}& & ={\text{σ}}_{\text{ν}}^2\{{\mathbf{U}}_1{\mathbf{R}}^{+}\left(\text{ρ}\right){\mathbf{U}}_1^{\prime }+\text{ζ}{\left(\text{ρ}\right)}^{-1}{\mathbf{V}}_{\xi }\},\hfill \end{array}$$

where ${\text{σ}}_{\text{ν}}^2{\mathbf{R}}^{+}\left(\text{ρ}\right)$ is a positive-definite $r\times r$ matrix, $\text{ζ}\left(\text{ρ}\right)>0$, by construction, and recall that ${\mathbf{V}}_{\xi }$ is known. Notice that all parameters in the SRE model’s covariance matrix (30) are expressed as functions of the SAR model’s parameters ${\text{σ}}_{\text{ν}}^2$ and ρ. Furthermore, notice that, while we have made $\mathbf{\Sigma }{(\text{ρ},{\text{σ}}_{\text{ν}}^2)}^{-1}$ close (in the Frobenius-norm sense) to $(\mathbf{I}-\text{ρ}{\mathbf{H}}^{\prime }){\mathbf{V}}_{\text{ν}}^{-1}(\mathbf{I}-\text{ρ}\mathbf{H})$, they are not equal; this difference is the only source of approximation in this article. The dominant (positive and negative) spatial dependencies in the SAR model are included in the calibrated model in a reduced-rank form, and residual variability is captured by setting ${\text{σ}}_{\xi }^2={\text{σ}}_{\text{ν}}^2/\text{ζ}$ (which is a function of ρ). In the next section, we fit the model,

$$\begin{array}{c}\hfill \mathbf{Z}\sim \text{Gau}(\mathbf{X}\mathbf{\beta },\mathbf{\Sigma }(\text{ρ},{\text{σ}}_{\text{ν}}^2)+{\text{σ}}_{\epsilon }^2{\mathbf{V}}_{\epsilon }),\end{array}$$

(31)

using likelihood-based methods, where $\mathbf{\Sigma }(\text{ρ},{\text{σ}}_{\text{ν}}^2)$ is given by (30).

5. Parameter Estimation

We use likelihood-based estimators of the parameters of the SAR-calibrated SRE model (31), because they have asymptotic optimality properties, and they are preserved under transformation. The SAR model’s parameters are β, ρ and ${\text{σ}}_{\text{ν}}^2$, which appear in the SRE likelihood defined by (30) and (31); we assume that the measurement-error component, ${\text{σ}}_{\epsilon }^2{\mathbf{V}}_{\epsilon }$, is known, which is often the case (Section 6.3).

We estimate the SAR covariance parameters ρ and ${\text{σ}}_{\text{ν}}^2$ using the log-restricted-likelihood function given by:

$$\begin{array}{cc}\hfill {\mathcal{L}}_R(\text{ρ},{\text{σ}}_{\text{ν}}^2)=\text{constant}-\frac{1}{2}& (log|\mathbf{\Phi }|+log|{\mathbf{X}}^{\prime }{\mathbf{\Phi }}^{-1}\mathbf{X}|+{\mathbf{Z}}^{\prime }\mathbf{\Pi }\mathbf{Z}),\hfill \end{array}$$

(32)

where $\mathbf{\Phi }\equiv {\text{σ}}_{\text{ν}}^2\left\{{\mathbf{U}}_1{\mathbf{R}}^{+}\left(\text{ρ}\right){\mathbf{U}}_1^{\prime }+\text{ζ}{\left(\text{ρ}\right)}^{-1}{\mathbf{V}}_{\xi }\right\}+{\text{σ}}_{\epsilon }^2{\mathbf{V}}_{\epsilon }$ and $\mathbf{\Pi }\equiv {\mathbf{\Phi }}^{-1}-{\mathbf{\Phi }}^{-1}\mathbf{X}{\left({\mathbf{X}}^{\prime }{\mathbf{\Phi }}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }{\mathbf{\Phi }}^{-1}$. We solve for (ρ, ${\text{σ}}_{\text{ν}}^2$) by minimising $-2{\mathcal{L}}_R(\text{ρ},{\text{σ}}_{\text{ν}}^2)$ in terms of ρ and ${\text{σ}}_{\text{ν}}^2\ge 0$ in two stages. First, for each value of ρ, we minimise $-2{\mathcal{L}}_R(\text{ρ},{\text{σ}}_{\text{ν}}^2)$ with respect to ${\text{σ}}_{\text{ν}}^2$, to obtain ${\widehat{\text{σ}}}_{\text{ν}}^2\left(\text{ρ}\right)$. Then, we minimise the profile, $-2{\mathcal{L}}_R(\text{ρ},{\widehat{\text{σ}}}_{\text{ν}}^2\left(\text{ρ}\right))$, with respect to ρ. The resulting $\widehat{\text{ρ}}$ and ${\widehat{\text{σ}}}_{\text{ν}}^2\equiv {\widehat{\text{σ}}}_{\text{ν}}^2\left(\widehat{\text{ρ}}\right)$ are the REML estimates. When making likelihood calculations, we use the Sherman–Morrison–Woodbury formulas [43] for both the inverse and the determinant of the SRE model’s covariance function.

This estimation technique requires a starting value for ${\text{σ}}_{\text{ν}}^2$ for each value of ρ; we use method-of-moments estimation to obtain the starting values. Specifically, define $\mathbf{W}\equiv {\mathbf{P}}^{\perp }\mathbf{Z}$, and hence:

$$\mathbf{W}\sim \text{Gau}(\mathbf{0},{\text{σ}}_{\text{ν}}^2\mathbf{L}\left(\text{ρ}\right)+\mathbf{M}),$$

(33)

where $\mathbf{L}\left(\text{ρ}\right)\equiv {\mathbf{P}}^{\perp }\{{\mathbf{U}}_1{\mathbf{R}}^{+}\left(\text{ρ}\right){\mathbf{U}}_1^{\prime }+\text{ζ}\left(\text{ρ}\right){\mathbf{V}}_{\xi }\}{\mathbf{P}}^{\perp }$ and $\mathbf{M}\equiv {\text{σ}}_{\epsilon }^2{\mathbf{P}}^{\perp }{\mathbf{V}}_{\epsilon }{\mathbf{P}}^{\perp }$. Consequently,

$$E\left({\mathbf{W}}^{\prime }\mathbf{N}\mathbf{W}\right)=\text{trace}\left(\mathbf{N}\{{\text{σ}}_{\text{ν}}^2\mathbf{L}\left(\text{ρ}\right)+\mathbf{M}\}\right),$$

where $\mathbf{N}$ is any symmetric, $n\times n$ matrix. For example, in Section 6, we choose $\mathbf{N}=\text{diag}(n_1,...,n_m)$, where $n_i$, $i=1,...,m$ is the number of people in area i. The resulting method-of-moments estimator is,

$${\text{σ}}_{\text{ν},0}^2\left(\text{ρ}\right)=\frac{{\mathbf{Z}}^{\prime }{\mathbf{P}}^{\perp }\mathbf{N}{\mathbf{P}}^{\perp }\mathbf{Z}-\text{trace}\left({\text{σ}}_{\epsilon }^2\mathbf{N}{\mathbf{P}}^{\perp }{\mathbf{V}}_{\epsilon }{\mathbf{P}}^{\perp }\right)}{\text{trace}\left(\mathbf{N}\{{\mathbf{P}}^{\perp }{\mathbf{U}}_1{\mathbf{R}}^{+}\left(\text{ρ}\right){\mathbf{U}}_1^{\prime }{\mathbf{P}}^{\perp }+\text{ζ}{\left(\text{ρ}\right)}^{-1}{\mathbf{P}}^{\perp }{\mathbf{V}}_{\xi }{\mathbf{P}}^{\perp }\}\right)},$$

(34)

which we use as a starting value for solving:

$${\widehat{\text{σ}}}_{\text{ν}}{\left(\text{ρ}\right)}^2=arg\underset{{\text{σ}}_{\text{ν}}{\left(\text{ρ}\right)}^2}{min}\left\{-2{\mathcal{L}}_R(\text{ρ},{\text{σ}}_{\text{ν}}^2)\right\},$$

(35)

for each ρ in its parameter space. The profile restricted log-likelihood for ρ is given by,

$${\mathcal{L}}_R^{*}\left(\text{ρ}\right)\equiv {\mathcal{L}}_R(\text{ρ},{\widehat{\text{σ}}}_{\text{ν}}^2\left(\text{ρ}\right)),$$

for ρ in its parameter space. The minimum of $-2{\mathcal{L}}_R^{*}\left(\text{ρ}\right)$ is found using a univariate grid-search algorithm.

A likelihood-based confidence interval (CI) for a model parameter, with a significance level of α, can be obtained by inverting the (restricted) likelihood-ratio test. Hence, an approximate $100(1-\alpha )$% CI for ρ is given by:

$$\{\text{ρ}:{\mathcal{L}}_R^{*}\left(\text{ρ}\right)>{\mathcal{L}}_R(\widehat{\text{ρ}},{\widehat{\text{σ}}}_{\text{ν}}^2\left(\widehat{\text{ρ}}\right))-{\chi }_1^2(1-\alpha )/2\},$$

(36)

where ${\chi }_1^2(1-\alpha )$ is the $(1-\alpha )100$ percentile of the chi-squared distribution on one degree of freedom.

Finally, the empirical generalised-least-squares estimate for β is,

$$\widehat{\mathbf{\beta }}\equiv {\left({\mathbf{X}}^{\prime }\mathbf{\Phi }{(\widehat{\text{ρ}},{\widehat{\text{σ}}}_{\text{ν}}^2\left(\widehat{\text{ρ}}\right))}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }\mathbf{\Phi }{(\widehat{\text{ρ}},{\widehat{\text{σ}}}_{\text{ν}}^2\left(\widehat{\text{ρ}}\right))}^{-1}\mathbf{Z},$$

(37)

and an estimate of its covariance matrix is ${\left({\mathbf{X}}^{\prime }\mathbf{\Phi }{(\widehat{\text{ρ}},{\widehat{\text{σ}}}_{\text{ν}}^2\left(\widehat{\text{ρ}}\right))}^{-1}\mathbf{X}\right)}^{-1}$.

6. Fitting a Spatial Model to Mean Usual Weekly Income

We illustrate the methodology developed in Section 2, Section 3, Section 4 and Section 5 using Australian Census data [44] from the 2011 Census conducted by the Australian Bureau of Statistics (ABS). Our interest lies in modelling the spatial distribution of mean usual weekly income (MWI) for small geographic areas in the state of New South Wales (NSW), Australia, and then fitting the spatial model by estimating its parameters. In the following sections, we describe the data, the fitting of an SAR-calibrated SRE model, our assumptions regarding measurement error and the statistical inferences that follow.

6.1. Imputed MWI

In NSW, 2011 Census summaries are available for four nested sets of geographic areas, called Statistical Areas (SA), where the most disaggregated level consists of SA1 areas. We fit a spatial model to MWI for the $16,850$ SA1 areas (Figure 1) with a population >100 persons, which is 97% of the total SA1 areas; populations ≤100 are too granular for the statistical analysis that we shall undertake. SA1 areas in NSW exhibit a wide range in both population and area, since the state’s population is concentrated along the coastline (in the east of the state). For the SA1 areas we consider, their populations range from 100 to 13,000 (with a median of 395).

In the Census, each person ≥15 years specifies their usual mean weekly income by selecting from 11 possible income ranges {≤$0, $1 to $199, $200 to $299, $300 to $399, $400 to $599, $600 to $799, $800 to $999, $1,000 to $1,249, $1,250 to $1,499, $1,500 to $1,999, ≥$2,000}. For the $i=1,...,16,850$ SA1 areas, the Census data are the number of persons, $n_{ij}$, in the j-th income range. Thus, the 2011 Census does not directly measure MWI for SA1s; denote the true MWI process as $Y(·)$. However, in the 2011 to 2012 Australian Survey of Income and Housing (SIH) [45], the mean usual weekly income of a large sample of individuals was measured, but with limited spatial information. Hence, we use the SIH microdata to impute an MWI value, $\overline{Z}\left({\mathbf{s}}_i\right)$, for the 2011 Census SA1 area $i=1,...,16,850$, where the “bar” on $\overline{Z}(·)$ is added to emphasise that MWI refers to mean income. In this analysis, we use only the SIH data and Census counts from the $J=9$ Census income ranges corresponding to $0 < MWI < $2,000, which includes most (approximately 94%) of the incomes in the Census. Negative, zero and very large incomes are set aside. These incomes, along with SA1 areas with populations less than 100, need special attention that is beyond the goal of this section.

Figure 1. (a) Chloropleth map of imputed mean usual weekly income (MWI) ($0 < MWI < $2,000) for Statistical Area 1 (SA1) areas in NSW, Australia; the resolution around Sydney is enhanced. (b) An excerpt of SA1 areas from (a) showing how asymmetry can occur in the spatial dependence matrix. The lines join the centroids of neighbouring areas, which are obtained using a K-nearest-neighbour algorithm, where here $K=8$.

Figure 1. (a) Chloropleth map of imputed mean usual weekly income (MWI) ($0 < MWI < $2,000) for Statistical Area 1 (SA1) areas in NSW, Australia; the resolution around Sydney is enhanced. (b) An excerpt of SA1 areas from (a) showing how asymmetry can occur in the spatial dependence matrix. The lines join the centroids of neighbouring areas, which are obtained using a K-nearest-neighbour algorithm, where here $K=8$.

To impute the MWI values, we first stratify the Australia-wide SIH data using the Census income ranges, resulting in ${\tilde{n}}_j$ SIH records $\{{\tilde{Z}}_{jl}:l=1,...,{\tilde{n}}_j\}$ for the $j=1,...,9$ income ranges. The SIH data include weights that are used to account for unequal selection probabilities when calculating population estimates. We use them to impute $n_{ij}$ MWI values, $\{Z_{ijk}:k=1,...,n_{ij}\}$, by sampling $n_{ij}$ SIH records, with replacement, from $\{{\tilde{Z}}_{jl}:l=1,...,{\tilde{n}}_j\}$, for each of $i=1,...,16,850$, $j=1,...,9$. As mentioned above, the MWI values in the samples are selected using the same probability distribution that was used to select individuals in the SIH survey. We average the imputed values for each area, to obtain an imputed MWI value:

$$\overline{Z}\left({\mathbf{s}}_i\right)\equiv \frac{{\sum }_{j=1}^9{\sum }_{k=1}^{n_{ij}}{Z}_{ijk}}{n_i},\text{for}i=1,...,16,850,$$

(38)

where $n_i\equiv {\sum }_{j=1}^9{n}_{ij}$. Hence, the spatial data that we shall analyse in this section are,

$$\mathbf{Z}\equiv {(\overline{Z}\left({\mathbf{s}}_1\right),...,\overline{Z}\left({\mathbf{s}}_{16850}\right))}^{\prime }.$$

(39)

6.2. A Simple Heteroskedastic Model for MWI

A basic heteroskedastic model for MWI was explored and found to fit reasonably well: assume that for each area i and stratum j, income values are a random sample from a distribution with mean ${\text{μ}}_j$ and variance ${\text{σ}}_j^2$. From (38) and for given $\left\{n_{ij}\right\}$,

$$\text{E}\left(\overline{Z}\left({\mathbf{s}}_i\right)\right)=\frac{{\sum }_{j=1}^9{n}_{ij}{\text{μ}}_j}{n_i},\text{and}\text{var}\left(\overline{Z}\left({\mathbf{s}}_i\right)\right)=\frac{{\sum }_{j=1}^9(n_{ij}/n_i){\text{σ}}_j^2}{n_i}.$$

(40)

Contrast this with an overly simple model that assumes a constant mean and variance for the elements of $\mathbf{Z}$.

In Figure 2, the left panel gives a histogram and Q-Q plot of unstandardised $\left\{\overline{Z}\left({\mathbf{s}}_i\right)\right\}$ motivated by the overly simple model that assumes constant mean and variance over the SA1 areas. The right panel gives a histogram and Q-Q plot of the standardised values from the basic heteroskedastic model:

$$\frac{\overline{Z}\left({\mathbf{s}}_i\right)-E\left(\overline{Z}\left({\mathbf{s}}_i\right)\right)}{\text{var}{\left(\overline{Z}\left({\mathbf{s}}_i\right)\right)}^{1/2}},\text{for}i=1,...,16,850,$$

(41)

where $\left\{{\text{μ}}_j\right\}$ and $\left\{{\text{σ}}_j^2\right\}$ in (40) are estimated straightforwardly from the SIH data. The fit of the histogram in the right panel to a $\text{Gau}(0,1)$ distribution is striking and indicates that any spatial model we fit should aim to capture the same types of heteroskedastic behaviour given by (40). Of course, $\left\{n_{ij}\right\}$ is an important part of the randomness in $\mathbf{Z}$ defined by (38) and (39), and in a formal (spatial) model, we cannot have them describing the first two moments (as we do informally in (40)).

To capture the heteroskedasticity in (40) and illustrated in Figure 2, recall from (5) that ${\mathbf{\Sigma }}_{\text{ν}}={\text{σ}}_{\text{ν}}^2{\mathbf{V}}_{\text{ν}}$. We estimate $V_{\text{ν}}\left({\mathbf{s}}_i\right)$ using the $m=1,...,100$ centiles of the Socio-Economic Indicator For Areas Index of Economic Resources (IER), a well-established socio-economic indicator produced for SA1 areas in Australia by the ABS [46]. There are 35 SA1 areas out of 16,850 that do not have an IER. For each of these, we obtained a value for the IER based on a regression of the 16,815 IER values that we do have, regressed on the covariates included in our model (see Section 6.5). Finally, then, we have 16,850 IER values whose centiles define 100 groups of SA1 areas, as follows:

$$I_m\equiv \{i:i\text{-th}\text{SA1}\text{area}\in m\text{-th}\text{decile}\},\text{for}m=1,...,100.$$

Define $q_m\equiv \left|I_m\right|$, ${\overline{n}}_{mj}\equiv {\sum }_{i\in I_m}{n}_{ij}/q_m$, ${\overline{n}}_m\equiv {\sum }_{i\in I_m}{n}_i/q_m$, and hence, we set:

$$V_{\text{ν}}\left({\mathbf{s}}_i\right)\equiv \sum _{j=1}^9\frac{({\overline{n}}_{mj}/{\overline{n}}_m){\text{σ}}_j^2}{n_i},\text{for}i\in {\text{I}}_m,m=1,...,100.$$

(42)

That is, in (42), we replace the random weights $\{(n_{ij}/n_i):j=1,...,9\}$ given in (40) with centile-smoothed weights. This leaves a single variance, ${\text{σ}}_{\text{ν}}^2$, to be estimated.

Figure 2. Histogram and Q-Q plot of $\mathbf{Z}$ for 16,850 SA1 areas (left panels). Histogram and Q-Q plot of standardised values (41) for the same SA1 areas (right panels).

Figure 2. Histogram and Q-Q plot of $\mathbf{Z}$ for 16,850 SA1 areas (left panels). Histogram and Q-Q plot of standardised values (41) for the same SA1 areas (right panels).

6.3. Measurement Error in MWI

Known sources of measurement error in the Census are due to non-response, invalid responses and random error purposely applied to the SA1 areas that have small populations (for confidentiality reasons; see [44] for further details). In the 2011 Census, item non-response for mean usual weekly income at the person level was approximately 7.9% [47].

Here, we model the non-response as “missing completely at random.” Then, for the i-th area, the observed $n_{ij}$ is the sum of Bernoulli random variables, $n_{ij}={\sum }_{k=1}^{N_{ij}}{J}_{ijk};j=1,...,9$, where $J_{ijk}$ is an indicator variable with $\text{Pr}(J_{ijk}=1)=p=0.921$, obtained from the item non-response, and $N_{ij}$ is the true number of people. From the randomness associated with $\{J_{ijk}:k=1,...,N_{ij}\}$, we calculate $\text{var}\left(\overline{Z}\left({\mathbf{s}}_i\right)\right|\left\{Z_{ijk}\right\})$ using the delta method. That is,

$$\begin{array}{cc}\hfill {\text{σ}}_{\epsilon }^2{V}_{\epsilon }\left({\mathbf{s}}_i\right)& ={\text{E}}_{\mathbf{n}}\left(\text{var}\left(\overline{Z}\left({\mathbf{s}}_i\right)\right|\left\{Z_{ijk}\right\})\right)+{\text{var}}_{\mathbf{n}}\left(\text{E}\left(\overline{Z}\left({\mathbf{s}}_i\right)\right|\left\{Z_{ijk}\right\})\right)\hfill \\ \mathrm{(43)}& & \simeq \left(\frac{N_{i}p+p-1}{N_i{(1-p+N_{i}p)}^2}+\frac{4p(1-p)}{{(1-p+N_{i}p)}^3}-\frac{1}{N_i^2}+\frac{2(1-p)}{N_i^{3}p}\right)\sum _{j=1}^9{N}_{ij}{\text{σ}}_j^2;i=1,...,16,850.\hfill \end{array}$$

Recall that $N_{ij}$ is the fully enumerated population total for area i and stratum j, and we define $N_i\equiv {\sum }_{j=1}^9{N}_{ij}$. To estimate (43), we use $n_{ij}{p}^{-1}$ for $N_{ij}$ and $n_i{p}^{-1}$ for $N_i$. Because the measurement-error calculations come from a separate study of item non-response, the values given by (43) are considered as known, and hence, ${\text{σ}}_{\epsilon }^2{\mathbf{V}}_{\epsilon }\equiv \text{diag}({\text{σ}}_{\epsilon }^2{V}_{\epsilon }\left({\mathbf{s}}_i\right):i=1,...,16,850)$ is known.

A summary statistic giving the measurement-error component of variability as a fraction of the total variability is:

$$\frac{{\text{σ}}_{\epsilon }^2{\sum }_{i=1}^{16850}{V}_{\epsilon }\left({\mathbf{s}}_i\right)}{{\sum }_{i=1}^{16850}{(\overline{Z}\left({\mathbf{s}}_i\right)-\overline{\overline{Z}})}^2}\text{,}$$

where $\overline{\overline{Z}}\equiv {\sum }_{i=1}^{16850}\overline{Z}\left({\mathbf{s}}_i\right)/16,850$. In our case, this is equal to 0.01%, which is very small.

6.4. Specifying the Spatial Dependence Matrix

The simple exploratory model we proposed in Section 6.2 does not account for spatial dependence between the SA1 areas. Now that we have captured the heteroskedasticity (Section 6.2) and measurement error (Section 6.3), we can build a SAR model (see Section 2.1) for the MWI values, $\left\{Y\left({\mathbf{s}}_i\right)\right\}$, defined for a fully enumerated population. That is, assume that $\mathbf{Y}\equiv {(Y\left({\mathbf{s}}_1\right),...,Y\left({\mathbf{s}}_{16850}\right))}^{\prime }$ has a covariance matrix given by (3), which is parametrised in terms of ρ and ${\text{σ}}_{\text{ν}}^2$.

Because there are $n=16,850$ small areas, fitting the SAR model is numerically challenging, so we turn to estimating the SAR parameters ρ and ${\text{σ}}_{\text{ν}}^2$ by fitting the SAR-calibrated SRE model given in Section 4. Recall that we choose ${\mathbf{U}}_1=({\mathbf{u}}_1,...,{\mathbf{u}}_r)$, to be the matrix containing the n-dimensional eigenvectors corresponding to the r largest absolute eigenvalues in (23), which are obtained from a spectral decomposition of the GMO (21) for a heteroskedastic SAR model.

The basis functions ${\mathbf{U}}_1$ and the SRE model parameters $\mathbf{K}$ and ${\text{σ}}_{\xi }^2$ all depend on the spatial dependence matrix $\mathbf{H}$ in (3). In NSW, there is substantial variation in the physical areas of the SA1s and in the spatial proximity of neighbouring areas. Along the eastern coastline and in the major urban centres, SA1s are reasonably uniformly distributed, whereas in other parts of the state, there is much greater asymmetry in the number and proximity of neighbouring areas.

A common choice for $\mathbf{H}$ is the symmetric adjacency matrix:

$$\mathbf{H}=\mathbf{A}\equiv \left(a_{i{i}^{\prime }}\right),$$

(44)

where $a_{i{i}^{\prime }}=1$ if the i-th area shares a common boundary with the $i^{\prime }$-th area; and $a_{i{i}^{\prime }}=0$, otherwise. However, an asymmetric $\mathbf{H}$ is more reasonable in this instance where distance or the shortest travel time between areas is more important than sharing a boundary. On a regular spatial lattice, the eight nearest neighbours are also the adjacent neighbours, but on an irregular lattice, they can be very different (see Figure 1b). Here, we model $\mathbf{H}$ based on the eight nearest neighbours of the irregular SA1 spatial lattice; that is,

$$\mathbf{H}=\mathbf{E}\equiv \left(e_{i{i}^{\prime }}\right),$$

(45)

where $e_{i{i}^{\prime }}=I(i^{\prime }\in E\left(i\right))$, for:

$$E\left(i\right)\equiv \{i^{\prime }:d_{i{i}^{\prime }}\le \text{eighth-largest}\text{distance}\text{from}i\}$$

and $d_{i{i}^{\prime }}$ is the distance between the respective centroids of the i-th and $i^{\prime }$-th areas,

Clearly, from Figure 1b, if $i^{\prime }\in E\left(i\right)$, it may not be true that $i\in E\left(i^{\prime }\right)$, and hence, $\mathbf{E}$ will generally be asymmetric. In our analysis, we fitted the SAR-calibrated SRE model using $\mathbf{H}$ equal to both (row-standardised) $\mathbf{A}$ and (row-standardised) $\mathbf{E}$, and we compared them.

6.5. Statistical Inference

In this section, we obtain parameter estimates for ρ, ${\text{σ}}_{\text{ν}}^2$, and β from the SAR-calibrated SRE model, using the methodology outlined in Section 5. Initially, we estimate the parameters for a small subset of the SA1 data obtained from the 1,230 SA1 areas that are located within the 2011 Census SA4 region of southwestern (SW) Sydney (Figure 3). We use these results to help select the number of basis functions, r, to use when modelling the full NSW dataset. With $n=1230$, the SAR model estimate ${\widehat{\text{ρ}}}_{SAR}$ for ρ can be computed, and hence, r is chosen so that $\widehat{\text{ρ}}$ (which depends on r) is close to the target value ${\widehat{\text{ρ}}}_{SAR}$.

We subsequently analyse the full NSW dataset using this selected value of r in a SAR-calibrated SRE model. Clearly, direct likelihood-based parameter estimates based on the SAR model are no longer possible for all of NSW, due to the size of the dataset ($n=16,850$).

Figure 3. SA1 level imputed MWI for SW Sydney.

Figure 3. SA1 level imputed MWI for SW Sydney.

While the basic heteroskedastic model for MWI in Section 6.2 fitted well, we cannot use it for inference, because its mean depends on the random counts $\left\{n_{ij}\right\}$ used to impute the data. Socio-economic and demographic covariates are generally used in an attempt to explain the variability in MWI. Using various selection criteria, we settled on the following explanatory variables: the proportion of males, out of the total population, in the age groups 15 to 24, 25 to 64 and ≥65 (three covariates), the proportion of females, out of the total population, in the same age groups (three covariates), the median mortgage repayment, the median rent, the average household size and the average number of persons per bedroom (four covariates) and the overall mean effect (one covariate). That is, we fitted a regression, $\mathbf{X}\mathbf{\beta }$, for the mean of the spatial model, where the number of regression parameters was $p=11$. The estimate of β is important and interpretable, and it may even be used by governments to make policy.

We specify the matrices that we use to represent heteroskedasticity, measurement error and spatial dependence in the model using (42), (43) and (44) or (45) from Section 6.2, Section 6.3 and Section 6.4, respectively. To complete our specification of the model, we use the term ${\text{σ}}_{\xi }^2{\mathbf{V}}_{\xi }$ to capture the additional variance not captured by the spatial basis functions included in the model. We expect the variance to depend on the number of observations in each area, and hence, we set $V_{\xi }\left({\mathbf{s}}_i\right)=c/n_i;i=1,...,16,850$, where $c=\text{trace}\left({\mathbf{V}}_{\text{ν}}\right)/({\sum }_{i=1}^{16850}1/n_i)$ is a scaling constant that makes ${\text{σ}}_{\xi }^2$ comparable to ${\text{σ}}_{\text{ν}}^2$.

We use the restricted-maximum-likelihood-estimation method described in Section 5 to obtain REML estimates, $\widehat{\text{ρ}}$ and ${\widehat{\text{σ}}}_{\text{ν}}^2$, for the SW Sydney and the NSW datasets, for the spatial dependence matrices $\mathbf{A}$ and $\mathbf{E}$ defined by (44) and (45), respectively. In each case, we used a grid-search algorithm to obtain $\widehat{\text{ρ}}$. For each value of ρ, we obtained a starting value, ${\text{σ}}_{\text{ν},0}^2\left(\text{ρ}\right)$ using (34) with $\mathbf{N}=\text{diag}(n_i:i=1,...,16,850)$, and then, we iteratively estimated ${\widehat{\text{σ}}}_{\text{ν}}^2\left(\text{ρ}\right)$. We calculated the profile restricted log-likelihood for each value of (ρ,${\widehat{\text{σ}}}_{\text{ν}}^2\left(\text{ρ}\right)$), from which we obtained $\widehat{\text{ρ}}\equiv argmin\left({\mathcal{L}}_R^{*}\left(\text{ρ}\right)\right)$, and hence, ${\widehat{\text{σ}}}_{\text{ν}}^2\equiv {\widehat{\text{σ}}}_{\text{ν}}^2\left(\widehat{\text{ρ}}\right)$. Finally, we used (37) to obtain generalised least-squares parameter estimates, $\widehat{\mathbf{\beta }}$.

For the SW Sydney dataset, we compared $\widehat{\text{ρ}}$ and ${\widehat{\text{σ}}}_{\text{ν}}^2$ with the target SAR model estimates. Hence, we use restricted maximum likelihood estimation and account for the measurement error to obtain SAR model parameter estimates, ${\widehat{\text{ρ}}}_{SAR}$ and ${\widehat{\text{σ}}}_{\text{ν},SAR}^2$.

Based on Figure 4 for SW Sydney, we see that for $r\ge 25$ spatial basis functions, the SAR-calibrated SRE model yields a 95% confidence interval for ρ that contains ${\widehat{\text{ρ}}}_{SAR}=0.8$ for $\mathbf{E}$ and ${\widehat{\text{ρ}}}_{SAR}=0.74$ for $\mathbf{A}$. From this and plots of SAR-calibrated estimates of ρ and ${\text{σ}}_{\text{ν}}^2$ versusr for all of NSW, we chose $r=25$ basis functions for our dimension-reduced analysis of MWI in the 16,850 SA1 areas of NSW shown in Figure 1. To choose which spatial dependence matrix to use for the NSW dataset, we compared the magnitude of the restricted log-likelihood functions for SW Sydney obtained using $\mathbf{A}$ and using $\mathbf{E}$, since the number of parameters is the same in both models. The asymmetric spatial dependence matrix $\mathbf{E}$ resulted in a smaller value for the negative restricted log-likelihood, so we used $\mathbf{E}$ to fit the model for NSW.

Figure 4. SAR-calibrated spatial random effects (SRE) model restricted (or residual) maximum likelihood (REML) estimates for $\widehat{\text{ρ}}$ and ${\widehat{\text{σ}}}_{\text{ν}}^2$ for SW Sydney as a function of r for $\mathbf{E}$ (left-hand panel) and for $\mathbf{A}$ (right-hand panel); the SAR model estimate ${\widehat{\text{ρ}}}_{SAR}$ and 95% profile log-likelihood confidence intervals (CI) for ρ are included for comparison (36).

Figure 4. SAR-calibrated spatial random effects (SRE) model restricted (or residual) maximum likelihood (REML) estimates for $\widehat{\text{ρ}}$ and ${\widehat{\text{σ}}}_{\text{ν}}^2$ for SW Sydney as a function of r for $\mathbf{E}$ (left-hand panel) and for $\mathbf{A}$ (right-hand panel); the SAR model estimate ${\widehat{\text{ρ}}}_{SAR}$ and 95% profile log-likelihood confidence intervals (CI) for ρ are included for comparison (36).

We fit the resulting SAR-calibrated SRE model given by (31) and (30) to the full NSW dataset, to obtain SAR-calibrated SRE model REML estimates, $\widehat{\text{ρ}}=0.74$ and ${\widehat{\text{σ}}}_{\text{ν}}^2\left(\widehat{\text{ρ}}\right)={\left(87\right)}^2$. We also obtain $\widehat{\text{ζ}}\equiv \text{ζ}\left(\widehat{\text{ρ}}\right)=0.87$, and hence, ${\widehat{\text{σ}}}_{\xi }^2={(100.3)}^2$. From these results, we conclude that MWI exhibits substantial spatial dependence in NSW, after allowing for the covariates included in the regression model.

Table 1 compares the SAR-calibrated SRE model estimates of β given by (37) with the ordinary least-squares (OLS) regression parameter estimates, ${\widehat{\mathbf{\beta }}}_{OLS}\equiv \mathbf{X}{\left({\mathbf{X}}^{\prime }\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }\mathbf{Z}$; the latter estimates ignore heteroskedasticity and spatial dependence. For all of the covariates except the intercept and “Female, 15 to 24 years”, the OLS estimates do not fall within the 95% CI for their respective β’s obtained using the SAR-calibrated SRE model. The OLS estimate ${\widehat{\mathbf{\beta }}}_{OLS}$ is based on the overly simple model, $\mathbf{Z}\sim \text{Gau}(\mathbf{X}\mathbf{\beta },{\text{σ}}^2\mathbf{I})$, and for this dataset, it is inadequate. This demonstrates the importance of capturing small-scale heteroskedastic spatial variation in our analysis of MWI.

Table 1. MWI for NSW. Ordinary least-squares regression-parameter estimates, generalised least-squares regression-parameter estimates obtained from fitting an SAR-calibrated SRE model and 95% confidence intervals (CIs). The SRE model is based on 25 basis functions and $\mathbf{H}=\mathbf{E}$ given by (45).

	${\text{β}}_{\mathbf{0}}$	${\text{β}}_{\mathbf{15}\mathbf{-}\mathbf{24}}^{\mathbf{Male}}$	${\text{β}}_{\mathbf{25}\mathbf{-}\mathbf{64}}^{\mathbf{Male}}$	${\text{β}}_{\mathbf{\ge }\mathbf{65}}^{\mathbf{Male}}$	${\text{β}}_{\mathbf{15}\mathbf{-}\mathbf{24}}^{\mathbf{Fem}\mathbf{.}}$	${\text{β}}_{\mathbf{25}\mathbf{-}\mathbf{64}}^{\mathbf{Fem}\mathbf{.}}$
${\widehat{\mathbf{\beta }}}_{SRE}$	463	47.5	773	−782	−201	675
${\text{CI}}_{SRE}$	(432,494)	(−6.7,102)	(732,814)	(−836,−729)	(−256,−146)	(622,728)
${\tilde{\mathbf{\beta }}}_{OLS}$	481	−95.6	719	−695	−177	789
	${\text{β}}_{\ge 65}^{Fem.}$	${\text{β}}_{Mort.}$	${\text{β}}_{Rent}$	${\text{β}}_{Size}$	${\text{β}}_{Beds}$
${\widehat{\mathbf{\beta }}}_{SRE}$	346	0.0443	0.298	−234	−15.1
${\text{CI}}_{SRE}$	(302,390)	(0.042,0.0466)	(0.287,0.309)	(−244,−224)	(−18.7,−11.3)
${\tilde{\mathbf{\beta }}}_{OLS}$	219	0.0472	0.197	−245	−7.84

Table 1. MWI for NSW. Ordinary least-squares regression-parameter estimates, generalised least-squares regression-parameter estimates obtained from fitting an SAR-calibrated SRE model and 95% confidence intervals (CIs). The SRE model is based on 25 basis functions and $\mathbf{H}=\mathbf{E}$ given by (45).

	${\text{β}}_{\mathbf{0}}$	${\text{β}}_{\mathbf{15}\mathbf{-}\mathbf{24}}^{\mathbf{Male}}$	${\text{β}}_{\mathbf{25}\mathbf{-}\mathbf{64}}^{\mathbf{Male}}$	${\text{β}}_{\mathbf{\ge }\mathbf{65}}^{\mathbf{Male}}$	${\text{β}}_{\mathbf{15}\mathbf{-}\mathbf{24}}^{\mathbf{Fem}\mathbf{.}}$	${\text{β}}_{\mathbf{25}\mathbf{-}\mathbf{64}}^{\mathbf{Fem}\mathbf{.}}$
${\widehat{\mathbf{\beta }}}_{SRE}$	463	47.5	773	−782	−201	675
${\text{CI}}_{SRE}$	(432,494)	(−6.7,102)	(732,814)	(−836,−729)	(−256,−146)	(622,728)
${\tilde{\mathbf{\beta }}}_{OLS}$	481	−95.6	719	−695	−177	789
	${\text{β}}_{\ge 65}^{Fem.}$	${\text{β}}_{Mort.}$	${\text{β}}_{Rent}$	${\text{β}}_{Size}$	${\text{β}}_{Beds}$
${\widehat{\mathbf{\beta }}}_{SRE}$	346	0.0443	0.298	−234	−15.1
${\text{CI}}_{SRE}$	(302,390)	(0.042,0.0466)	(0.287,0.309)	(−244,−224)	(−18.7,−11.3)
${\tilde{\mathbf{\beta }}}_{OLS}$	219	0.0472	0.197	−245	−7.84

7. Discussion and Conclusions

This paper has considered the specification and estimation of a hierarchical model that includes a computationally-efficient spatial random effects (SRE) model for a latent spatial process and an additive measurement error model for the observed data. Using the generalised Moran operator, we calibrate this model to an SAR model, a model that is widely used in spatial econometrics, but that is computationally prohibitive for large datasets. We use restricted maximum likelihood on a computationally-efficient SRE model to estimate the SAR model’s covariance parameters.

Implementation of our model is demonstrated using Australian Census data for mean income. It is clear from the SAR-calibrated SRE model fitted to NSW’s 16,850 SA1 areas that a regression model with the 11 covariates and independent residuals does not properly capture the small-scale variability. This (heteroskedastic) spatial variability is important, which is reflected in the estimate for ρ, namely $\widehat{\text{ρ}}=0.74$ for 25 basis functions. It is also reflected in the change of regression coefficient estimates when the spatial variability is used to improve both the estimates and their estimation variance.

Notice the role of model calibration here, which is different from model fitting. The SRE model is calibrated to be close to the SAR model of interest, and then, the computationally-efficient SAR-calibrated SRE model is fitted. Importantly, the parameters that are fitted are those of the original, interpretable SAR model. This enables SAR model parameters to be estimated for much larger datasets than is currently possible.

It is easy to see that our SAR-calibrated SRE model can be adapted to obtain a CAR-calibrated SRE model. In some sense, this represents the converse of the calibration undertaken by Rue and Tjelmeland [48], who find the Gaussian–Markov random field (GMRF) that comes closest to a geostatistical model; note that the CAR model is a GMRF and that the SRE model is a geostatistical model.

Besides the advantages of the SAR-calibrated SRE model in producing parameter estimates, the r spatial basis functions in $\mathbf{S}$ could be used to identify hidden, possibly overlapping, regions or classes of areas. The associated random effects can be easily predicted, and hence, regions or classes having higher or lower values of MWI could also be predicted. Moreover, as the covariance matrix $\mathbf{K}$ of the random effects is not constrained to be diagonal, the SRE model captures additional dependence associated with these classes.