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 standardisation 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 down-weight 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 biased 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 earliest 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}, \ \ {\mathbf y} = {\mathbf X}{\mathbf \beta} + ({\mathbf I} - \rho_{\mathrm{Err}} {\mathbf W})^{-1} {\mathbf \varepsilon}, \ \ ({\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 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 SEM with \({\mathbf \gamma} = 0\) and \(\rho_{\mathrm{Lag}} = 0\), and the 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 taken 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 dataset. 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 a 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 Raymond J. G. M. Florax, Folmer, and Rey (2003) (see also comments by Hendry (2006) and Raymond J. G. M. 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 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, by the new package spldv (Sarrias and Piras 2022), 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 generalised method of moments (GMM) estimators in spatialreg (Bivand and Piras 2022) 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 Bivand, Millo, and Piras (2021) covers many of the features of R packages for spatial econometrics, updating 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 Bivand (2017) using the Boston house value dataset.

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 datasets. 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 datasets 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 optimisation, using a heuristic to choose starting values. Like lagsarlm(), this function does not take a weights= argument.

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

Boston house value dataset examples

The examples use the objects read and created in Chapter 16, based on 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 dataset 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))
# Warning in RET$pfunction("adjusted", ...): Completion with error >
# abseps

# Warning in RET$pfunction("adjusted", ...): Completion with error >
# abseps

# Warning in RET$pfunction("adjusted", ...): Completion with error >
# abseps

# Warning in RET$pfunction("adjusted", ...): Completion with error >
# abseps

# Warning in RET$pfunction("adjusted", ...): Completion with error >
# abseps
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:

o <- lmtest::lrtest(SEM_489, SDEM_489)
attr(o, "heading")[2] <- "Model 1: SEM_489\nModel 2: SDEM_489"
o
# Likelihood ratio test
# 
# Model 1: SEM_489
# Model 2: SDEM_489
#   #Df LogLik Df Chisq Pr(>Chisq)    
# 1  16    274                        
# 2  29    311 13  74.4    1.2e-10 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
o <- lmtest::lrtest(SEM_94, SDEM_94)
attr(o, "heading")[2] <- "Model 1: SEM_94\nModel 2: SDEM_94"
o
# Likelihood ratio test
# 
# Model 1: SEM_94
# Model 2: SDEM_94
#   #Df LogLik Df Chisq Pr(>Chisq)    
# 1  16   59.7                        
# 2  29   81.3 13  43.2    4.2e-05 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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 dataset, we obviously do:

SLX_489 <- lmSLX(form, data = boston_489, listw = lw_q_489,
                 zero.policy = TRUE)
o <- lmtest::lrtest(SLX_489, SDEM_489)
attr(o, "heading")[2] <- "Model 1: SLX_489\nModel 2: SDEM_489"
o
# Likelihood ratio test
# 
# Model 1: SLX_489
# Model 2: SDEM_489
#   #Df LogLik Df Chisq Pr(>Chisq)    
# 1  28    231                        
# 2  29    311  1   159     <2e-16 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

but in the output zone case, we do not.

SLX_94 <- lmSLX(form, data = boston_94, listw = lw_q_94)
o <- lmtest::lrtest(SLX_94, SDEM_94)
attr(o, "heading")[2] <- "Model 1: SLX_94\nModel 2: SDEM_94"
o
# Likelihood ratio test
# 
# Model 1: SLX_94
# Model 2: SDEM_94
#   #Df LogLik Df Chisq Pr(>Chisq)
# 1  28   81.2                    
# 2  29   81.3  1  0.22       0.64

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))
o <- lmtest::lrtest(SLX_489w, SDEM_489w)
attr(o, "heading")[2] <- "Model 1: SLX_489w\nModel 2: SDEM_489w"
o
# Likelihood ratio test
# 
# Model 1: SLX_489w
# Model 2: SDEM_489w
#   #Df LogLik Df Chisq Pr(>Chisq)    
# 1  28    311                        
# 2  29    379  1   136     <2e-16 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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))
o <- lmtest::lrtest(SLX_94w, SDEM_94w)
attr(o, "heading")[2] <- "Model 1: SLX_94w\nModel 2: SDEM_94w"
o
# Likelihood ratio test
# 
# Model 1: SLX_94w
# Model 2: SDEM_94w
#   #Df LogLik Df Chisq Pr(>Chisq)
# 1  28   97.5                    
# 2  29   98.0  1  0.92       0.34

In this case and based on arguments advanced in 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 we 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 10 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 misspecification issues, and as Bivand (2017) discusses in more detail, it gives a willingness to pay for pollution abatement that is much larger than misspecified 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 dataset, 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 rewritten 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 datasets drives model specification choice. The choice of training and other datasets 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)
# Loading required package: foreach
# Loading required package: parallel
# Loading required package: sp
# This is INLA_23.04.24 built 2023-04-24 19:15:35 UTC.
#  - See www.r-inla.org/contact-us for how to get help.
#  - To enable PARDISO sparse library; see inla.pardiso()
# 
# Attaching package: 'INLA'
# The following object is masked _by_ '.GlobalEnv':
# 
#     f
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))
mv_mean <- exp(SDEM_506_slm$summary.fitted.values$mean[
               which(is.na(boston_506$median))])

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, censored = 
      boston_506$censored[as.integer(attr(t0, "region.id"))])
#     fit_TS fit_KP2 fit_KP5 INLA_slm censored
# 13   23912   29477   28147    30638    right
# 14   28126   27001   28516    31260    right
# 15   30553   36184   32476    39394    right
# 17   18518   19621   18878    20490    right
# 43    9564    6817    7561     7169     left
# 50    8371    7196    7383     6952     left
# 312  51477   53301   54173    56468    right
# 313  45921   45823   47095    46689    right
# 314  44196   44586   45361    42808    right
# 317  43427   45707   45442    47503    right
# 337  39879   42072   41127    41215    right
# 346  44708   46694   46108    45363    right
# 355  48188   49068   48911    48908    right
# 376  42881   45883   44966    46894    right
# 408  44294   44615   45670    46538    right
# 418  38211   43375   41914    42808    right
# 434  41647   41690   42398    41257    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 Raymond J. G. M. Florax, Folmer, and Rey (2003), if we choose to use a pre-test strategy, do linear models of the properties-only dataset 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 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?