Applied Spatial Data Analysis

Chapter 3 - Spatial Lattice Data I

Author

Affiliation

Dr. Şebnem Er

Department of Statistical Sciences, University of Cape Town

Modelling Lattice Data

Exploratory spatial data analysis is often a preliminary step to more formal modelling approaches. In spatial lattice data, the goal is often to explain variation in a variable observed for each areal unit using other variables recorded for those same areal units.

In this lecture, the focus moves from exploratory spatial dependence to spatial regression models in a simple cross-sectional setting.

The key question is:

what happens to regression modelling when observations are spatially dependent?

Tip — click for more

Exploratory maps and Moran-type statistics tell us whether spatial structure appears to be present. Spatial regression asks how that structure should enter a formal model: through the dependent variable, through the error term, or both.

Exploratory spatial data analysis is often a preliminary step to more formal modelling approaches that seek to establish relationships between the observations of a variable and the observations of other variables, recorded for each areal unit. The focus in this chapter is on spatial regression models in a simple crosssectional setting

Spatial Regression Models

Our starting point is the ordinary linear regression model. For each lattice observation or areal unit (i=1,,N), let

[ y_i = x_i^+ _i. ]

Equivalently, in matrix form,

[ y = X+ , ]

where (y) is an ((N)) vector of observations, (X) is an ((NF)) matrix of explanatory variables, () is the vector of regression coefficients, and () is the error vector.

In the classical regression specification:

[ E(_i)=0, ]

[ (_i)=^2, ]

and, for (ij),

[ (_i,_j)=0. ]

Thus the errors are assumed to be independent and identically distributed, with constant variance and no correlation across observations.

The assumption of independent observations greatly simplifies the model. However, for areal data this is often unrealistic because neighbouring areas may influence one another or may share omitted spatially structured variables.

Original lecture visual for slide 3

Tip — click for more

In spatial regression, the problem is not simply that observations have coordinates. The problem is that the model errors, the dependent variable, or the explanatory variables may be systematically related across neighbouring areas. If that dependence is ignored, standard OLS interpretation and inference may fail.

In this model, y is an (N x 1) vector of observations on a dependent variable recorded at each of N locations, X is an (N x F) matrix of exogenous variables with f= 1 ; . . . ;F (including a constant, or one), b eta(f) the matching regression coefficient, and epsilon i the error term. If the regressors, residuals or the dependent variable are spatially dependent, the model suffers from a misspecification problem and the results of the model are biased or inconsistent.

Columbus dataset

We begin with the Columbus data and model crime as a function of income and housing value using OLS:

[ _i = _0+ _1_i+ _2_i+ _i. ]

A typical R workflow is:

library(spData)
library(spdep)
library(spatialreg)

data("columbus", package = "spData")

OLScolumbus <- lm(CRIME ~ INC + HOVAL, data = columbus)
summary(OLScolumbus)

Original lecture visual for slide 4

Tip — click for more

The OLS model is the baseline. We first estimate a non-spatial regression, then examine whether the residuals retain spatial patterning. If residuals are spatially autocorrelated, the model has not fully captured the spatial structure in the data.

All estimated coefficients are strongly significant. The fit is reasonable (0.552 in terms of rsquared ) We are going to obviously look at the residuals.

OLS – Residuals vs Fitted Plot & Q-Q Plot

The first diagnostic step is to inspect standard OLS residual plots.

The residuals versus fitted plot checks patterns in the residuals, non-linearity, and changing variance.
The Q-Q plot checks whether residuals are approximately normally distributed.

par(mfrow = c(1, 2))
plot(OLScolumbus, which = 1)
plot(OLScolumbus, which = 2)
par(mfrow = c(1, 1))

Original lecture visual for slide 5

Tip — click for more

These are ordinary regression diagnostics. They do not directly test spatial dependence, but they help identify model problems such as non-linearity, outliers, heteroskedasticity, or non-normal residuals that may complicate spatial modelling.

Predicted and Residual Maps

