options(repr.plot.width=8.2, repr.plot.height=3.4)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(sf); library(spdep); library(spatialreg); library(spData); library(ggplot2); library(patchwork)})
col <- st_read(system.file("shapes/columbus.gpkg", package="spData"), quiet=TRUE)
st_crs(col) <- NA                                              # Columbus coords are arbitrary planar units; drop CRS so ggplot maps don't attempt a PROJ transformation
nb <- poly2nb(col); lw <- nb2listw(nb, style="W")
cat(nrow(col), "neighbourhoods; crime ~ income (INC) + housing value (HOVAL)\n")

Warning message:
"package 'spatialreg' was built under R version 4.6.1"

49 neighbourhoods; crime ~ income (INC) + housing value (HOVAL)

ols <- lm(CRIME ~ INC + HOVAL, data=col); col$resid <- resid(ols)
cat("OLS coefficients:\n"); print(round(coef(ols),2))
cat(sprintf("\nMoran's I on OLS residuals: p = %.4f (spatial autocorrelation present)\n", lm.morantest(ols, lw)$p.value))
lmt <- tryCatch(lm.RStests(ols, lw, test="all"), error=function(e) lm.LMtests(ols, lw, test="all"))
cat("\nLagrange-multiplier tests (small p = that alternative is needed):\n")
print(sapply(lmt, function(t) round(t$p.value,4)))
p1<-ggplot(col)+geom_sf(aes(fill=CRIME),color="white",linewidth=0.2)+scale_fill_distiller(palette="OrRd",direction=1)+labs(title="Crime",fill="")+theme_void()
p2<-ggplot(col)+geom_sf(aes(fill=resid),color="white",linewidth=0.2)+scale_fill_distiller(palette="RdBu")+labs(title="OLS residuals (clustered)",fill="")+theme_void()
print(p1 + p2)
cat("\nThe residuals are spatially autocorrelated, so OLS is mis-specified; the robust LM tests indicate which\n")
cat("spatial specification (lag or error) to prefer.\n")

OLS coefficients:

(Intercept)         INC       HOVAL 
      68.62       -1.60       -0.27

Moran's I on OLS residuals: p = 0.0023 (spatial autocorrelation present)

Lagrange-multiplier tests (small p = that alternative is needed):

   RSerr    RSlag adjRSerr adjRSlag    SARMA 
  0.0225   0.0029   0.8340   0.0533   0.0114

The residuals are spatially autocorrelated, so OLS is mis-specified; the robust LM tests indicate which

spatial specification (lag or error) to prefer.

m_lag <- lagsarlm(CRIME ~ INC + HOVAL, data=col, listw=lw)
cat(sprintf("spatial-lag rho = %.2f  (p = %.4f)\n", m_lag$rho, summary(m_lag)$LR1$p.value))
imp <- impacts(m_lag, listw=lw, R=2000)
cat("\ndirect / indirect / total effects:\n"); print(round(imp$res$direct,2));
eff <- data.frame(direct=imp$res$direct, indirect=imp$res$indirect, total=imp$res$total, row.names=c("INC","HOVAL"))
print(round(eff,2))
cat("\nAs in the Python notebook, income's TOTAL effect on crime exceeds its coefficient once the spillover\n")
cat("(indirect) part is included -- the effect, not the coefficient, is the policy-relevant quantity.\n")

spatial-lag rho = 0.42  (p = 0.0022)

direct / indirect / total effects:

      [,1]
[1,] -1.10
[2,] -0.28

      direct indirect total
INC    -1.10    -0.72 -1.82
HOVAL  -0.28    -0.18 -0.46

As in the Python notebook, income's TOTAL effect on crime exceeds its coefficient once the spillover

(indirect) part is included -- the effect, not the coefficient, is the policy-relevant quantity.

m_err <- errorsarlm(CRIME ~ INC + HOVAL, data=col, listw=lw)
cat(sprintf("spatial-error lambda = %.2f  (p = %.4f)\n\n", m_err$lambda, summary(m_err)$LR1$p.value))
comp <- rbind(OLS=coef(ols)[c("INC","HOVAL")],
              SAR=coef(m_lag)[c("INC","HOVAL")],
              SEM=coef(m_err)[c("INC","HOVAL")])
print(round(comp,2))
cat(sprintf("\nAIC: OLS %.0f, SAR %.0f, SEM %.0f  (lower = better)\n", AIC(ols), AIC(m_lag), AIC(m_err)))
cat(sprintf("Both improve on OLS, and both criteria point the same way. The ROBUST LM tests are the ones to\n"))
cat(sprintf("read: adjRSerr p = %.2f says no error term is needed once a lag is allowed, while adjRSlag p = %.3f\n",
            lmt[["adjRSerr"]]$p.value, lmt[["adjRSlag"]]$p.value))
cat(sprintf("says a lag is still needed once an error term is allowed. AIC agrees: SAR %.0f < SEM %.0f < OLS %.0f.\n",
            AIC(m_lag), AIC(m_err), AIC(ols)))
cat("So on Columbus the evidence favours the LAG specification -- neighbours' crime behaving as a genuine\n")
cat("spillover rather than a correlated omitted factor. That is a statistical verdict, not a causal one:\n")
cat("the tests choose between two error structures, and only theory can license the causal reading.\n")

spatial-error lambda = 0.55  (p = 0.0071)

      INC HOVAL
OLS -1.60 -0.27
SAR -1.05 -0.27
SEM -0.96 -0.30

AIC: OLS 383, SAR 375, SEM 377  (lower = better)

Both improve on OLS, and both criteria point the same way. The ROBUST LM tests are the ones to

read: adjRSerr p = 0.83 says no error term is needed once a lag is allowed, while adjRSlag p = 0.053

says a lag is still needed once an error term is allowed. AIC agrees: SAR 375 < SEM 377 < OLS 383.

So on Columbus the evidence favours the LAG specification -- neighbours' crime behaving as a genuine

spillover rather than a correlated omitted factor. That is a statistical verdict, not a causal one:

the tests choose between two error structures, and only theory can license the causal reading.

Spatial Econometrics (R)¶

`spatialreg`: lag, error, LM tests, and impacts¶

1. OLS, residual autocorrelation, and the lag-vs-error tests¶

2. Spatial lag (SAR) and its impacts¶

3. Spatial error (SEM), and the comparison¶

4. Summary¶

Spatial Econometrics (R)¶

spatialreg: lag, error, LM tests, and impacts¶

1. OLS, residual autocorrelation, and the lag-vs-error tests¶

2. Spatial lag (SAR) and its impacts¶

3. Spatial error (SEM), and the comparison¶

4. Summary¶

`spatialreg`: lag, error, LM tests, and impacts¶