Banner

3.4 Spatial Regression Model Families¶

Part: 3 — Spatial Regression (Chapter 3.4 of 7)

Dataset: diabetes_atlantic.shp (666 counties)

Focus: Fitting the spatial lag, spatial error, and Durbin-family models, and what each model's coefficients actually mean

Library: spatialstats.regression (PySpatialStats)

Table of Contents¶

  1. Overview

  2. Python Setup

  3. The Spatial Lag Model (SAR / SLM))

    • 3.1 Interpreting coefficients: direct, indirect, and total impacts

  4. The Spatial Error Model (SEM))

    • 4.1 Filtered vs. raw residuals

  5. Extensions: SLX, Spatial Durbin, and Durbin Error

    • 5.1 The common-factor idea: does SDM reduce to SEM?

  6. What This Library Does Not Implement

  7. Summary and Conclusions

  8. Resources

1. Overview¶

Chapter 3.3 recommended the spatial error model. This chapter fits it — and, for comparison and because the robust LM tests left room for a joint interpretation, the spatial lag model and both Durbin extensions as well — with the emphasis on what each model's fitted parameters actually mean, which is where the two families diverge most from OLS intuition.

Step What we do
SLM Fit $y = \rho Wy + X\beta + \varepsilon$ by maximum likelihood, and unpack why $\beta$ alone is not a marginal effect
Impacts Direct, indirect (spillover), and total effects, with simulation-based intervals
SEM Fit $y = X\beta + u,\ u = \lambda Wu + \varepsilon$, and show its spatial filter directly removing residual autocorrelation
Extensions SLX, Spatial Durbin (SDM), and Durbin Error (SDEM), plus the "common factor" relationship between SDM and SEM
Scope What this library's regression module does not implement, and why that matters for model choice

2. Python Setup¶

In [1]:
import sys, warnings
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt

warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=DeprecationWarning)


def find_repo_root(start: Path = Path.cwd()) -> Path:
    for p in [start, *start.parents]:
        if (p / "pyproject.toml").exists() and (p / "data").is_dir():
            return p
    raise FileNotFoundError("Run this notebook from inside the PySpatialStats repository.")


ROOT = find_repo_root()
DATA = ROOT / "data"
try:
    import spatialstats
except ImportError:
    sys.path.insert(0, str(ROOT))
    import spatialstats

from spatialstats.weights import Queen
from spatialstats.esda import Moran
from spatialstats.regression import OLS, SpatialLag, SpatialError, SpatialDurbin, SLX

plt.rcParams.update({
    "figure.dpi": 100, "savefig.dpi": 100, "font.size": 10, "axes.titlesize": 11,
    "axes.titleweight": "bold", "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25, "legend.frameon": False,
})
C = {"blue": "#2b6cb0", "red": "#c0392b", "green": "#2f855a", "orange": "#dd6b20",
     "grey": "#718096", "purple": "#6b46c1", "teal": "#2c7a7b"}

counties = gpd.read_file(DATA / "diabetes_atlantic.shp")
covariates = ["Obesity", "PhysInact", "Acc_Exer", "PCP_rate", "FoodEnvIdx", "SVI"]
X = counties[covariates].copy()
X["Urban"] = (counties["UrbRural"] == "Urban").astype(int)
y = counties["Dia_pct"]
w = Queen(counties)
print(f"spatialstats : version {spatialstats.__version__}")

spatialstats : version 0.2.0

3. The Spatial Lag Model (SAR / SLM)¶

$$ y = \rho W y + X\beta + \varepsilon, \qquad \varepsilon \sim N(0, \sigma^2 I). $$

$Wy$ is endogenous — it is correlated with $\varepsilon$ through the same simultaneity that makes $y$ depend on its own neighbours — so OLS on this equation would be inconsistent. SpatialLag instead maximises the exact Gaussian log-likelihood, which includes a log-Jacobian term $\ln|I - \rho W|$ that OLS's loss function does not have; $\rho$'s feasible range comes from the extreme eigenvalues of $W$, and asymptotic standard errors come from the (numerically inverted) information matrix.

In [2]:
slm = SpatialLag(y, X, w)
print(slm.summary())