Predicted values and residuals can be mapped to see whether model performance varies spatially.

If areas with similar residuals cluster together, this suggests that the OLS model has left spatial structure unexplained.

A typical mapping workflow is:

columbus$ols_fitted <- fitted(OLScolumbus)
columbus$ols_resid  <- residuals(OLScolumbus)

plot(columbus["ols_fitted"], main = "OLS fitted values")
plot(columbus["ols_resid"],  main = "OLS residuals")

Original lecture visual for slide 6

Tip — click for more

A residual map is often more informative than a residual vector. Spatial clustering of positive residuals means the model is underpredicting in neighbouring places. Spatial clustering of negative residuals means the model is overpredicting in neighbouring places.

This map does suggest that similarly colored areas tend to be in similar locations, which could indicate positive spatial autocorrelation (a Moran’s I test for residual spatial autocorrelation is positive and highly significant). Also, it indicates a tendency to overpredict (negative residuals) in the outlying areas and a tendency to underpredict (positive residuals) in the core, suggesting the possible presence of spatial heterogeneity in the form of spatial regimes.

Moran’s I Test for Residual Spatial Autocorrelation

The next step is to test for residual spatial autocorrelation using Moran’s (I) applied to the OLS model residuals.

In R:

col.moran <- lm.morantest(
  OLScolumbus,
  dnbTresh1.listw
)

summary(col.moran)

Original lecture visual for slide 7

Tip — click for more

lm.morantest() is designed for regression residuals. It adjusts the reference distribution for the fact that residuals are estimated from a fitted model. This is why we do not simply treat OLS residuals as an ordinary observed variable without considering the fitted regression model.

The LM test for spatial autocorrelation in the error term has a slightly lower significance than that indicated by a traditional Moran test, although this is a general statement we can make. Why not simply apply Moran’s I test on the residuals? (Ref: anselinintrospatregres.pdf page 11) Also: Roger Bivand’s book page: 281: Since the Moran test is intended to detect spatial autocorrelation, we can try to fit a model taking this into account.We should not, however, forget that the misspecifications detected by Moran’s I can have a range of causes (see Sect. 9.4). It is also the case that if the fitted model exhibits multi-collinearity, the results of the test may be affected because of the numerical consequences of the model matrix not being of full rank for the expectation and variance of the statistic.

Spatial Models

Spatial regression models introduce spatial dependence explicitly into the regression framework.

The main specifications in this lecture are:

spatial lag / spatial autoregressive models;
spatial error models;
combined spatial lag and spatial error models;
diagnostic tests that help choose between them.

Original lecture visual for slide 8

Tip — click for more

The purpose of spatial regression is not just to obtain a significant spatial parameter. It is to represent the mechanism that generated the dependence: interaction in the dependent variable, omitted spatially structured errors, or both.

General Spatial Autoregressive (Lag) Models

A general spatial autoregressive model for cross-sectional lattice data can be written as a SARAR or spatial lag-error model:

[ y = W_1 y + X+ u, ]

[ u = W_2 u + , ]

with

[ N(0,^2 I). ]

Here:

(y) is the dependent variable;
(X) is the matrix of explanatory variables;
(W_1) is the spatial weights matrix for the lagged dependent variable;
(W_2) is the spatial weights matrix for the error process;
() is the spatial autoregressive lag coefficient;
() is the spatial autoregressive error coefficient;
(u) is the spatially structured disturbance;
() is the innovation error term.

The model can also be written as:

[ (I-W_1)y = X+ u, ]

and

[ (I-W_2)u = . ]

Special models are obtained by imposing restrictions such as (), (), or removing (X).

Original lecture visual for slide 9

Tip — click for more

The two spatial parameters play different roles. The parameter ( ho) says that neighbouring values of the dependent variable directly influence each other. The parameter () says that unexplained shocks or omitted variables are spatially correlated.

