R. Kelley Pace

LREC Endowed Chair of Real Estate

Department of Finance

E.J. Ourso College of Business Administration

Louisiana State University

Baton Rouge, LA 70803-6308

(225)-388-6256

FAX: (225)-334-1227

www.spatial-statistics.com

kelley@pace.am

and

Dongya Zou

M.S. Candidate

Department of Computer Science

Louisiana State University

Baton Rouge, LA 70803

(225)-388-6238

FAX: (225)-334-1227

Abstract

Models of n² potential spatial dependencies among n observations spread irregularly over space seem unlikely to yield simple structure. However, the use of the nearest neighbor leads to a very parsimonious eigenstructure of the associated adjacency matrix which results in an extremely simple closed-form for the log-determinant. In turn, this leads to a closed form solution for the maximum likelihood estimates of the spatially autoregressive and mixed regressive spatially autoregressive models. With the closed-form solution, one can find the neighbors and compute maximum likelihood estimates for 100,000 observations in under one minute. The model has theoretical, pedagogical, diagnostic, modeling, and methodological applications. For example, the model could serve as a more enlightened null hypothesis for geographic data than spatial independence.

Keywords: spatial statistics, spatial autoregression, nearest neighbor, maximum likelihood, sparse matrices, log-determinants

I. Introduction

The task of modeling the n² potential dependencies among n spatially scattered observations has generated a substantial amount of intellectual activity. Naturally, many possible ways exist to specify the spatial interactions. Virtually all distance-based weighting methods such as variograms or spatial weight matrices assign some positive influence to the nearest observation, the nearest neighbor. In this sense, the nearest neighbor is common to most methods.

Given the ubiquitous nature of nearest neighbors, it seems worthwhile to examine a limiting model of spatial dependence, a model where observations depend only upon their nearest neighbor. Such special cases often help illustrate the more general case.

The nearest neighbor graph exhibits a number of important regularities. These regularities ensure the associated adjacency matrix (spatial weight matrix in this case) has a remarkable structure — the eigenvalues are all –1,0, or 1 and the number of 1,-1 pairs equals the number of pairs of symmetric elements in the adjacency matrix. This leads to a simple closed form for , where represents a spatial autoregressive parameter and D represents the n by n nearest neighbor adjacency matrix. Specifically, , where K represents the number of symmetric pairs of elements in the adjacency matrix. While such simple forms for the exact eigenvalues and log-determinants have been obtained before (e.g., Ord (1975), Gasim (1988), Griffith and Sone (1995), Martin (1993)), such results pertain to various patterns based upon regular sets of points or to very small numbers of points. The nearest neighbor eigenvalue results apply to arbitrary sets of points scattered across one, two, or more dimensions.

This simple form for the log-determinant leads to a closed-form solution for the spatially autoregressive and mixed regressive spatially autoregressive maximum likelihood estimators. In fact, we obtain a closed form not just for the log-determinant of the nearest neighbor adjacency matrix but also the log-determinant of any linear combination of the odd powers of this matrix. Finding the nearest neighbors and using them in estimation requires approximately O(nlog(n)) operations.

We illustrate the operation of the nearest neighbor model with an application to the data set analyzed by Pace and Barry (1997). By modeling just one neighbor, we obtain estimates whose log-likelihood is approximately halfway between those of maximum likelihood using zero neighbors (OLS) and maximum likelihood using four neighbors.

In addition, we conduct a Monte Carlo investigation of the sensitivity of the estimator to a simple form of zonal anisotropy. As expected, the anisotropy reduces the benefits from using the estimator. Nevertheless, the estimator still yields substantial improvements over the use of OLS even in the presence of this form of isotropy.

The simple form of the estimator can help clarify more complicated spatial dependence, aid in theoretical investigations, and serve as a useful pedagogical tool. In addition, the nearest neighbor maximum likelihood estimate of spatial dependence could have diagnostic value. Finally, it can provide a naïve model useful for gauging the contribution of more complicated approaches while still keeping within the useful maximum likelihood framework.

Section II discusses the regularities in the nearest neighbor graph, section III uncovers the eigenvalues and log-determinant of the associated adjacency matrix, section IV shows how these lead to a closed-form solution to the maximum likelihood first order conditions, section V provides a more general model with a closed-form maximum likelihood solution, section VI applies the model to an actual data set, section VII examines the performance of the estimator in the presence of misspecified dependence and anisotropy, section VIII discusses the time required for computing the estimates, and section IX concludes with a discussion of the key results.

