# 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:

broom::tidy(lmtest::lrtest(SEM_489, SDEM_489))
# # 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
broom::tidy(lmtest::lrtest(SEM_94, SDEM_94))
# # 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

1. 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?
2. Could the inclusion of municipality department dummies, or a municipality department regimes model assist in reducing residual spatial dependence?
3. 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?
4. 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?
5. 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.