Where rho represents a regression parameter to be estimated and epsilon denotes the stochastic disturbance in the relationship. The parameter rho would reflect the spatial dependence inherent in our sample data, measuring the average influence of neighboring or contiguous observations on observations in the vector y . If we posit spatial dependence between the individual observations in the data sample y , some part of the total variation in y across the spatial sample would be explained by each observation’s dependence on its neighbors . The parameter would reflect this in the typical sense of regression. In addition, we could calculate the proportion of the total variation in y that is explained by spatial dependence.

Special Forms of Spatial Autoregressive (Lag) Models

From the general form, we obtain several special cases.

First-order spatial autoregressive model

If there are no explanatory variables and no spatial error process:

[ y = Wy + . ]

This model explains variation in (y) through neighbouring values of (y).

Spatial lag model / SAR model

If explanatory variables are included but the error term is not spatially autocorrelated:

[ y = Wy + X+ . ]

Spatial error model

If spatial dependence enters through the disturbances:

[ y = X+ u, ]

[ u = Wu + . ]

Spatial lag-error model / SARAR model

If both processes are present:

[ y = W_1y + X+ u, ]

[ u = W_2u + . ]

Moran’s (I) for the OLS residuals may indicate spatial dependence, but by itself it does not tell us whether the dependence should be modelled as a spatial lag, a spatial error process, or both.

Original lecture visual for slide 10

Tip — click for more

Spatial lag and spatial error models have different interpretations. A spatial lag model implies substantive spillovers or interaction in the dependent variable. A spatial error model often points to omitted spatially structured variables, measurement problems, or spatially correlated shocks.

From the general form we can derive special models by imposing restrictions. First order spatial autoregressive model attempts to explain variation in y as a linear combination of contiguous or neighboring units with no other explanatory variables S patial error dependence : spatial correlation in the error term Morans I for the OLS residuals indicate a spatial dependence but which way around is it? Is it a spatial lag, error or both? Hence, it has become convenient to distinguish between spatial lag and spatial error model specifications. These models require the application of special estimation procedures. Why not OLS?

Problems with OLS

OLS is problematic when the true model contains spatial dependence.

The consequences depend on where the spatial dependence appears:

if the dependent variable is spatially lagged, the spatial lag (Wy) is endogenous;
if the error term is spatially autocorrelated, the OLS coefficients may be inefficient and the usual standard errors are invalid;
if the spatial process is misspecified, coefficient estimates and inference can be biased or inconsistent.

This motivates special estimation procedures such as maximum likelihood and instrumental variables.

Original lecture visual for slide 11

Tip — click for more

The central issue is that spatial dependence violates the assumptions used to justify ordinary least squares. In a spatial lag model, the regressor (Wy) contains neighbouring values of (y), which are themselves related to the error process through the spatial system.

First Order Spatial Lag Model with OLS

Consider the simplest first-order spatial lag model without explanatory variables:

[ y = Wy + . ]

The SAR specification uses neighbouring values of (y) to account for spatial dependence.

The spatially lagged variable is

[ Wy, ]

which is usually interpreted as a weighted average of neighbouring values when (W) is row-standardised.

If (y) is mean standardised, the model focuses on deviations from the overall mean.

Original lecture visual for slide 12

Tip — click for more

This model is useful pedagogically because it isolates the spatial lag mechanism. However, in applied regression we usually include explanatory variables as well, leading to (y=Wy+X+).

Here show the OLS estimate calculation. Then the dependency. Let us take the simple form with no X variables. Y variable is mean standardized

First Order Spatial Lag Model with OLS

For the simple model

[ y=Wy+, ]

OLS treats (Wy) as an explanatory variable. The OLS estimator of () is

[ _{OLS} = ^{-1}(Wy)’y. ]

Equivalently,

[ _{OLS} = . ]

If (W) is symmetric, this is often written as

[ _{OLS} = . ]

The problem is that (Wy) is not exogenous. Because (y) itself is generated through a spatial feedback process, (Wy) is correlated with the error structure.

Original lecture visual for slide 13

Tip — click for more