II. The Nearest Neighbor Graph

Given n points or nodes, we consider their nearest neighbor graph. We make the following assumption:

This Assumption eliminates the situation depicted in Diagram 1, where the arrows point to the nearest neighbor. Hence, nodes c and a cannot simultaneously be the nearest neighbors of node b. Given Assumption 1, the results developed below apply only to irregularly spaced sets of nodes.

(Diagram 1)

(Diagram 2)

(Diagram 3)

(Diagram 4)

(Diagram 5)

(Diagram 6)

(Diagram 7)

As an additional illustration, Diagram 8 presents a more complicated individual subgraph. The diagram shows how multiple nodes can have an individual node as nearest neighbors but not vice-versa. Note, the existence of only one 2-node-cycle between nodes a and b.

(Diagram 8)

To give a better idea of the effects of the three propositions, we present a diagram of the United States counties and their nearest neighbor graph (without arrows) as depicted in Diagram 9. As expected from the Propositions stated above, all the subgraphs are disconnected, and each one has exactly one 2-node-cycle and no simple cycle (cycles containing three or more nodes). Intriguingly, the nearest neighbor graph leads to a potential mathematical description of regions.

III. The Nearest Neighbor Adjacency Matrix and its Eigenvalues

Given the nearest neighbor graph discussed in section II, we now focus upon the associated adjacency matrix and its properties. We show the adjacency matrix has particularly simple eigenvalues and an extremely simple closed form determinant formula. Insofar as the matrix has no negative values, has only zeros on the main diagonal, and each row contains a single 1, this adjacency matrix is also a row-stochastic spatial weight matrix and hence the log-determinant of this has immediate relevance in computing maximum likelihood estimates, a topic explored in section IV.

We begin by examining Proposition 4 which clarifies the structure of the adjacency matrix.

Theorem 1 has many matrix function applications. The most obvious application leads to a closed-form solution for the determinant of functions involving the nearest neighbor adjacency matrix.

Other applications of Theorem 1 are possible. For example, Bavaud (1998) discussed the possible importance of matrix functions of row-stochastic weight matrices. Some of these functions of the nearest neighbor weight matrix also have simple eigenvalues or determinants. As a limitation, attempts to create symmetry by adding the nearest neighbor adjacency matrix and its transpose destroy the simplifying structure imposed by using the nearest neighbor and do not yield simple eigenstructure.

IV. A Spatial Maximum Likelihood Estimator with Closed Form Solution

Computation of the spatial maximum likelihood estimators has traditionally required iterative, sometimes arduous, computations. The foregoing closed-form expression for the log-determinant enables the spatially autoregressive model to have a closed-form solution as well.

Let Y represent an n element vector of observations on the dependent variable, X represent an n by p matrix of observations on the independent variable (X may include functions of the spatial weight matrix), let D represent the n by n nearest neighbor adjacency matrix as specified above, let represent a p element vector of parameters associated with the independent variables, and let represent a scalar autoregressive parameter. We fit the model . The profile log-likelihood function is,

where C represents a scalar constant. Let , ,, . Following Anselin (1988, p. 181) or Pace and Barry (1997, p. 235-236), we can write the sum-of-squared errors as a function of the autoregressive parameter as in .

Substituting the expressions for from and from into the profile log-likelihood equation yields ,

where K represents the number of symmetric pairs in the adjacency matrix. We wish to find its maximum point with respect to the a parameter in (-1,1). We assume SSE>0, and given this we calculate the first derivative of L which appears in .

Setting and substituting the representation for SSE from , we obtain .

We assume that not all points are symmetric (i.e., ). If , the cubic term vanishes and the first order condition in becomes a simpler quadratic equation. Assuming , we can rewrite the coefficients from as in .

and rewrite the cubic equation from using the coefficients from as .

We define in the following constants used for solving the cubic equation.

Finally, the solutions appear in .

By the properties of cubic equations, at least one of these solutions is real.

Exactly one of the solutions from falls in (-1, 1). To see this, in we take the second derivative of L with respect to a .

