Estimation of spatial autoregressive models with measurement error for large data sets

Abstract

Maximum likelihood (ML) estimation of spatial autoregressive models for large spatial data sets is well established by making use of the commonly sparse nature of the contiguity matrix on which spatial dependence is built. Adding a measurement error that naturally separates the spatial process from the measurement error process are not well established in the literature, however, and ML estimation of such models to large data sets is challenging. Recently a reduced rank approach was suggested which re-expresses and approximates such a model as a spatial random effects model (SRE) in order to achieve fast fitting of large data sets by fitting the corresponding SRE. In this paper we propose a fast and exact method to accomplish ML estimation and restricted ML estimation of complexity of $O(n^{3/2})$ operations when the contiguity matrix is based on a local neighbourhood. The methods are illustrated using the well known data set on house prices in Lucas County in Ohio.

Appendices

A. Simulation results for the SEM-M

See Tables 7, 8, 9, and 10.

Table 7 Empirical mean of $\rho$, $\sigma ^2_y$, $\sigma ^2_{\epsilon }$ and $\beta$ estimates over 10,000 simulations for the SEM-M model with $n=506$, $\sigma ^2_y=1$, $\sigma ^2_{\epsilon }=2$ and $\beta _0=1$ and $\beta _1=5$

Table 8 Empirical mean of $\rho$, $\sigma ^2_y$, $\sigma ^2_{\epsilon }$ and $\beta$ estimates over 10,000 simulations excluding those for which $\hat{\sigma }^2_{\epsilon }=0$ for the SEM-M model with $n=506$, $\sigma _y=1$, $\sigma ^2_{\epsilon }=2$ and $\beta _0=1$ and $\beta _1=5$

Table 9 Coverage of Wald-type CI for $\rho$, $\sigma ^2_y$, $\sigma ^2_{\epsilon }$ and $\beta$ estimates over 10,000 simulations for the SEM-M model with $n=506$, $\sigma ^2_y=1$, $\sigma ^2_{\epsilon }=2$ and $\beta _0=1$ and $\beta _1=5$

Table 10 Coverage of Wald-type CI for $\rho$, $\sigma ^2_y$, $\sigma ^2_{\epsilon }$ and $\beta$ estimates over 10,000 simulations excluding those for which $\hat{\sigma }^2_{\epsilon }=0$ for the SEM-M model with $n=506$, $\sigma ^2_y=1$, $\sigma ^2_{\epsilon }=2$ and $\beta _0=1$ and $\beta _1=5$

B. Reduced rank approach

To fit the SEM-M, Burden et al. (2015) considered a reduced rank approach that relies on a spectral decomposition. Once this has been accomplished the computational complexity is reduced to that of a SRE model. The spectral decomposition of the matrix

$$\begin{aligned} \mathbf {B}= \frac{ \mathbf {P}\left( \mathbf {W} + \mathbf {W}^{\top }\right) \mathbf {P}}{2}= \mathbf {U}\mathbf {\Omega }\mathbf {U}, \end{aligned}$$

(12)

see formula (22) in Burden et al. (2015), is computationally intensive. The matrix $\mathbf {P}$ is the projection matrix defined in Sect. 3.2. The matrix $\mathbf {U}$ contains the eigenvectors of $\mathbf {B}$ and $\mathbf {\Omega }$ is a diagonal matrix with the eigenvalues of $\mathbf {B}$ on its diagonal. Notice that the matrix $\mathbf {B}$ is a first order approximation of the (co)variance matrix that is obtained from applying REML, ignoring second order terms. The decomposition can be written as

$$\begin{aligned} \mathbf {B}=\mathbf {U}_1\mathbf {\Omega }_1\mathbf {U}_1 + \mathbf {U}_2\mathbf {\Omega }_2\mathbf {U}_2, \end{aligned}$$

where now $\mathbf {\Omega }_1$ contains the largest r absolute eigenvalues and $\mathbf {\Omega }_2$ the remaining $n-r$, and where $\mathbf {U}_j$ is the corresponding matrix containing the eigenvectors. The authors suggested then to use $\mathbf {U}_1\mathbf {\Omega }_1\mathbf {U}_1$ as an approximation to $\mathbf {B}$ and used this approximation to base estimation of the SEM-M on estimation of SREs.

Fig. 4

Boxplot showing the distribution of the eigenvalues of $\mathbf {B}$ for the house data set

We use the house data set with $n=25{,}357$, see Sect. 5, to illustrate the preprocessing times of the reduced rank approach. Instead of calculating all or a large number of eigenvalues we used Matlab and the R interface of the ARPACK software (Lehoucq et al. 2015) to calculate only the largest absolute $r=25$ eigenvalues and corresponding eigenvectors. We would expect this process to take only a fraction of the time compared to calculating all $n=25{,}357$ eigenvalues, presented in Fig. 4, which took 1 h and 36 minutes (min) using Matlab. However this process of calculating 25 eigenvalues took 54 min 40 s (s) with Matlab (using command ‘eigs’ and option ‘largest magnitude’) and 54 min 32 s with ARPACK (also option ‘largest magnitude’), approximately two thirds of the time even though only $\approx 1/1000$ (25 / 25, 357) of the eigenvalues and vectors were calculated. Also the matrix $\mathbf {B}$ requires storage. In fact 4.9 GB of memory is required with R indicating that the calculation of $\mathbf {B}$ is limited by the available hardware. The objects of the proposed approach of this paper do not require much memory, for example $\mathbf {W}$ or the Cholesky decomposition of $\mathbf {A}\mathbf {A}^{\top }$ only require a couple of (up to 6) Megabytes.

The 54 min serves as a very optimistic lower bound for the time needed to fit this data set with the reduced rank approach. The SRE also has to be fitted, meaning that the total calculation time will clearly exceed this value. Figure 4 in Burden et al. (2015) demonstrates the effect of choosing different values of $r\le 100$ on $\hat{\rho }$, showing the approximate nature of the method for various values of r.