Spatial Lag Model (ML)   (dependent variable: Dia_pct)
======================================================
  Observations      666    Log likelihood    -767.7914
  Parameters          9    AIC               1553.5828
  Pseudo R²      0.7607    BIC               1594.0944
  rho            0.1784    std err   0.0312    p = 1.093e-08

  variable                 coef    std err        z         p
  const                  1.7706     0.4729    3.744    0.0002 ***
  Obesity                0.1368     0.0110   12.454    0.0000 ***
  PhysInact              0.1415     0.0143    9.881    0.0000 ***
  Acc_Exer              -0.0039     0.0016   -2.463    0.0138 *
  PCP_rate              -0.0005     0.0009   -0.548    0.5838 
  FoodEnvIdx            -0.1818     0.0443   -4.100    0.0000 ***
  SVI                    0.9709     0.1631    5.953    0.0000 ***
  Urban                 -0.0159     0.0724   -0.219    0.8264

$\rho = 0.178$ (s.e. 0.031, $p<0.001$): each county's diabetes prevalence is, net of its own covariates, pulled toward its neighbours' prevalence. Compare the log-likelihood to Chapter 3.2's OLS ($-784.2 \to -767.8$, an improvement Chapter 3.5 formalises with a likelihood-ratio test) and the reduction in residual autocorrelation.

In [3]:
print(f"OLS  residual Moran's I : {OLS(y, X, w=w).residual_moran(w)['I']:.4f}")
print(f"SLM  residual Moran's I : {slm.residual_moran(w)['I']:.4f}")

OLS  residual Moran's I : 0.1812
SLM  residual Moran's I : 0.0753

3.1 Interpreting coefficients: direct, indirect, and total impacts¶

