# Chapter 17 Spatial econometrics models

Spatial autoregression models using spatial weights matrices were described in some detail using maximum likelihood estimation some time ago (Cliff and Ord 1973, 1981). A family of models was elaborated in spatial econometric terms extending earlier work, and in many cases using the simultaneous autoregressive framework and row standardization of spatial weights (Anselin 1988). The simultaneous and conditional autoregressive frameworks can be compared, and both can be supplemented using case weights to reflect the relative importance of different observations (Waller and Gotway 2004).

Before moving to presentations of issues raised in fitting spatial regression models, it is worth making a few further points. A recent review of spatial regression in a spatial econometrics setting is given by Kelejian and Piras (2017); note that their usage is to call the spatial coefficient of the lagged response \(\lambda\) and that of the lagged residuals \(\rho\), the reverse of other usage (Anselin 1988; LeSage and Pace 2009); here we use \(\rho_{\mathrm{Lag}}\) for the spatial coefficient in the spatial lag model, and \(\rho_{\mathrm{Err}}\) for the spatial error model. One interesting finding is that relatively dense spatial weights matrices may downweight model estimates, suggesting that sparser weights are preferable (Smith 2009). Another useful finding is that the presence of residual spatial autocorrelation need not bias the estimates of variance of regression coefficients, provided that the covariates themselves do not exhibit spatial autocorrelation (Smith and Lee 2012). In general, however, the footprints of the spatial processes of the response and covariates may not be aligned, and if covariates and the residual are autocorrelated, it is likely that the estimates of variance of regression coefficients will be biassed downwards if attempts are not made to model the spatial processes.

## 17.1 Spatial econometric models: definitions

In trying to model spatial processes, one of the earliiest spatial econometric representations is to model the spatial autocorrelation in the residual (spatial error model, SEM):

\[ {\mathbf y} = {\mathbf X}{\mathbf \beta} + {\mathbf u}, \qquad {\mathbf u} = \rho_{\mathrm{Err}} {\mathbf W} {\mathbf u} + {\mathbf \varepsilon}, \] where \({\mathbf y}\) is an \((N \times 1)\) vector of observations on a response variable taken at each of \(N\) locations, \({\mathbf X}\) is an \((N \times k)\) matrix of covariates, \({\mathbf \beta}\) is a \((k \times 1)\) vector of parameters, \({\mathbf u}\) is an \((N \times 1)\) spatially autocorrelated disturbance vector, \({\mathbf \varepsilon}\) is an \((N \times 1)\) vector of independent and identically distributed disturbances and \(\rho_{\mathrm{Err}}\) is a scalar spatial parameter.

This model, and other spatial econometric models, do not fit into the mixed models framework. Here the modelled spatial process interacts directly with the response, covariates, and their coefficients. This modelling framework appears to draw on an older tradition extending time series to two dimensions:

\[ {\mathbf u} = ({\mathbf I} - \rho_{\mathrm{Err}} {\mathbf W})^{-1} {\mathbf \varepsilon}, \qquad {\mathbf y} = {\mathbf X}{\mathbf \beta} + ({\mathbf I} - \rho_{\mathrm{Err}} {\mathbf W})^{-1} {\mathbf \varepsilon}, \qquad ({\mathbf I} - \rho_{\mathrm{Err}} {\mathbf W}) {\mathbf y} = ({\mathbf I} - \rho_{\mathrm{Err}} {\mathbf W}) {\mathbf X}{\mathbf \beta} + {\mathbf \varepsilon}. \]

If the processes in the covariates and the response match, we should find little difference between the coefficients of a least squares and a SEM, but very often they diverge, suggesting that a Hausman test for this condition should be employed (Pace and LeSage 2008). This may be related to earlier discussions of a spatial equivalent to the unit root and cointegration where spatial processes match (Fingleton 1999).

A model with a spatial process in the response only is termed a spatial lag model (SLM, often SAR - spatial autoregressive) (LeSage and Pace 2009). Durbin models add the spatially lagged covariates to the covariates included in the spatial model; spatial Durbin models are reviewed by Mur and Angulo (2006). If it is chosen to admit a spatial process in the residuals in addition to a spatial process in the response, again two models are formed, a general nested model (GNM) nesting all the others, and a model without spatially lagged covariates (SAC, also known as SARAR - Spatial AutoRegressive-AutoRegressive model). If neither the residuals nor the response are modelled with spatial processes, spatially lagged covariates may be added to a linear model, as a spatially lagged X model (SLX) (Elhorst 2010; Bivand 2012; LeSage 2014; Halleck Vega and Elhorst 2015). We can write the general nested model (GNM) as:

