REML estimation with intrinsic Matérn dependence in the spatial linear mixed model

Abstract

We present a new matrix-free residual maximum likelihood (REML) analysis for irregularly spaced spatial data, where observations usually represent average values over very small regions that are interpreted as points. The REML analysis is obtained after embedding the sampling locations in a fine scale rectangular lattice, treating unobserved sites as missing data. The spatial random fields considered here are based on fractional Laplacian differencing on the lattice and they are unique in approximating continuum intrinsic Matérn dependence. Here, using the h-likelihood method, we derive REML estimating equations that allow for singular precision matrices, estimation of covariate effects, prediction of unobserved spatial effects and REML estimation of precision parameters as a solution to an explicit gamma non-linear model. Furthermore, we devise a sophisticated computational algorithm that enables us to achieve scalable matrix-free statistical computations. In particular, these matrix-free computations include the use of (1) the two-dimensional discrete cosine transformation that arises in the spectral decomposition of the precision matrix of our spatial random fields and that allows fast matrix-free matrix-vector multiplication, (2) a matrix-free pre-conditioned Lanczos algorithm that solves non-sparse matrix equations with linear constraints, (3) a matrix-free Hutchinson’s trace estimator that stochastically approximates the trace of a matrix, (4) a robust trust region method that always finds a local maximum of the non-concave residual log-likelihood function and (5) some preliminary computations of the log REML likelihood function based on Taylor series approximation. Using computer experiments, we provide further understanding on not just the number and values but also the basins of attraction of the local and global maxima of the REML function. This understanding significantly simplifies the problem of finding global maxima. We further demonstrate through computer experiments that our matrix-free REML estimators attain both efficiency and geostatistical inference, and surpass the widely used INLA methods in computational times. We provide an extensive application on mapping ground water arsenic concentration in Bangladesh, indicating numeric consistency of results and robustness of inference to changes of lattice spacing. The paper closes with some discussions that include computations in the stationary case, conditional simulations and matrix-free MCMC computations.