Likelihood estimation and inference for the autologistic model.

The autologistic model is commonly used to model spatial binary data on the lattice. However, if the lattice size is too large, then exact calculation of its normalizing constant poses a major difficulty. Various different methods for estimation of model parameters, such as pseudo-likelihood, have been proposed to overcome this problem. This article presents a method to estimate the normalizing constant in an efficient manner. In particular, this allows tasks such as maximum likelihood estimation and inference for model parameters. We also consider the true likelihood approximated by the product of likelihoods for which the normalizing constant can be found by an analytic computational method by wrapping the lattice on the cylinder. This gives a simulation-free method of inference. We compare estimates of model parameters based on our new methods with the commonly used pseudo-likelihood approach. Although we have not considered Bayesian inferences here, the method can be straightforwardly extended to find posterior distributions. We apply our methods to the well-known endive data and to simulated data and find that our methods give substantially increased accuracy of estimation of model parameters.

Key Words: Autologistic distribution; Gibbs distribution; Ising model; Maximum likelihood; Markov chain Monte Carlo; Pseudo-likelihood.

1. INTRODUCTION

The autologistic model was first proposed by Besag (1972, 1974) and is widely used to model binary spatial data. A major drawback with its use, however, is that the normalizing constant is generally unknown analytically. In particular, this makes tasks such as maximum likelihood estimation impossible to carry out; see, for example, Huang and Ogata (2001). Many different methods have been introduced to overcome this problem by introducing approximate methods, for example, based on estimating equations, pseudo-likelihood or coding (Besag 1986). In fact pseudo-likelihood has been applied to a wide variety of spatial data (see, e.g., Augustin, Mugglestone, and Buckland 1996; Wu and Huffer 1997). Geyer and Thompson (1992) provided a method for finding maximum likelihood estimates--Monte Carlo maximum likelihood--but this method, based on importance sampling ideas, is potentially very slow numerically and unstable due to calculation of the density and therefore exponentiation of large-valued statistics. The log-density is used in our methods, which are therefore numerically stable.

Pettitt and Friel (2001) showed how the normalizing constant for the autologistic model can be found by a computational analytic method if the lattice is wrapped onto the cylinder. But this method is restricted to lattices of size m x n, (m [less than or equal to] n), with m [less than or equal to] 10, as the method depends on finding eigenvalues of a [2.sup.m] x [2.sup.m] matrix. Pettitt and Friel (2001) showed how it is possible to extend this result on the cylinder lattice to find the normalizing constant for the free boundary lattice using a Markov chain Monte Carlo scheme known as path sampling (Gelman and Meng 1998). This article shows that it is possible to extend these results to larger, arbitrary sized lattices. Further, we use this methodology to find the maximum likelihood estimate or posterior mode of model parameters for larger lattices and illustrate the method with the well known 14 x 179 lattice involving the endive data (Besag 1978) and simulated data.

For the endive data, we additionally find estimates and confidence intervals based on an approximation of the likelihood by approximating the normalizing constant as the product of two normalizing constants based on two 7 x 179 lattices found by splitting the original data into two sub-lattices. For each sub-lattice the normalizing constant is approximated by the corresponding normalizing constant for the cylinder. For the normalizing constants calculation we effectively treat the two sub-lattices as independent whereas the unnormalized distribution is taken to be its true value. We compare our two new methods of inference for the parameters with the well-known pseudo-likelihood approach. For the endive data we find good agreement for values of estimates and confidence intervals for our new methods, while the value of the pseudo-likelihood estimate is statistically significantly different from the maximum likelihood estimate. For simulated autologistic data there is a substantial difference between our new methods and the pseudo-likelihood approach. Our simulation results suggest that the true likelihood method and our good approximation to it are substantially more efficient than the pseudo-likelihood method. We find that the efficiency (defined as the ratio of the mean square error of the method and that of the true likelihood estimate) of the pseudo-likelihood estimate is typically 15% and our new approximate estimate about 60% efficient for the range of cases considered. These findings are also supported by Gu and Zhu (2001), who derived maximum likelihood estimates for the Ising model using an MCMC based Newton-Raphson algorithm. Additionally, our approximate methods can be used to find posterior distributions of parameters without recourse to simulation and only using MCMC to refine the normalizing constant approximation. However, here we illustrate the technique by finding confidence intervals for parameters based on the profile likelihood.

Section 2 briefly introduces the autologistic model and pseudo-likelihood estimation. Section 3 outlines the approach to calculate the normalizing constant for the lattice wrapped on the cylinder, and then shows how this result, when used with an MCMC scheme, may be applied to calculate the normalizing constant for the free boundary lattice. Section 4 presents results of our maximum likelihood estimation procedures, together with confidence intervals based on profile likelihoods. Finally, we present conclusions to this work in Section 5.

2. AUTOLOGISTIC MODEL

Consider a binary random variable [x.sub.ij] taking the values {0, 1} at each site (i, j) on a regular m x n lattice. The unnormalized autologistic distribution on the lattice may be written in exponential form as

(2.1) q(x|[THETA])=exp{[[THETA].sup.T]V(x)} =exp{[[theta].sub.0][V.sub.0](x)+[[theta].sub.f][V.sub.f](x)}.

Here [THETA]=([[theta].sub.0], [[theta].sub.f]) and V(x)=([V.sub.0](x), [V.sub.f](x)), with

(2.2) [V.sub.0](x)=[m.summation over i=1][n.summation over j=1][x.sub.ij]

(2.3) [V.sub.f](x)=[m-1.summation over i=1][n.summation over j=1]I [[x.sub.ij]=[x.sub.i+1,j]]+[m.summation over i=1][n-1.summation over j=1] I[[x.sub.ij]=[x.sub.i,j+1]].

Here I[[x.sub.ij]=[x.sub.i+1,j]] is the usual indicator function, taking the value 1 when [x.sub.ij]=[x.sub.i+1,j], and 0 otherwise. It is seen that [V.sub.0](x) is the number of 1's in the lattice. The first term of [V.sub.f](x) counts the number of like-valued nearest neighbors within each column. The second term counts a similar …