\[ {\mathbf y} = \rho_{\mathrm{Lag}} {\mathbf W}{\mathbf y} + {\mathbf X}{\mathbf \beta} + {\mathbf W}{\mathbf X}{\mathbf \gamma} + {\mathbf u}, \qquad {\mathbf u} = \rho_{\mathrm{Err}} {\mathbf W} {\mathbf u} + {\mathbf \varepsilon}, \] where \({\mathbf \gamma}\) is a \((k' \times 1)\) vector of parameters. \(k'\) defines the subset of the intercept and covariates, often \(k' = k-1\) when using row standardised spatial weights and omitting the spatially lagged intercept.

This may be constrained to the double spatial coefficient model SAC/SARAR by setting \({\mathbf \gamma} = 0\), to the spatial Durbin (SDM) by setting \(\rho_{\mathrm{Err}} = 0\), and to the error Durbin model (SDEM) by setting \(\rho_{\mathrm{Lag}} = 0\). Imposing more conditions gives the spatial lag model (SLM) with \({\mathbf \gamma} = 0\) and \(\rho_{\mathrm{Err}} = 0\), the spatial error model (SEM) with \({\mathbf \gamma} = 0\) and \(\rho_{\mathrm{Lag}} = 0\), and the spatially lagged X model (SLX) with \(\rho_{\mathrm{Lag}} = 0\) and \(\rho_{\mathrm{Err}} = 0\).

Although making predictions for new locations for which covariates are observed was raised as an issue some time ago, it has many years to make progress in reviewing the possibilities (Bivand 2002; Goulard, Laurent, and Thomas-Agnan 2017; Laurent and Margaretic 2021). The prediction methods for SLM, SDM, SEM, SDEM, SAC and GNM models fitted with maximum likelihood were contributed as a Google Summer of Coding project by Martin Gubri. This work, and work on similar models with missing data (Suesse 2018) is also relevant for exploring censored median house values in the Boston data set. Work on prediction also exposed the importance of the reduced form of these models, in which the spatial process in the response interacts with the regression coefficients in the SLM, SDM, SAC and GNM models.

The consequence of these interactions is that a unit change in a covariate will only impact the response as the value of the regression coefficient if the spatial coefficient of the lagged response is zero. Where it is non-zero, global spillovers, impacts, come into play, and these impacts should be reported rather than the regression coefficients (LeSage and Pace 2009; Elhorst 2010; Bivand 2012; LeSage 2014; Halleck Vega and Elhorst 2015). Local impacts may be reported for SDEM and SLX models, using linear combination to calculate standard errors for the total impacts of each covariate (sums of coefficients on the covariates and their spatial lags).

This can be seen from the GNM data generation process:

\[ ({\mathbf I} - \rho_{\mathrm{Err}} {\mathbf W})({\mathbf I} - \rho_{\mathrm{Lag}} {\mathbf W}){\mathbf y} = ({\mathbf I} - \rho_{\mathrm{Err}} {\mathbf W})({\mathbf X}{\mathbf \beta} + {\mathbf W}{\mathbf X}{\mathbf \gamma}) + {\mathbf \varepsilon}, \] re-writing:

\[ {\mathbf y} = ({\mathbf I} - \rho_{\mathrm{Lag}} {\mathbf W})^{-1}({\mathbf X}{\mathbf \beta} + {\mathbf W}{\mathbf X}{\mathbf \gamma}) + ({\mathbf I} - \rho_{\mathrm{Lag}} {\mathbf W})^{-1}({\mathbf I} - \rho_{\mathrm{Err}} {\mathbf W})^{-1}{\mathbf \varepsilon}. \] There is interaction between the \(\rho_{\mathrm{Lag}}\) and \({\mathbf \beta}\) (and \({\mathbf \gamma}\) if present) coefficients. This can be seen from the partial derivatives: \(\partial y_i / \partial x_{jr} = (({\mathbf I} - \rho_{\mathrm{Lag}} {\mathbf W})^{-1} ({\mathbf I} \beta_r + {\mathbf W} \gamma_r))_{ij}\). This dense matrix \(S_r({\mathbf W}) = (({\mathbf I} - \rho_{\mathrm{Lag}} {\mathbf W})^{-1} ({\mathbf I} \beta_r + {\mathbf W} \gamma_r))\) expresses the direct impacts (effects) on its principal diagonal, and indirect impacts in off-diagonal elements.

Piras and Prucha (2014) revisit and correct Florax, Folmer, and Rey (2003) (see also comments by Hendry (2006) and Florax, Folmer, and Rey (2006)), finding that the common use of pre-test strategies for model selection probably ought to be replaced by the estimation of the most general model appropriate for the relationships being modelled. In the light of this finding, pre-test model selection will not be used here.

Current work in the **spatialreg** package is focused on refining the handling of spatially lagged covariates using a consistent `Durbin=`

argument taking either a logical value or a formula giving the subset of covariates to add in spatially lagged form. There is a speculation that some covariates, for example some dummy variables, should not be added in spatially lagged form. This then extends to handling these included spatially lagged covariates appropriately in calculating impacts. This work applies to cross-sectional models fitted using MCMC or maximum likelihood, and will offer facilities to spatial panel models.

It is worth mentioning the almost unexplored issues of functional form assumptions, for which flexible structures are useful, including spatial quantile regression presented in the **McSpatial** package (McMillen 2013). There are further issues with discrete response variables, covered by some functions in **McSpatial**, and in the **spatialprobit** and **ProbitSpatial** packages (Wilhelm and Matos 2013; Martinetti and Geniaux 2017); the MCMC implementations of the former are based on LeSage and Pace (2009). Finally, Wagner and Zeileis (2019) show how an SLM model may be used in the setting of recursive partitioning, with an implementation using `spatialreg::lagsarlm()`

in the **lagsarlmtree** package.

The review of cross-sectional maximum likelihood and generalized method of moments (GMM) estimators in **spatialreg** (R. Bivand and Piras 2021) and **sphet** for spatial econometrics style spatial regression models by Bivand and Piras (2015) is still largely valid. In the review, estimators in these R packages were compared with alternative implementations available in other programming languages elsewhere. The review did not cover Bayesian spatial econometrics style spatial regression. More has changed with respect to spatial panel estimators described in Millo and Piras (2012), but will not be covered here.

Because R. Bivand, Millo, and Piras (2021) covers many of the features of R packages for spatial econometrics, updating Roger S. Bivand and Piras (2015), and including recent advances in General Method of Moments and spatial panel modelling, this chapter will be restricted to a small number of examples drawing on R. Bivand (2017) using the Boston house value data set.

## 17.2 Maximum likelihood estimation in **spatialreg**

For models with single spatial coefficients (SEM and SDEM using `errorsarlm()`

, SLM and SDM using `lagsarlm()`

), the methods initially described by Ord (1975) are used. The following table shows the functions that can be used to estimate the models described above using maximum likelihood.

model | model name | maximum likelihood estimation function |
---|---|---|

SEM | spatial error | `errorsarlm(..., Durbin=FALSE, ...)` |

SEM | spatial error | `spautolm(..., family="SAR", ...)` |

SDEM | spatial Durbin error | `errorsarlm(..., Durbin=TRUE, ...)` |

SLM | spatial lag | `lagsarlm(..., Durbin=FALSE, ...)` |

SDM | spatial Durbin | `lagsarlm(..., Durbin=TRUE, ...)` |

SAC | spatial autoregressive combined | `sacsarlm(..., Durbin=FALSE, ...)` |

GNM | general nested | `sacsarlm(..., Durbin=TRUE, ...)` |

The estimating functions `errorsarlm()`

and `lagsarlm()`

take similar arguments, where the first two, `formula=`

and `data=`

are shared by most model estimating functions. The third argument is a `listw`

spatial weights object, while `na.action=`

behaves as in other model estimating functions if the spatial weights can reasonably be subsetted to avoid observations with missing values. The `weights=`

argument may be used to provide weights indicating the known degree of per-observation variability in the variance term - this is not available for `lagsarlm()`

.

The `Durbin=`

argument replaces the earlier `type=`

and `etype=`

arguments, and if not given is taken as `FALSE`

. If given, it may be `FALSE`

, `TRUE`

in which case all spatially lagged covariates are included, or a one-sided formula specifying which spatially lagged covariates should be included. The `method=`

argument gives the method for calculating the log determinant term in the log likelihood function, and defaults to `"eigen"`

, suitable for moderately sized data sets. The `interval=`

argument gives the bounds of the domain for the line search using `stats::optimize()`

used for finding the spatial coefficient. The `tol.solve()`

argument, passed through to `base::solve()`

, was needed to handle data sets with differing numerical scales among the coefficients which hindered inversion of the variance-covariance matrix; the default value in `base::solve()`

used to be much larger. The `control=`

argument takes a list of control values to permit more careful adjustment of the running of the estimation function.

The `sacsarlm()`

function may take second spatial weights and interval arguments if the spatial weights used to model the two spatial processes in the SAC and GNM specifications differ. By default, the same spatial weights are used. By default, `stats::nlminb()`

is used for numerical optimization, using a heuristic to choose starting values. Like `lagsarlm()`

, this function does not take a `weights=`

argument.

Where larger data sets are used, a numerical Hessian approach is used to calculate the variance-covariance matrix of coefficients, rather than an analytical asymptotic approach.

### 17.2.1 Boston house value data set examples

The examples use the objects read and created in chapter 16, based on R. Bivand (2017).

```
library(spatialreg)
eigs_489 <- eigenw(lw_q_489)
SDEM_489 <- errorsarlm(form, data = boston_489, listw = lw_q_489, Durbin = TRUE,
zero.policy = TRUE, control = list(pre_eig = eigs_489))
SEM_489 <- errorsarlm(form, data = boston_489, listw = lw_q_489,
zero.policy = TRUE, control = list(pre_eig = eigs_489))
```

Here we are using the `control=`

list argument to pass through pre-computed eigenvalues for the default `"eigen"`

method.

```
cbind(data.frame(model=c("SEM", "SDEM")),
rbind(broom::tidy(Hausman.test(SEM_489)),
broom::tidy(Hausman.test(SDEM_489))))[,1:4]
```

```
# model statistic p.value parameter
# 1 SEM 52.0 2.83e-06 14
# 2 SDEM 48.7 6.48e-03 27
```

Both Hausman test results for the 489 tract data set suggest that the regression coefficients do differ from their non-spatial counterparts, perhaps indicating that the footprints of the spatial processes do not match.

```
eigs_94 <- eigenw(lw_q_94)
SDEM_94 <- errorsarlm(form, data=boston_94, listw=lw_q_94, Durbin=TRUE,
control=list(pre_eig=eigs_94))
SEM_94 <- errorsarlm(form, data=boston_94, listw=lw_q_94, control=list(pre_eig=eigs_94))
```

For the 94 air pollution model output zones, the Hausman tests find little difference between coefficients:

```
cbind(data.frame(model=c("SEM", "SDEM")),
rbind(broom::tidy(Hausman.test(SEM_94)),
broom::tidy(Hausman.test(SDEM_94))))[, 1:4]
```

```
# model statistic p.value parameter
# 1 SEM 15.66 0.335 14
# 2 SDEM 9.21 0.999 27
```

This is related to the fact that the SEM and SDEM models add little to least squares or SLX at the air pollution model output zone level, using likelihood ratio tests:

```
cbind(data.frame(model=c("SEM", "SDEM")), rbind(broom::tidy(LR1.Sarlm(SEM_94)),
broom::tidy(LR1.Sarlm(SDEM_94))))[,c(1, 4:6)]
```

```
# model statistic p.value parameter
# 1 SEM 2.593 0.107 1
# 2 SDEM 0.216 0.642 1
```

We can use `spatialreg::LR.Sarlm()`

to apply a likelihood ratio test between nested models, but here choose `lmtest::lrtest()`

, which gives the same results, preferring models including spatially lagged covariates both for tracts and model output zones:

```
# # A tibble: 2 × 5
# X.Df LogLik df statistic p.value
# <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 16 273. NA NA NA
# 2 29 311. 13 74.4 1.23e-10
```

```
# # A tibble: 2 × 5
# X.Df LogLik df statistic p.value
# <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 16 59.7 NA NA NA
# 2 29 81.3 13 43.2 0.0000421
```

The SLX model is fitted using least squares, and also returns a log likelihood value, letting us test whether we need a spatial process in the residuals. In the tract data set we obviously do:

```
SLX_489 <- lmSLX(form, data = boston_489, listw = lw_q_489, zero.policy = TRUE)
broom::tidy(lmtest::lrtest(SLX_489, SDEM_489))
```

```
# # A tibble: 2 × 5
# X.Df LogLik df statistic p.value
# <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 28 231. NA NA NA
# 2 29 311. 1 159. 1.55e-36
```

but in the output zone case we do not.

```
SLX_94 <- lmSLX(form, data = boston_94, listw = lw_q_94)
broom::tidy(lmtest::lrtest(SLX_94, SDEM_94))
```

```
# # A tibble: 2 × 5
# X.Df LogLik df statistic p.value
# <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 28 81.2 NA NA NA
# 2 29 81.3 1 0.216 0.642
```

These outcomes are sustained also when we use the counts of house units by tract and output zones as case weights:

```
SLX_489w <- lmSLX(form, data = boston_489, listw = lw_q_489, weights = units, zero.policy = TRUE)
SDEM_489w <- errorsarlm(form, data = boston_489, listw = lw_q_489, Durbin =TRUE,
weights = units, zero.policy = TRUE, control = list(pre_eig = eigs_489))
broom::tidy(lmtest::lrtest(SLX_489w, SDEM_489w))
```

```
# # A tibble: 2 × 5
# X.Df LogLik df statistic p.value
# <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 28 311. NA NA NA
# 2 29 379. 1 136. 1.70e-31
```

```
SLX_94w <- lmSLX(form, data = boston_94, listw = lw_q_94, weights = units)
SDEM_94w <- errorsarlm(form, data = boston_94, listw = lw_q_94, Durbin = TRUE, weights = units,
control = list(pre_eig = eigs_94))
broom::tidy(lmtest::lrtest(SLX_94w, SDEM_94w))
```

```
# # A tibble: 2 × 5
# X.Df LogLik df statistic p.value
# <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 28 97.5 NA NA NA
# 2 29 98.0 1 0.917 0.338
```

In this case and based on arguments advanced in R. Bivand (2017), the use of weights is justified because tract counts of reported housing units underlying the weighted median values vary from 5 to 3,031, and air pollution model output zone counts vary from 25 to 12,411. Because of this, and because a weighted general nested model has not been developed, we cannot take the GNM as the starting point for general-to-simpler testing, but start rather from the SDEM model, and use the Hausman test to guide the choice of units of observation.

## 17.3 Impacts

Global impacts have been seen as crucial for reporting results from fitting models including the spatially lagged response (SLM, SDM, SAC. GNM) for over ten years (LeSage and Pace 2009). Extension to other models including spatially lagged covariates (SLX, SDEM) has followed (Elhorst 2010; Bivand 2012; Halleck Vega and Elhorst 2015). For SLM, SDM, SAC and GNM models fitted with maximum likelihood or GMM, the variance-covariance matrix of the coefficients is available, and can be used to make random draws from a multivariate Normal distribution with mean set to coefficient values and variance to the estimated variance-covariance matrix. For these models fitted using Bayesian methods, draws are already available. In the SDEM case, the draws on the regression coefficients of the unlagged covariates represent direct impacts, and draws on the coefficients of the spatially lagged covariates represent indirect impacts, and their by-draw sums the total impacts.

Since sampling is not required for inference for SLX and SDEM models, linear combination is used for models fitted using maximum likelihood; results are shown here for the air pollution variable only. The literature has not yet resolved the question of how to report model output, as each covariate is now represented by three impacts. Where spatially lagged covariates are included, two coefficients are replaced by three impacts, here for the air pollution variable of interest.

```
sum_imp_94_SDEM <- summary(impacts(SDEM_94))
rbind(Impacts = sum_imp_94_SDEM$mat[5,], SE = sum_imp_94_SDEM$semat[5,])
```

```
# Direct Indirect Total
# Impacts -0.01276 -0.01845 -0.0312
# SE 0.00235 0.00472 0.0053
```

In the SLX and SDEM models, the direct impacts are the consequences for the response of changes in air pollution in the same observational entity, and the indirect (local) impacts are the consequences for the response of changes in air pollution in neighbouring observational entities.

```
sum_imp_94_SLX <- summary(impacts(SLX_94))
rbind(Impacts = sum_imp_94_SLX$mat[5,], SE = sum_imp_94_SLX$semat[5,])
```

```
# Direct Indirect Total
# Impacts -0.0128 -0.01874 -0.03151
# SE 0.0028 0.00556 0.00611
```

Applying the same approaches to the weighted spatial regressions, the total impacts of air pollution on house values are reduced, but remain significant:

```
sum_imp_94_SDEMw <- summary(impacts(SDEM_94w))
rbind(Impacts = sum_imp_94_SDEMw$mat[5,], SE = sum_imp_94_SDEMw$semat[5,])
```

```
# Direct Indirect Total
# Impacts -0.00592 -0.01076 -0.01668
# SE 0.00269 0.00531 0.00559
```

On balance, using a weighted spatial regression representation including only the spatially lagged covariates aggregated to the air pollution model output zone level seems to clear most of the mis-specification issues, and as R. Bivand (2017) discusses in more detail, gives a willingness to pay for pollution abatement that is much larger than mis-specified alternative models:

```
sum_imp_94_SLXw <- summary(impacts(SLX_94w))
rbind(Impacts = sum_imp_94_SLXw$mat[5,], SE = sum_imp_94_SLXw$semat[5,])
```

```
# Direct Indirect Total
# Impacts -0.00620 -0.01221 -0.01842
# SE 0.00326 0.00628 0.00629
```

## 17.4 Predictions

In the Boston tracts data set, 17 observations of median house values, the response, are censored. We will use the `predict()`

method for `"Sarlm"`

objects to fill in these values; the method was re-written by Martin Gubri based on Goulard, Laurent, and Thomas-Agnan (2017; see also Laurent and Margaretic 2021). The `pred.type=`

argument specifies the prediction strategy among those presented in the article.

Using these as an example and comparing some `pred.type=`

variants for the SDEM model and predicting out-of-sample, we can see that there are differences, suggesting that this is a fruitful area for study. There have been a number of alternative proposals for handling missing variables (Gómez-Rubio, Bivand, and Rue 2015; Suesse 2018). Another reason for increasing attention on prediction is that it is fundamental for machine learning approaches, in which prediction for validation and test data sets drives model specification choice. The choice of training and other data sets with dependent spatial data remains an open question, and is certainly not as simple as with independent data.

Here, we’ll list the predictions for the censored tract observations using three different prediction types, taking the exponent to get back to the USD median house values. Note that the `row.names()`

of the `newdata=`

object are matched with the whole-data spatial weights matrix `"region.id"`

attribute to make out-of-sample prediction possible:

```
nd <- boston_506[is.na(boston_506$median),]
t0 <- exp(predict(SDEM_489, newdata = nd, listw = lw_q, pred.type = "TS", zero.policy =TRUE))
suppressWarnings(t1 <- exp(predict(SDEM_489, newdata = nd, listw = lw_q, pred.type = "KP2",
zero.policy = TRUE)))
suppressWarnings(t2 <- exp(predict(SDEM_489, newdata = nd, listw = lw_q, pred.type = "KP5",
zero.policy = TRUE)))
```

We can also use the `"slm"`

model in INLA to predict missing response values as part of the model fitting function call. A certain amount of set-up code is required as the `"slm"`

model is still experimental:

```
library(INLA)
W <- as(lw_q, "CsparseMatrix")
n <- nrow(W)
e <- eigenw(lw_q)
re.idx <- which(abs(Im(e)) < 1e-6)
rho.max <- 1/max(Re(e[re.idx]))
rho.min <- 1/min(Re(e[re.idx]))
rho <- mean(c(rho.min, rho.max))
boston_506$idx <- 1:n
zero.variance = list(prec=list(initial = 25, fixed = TRUE))
args.slm <- list(rho.min = rho.min, rho.max = rho.max, W = W, X = matrix(0, n, 0),
Q.beta = matrix(1,0,0))
hyper.slm <- list(prec = list(prior = "loggamma", param = c(0.01, 0.01)),
rho = list(initial = 0, prior = "logitbeta", param = c(1,1)))
WX <- create_WX(model.matrix(update(form, CMEDV ~ .), data=boston_506), lw_q)
SDEM_506_slm <- inla(update(form, . ~ . + WX + f(idx, model = "slm", args.slm = args.slm,
hyper = hyper.slm)), data = boston_506, family = "gaussian", control.family =
list(hyper = zero.variance), control.compute = list(dic = TRUE, cpo = TRUE))
mvs <- SDEM_506_slm$marginals.fitted.values[which(is.na(boston_506$median))]
mv_mean <- sapply(mvs, function(x) mean(exp(x[, 1]), ))
mv_sd <- sapply(mvs, function(x) sd(exp(x[, 1])))
```

INLA also provide gridded estimates of the marginal distributions of the predictions, offering a way to assess the uncertainty associated with the predicted values:

```
data.frame(fit_TS = t0[,1], fit_KP2 = c(t1), fit_KP5 = c(t2), INLA_slm = mv_mean,
INLA_slm_sd = mv_sd, censored = boston_506$censored[as.integer(attr(t0, "region.id"))])
```

```
# fit_TS fit_KP2 fit_KP5 INLA_slm INLA_slm_sd censored
# 13 23912 29477 28147 33967 15983 right
# 14 28126 27001 28516 34537 17242 right
# 15 30553 36184 32476 45116 21699 right
# 17 18518 19621 18878 22861 10094 right
# 43 9564 6817 7561 7420 3348 left
# 50 8371 7196 7383 7512 3535 left
# 312 51477 53301 54173 62148 31646 right
# 313 45921 45823 47095 51156 25466 right
# 314 44196 44586 45361 46806 22276 right
# 317 43427 45707 45442 52293 24265 right
# 337 39879 42072 41127 44790 19660 right
# 346 44708 46694 46108 49180 20341 right
# 355 48188 49068 48911 53208 23714 right
# 376 42881 45883 44966 51724 23155 right
# 408 44294 44615 45670 50541 24077 right
# 418 38211 43375 41914 47651 21600 right
# 434 41647 41690 42398 45074 20335 right
```

The spatial regression toolbox remains incomplete, and it will take time to fill in blanks. It remains unfortunate that the several traditions in spatial regression seldom seem to draw on each others’ understandings and advances.

## 17.5 Exercises

- Referring to Piras and Prucha (2014) and Florax, Folmer, and Rey (2003), if we choose to use a pre-test strategy, do linear models of the properties-only data set and the properties with added municipality department variables show residual spatial dependence? Which model specifications might the pre-tests indicate?
- Could the inclusion of municipality department dummies, or a municipality department regimes model assist in reducing residual spatial dependence?
- Attempt to fit a SEM specification by maximum likelihood (see R. Bivand, Millo, and Piras (2021) for GMM code examples) to the properties-only and the properties with added municipality department variables models; extend to an SDEM model. Repeat with SLX models; how might the changes in the tests of residual autocorrelation in the SLX models be interpreted? How might you interpret the highly significant outcomes of Hausman tests on the SEM and SDEM models?
- Fit GNM specifications to the properties-only and the properties with added municipality department variables models; can these models be simplified to say SDM or SDEM representations?
- Do the model estimates reached in the chapter 16 exercises provide more clarity than those in this chapter?

### References

Anselin, L. 1988. *Spatial Econometrics: Methods and Models*. Kluwer Academic Publishers.

Bivand, Roger. 2017. “Revisiting the Boston Data Set — Changing the Units of Observation Affects Estimated Willingness to Pay for Clean Air.” *REGION* 4 (1): 109–27. https://doi.org/10.18335/region.v4i1.107.

Bivand, Roger, Giovanni Millo, and Gianfranco Piras. 2021. “A Review of Software for Spatial Econometrics in R.” *Mathematics* 9 (11). https://doi.org/10.3390/math9111276.

Bivand, Roger, and Gianfranco Piras. 2021. *Spatialreg: Spatial Regression Analysis*. https://CRAN.R-project.org/package=spatialreg.

Bivand, Roger S. 2012. “After ’Raising the Bar’: Applied Maximum Likelihood Estimation of Families of Models in Spatial Econometrics.” *Estadística Española* 54: 71–88.

Bivand, Roger S., and Gianfranco Piras. 2015. “Comparing Implementations of Estimation Methods for Spatial Econometrics.” *Journal of Statistical Software* 63 (1): 1–36. https://doi.org/10.18637/jss.v063.i18.

Bivand, R. S. 2002. “Spatial Econometrics Functions in R: Classes and Methods.” *Journal of Geographical Systems* 4: 405–21.

Cliff, A. D., and J. K. Ord. 1973. *Spatial Autocorrelation*. London: Pion.

Cliff, A. 1981. *Spatial Processes*. London: Pion.

Elhorst, J. Paul. 2010. “Applied Spatial Econometrics: Raising the Bar.” *Spatial Economic Analysis* 5: 9–28.

Fingleton, B. 1999. “Spurious spatial regression: Some Monte Carlo results with a spatial unit root and spatial cointegration.” *Journal of Regional Science* 9: 1–19.

Florax, Raymond J.G.M, Hendrik Folmer, and Sergio J Rey. 2003. “Specification Searches in Spatial Econometrics: The Relevance of Hendry’s Methodology.” *Regional Science and Urban Economics* 33 (5): 557–79. https://doi.org/10.1016/S0166-0462(03)00002-4.

Florax, Raymond J.G.M., Hendrik Folmer, and Sergio J. Rey. 2006. “A Comment on Specification Searches in Spatial Econometrics: The Relevance of Hendry’s Methodology: A Reply.” *Regional Science and Urban Economics* 36 (2): 300–308. https://doi.org/10.1016/j.regsciurbeco.2005.10.002.

Gómez-Rubio, Virgilio, Roger Bivand, and Håvard Rue. 2015. “A New Latent Class to Fit Spatial Econometrics Models with Integrated Nested Laplace Approximations.” *Procedia Environmental Sciences* 27: 116–18. https://doi.org/https://doi.org/10.1016/j.proenv.2015.07.119.

Goulard, Michel, Thibault Laurent, and Christine Thomas-Agnan. 2017. “About Predictions in Spatial Autoregressive Models: Optimal and Almost Optimal Strategies.” *Spatial Economic Analysis* 12 (2-3): 304–25. https://doi.org/10.1080/17421772.2017.1300679.

Halleck Vega, Solmaria, and J. Paul Elhorst. 2015. “The SLX Model.” *Journal of Regional Science* 55 (3): 339–63. https://doi.org/10.1111/jors.12188.

Hendry, David F. 2006. “A Comment on ‘Specification Searches in Spatial Econometrics: The Relevance of Hendry’s Methodology’.” *Regional Science and Urban Economics* 36 (2): 309–12. https://doi.org/10.1016/j.regsciurbeco.2005.10.001.

Kelejian, Harry, and Gianfranco Piras. 2017. *Spatial Econometrics*. London: Academic Press.

Laurent, Thibault, and Paula Margaretic. 2021. “Predictions in Spatial Econometric Models: Application to Unemployment Data.” In *Advances in Contemporary Statistics and Econometrics: Festschrift in Honor of Christine Thomas-Agnan*, edited by Abdelaati Daouia and Anne Ruiz-Gazen, 409–26. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-030-73249-3_21.

LeSage, James P., and Kelley R. Pace. 2009. *Introduction to Spatial Econometrics*. Boca Raton, FL: CRC Press.

LeSage, J. P. 2014. “What Regional Scientists Need to Know About Spatial Econometrics.” *Review of Regional Studies* 44: 13–32. https://journal.srsa.org/ojs/index.php/RRS/article/view/44.1.2.

Martinetti, Davide, and Ghislain Geniaux. 2017. “Approximate Likelihood Estimation of Spatial Probit Models.” *Regional Science and Urban Economics* 64: 30–45. https://doi.org/https://doi.org/10.1016/j.regsciurbeco.2017.02.002.

McMillen, D. P. 2013. *Quantile Regression for Spatial Data*. Heidelberg: Springer-Verlag.

Millo, Giovanni, and Gianfranco Piras. 2012. “splm: Spatial Panel Data Models in R.” *Journal of Statistical Software* 47 (1): 1–38.

Mur, Jesús, and Ana Angulo. 2006. “The Spatial Durbin Model and the Common Factor Tests.” *Spatial Economic Analysis* 1 (2): 207–26. https://doi.org/10.1080/17421770601009841.

Ord, J. K. 1975. “Estimation Methods for Models of Spatial Interaction.” *Journal of the American Statistical Association* 70 (349): 120–26.

Pace, RK, and JP LeSage. 2008. “A Spatial Hausman Test.” *Economics Letters* 101: 282–84.

Piras, Gianfranco, and Ingmar R. Prucha. 2014. “On the Finite Sample Properties of Pre-Test Estimators of Spatial Models.” *Regional Science and Urban Economics* 46: 103–15. https://doi.org/10.1016/j.regsciurbeco.2014.03.002.

Smith, T. E., and K. L. Lee. 2012. “The effects of spatial autoregressive dependencies on inference in ordinary least squares: a geometric approach.” *Journal of Geographical Systems* 14 (January): 91–124. https://doi.org/10.1007/s10109-011-0152-x.

Smith, Tony E. 2009. “Estimation Bias in Spatial Models with Strongly Connected Weight Matrices.” *Geographical Analysis* 41 (3): 307–32. https://doi.org/10.1111/j.1538-4632.2009.00758.x.

Suesse, Thomas. 2018. “Marginal Maximum Likelihood Estimation of Sar Models with Missing Data.” *Computational Statistics & Data Analysis* 120: 98–110. https://doi.org/https://doi.org/10.1016/j.csda.2017.11.004.

Wagner, Martin, and Achim Zeileis. 2019. “Heterogeneity and Spatial Dependence of Regional Growth in the EU: A Recursive Partitioning Approach.” *German Economic Review* 20 (1): 67–82. https://doi.org/10.1111/geer.12146.

Waller, Lance A., and Carol A. Gotway. 2004. *Applied Spatial Statistics for Public Health Data*. Hoboken, NJ: John Wiley & Sons.

Wilhelm, Stefan, and Miguel Godinho de Matos. 2013. “Estimating Spatial Probit Models in R.” *The R Journal* 5 (1): 130–43. https://doi.org/10.32614/RJ-2013-013.