This is the spatial analogue of an endogeneity problem. OLS assumes the explanatory variable is uncorrelated with the disturbance term. But in (y=Wy+), the variable (Wy) is built from (y), and (y) contains the disturbances propagated through the spatial system.

Let us take the simple form with no X variables. Y variable is mean standardized

First Order Spatial Autoregressive Model with OLS

For OLS to be consistent, the explanatory variable must be asymptotically uncorrelated with the error term. In the simple spatial lag model this requires a condition of the form

[ (Wy)’. ]

However, in the spatial case this condition generally fails. Since

[ y = (I-W)^{-1}, ]

we have

[ Wy = W(I-W)^{-1}. ]

Therefore,

[ (Wy)‘ = ’‘W’. ]

This is a quadratic form in the error terms, and it is not generally equal to zero.

Consequently, the OLS estimator is biased and inconsistent for the spatial lag parameter.

Original lecture visual for slide 14

Tip — click for more

This slide is the formal reason why OLS fails for a spatial lag model. The spatial multiplier ((I- ho W)^{-1}) propagates shocks through the network of neighbours, so the spatial lag regressor inherits dependence on the error term.

Refer to the hand written notes here.

Maximum Likelihood Estimation

Because OLS is inappropriate for models with spatial dependence, maximum likelihood estimation becomes a natural alternative.

Maximum likelihood estimates the model parameters by choosing the values that make the observed vector (y) most likely under the assumed spatial data-generating process.

The likelihood must account for two important transformations:

the transformation from (y) to the spatially filtered response ((I-W)y);
the transformation from the spatially structured disturbance (u) to the innovation ((I-W)u).

These transformations introduce Jacobian terms into the likelihood.

Original lecture visual for slide 15

Tip — click for more

The Jacobian term is the part that makes spatial ML different from ordinary regression ML. It adjusts the density when we transform the dependent variable through matrices such as (I- ho W).

Spatial Lag + Error Model

The spatial lag-error model can be written as

[ y = W_1^s y + Z+ u, ]

[ u = W_2^s u + . ]

Equivalently,

[ (I-W_1^s)y = Z+ u, ]

and

[ (I-W_2^s)u = . ]

Solving for the disturbance gives

[ u=(I-W_1^s)y-Z. ]

The innovation term is therefore

[ = (I-W_2^s) . ]

Original lecture visual for slide 16

Tip — click for more

The superscript (s) indicates row-standardised weights. Row-standardisation is common in applied work because the spatial lag becomes a neighbour average, but the interpretation of parameters always depends on the chosen weights matrix.

The Likelihood Function and the Jacobian for the General Model

For the general spatial lag-error model, the likelihood is built from the transformation between (y) and the normally distributed innovation term.

Define

[ B()=I-W_1^s, ]

and

[ A()=I-W_2^s. ]

Then

[ u=B()y-Z, ]

and

[ =A()u = A(){B()y-Z}. ]

The change of variables from (y) to () contributes the Jacobian terms

[ |B()| = |I-W_1^s| ]

and

[ |A()| = |I-W_2^s|. ]

Original lecture visual for slide 17

Tip — click for more

A determinant tells us how a matrix transformation changes volume. In likelihood terms, when (y) is transformed into spatially filtered errors, the determinant of the transformation matrix must be included.

The Likelihood Function and the Jacobian for the General Model

The log-likelihood for the general model contains the normal density component plus the spatial Jacobian terms.

Using the notation from the lecture, one form is

[ = -() - (||) + (|I-W_1^s|) + (|I-W_2^s|) - v’v. ]

Here (v) is the transformed, standardised residual vector implied by the spatial model.

The maximum likelihood estimator chooses the values of the parameters, including (), (), (), and variance parameters, that maximise this log-likelihood.

Ref: Anselin, p. 63; see also Anselin et al. (1996).

Original lecture visual for slide 18

Tip — click for more

The two log-determinant terms are the spatial Jacobians. They are computationally important because calculating (|I- ho W|) can become expensive for large spatial data sets.