The first and second terms within the brackets are unambiguously positive. The third term (the term in braces) is identical to the second term of the first order conditions in . As the first term in the first order conditions is negative, it follows the only other term (the second term) must be non-negative at the optimum. The first two terms in brackets are positive and the third term is non-negative and thus the expression in brackets is positive. The positive expression within brackets is multiplied by a negative term in parentheses. Thus, for and hence is strictly concave. Since , there exists at least one interior maxima. By the strict concavity of , this maxima will be the maximum. Hence, a value of from (-1, 1) uniquely maximizes the profile log-likelihood.

V. A More General Spatial Maximum Likelihood Estimator

We consider a generalization of the model in section IV. We now find the closed-form solution for the autoregressive parameters associated with the odd powers of the nearest neighbor adjacency matrix. The odd powers of the nearest neighbor adjacency matrix have the desirable feature of zeros along the diagonals (unlike the even powers). Provided the linear combination of the odd powers have coefficients summing to unity, the resulting linear combination will have the same set of symmetric elements as the first-order nearest neighbor adjacency matrix, but will have different weights for the asymmetric elements. Because they share the same set of symmetric elements (same set of 2-node-cycles), they will have the same set of eigenvalues, as shown in Application 2 of Theorem 1. Hence, the linear combination provides a way of incorporating the effects of more distant observations without complicating the computation of the log-determinant. Thus, this generalization can improve the richness of the nearest neighbor model of spatial dependence at little additional computational cost.

Now we would like to consider the maximum likelihood function in ,

where ,, , , and . We assume E has full rank. Substituting the relation from into yields .

To find the maximum of L we proceed in two stages. First, we maximize L with respect to subject to holding constant where . Holding constant makes maximizing the log-likelihood equivalent to minimizing a quadratic form subject to linear equality constraints. We substitute the optimal into L to form the profile . Second, we maximize the resulting profile log-likelihood with respect to .

To begin the first stage of optimization, let

be the Lagrange function, where. Taking first derivatives yields .

From the second equation, we get ; from this and the first one we get . Substituting this back into the first equation yields the optimal solution as a function of in .

The matrix is a variance-covariance matrix and is positive definite (we assumed E had full rank). Hence, this solution should minimize the sum-of-squared errors subject to the constraints and thus maximize the log-likelihood subject to the equality constraint for a given value of . Substituting the optimal into the log-likelihood in yields the profile log-likelihood in .

Let , a 2 by 2 matrix. We can rewrite the profile log-likelihood in as .

This is exactly the same problem as analyzed in section IV with serving the same role as , serving the same role as , and serving these same role as . Hence, we can use the solution developed there for the second stage of the optimization of the profile log-likelihood.

VI. Application of the Estimator to Geographical Data

To illustrate the effects of modeling spatial data with the nearest neighbor as opposed to multiple neighbors, we fitted the model and data introduced by Pace and Barry (1997). They examined voting behavior for the 1980 presidential election across 3,107 counties in the continental US using a mixed regressive spatially autoregressive estimator (Ord (1975), Anselin (1988)). Table 1 reproduces their OLS and maximum likelihood estimates using the four nearest neighbors () and augments these with the estimates using the nearest neighbor () and the nearest neighbor with the addition of a third-order nearest neighbor (). We employed the Euclidean metric using unprojected latitude and longitude to determine the neighbors. Table 1 only shows the estimates for the independent variables common to all estimators. For brevity, it does not show the results for the numerous spatially lagged independent variables. We report likelihood ratio statistics (i.e., twice the difference between the restricted and unrestricted log-likelihoods) for the null hypothesis that a variable has no effect or in many cases that a variable and its spatial lag have no effect.

The third-order neighbor corresponds to using the nearest neighbor of the nearest neighbor of the nearest neighbor (i.e., ). For symmetric elements, the third-order neighbor adjacency matrix has the same entries as the first-order nearest neighbor adjacency matrix (in Diagram 8 the nearest neighbor to node b is node a, its nearest neighbor is node b, and finally its nearest neighbor is back to node a). For asymmetric elements, the third-order neighbor adjacency matrix elements are not the same as the first-order neighbor adjacency matrix elements (in Diagram 8 the nearest neighbor to node e is node d, its nearest neighbor is node c, and the third-order nearest neighbor is node b). For graphs with many asymmetric nodes, the use of the third-order or higher odd-order adjacency matrices could capture some of the more complicated spatial dependence. See Anselin and Smirnov (1996) for more on higher-order weight matrices.