Because $y = (I-\rho W)^{-1}(X\beta + \varepsilon)$, changing $x_k$ in county $i$ affects $y_i$ and, through $\rho$, every county connected to $i$ through the neighbour graph. $\beta_k$ alone is not the effect of $x_k$; impacts() decomposes it into:

  • direct — the average effect of a change in $x_k$ at a place on its own $y$ (including the feedback loop where $i$'s change affects its neighbours, who affect $i$ back);
  • indirect — the average effect on other places' $y$ from the same change (the spillover);
  • total = direct + indirect.
In [4]:
impacts = slm.impacts(seed=1)
impacts.round(4)
Out[4]:
variable direct direct_lo direct_hi indirect indirect_lo indirect_hi total total_lo total_hi
0 Obesity 0.1376 0.1125 0.1574 0.0288 0.0178 0.0408 0.1665 0.1372 0.1938
1 PhysInact 0.1424 0.1185 0.1712 0.0298 0.0183 0.0445 0.1722 0.1405 0.2094
2 Acc_Exer -0.0040 -0.0071 -0.0010 -0.0008 -0.0016 -0.0002 -0.0048 -0.0084 -0.0012
3 PCP_rate -0.0005 -0.0022 0.0013 -0.0001 -0.0005 0.0003 -0.0006 -0.0028 0.0016
4 FoodEnvIdx -0.1830 -0.2626 -0.1026 -0.0383 -0.0626 -0.0186 -0.2213 -0.3214 -0.1240
5 SVI 0.9771 0.6868 1.2617 0.2046 0.1155 0.3244 1.1817 0.8192 1.5363
6 Urban -0.0160 -0.1570 0.1260 -0.0033 -0.0337 0.0278 -0.0193 -0.1920 0.1492

For Obesity, the direct effect (0.138) is barely distinguishable from the SLM's own raw coefficient $\beta$ (0.137, from the summary() table above) — the feedback loop where a county's own change echoes back to itself through its neighbours is a tiny correction, because Queen weights are sparse and only mildly clustered. The indirect effect (0.029), by contrast, is not part of $\beta$ at all: it is the average spillover onto every other county from a one-unit change in one county's obesity rate, propagated through $\rho$. For a row-standardised $W$, the exact identity total $= \beta / (1-\rho)$ holds:

In [5]:
beta_series = pd.Series(slm.beta, index=slm.names)
check = impacts.assign(beta=lambda d: d["variable"].map(beta_series),
                        beta_over_1_minus_rho=lambda d: d["variable"].map(beta_series) / (1 - slm.rho))
check[["variable", "beta", "total", "beta_over_1_minus_rho"]].round(6)
Out[5]:
variable beta total beta_over_1_minus_rho
0 Obesity 0.136762 0.166452 0.166452
1 PhysInact 0.141505 0.172226 0.172226
2 Acc_Exer -0.003931 -0.004784 -0.004784
3 PCP_rate -0.000480 -0.000584 -0.000584
4 FoodEnvIdx -0.181827 -0.221301 -0.221301
5 SVI 0.970915 1.181695 1.181695
6 Urban -0.015871 -0.019317 -0.019317

Exact agreement, as the identity predicts. The gap between $\beta$ and the total impact — here, a multiplier of $1/(1-0.178) \approx 1.217$ — is the "spatial multiplier": a spillover process amplifies every regressor's effect uniformly by the same factor, on top of whatever direct effect it already had.

4. The Spatial Error Model (SEM)¶

$$ y = X\beta + u, \qquad u = \lambda W u + \varepsilon, \qquad \varepsilon \sim N(0, \sigma^2 I). $$

Here spatial dependence lives entirely in the error, not in $y$ itself. $\beta$ keeps its ordinary OLS interpretation — a marginal effect, no spillover, no impacts() needed — and $\lambda \ne 0$ is read as "some spatially structured influence is missing from $X$", not as evidence that $\beta$ is wrong.

In [6]:
sem = SpatialError(y, X, w)
print(sem.summary())

Spatial Error Model (ML)   (dependent variable: Dia_pct)
========================================================
  Observations      666    Log likelihood    -761.8263
  Parameters          9    AIC               1541.6527
  Pseudo R²      0.7454    BIC               1582.1643
  lam            0.3521    std err   0.0510    p = 5.244e-12

  variable                 coef    std err        z         p
  const                  2.8805     0.4631    6.220    0.0000 ***
  Obesity                0.1478     0.0117   12.627    0.0000 ***
  PhysInact              0.1501     0.0150    9.999    0.0000 ***
  Acc_Exer              -0.0047     0.0017   -2.785    0.0053 **
  PCP_rate              -0.0001     0.0009   -0.113    0.9097 
  FoodEnvIdx            -0.1669     0.0470   -3.553    0.0004 ***
  SVI                    0.9464     0.1731    5.469    0.0000 ***
  Urban                 -0.0657     0.0745   -0.882    0.3775

$\lambda = 0.352$ (s.e. 0.051): stronger than the SLM's $\rho$, consistent with Chapter 3.3's robust-test comparison favouring the error interpretation. Compare $\beta$ coefficients directly against OLS — this is the practical signature of a genuine SEM situation (Chapter 3.7, Exercise 1 asks you to verify this explicitly):

In [7]:
ols = OLS(y, X, w=w)
coef_compare = pd.DataFrame({
    "OLS": pd.Series(ols.beta, index=ols.names),
    "SEM": pd.Series(sem.beta, index=sem.names),
    "SLM": pd.Series(slm.beta, index=slm.names),
}).round(4)
coef_compare
Out[7]:
OLS SEM SLM
const 3.1546 2.8805 1.7706
Obesity 0.1462 0.1478 0.1368
PhysInact 0.1520 0.1501 0.1415
Acc_Exer -0.0060 -0.0047 -0.0039
PCP_rate -0.0003 -0.0001 -0.0005
FoodEnvIdx -0.2002 -0.1669 -0.1818
SVI 1.0327 0.9464 0.9709
Urban -0.0144 -0.0657 -0.0159

OLS and SEM coefficients are close throughout (e.g. Obesity: 0.146 vs. 0.148); SLM's are somewhat smaller across the board (0.137) because part of what OLS attributed to the covariates, SLM instead attributes to the $\rho Wy$ term. This pattern — SEM close to OLS, SLM systematically different — is exactly what Chapter 3.3's diagnosis (error, not lag) predicts.

4.1 Filtered vs. raw residuals¶

SpatialError exposes both residuals_raw ($y - X\beta$, before any spatial adjustment) and residuals (the spatially filtered residual $\hat\varepsilon = (I-\hat\lambda W)(y - X\beta)$, which the model claims is i.i.d.). Comparing their autocorrelation is the most direct demonstration of what $\lambda$ is doing.

In [8]:
moran_raw = Moran(sem.residuals_raw, w, permutations=0)
moran_filtered = Moran(sem.residuals, w, permutations=0)
print(f"raw residuals (y - X*beta)         Moran's I = {moran_raw.I:.4f}")
print(f"filtered residuals (I - lam*W)(..) Moran's I = {moran_filtered.I:.4f}")

raw residuals (y - X*beta)         Moran's I = 0.2011
filtered residuals (I - lam*W)(..) Moran's I = -0.0216

The raw residuals are about as autocorrelated as OLS's (I=0.201, essentially the same problem OLS has), but the spatially filtered residuals are indistinguishable from zero (I=−0.022) — a direct, visual confirmation that the SEM's error structure is doing exactly what it claims to.

5. Extensions: SLX, Spatial Durbin, and Durbin Error¶

Three further specifications add spatially lagged covariates ($WX$) — a different, additive way of letting neighbours' characteristics matter, without the endogeneity of a lagged $y$:

  • SLX: $y = X\beta + WX\theta + \varepsilon$ — plain OLS, since $WX$ is a (pre-computed, exogenous) regressor, not an endogenous term.
  • Spatial Durbin (SDM): SLM plus $WX\theta$ — SpatialLag(slx=True) / SpatialDurbin(...).
  • Durbin Error (SDEM): SEM plus $WX\theta$ — SpatialError(slx=True).
In [9]:
slx = SLX(y, X, w=w)
sdm = SpatialDurbin(y, X, w)
sdem = SpatialError(y, X, w, slx=True)

for label, m in [("SLX", slx), ("SDM", sdm), ("SDEM", sdem)]:
    sp = f", rho={m.rho:.3f}" if hasattr(m, "rho") else (f", lam={m.lam:.3f}" if hasattr(m, "lam") else "")
    print(f"{label:5s}  loglik={m.loglik:9.2f}  AIC={m.aic:9.2f}  n_params={m.n_params:3d}{sp}")

SLX    loglik=  -772.81  AIC=  1575.62  n_params= 15
SDM    loglik=  -753.11  AIC=  1538.22  n_params= 16, rho=0.321
SDEM   loglik=  -753.87  AIC=  1539.74  n_params= 16, lam=0.320

Every Durbin-family model roughly doubles the parameter count (one $W$-lagged term per original covariate: 16 parameters for SDM/SDEM against SLM/SEM's 9), which the AIC/BIC comparability discussed in Chapter 3.2 accounts for automatically — but with 7 extra $WX$ terms, individual coefficients are rarely all significant (Chapter 3.5's compare_models table quantifies the trade-off between the improved fit and the larger, harder-to-interpret model).

5.1 The common-factor idea: does SDM reduce to SEM?¶

A classical result connects SDM back to SEM: if the true model is SEM, an SDM fit on the same data should satisfy the common factor restriction $\theta \approx -\rho\beta$ — the $WX$ coefficients are approximately the negative of $\rho$ times the original coefficients. Checking this informally, term by term:

In [10]:
names = sdm.names
p = (len(names) - 1) // 2  # excluding the constant
beta_sdm = sdm.beta[1:1 + p]
theta_sdm = sdm.beta[1 + p:1 + 2 * p]
common_factor = pd.DataFrame({
    "variable": names[1:1 + p],
    "beta": beta_sdm.round(4),
    "theta (fitted)": theta_sdm.round(4),
    "-rho*beta (predicted under common factor)": (-sdm.rho * beta_sdm).round(4),
})
common_factor
Out[10]:
variable beta theta (fitted) -rho*beta (predicted under common factor)
0 Obesity 0.1498 -0.0457 -0.0481
1 PhysInact 0.1492 -0.0457 -0.0479
2 Acc_Exer -0.0017 -0.0060 0.0005
3 PCP_rate -0.0002 0.0004 0.0001
4 FoodEnvIdx -0.1221 0.0170 0.0392
5 SVI 0.8196 0.1141 -0.2630
6 Urban -0.0638 0.1480 0.0205

The restriction holds reasonably well for Obesity and PhysInact (fitted $\theta$ close to the common-factor prediction) but fails clearly for SVI (fitted $\theta=0.114$ vs. a predicted $-0.263$ — opposite sign) and several others. A formal common-factor likelihood-ratio test (comparing SDM's log-likelihood to a version constrained to satisfy $\theta=-\rho\beta$ exactly) is not implemented in this library, but the informal pattern here already says the same thing Chapter 3.5's AIC/BIC comparison will confirm: SDM is not simply SEM in disguise for this dataset — it is capturing genuine additional structure, particularly through SVI, that a pure SEM specification would miss.

6. What This Library Does Not Implement¶

Two entries from Chapter 3.1's taxonomy have no estimator here:

  • SARAR ($y=\rho Wy + X\beta + u,\ u=\lambda Wu+\varepsilon$) — both mechanisms simultaneously. When Chapter 3.3's LM-SARMA test is significant (it is: 59.4, $p\approx1.3\times10^{-13}$) but the robust tests disagree on which single mechanism dominates, a SARAR model is the textbook next step; here, the Durbin models (Section 5) are the available substitute, since they can absorb some of the same joint structure through $WX$.
  • 2SLS/GMM estimation — this library's SLM and SEM are both maximum-likelihood only. GMM-based spatial estimators (Kelejian & Prucha) are a common alternative in the wider literature, particularly favoured for large $n$ where ML's dense-matrix operations (Chapter 3.2's dense_inverse_ok guard, $n \le 6000$) become costly.

Knowing what is not available is as important as knowing what is: a workflow that needs SARAR should treat this library's Durbin-model fits as an approximation, and cross-check important conclusions against a package that implements it directly (Part 2's to_libpysal() interoperability is the natural bridge).

7. Summary and Conclusions¶

This chapter fitted every model family in this library's spatial regression taxonomy and unpacked what their parameters mean.

What we did¶

  1. Fit the spatial lag model: $\rho=0.178$ (s.e. 0.031), residual Moran's I falling from OLS's 0.181 to 0.075.
  2. Computed direct/indirect/total impacts, and verified the exact identity total $=\beta/(1-\rho)$ for every covariate under row-standardised weights.
  3. Fit the spatial error model: $\lambda=0.352$ (s.e. 0.051), and showed its spatially filtered residuals (I=−0.022) are indistinguishable from zero while the raw, unfiltered residuals (I=0.201) are not.
  4. Compared OLS, SEM, and SLM coefficients directly: SEM tracks OLS closely (as SEM-appropriate data should), SLM's coefficients differ systematically.
  5. Fit SLX, SDM, and SDEM, and checked the common-factor restriction connecting SDM to SEM informally — it holds for some covariates (Obesity, PhysInact) but clearly fails for others (SVI), indicating SDM captures real structure beyond a disguised SEM.
  6. Noted what this library does not implement (SARAR, GMM estimation) and when that gap matters.

Practical takeaways¶

  • Never report a spatial lag model's raw $\beta$ as "the effect of $x$" — use impacts(), and report direct, indirect, and total separately; they answer different questions.
  • SEM coefficients are directly comparable to OLS; SLM coefficients are not, by construction.
  • Comparing filtered vs. raw SEM residuals is a fast, concrete way to confirm a fitted spatial error model is actually doing its job, beyond trusting the summary table.

Next steps¶

Chapter 3.5 (Comparison, Reporting and Diagnostics) puts every model fitted in this chapter side by side — compare_models, lr_test, and the acceptance criterion that residual Moran's I should be close to zero after fitting.

8. Resources¶

Textbooks¶

Reference Notes
LeSage, J. & Pace, R. K. (2009). Introduction to Spatial Econometrics. CRC Press. The definitive reference for impacts and the direct/indirect/total decomposition
Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer Academic. Maximum-likelihood estimation of SAR and SEM, including the log-Jacobian term
Elhorst, J. P. (2014). Spatial Econometrics: From Cross-Sectional Data to Spatial Panels. Springer. The common-factor restriction connecting SDM and SEM

Key papers¶

Paper Contribution
LeSage, J. P. & Pace, R. K. (2009). Chapter 2. Interpreting spatial models. In Introduction to Spatial Econometrics. Origin of the impacts framework used in Section 3.1
Burridge, P. (1981). Testing for a common factor in a spatial autoregression model. Environment and Planning A, 13, 795–800. Formal common-factor test referenced in Section 5.1
Kelejian, H. H. & Prucha, I. R. (1998). A generalized spatial two-stage least squares procedure. Journal of Real Estate Finance and Economics, 17, 99–121. GMM/2SLS alternative not implemented in this library (Section 6)

In this library¶

Object Role
spatialstats.regression.SpatialLag SAR/SLM by ML; .impacts() for direct/indirect/total effects
spatialstats.regression.SpatialError SEM by ML; .residuals (filtered) vs. .residuals_raw
spatialstats.regression.SpatialDurbin SDM = SpatialLag(slx=True)
spatialstats.regression.SLX OLS with WX; SpatialError(slx=True) gives SDEM
Chapter 3.5 compare_models, lr_test — the side-by-side comparison across every model fitted here