We can either use, Maximum Likelihood estimation Instrumental Variables We will look at MLE. Y variable is mean standardized

Calculation of the Jacobian Term

There are several ways to calculate the Jacobian term, especially the log determinant

[ |I-W|. ]

Exact approaches

eigenvalues, Ord (1975), used in R with method = "eigen";
LU decomposition, method = "LU";
Cholesky decomposition, method = "Matrix";
sparse Cholesky decomposition, method = "spam".

Approximate approaches

Chebyshev approximation, Pace and LeSage (2004), method = "Chebyshev";
Monte Carlo approximation, Barry and Pace (1999), method = "MC".

Exact methods should agree up to numerical precision. Approximate methods should give similar values when the approximation is sufficiently accurate.

Original lecture visual for slide 19

Tip — click for more

Eigenvalue methods are convenient for small data sets, but they become infeasible when (N) is large. Sparse matrix methods are usually more scalable because spatial weights matrices are often mostly zeros.

The default method of method=“full” uses eigenvalues, and can thus also set the lower and upper bounds for the line search for λ accurately (as [1 / min i ( ζ i ) , 1 / max i ( ζ i )]), but is not feasible for large n . It should also be noted that although eigenvalues are computed for intrinsically asymmetric spatial weights matrices, their imaginary parts are discarded, so that even for method=“full” , the consequences of using such asymmetric weights matrices are unknown. Alternative approaches involve finding the log determinant of a Cholesky decomposition of the sparse matrix ( I − λW ) directly. Here it is not possible to pre-compute eigenvalues, so one log determinant is computed for each value of λ used, but the number needed is in general not excessive, and much larger n become feasible on ordinary computers. A number of different sparse matrix approaches have been tried, with the use of Matrix and method=“Matrix” , the one suggested currently. All of the sparse matrix approaches to computing the Jacobian require that matrix W be symmetric or at least similar to symmetric, thus providing for weights with “W” and “S” styles based on symmetric neighbour lists and symmetric general spatial weights, such as inverse distance. Matrices that are similar to symmetric have the same eigenvalues, so that the eigenvalues of symmetric W ∗ = D 1 / 2 WD 1 / 2 and row-standardised W = DB are the same, for symmetric binary or general weights matrix B , and D a diagonal matrix of inverse row sums of B , d ii = 1 / n j =1 b ij (Ord, 1975, p. 125).

Hypothesis Tests Based on the ML Principles

Maximum likelihood estimation naturally leads to several hypothesis-testing frameworks.

Wald test

Tests restrictions using the estimated parameter and its estimated variance.

Likelihood ratio test

Compares the maximised likelihoods of a restricted and an unrestricted model.

Lagrange multiplier test

Evaluates whether adding a spatial parameter would improve the model, using estimates from the restricted model, usually OLS.

Tip — click for more

The LM test is especially useful as a diagnostic because it can be computed after fitting the simpler OLS model. It asks whether the omitted spatial lag or spatial error component appears to be needed.

Example: Columbus Dataset

The Columbus data frame has 49 rows and 22 columns.

The main variables used here are:

CRIME: residential burglaries and vehicle thefts per thousand households in the neighbourhood;
INC: household income;
HOVAL: housing value.

library(spData)
data("columbus", package = "spData")

nrow(columbus)
ncol(columbus)
names(columbus)
summary(columbus[, c("CRIME", "INC", "HOVAL")])

Tip — click for more

The Columbus example is small, which makes it useful for teaching. We can estimate OLS, spatial lag, spatial error, and combined models without computational difficulty.

Model 0: OLS with INC and HOVAL

Model 0 is the OLS model with income and housing value:

[ _i = _0+ _1_i+ _2_i+ _i. ]

In R:

OLScolumbus <- lm(
  CRIME ~ INC + HOVAL,
  data = columbus
)

summary(OLScolumbus)

Original lecture visual for slide 22

Tip — click for more

This model treats the 49 neighbourhoods as independent observations after controlling for income and housing value. The rest of the lecture examines whether that assumption is adequate.