As one might expect, the nearest neighbor maximum likelihood estimates usually fell between the OLS and maximum likelihood using the four nearest neighbors. The log-likelihood was –6307.1 for OLS without lagged independent variables, -5964.6 for maximum likelihood with the nearest neighbor, -5954.9 for maximum likelihood with the nearest neighbor and the third-order neighbor, and –5678.8 for maximum likelihood with four nearest neighbors. As expected, maximum likelihood with four neighbors significantly outperformed the other specifications. However, the nearest neighbor maximum likelihood estimator clearly outperformed OLS. Its log-likelihood fell about halfway between that of OLS and the four nearest neighbor maximum likelihood estimators. The addition of the third-order neighbor increased the log-likelihood significantly, but did not materially change the estimates or the goodness-of-fit.

Figure 1 displays the profile log-likelihood (actually twice the log-likelihood) as a function of the autoregressive parameter for both the unrestricted model and for the relevant restricted models. The restricted models include one without an intercept (1 restriction), one without the population over 18 years of age and its spatial lag (2 restrictions), one without the population with more than 12 years of education and its spatial lag (2 restrictions), one without the owner-occupied housing units and its spatial lag (2 restrictions), as well as one without aggregate income and its spatial lag (2 restrictions). For each profile log-likelihood curve, an asterisk denotes the producing the maximum log-likelihood for that model. As expected, all of the curves have clear concave shapes without any visible irregularities.

The profile log-likelihood curves provide a wealth of information on the importance of the variables, the autoregressive parameter, and their interactions. First, variables whose deletion (restriction to zero effect) results in a curve substantially below the unrestricted curve are clearly statistically significant. Under the null hypothesis of no effect, single restrictions (i.e., no intercept) should have a likelihood ratio statistic (twice the difference between the restricted and unrestricted log-likelihood) of less than 4 at approximately the 5% significance level and of less than 6 at approximately the 1% significance level. Similarly under the null hypothesis of no effect, dual restrictions (e.g., no effect of income and spatially lagged income) should have a likelihood ratio statistic (twice the difference between the restricted and unrestricted log-likelihood) of less than 6 at approximately the 5% significance level and of less than 10 at approximately the 1% significance level. All of the variables in the model have likelihood ratio statistics exceeding these critical values. Second, the curvature of each curve conveys information concerning the estimation error associated with the autoregressive parameter. Third, the shift of the autoregressive parameter across the curves reveals the extent to which the autocorrelation adjustment can substitute for the deletion of a variable. For example, the autoregressive parameter rises from 0.34 to 0.49 in response to the deletion of the education variable and its spatial lag.

The optimal value of for the nearest neighbor estimator was 0.3352 which contrasts with the optimal value for the four neighbor estimator of 0.6150. The optimal value for was 0.341 ( was 0.733, 0.267), a small increase over the nearest neighbor result. Median errors and other statistics show the nearest neighbor estimator results continue to fall between the results for OLS and the four neighbor maximum likelihood estimator.

Undoubtedly, other ways of specifying the model and the form of spatial dependence would result in superior performance to the simple four neighbor model used as a standard of comparison. In addition, other data and models may yield different degrees of improvement. To partially address this and to address the important problem of anisotropy, we conduct a Monte Carlo study in the following section.

VII. Simulated Performance With Misspecified Dependence and Anisotropy

Real world data exhibit far more complicated spatial patterns than nearest neighbor spatial dependency. Actual spatial processes almost certainly involve multiple neighbors and may vary with direction or by zone (directional or zonal anisotropy). To better understand the roles of the various estimators, we devised a Monte Carlo experiment which employs random variables with multiple neighbor dependency and with a particularly simple form of zonal anisotropy.

We begin by examining the point distribution of the spatially correlated random variates. Points were uniformly distributed over the regions depicted in Diagram 10. We begin with the pair of boxes on the left separated by a heavy line. Each of these boxes has an aspect ratio (length/width) of 2. In contrast, the pair of boxes on the right have an aspect ratio of 10. As the aspect ratio becomes very large, the regions approach lines.

In the experiment, we will generate pairs of correlated random variables (with spatially random point distributions). One of each pair will be distributed in the top box and the other will be distributed in the bottom box. Each set of variates in the pair will be statistically independent of the other. Within each box the random variates are spatially correlated, but between boxes they are independent. Hence, this provides a simple way to simulate a form of zonal anisotropy.

