Marginal maximum likelihood estimation of SAR models with missing data

Abstract

Maximum likelihood (ML) estimation of simultaneous autocorrelation models is well known. Under the presence of missing data, estimation is not straightforward, due to the implied dependence of all units. The EM algorithm is the standard approach to accomplish ML estimation in this case. An alternative approach is considered, the method of maximising the marginal likelihood. At first glance the method is computationally complex due to inversion of large matrices that are of the same size as the complete data, but these can be avoided, leading to an algorithm that is usually much faster than the EM algorithm and without typical EM convergence issues. Another approximate method is also proposed that serves as an alternative, for example when the contiguity matrix is dense. The methods are illustrated using a well known data set on house prices with 25,357 units.

Introduction

Simultaneous autoregressive models (SAR) are popular linear regression models for spatially distributed data that take the dependence of the response variable of neighbouring units into account. Maximum likelihood (ML) estimation for SAR models is well established (Ord, 1975). The dependence of neighbouring units is represented by the contiguity matrix $\mathbf{W}$, which has non-zero entries $W_{ij}$ if $i$ is close to $j$. Various choices of how to define $\mathbf{W}$ are available. The $\mathbf{W}$ of size $n\times n$ matrix refers to the $n$ units of interest. If complete data are observed for a response variable $Y$ of interest, then ML estimation can be accomplished relatively effectively following the methods proposed by Ord (1975). However when missing data are present, i.e. the response variable $Y$ is only observed for a subset of $n_s$ units with $n_s<n$, estimation becomes more difficult as the contiguity matrix refers to $n$ units while $Y$ refers to $n_s$ units. ML estimation for this situation has been investigated by Lesage and Pace (2004) using the EM algorithm.

Kato (2013) considered estimation and prediction of several spatial covariance models, including SAR models, using the EM algorithm and also a quasi-likelihood method for estimation. Goulard et al. (2017) used the EM algorithm to investigate the performance of various in-sample and out-of-sample predictors, but also used an ML algorithm. Computational improvements of the EM algorithm based on Lesage and Pace (2004) have been considered by Suesse and Zammit Mangion (2017).

Other estimation methods for SAR models under the presence of missing data have been extensively studied using, for example, the generalised method of moments and least squares (Wang and Lee, 2013) and approximate Bayesian methods using Integrated Laplace (INLA) approximation Bivand et al. (2014), Gomez-Rubio et al. (2015), Gomez-Rubio et al. (2017).

In this paper I follow a different approach for ML estimation, by maximising the marginal likelihood directly instead of maximising it indirectly via the EM algorithm, leading to a generally much faster algorithm by applying some computational tricks. The proposed method is a true marginal ML method and requires out-of-sample information, which is different from the ML methods used by Kato (2013) and Goulard et al. (2017) that only use in-sample information. The proposed algorithm also does not suffer from convergence problems, that are common and well known for the EM algorithm (McLachlan and Krishnan, 2008), in particular when the sample size is very small.

Both the EM algorithm and the marginal ML approach generally lead to the same ML estimates, provided convergence to a global maximum has been achieved. In analogy Laird and Ware (1982) considered ML estimation via the EM algorithm for linear mixed models, whereas Lindstrom and Bates (1988) considered maximising the marginal likelihood. Both are also common techniques for restricted ML estimation, a method to achieve unbiased variance estimates.

In Section 2, the spatial models under consideration are introduced. Then in Section 3 standard ML estimation of complete data is reviewed. In Section 4 I introduce an exact method of maximising the marginal likelihood and an alternative approximate method. Section 4 presents a simulation study to investigate the performance of the various methods and Section 5 applies the methods to a well known data set on house prices by also illustrating typical computation times and convergence issues of the EM algorithm for small sampling ratios. The article concludes with a discussion.