All estimated coefficients are strongly significant. The fit is reasonable (0.552 in terms of rsquared ) We are going to obviously look at the residuals.

Moran’s I Test for Residual Spatial Autocorrelation

The OLS residuals are tested for spatial autocorrelation.

col.moran <- lm.morantest(
  OLScolumbus,
  dnbTresh1.listw
)

summary(col.moran)

Original lecture visual for slide 23

Tip — click for more

A significant residual Moran test tells us that the OLS residuals still contain spatial structure. It is a diagnostic signal, not yet a complete model selection rule.

Spatial Lag Model (SLag) estimation in R with MLE

The spatial lag model introduces spatial correlation in the dependent variable:

[ y = Wy + X+ . ]

For the Columbus example, two models are considered:

Model 1: no independent variables;
Model 2: INC and HOVAL included as independent variables.

In R, maximum likelihood estimation can be performed with:

lagsarlm(
  formula,
  data = columbus,
  listw = dnbTresh1.listw,
  method = "eigen"
)

Original lecture visual for slide 24

Tip — click for more

The spatial lag model is appropriate when the outcome in one area is directly related to outcomes in neighbouring areas. For crime, this may represent diffusion, displacement, or interaction across neighbourhood boundaries.

Model 1: First order SLag Model

Model 1 is the first-order spatial lag model without independent variables:

[ y = Wy + . ]

In R:

SLag_model1 <- lagsarlm(
  CRIME ~ 1,
  data = columbus,
  listw = dnbTresh1.listw,
  method = "eigen"
)

summary(SLag_model1)

Original lecture visual for slide 25

Tip — click for more

This model is intentionally simple. It asks how much of the spatial pattern in crime can be represented solely through the spatial lag of crime, without additional explanatory variables.

Model 2: SLag Model with INC and HOVAL

Model 2 is the spatial lag model with INC and HOVAL:

[ = W + _0 + _1 + _2 + . ]

In R:

SLag_model2 <- lagsarlm(
  CRIME ~ INC + HOVAL,
  data = columbus,
  listw = dnbTresh1.listw,
  method = "eigen"
)

summary(SLag_model2)

Original lecture visual for slide 26

Tip — click for more

After adding (W y), coefficients of INC and HOVAL no longer have exactly the same interpretation as in OLS. In spatial lag models, impacts may include direct and indirect/spillover effects.

We can either use, Maximum Likelihood estimation Instrumental Variables We will look at MLE. Y variable is mean standardized

Spatial Error Model (SError) Estimation in R with MLE

A spatial error model is a regression with a non-spherical error term. Spatial dependence appears in the disturbances:

[ y = X+ u, ]

[ u = Wu + . ]

Equivalently,

[ u = (I-W)^{-1}. ]

The off-diagonal elements of the implied covariance matrix express spatial dependence between errors.

In R:

errorsarlm(
  formula,
  data = columbus,
  listw = dnbTresh1.listw,
  method = "eigen"
)

Original lecture visual for slide 27

Tip — click for more

A spatial error model is often interpreted as evidence of omitted spatially structured variables or spatially correlated shocks. The spatial parameter () does not represent direct interaction in the dependent variable in the same way as ( ho) does.

Model 3: SError Model with INC and HOVAL

Model 3 is the spatial error model with income and housing value:

[ = _0+ _1+_2+u, ]

[ u=Wu+. ]

In R:

SError_model3 <- errorsarlm(
  CRIME ~ INC + HOVAL,
  data = columbus,
  listw = dnbTresh1.listw,
  method = "eigen"
)

summary(SError_model3)

Original lecture visual for slide 28

Tip — click for more

Compare the spatial error model with the OLS model. If the main change is in standard errors and residual spatial autocorrelation, the error specification may be capturing spatially structured omitted information.

Model 4: SLag and SError Model

Model 4 includes both spatial lag and spatial error effects:

[ y = Wy + X+ u, ]

[ u = Wu + . ]

This is a spatial lag-error model, also called a SARAR-type specification.