The anisotropy in this case means that movements in one direction (along the longest dimension) can have different dependence than movements in the other direction (along the shortest dimension). One could imagine the heavy line represents a closed border, for example, and economic activity on each side is correlated, but that the economic activities of each side are independent of the other side due to lack of trade.

(Diagram 10)

This situation presents a way of testing the performance of the nearest neighbor estimator in the presence of a form of anisotropy. Since estimator relies upon the nearest neighbor, selection of the neighbors disregarding the closed border will sometimes result in the selection of an across-the-border neighbor whose value is statistically independent of the dependent variable. Obviously, this should degrade the performance of the nearest neighbor estimator. As the aspect ratio rises (holding n constant), the chance of the nearest neighbor appearing across the border rises. Hence, one would expect a larger aspect ratio would hurt the performance of the nearest neighbor estimator. By the same argument, other estimators using a weight matrix based upon isotropic dependencies should also suffer a loss in performance.

Given the point distributions of the random variates, we simulated a mixed-regressive spatially autoregressive process where Y is an n vector of the correlated dependent variable, is the autoregressive parameter, W is an n by n Delaunay spatial weight matrix formed from the spatially random point distribution, is an n element N(0,1) vector, , and is an n element vector of N(0,1) random errors.

(30)

For each of the six values of the autoregressive parameter (), we ran a pair of separate realizations of 5,000 observations each. We repeated this 100 times. Hence, we ran 100 pairs of realizations of the 5,000 element vector of correlated random variates for each of the six values of the autoregressive parameter.

To place the results in perspective, we compared the nearest neighbor estimates based upon the pooled observations to (1) estimates from using the correct model (Delaunay weights) for each one of the series comprising the pair, (2) estimates from using the Delaunay weights based upon the pooled observations, and (3) the OLS estimates based upon the pooled observations. Table 2 presents the results of the experiment. As an illustration of the accuracy of the simulation and estimation, the average estimated autoregressive parameter was on average within two places of the true autoregressive parameter.

The root mean-squared error of the nearest neighbor estimates relative to the root mean-squared error of the estimates based upon the true model appear in the column , the root mean-squared error of the maximum likelihood estimates using Delaunay weights for the pooled observations (ignoring the boundary) relative to the root mean-squared error of the estimates based upon the true model appear in the column , while the root mean-squared error of the OLS estimates relative to the estimates based upon the true model appear in the column .

The performance of the nearest neighbor (ML-1), the Delaunay model applied to the pooled observations, and OLS degrade as the degree of autocorrelation increases. However, the nearest neighbor performs relatively better than OLS for high levels of autocorrelation. For an aspect ratio of 2 and the closest neighbor estimator does 67% worse than the estimator based upon the true model but the OLS estimator performs 166% worse and the estimator using the Delaunay weights based upon the pooled observations performs only 2% worse.

The increase in the aspect ratio adversely affects the performance of all the estimators. However, OLS suffers relatively smaller losses compared to the nearest neighbor. For an aspect ratio of 200,000 and the closest neighbor estimator does 171% worse than the estimator based upon the true model, the OLS estimate performs 222% worse, and the estimator using the Delaunay weights based upon the pooled observations performs 126% worse. Interestingly, the anistropy in this case produced a greater deleterious effect upon the more complicated model (Delaunay on the pooled observations) than it did upon the simpler nearest neighbor and a fortiori OLS. The point geometry (e.g., aspect ratio) changes the nature of the spatial dependence and hence the performance of the estimators vary.

Note, the Delaunay W is symmetric and for this two-dimensional example will on average have six neighbors for each observation. Although the true process is symmetric, two-dimensional, and anisotropic, the nearest neighbor method outperforms OLS. For this experiment, an imperfect but limited model of spatial dependence can still outperform OLS, which assumes no dependence.

VIII. Computational Considerations

To find the nearest neighbors one usually employs a quadtree algorithm or a Delaunay triangularization algorithm. Optimal versions of these algorithms require O(nlog(n)) operations (Eppstein, Paterson, and Yao (1997)). The OLS computations used by the estimator require O(nk³) operations.

To examine this further, we generated some random points for the dependent variable, the independent variables (k=21), and their two-dimensional locational coordinates. Table 3 contains the timings. The timings coincide well with the predictions of O(nlog(n)) performance.