In R, depending on package version and function availability, this class of model is commonly estimated with a SAC/SARAR function such as:

SLag_SError_model4 <- sacsarlm(
  CRIME ~ INC + HOVAL,
  data = columbus,
  listw = dnbTresh1.listw,
  method = "eigen"
)

summary(SLag_SError_model4)

Original lecture visual for slide 29

Tip — click for more

The combined model is more flexible, but also more demanding. It should not be chosen just because it is larger; the spatial lag and spatial error parameters should make substantive and statistical sense.

Lagrange Multiplier Test Statistics for Spatial Autocorrelation

Burridge (1980) and Anselin (1988) developed Lagrange Multiplier tests for spatial dependence.

The null hypothesis for both basic tests is the OLS model:

[ H_0:,. ]

The alternatives are:

spatial lag dependence;
spatial error dependence.

The LM test statistics are asymptotically distributed as chi-square random variables under the null.

If both the basic LM lag and LM error tests are significant, this does not by itself identify which model should be used. The two tests are not isolated from one another. Anselin et al. (1996) therefore proposed robust versions of these tests.

The robust test that remains significant indicates which spatial alternative is more strongly supported after accounting for the other possibility.

Original lecture visual for slide 30

Tip — click for more

The practical rule often taught is: if only one simple LM test is significant, consider that model; if both are significant, compare the robust LM tests. If robust LM lag is significant but robust LM error is not, the spatial lag model is favoured, and vice versa.

The robust LM tests take into account the alternative possibility, that is the LMerr test will respond to both an omitted spatially lagged dependent variable and spatially autocorrelated residuals, while the robust RLMerr is designed to test for spatially autocorrelated residuals in the possible presence of an omitted spatially lagged dependent variable. The lm.LMtests function here returns a list of five LM tests, which seem to point to a spatial lag specification. Again, it is the case that if the fitted model exhibits multicollinearity, the results of the tests will be affected. # LMerror , LMlag are both significant here then check the robust values and choose the higher robust value. Here SARMA is significant as well but forget about it. robust version means H0: no spatial autocorr H1: spatial error but taking into account that it could be spatial lag.

Lagrange Multiplier Test Statistics for Spatial Autocorrelation

There is also a joint test for spatial dependence.

The SARMA statistic combines evidence from lag and error alternatives. It can be constructed as a sum of LM components and is compared to a chi-square distribution with two degrees of freedom:

[ ^2_2. ]

The joint test asks whether some spatial dependence is present, but it does not by itself tell us whether the preferred specification is spatial lag, spatial error, or both.

Original lecture visual for slide 31

Tip — click for more

A significant SARMA test is useful evidence that the OLS model is inadequate, but it is less specific than the robust LM tests for choosing between lag and error alternatives.

Lagrange Multiplier Test Statistics for Spatial Autocorrelation

The lecture uses the lm.LMtests() function to perform the LM diagnostics.

The reported output is:

         Statistics df p-value
LMerr     6.3670738  1 0.0116257146
LMlag    13.6904198  1 0.0002155513
RLMerr    0.1016064  1 0.7499103168
RLMlag    7.4249524  1 0.0064325528
SARMA    13.7920262  2 0.0010118114

The simple LM error and simple LM lag tests are both significant. However, among the robust tests, the robust lag test is significant while the robust error test is not. This points to a spatial lag specification.

In R:

lm.LMtests(
  OLScolumbus,
  dnbTresh1.listw,
  test = c("LMerr", "LMlag", "RLMerr", "RLMlag", "SARMA")
)

Original lecture visual for slide 32

Tip — click for more

The robust LM lag result is the key line here. RLMlag remains significant after accounting for possible spatial error dependence, whereas RLMerr does not remain significant after accounting for a possible spatial lag.

GeoDa Workbook Reference

The GeoDa workbook provides a worked treatment of these spatial regression diagnostics and model-selection ideas.

Reference:

Anselin (2005), Exploring Spatial Data with GeoDa: A Workbook, p. 217.

https://geodacenter.github.io/docs/geodaworkbook.pdf

Original lecture visual for slide 33

Tip — click for more

This slide is a useful bridge between R and GeoDa. The same conceptual diagnostics—OLS residual spatial autocorrelation, LM lag, LM error, robust LM tests—are central in both workflows.

In Python

Python spatial analysis resources are available through PySAL:

https://pysal.org/

PySAL provides tools for spatial weights, exploratory spatial data analysis, spatial econometrics, and regional/spatial modelling.

Tip — click for more

PySAL is the major Python ecosystem for spatial analysis. In R, analogous workflows often use packages such as sf, spdep, spatialreg, tmap, and spData.

References

Selected references from the lecture:

http://www.eurojournals.com/ejss_18_1_02.pdf
http://www.ecomod.org/files/papers/486.pdf
http://www.people.fas.harvard.edu/~zhukov/Spatial6.pdf
http://www.drs.wisc.edu/documents/articles/curtis/cesoc977-12/W3_W6_W9_GeodaWorkbook.pdf
http://dae.unizar.es/docencia/regional/spacestat%20Tutorial.pdf
http://www.people.fas.harvard.edu/~zhukov/spatial.html
https://geodacenter.asu.edu/system/files/rex1.pdf
http://scc.stat.ucla.edu/page_attachments/0000/0094/spatial_R_1_09S.pdf
Giuseppe Arbia, Spatial Econometrics, Springer, 2005.
Roger Bivand et al., Applied Spatial Data Analysis with R, Springer, 2008.
http://www.jstatsoft.org/v47/i01/paper
http://cran.r-project.org/web/packages/splm/splm.pdf
http://www.r-project.org/conferences/useR-2009/slides/Millo+Piras.pdf
http://cran.r-project.org/web/packages/plm/plm.pdf
http://www.jstatsoft.org/v27/i02/paper
http://www.spatial-econometrics.com/html/wbook.pdf
http://geodacenter.asu.edu/

Other Resources

Additional resources:

SEAI, May 2012 Rome Lectures: G. Arbia, L. Anselin, I. Prucha, J. H. P. Paelinck, D. Arribas-Bel, G. Piras.
Giuseppe Arbia, Spatial Econometrics, Stellenbosch University, Fall 2017.

Modelling Lattice Data

Spatial Regression Models

Columbus dataset

OLS – Residuals vs Fitted Plot & Q-Q Plot

Predicted and Residual Maps

Moran’s I Test for Residual Spatial Autocorrelation

Spatial Models

General Spatial Autoregressive (Lag) Models

Special Forms of Spatial Autoregressive (Lag) Models

First-order spatial autoregressive model

Spatial lag model / SAR model

Spatial error model

Spatial lag-error model / SARAR model

Problems with OLS

First Order Spatial Lag Model with OLS

First Order Spatial Lag Model with OLS

First Order Spatial Autoregressive Model with OLS

Maximum Likelihood Estimation

Spatial Lag + Error Model

The Likelihood Function and the Jacobian for the General Model

The Likelihood Function and the Jacobian for the General Model

Calculation of the Jacobian Term

Exact approaches

Approximate approaches

Hypothesis Tests Based on the ML Principles

Wald test

Likelihood ratio test

Lagrange multiplier test

Example: Columbus Dataset

Model 0: OLS with INC and HOVAL

Moran’s I Test for Residual Spatial Autocorrelation

Spatial Lag Model (SLag) estimation in R with MLE

Model 1: First order SLag Model

Model 2: SLag Model with INC and HOVAL

Spatial Error Model (SError) Estimation in R with MLE

Model 3: SError Model with INC and HOVAL

Model 4: SLag and SError Model

Lagrange Multiplier Test Statistics for Spatial Autocorrelation

Lagrange Multiplier Test Statistics for Spatial Autocorrelation

Lagrange Multiplier Test Statistics for Spatial Autocorrelation

GeoDa Workbook Reference

Other Topics

In Python

References

Other Resources