The nearest neighbor spatial autoregressions require almost no time for the sample sizes used in practice, require less than one minute for 100,000 observations, and under four minutes for 500,000 observations. We measured the times on a dual 500 megahertz Pentium III with one gigabyte of memory running Matlab 5.3 under Windows NT 4.0.

We used a Delaunay triangle based method for finding the nearest neighbor. The routines do not require the use of sparse matrices and hence could be implemented in other computer languages. The data and software is available at www.finance.lsu.edu/re.

IX. Conclusion

Even from a set of irregularly arranged points, nearest neighbor graphs display several important features which lead to a very simple eigenstructure of the associated adjacency matrix. In turn, this simple eigenstructure leads to an extremely simple closed-form for the log-determinant. Finally, this leads to a simple closed-form solution to the spatially autoregressive or mixed regressive spatially autoregressive maximum likelihood first-order conditions.

The results for nearest neighbor dependency have a number of applications. First, its mathematical simplicity could aid theoretical investigations of more complicated situations where space plays a role. The nearest neighbor model combines certain aspects of spatial statistical and time series dependence. Some observations have the simultaneous dependence (i.e., the 2-node-cycles) the hallmark of spatial statistics while other observations depend strictly upon the previous node in the graph. As demonstrated by the Monte Carlo experiment, the low cost of computation facilitates the execution of simulation experiments. Hence, one could investigate other estimators than the ones studied here. For example, one could easily compute the maximum likelihood SAR estimator. The asymmetry of the weights makes it difficult, however, to develop a traditional maximum likelihood CAR version of the nearest neighbor model (Besag (1974), Besag (1975)).

Spatial robustness is one area where the simplicity of the nearest neighbor estimator might play a role. Intuitively, the nearest neighbor model could suffer more severely from spatial outliers than more complicated models of spatial dependence. However, its simple structure might permit robustification without making the resulting estimator overly difficult to use. In fact, the nearest neighbor model might provide a means of detecting such spatial outliers.

Second, the easy computation of the estimates combined with the simplicity of the method may confer a pedagogical advantage to the nearest neighbor model of spatial dependence. Since calculating the log-determinant only requires counting the number of symmetric elements in the weight matrix, students could compute simple examples by hand.

Third, since more general spatial dependence models usually include the nearest neighbor, focusing on the nearest neighbor has a certain generality with respect to the different approaches. As the nearest neighbor estimator should underfit the overall spatial dependence, it provides a conservative estimate of the spatial structure and of the possible improvement one could obtain with spatial methods. Both the empirical example and the Monte Carlo study provided herein corroborate this assertion. Indeed, underfitting may be preferable to overfitting. See Griffith (1996, p. 80)) for a discussion of the whether to underfit or overfit. Intuitively, underfitting may serve as a robust strategy in the presence of uncertain anisotropy. The Monte Carlo experiment conducted in this paper tends to support this view. In the experiment, anisotropy actually degraded the performance of the more complicated estimator (using Delaunay weights based upon all the points ignoring the boundary) relative to the simpler nearest neighbor estimator.

Overfitting may also play a role. A liberal estimate of improvement obtained by overfitting a spatial model (e.g., placing many variants of spatially lagged dependent and independent variables on the right-hand-side of the equation and estimating by OLS) could approximately bound the possible performance gains obtainable through spatial estimation.

Fourth, the generality of nearest neighbor dependence and its ease of computation suggest the use of the nearest neighbor maximum likelihood estimate of spatial dependence as a possible diagnostic. Again, it should underestimate the total degree of spatial dependence and serve as a useful foil to upwardly biased EGLS estimates of spatial dependence. In addition, decomposable or local spatial autocorrelation statistics have become quite popular. One could perform a similar decomposition (by subgraph) on the nearest neighbor estimator.

Finally, the nearest neighbor maximum likelihood estimator could provide a more enlightened null hypothesis than spatial independence. In many cases, such as with real estate pricing models, assuming the null hypothesis of spatial error independence is a "straw man," since such data virtually always exhibit strong spatial dependence. Naïve models can serve as worthy adversaries, forcing more complicated models to justify their existence. For example, in the forecasting arena, many sophisticated models do not dominate their simpler counterparts in ex-sample prediction (Fildes and Makridakis (1995)). Possibly, the same phenomenon could occur with spatial models.