3.7 Code Skeleton and Exercises

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

Dataset: diabetes_atlantic.shp (666 counties)

Focus: A compact, reusable workflow skeleton, and five solved exercises spanning every chapter in this Part

Library: spatialstats.regression, spatialstats.weights (PySpatialStats)

---

Table of Contents

1. Overview

2. Python Setup

3. The Code Skeleton

4. Exercises

   • 4.1 SEM vs. SLM: which stays close to OLS?

   • 4.2 LM tests vs. AIC: do they agree between SLM and SEM?

   • 4.3 The impacts identity under a non-row-standardised W

   • 4.4 Refitting with KNN(6): do conclusions change?:-do-conclusions-change?)

   • 4.5 Why AIC and BIC disagree about the SDM

5. Summary and Conclusions

6. Resources

---

1. Overview

This closing chapter has two parts: a minimal, reusable skeleton that reproduces the entire Chapters 3.2–3.5 workflow in about a dozen lines, and five solved exercises that revisit specific claims from earlier chapters with fresh checks — including one exercise whose answer complicates the textbook expectation in an honest, data-driven way (Section 4.1).

Step	What we do
Skeleton	The complete fit → diagnose → compare → interpret pipeline, condensed
Exercises	Five solved problems spanning Chapters 3.2 (OLS), 3.3 (LM tests), 3.4 (impacts), and 3.5 (comparison)

---

2. Python Setup

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, KNN
from spatialstats.regression import OLS, SpatialLag, SpatialError, SpatialDurbin, compare_models, lr_test

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")
print(f"spatialstats : version {spatialstats.__version__}")

spatialstats : version 0.2.0

---

3. The Code Skeleton

Everything from Chapters 3.2–3.5, condensed to the essential calls:

from spatialstats.weights import Queen
from spatialstats.regression import OLS, SpatialLag, SpatialError, SpatialDurbin, compare_models, lr_test

covariates = ["Obesity", "PhysInact", "Acc_Exer", "PCP_rate", "FoodEnvIdx", "SVI"]
X = counties[covariates].copy()
X["Urban"] = (counties["UrbRural"] == "Urban").astype(int)   # code categoricals before fitting
y = counties["Dia_pct"]

w = Queen(counties)                                # row-standardised by default
ols = OLS(y, X, w=w)
print(ols.lm_tests().recommendation())              # 3.3: which spatial model does the data ask for?

slm, sem = SpatialLag(y, X, w), SpatialError(y, X, w)
sdm = SpatialDurbin(y, X, w)
print(compare_models({"OLS": ols, "SLM": slm, "SEM": sem, "SDM": sdm}))   # 3.5: side-by-side

print(slm.impacts(seed=1))                          # 3.4: direct / indirect / total, with intervals
print(lr_test(ols, sem))                             # 3.5: formal nested-model test

Spatial error model (both robust tests significant; larger statistic chosen — consider SARMA/Durbin)

       n_params      loglik          aic          bic  pseudo_r2       rho  \
model                                                                        
OLS           8 -784.201471  1584.402942  1620.413259   0.747023       NaN   
SLM           9 -767.791416  1553.582832  1594.094439   0.760679  0.178371   
SEM           9 -761.826339  1541.652677  1582.164285   0.745435       NaN   
SDM          16 -753.108919  1538.217838  1610.238472   0.774397  0.320848   

         lambda  resid_moran_I  resid_moran_p  
model                                          
OLS         NaN       0.181223   2.512732e-14  
SLM         NaN       0.075346   1.769375e-03  
SEM    0.352145      -0.021580   4.139998e-01  
SDM         NaN      -0.020798   4.324358e-01  

     variable    direct  direct_lo  direct_hi  indirect  indirect_lo  \
0     Obesity  0.137632   0.112471   0.157394  0.028821     0.017773   
1   PhysInact  0.142405   0.118459   0.171209  0.029821     0.018326   
2    Acc_Exer -0.003956  -0.007106  -0.000951 -0.000828    -0.001632   
3    PCP_rate -0.000483  -0.002235   0.001310 -0.000101    -0.000491   
4  FoodEnvIdx -0.182983  -0.262619  -0.102632 -0.038318    -0.062556   
5         SVI  0.977087   0.686752   1.261708  0.204608     0.115508   
6       Urban -0.015972  -0.157046   0.125970 -0.003345    -0.033701   

   indirect_hi     total  total_lo  total_hi  
0     0.040759  0.166452  0.137158  0.193788  
1     0.044504  0.172226  0.140512  0.209380  
2    -0.000192 -0.004784 -0.008449 -0.001163  
3     0.000284 -0.000584 -0.002771  0.001595  
4    -0.018552 -0.221301 -0.321443 -0.124028  
5     0.324449  1.181695  0.819194  1.536268  
6     0.027766 -0.019317 -0.192048  0.149162  
{'statistic': 44.750264218103894, 'df': 1, 'p_value': 2.2383673975788455e-11}

Running this cell reproduces every headline number from Chapters 3.2–3.5 in one place — a useful starting point to adapt for a new dataset or a new set of covariates.

---

4. Exercises

4.1 SEM vs. SLM: which stays close to OLS?

Theory (Chapter 3.4) predicts SEM coefficients should track OLS closely (spatial dependence is a nuisance, not a structural change), while SLM's should differ more systematically (part of the covariates' apparent effect is reattributed to $\rho Wy$). Check this claim coefficient by coefficient, not just by eyeballing the summary tables.

w = Queen(counties)
X_ex = counties[covariates].copy()
X_ex["Urban"] = (counties["UrbRural"] == "Urban").astype(int)
y_ex = counties["Dia_pct"]

ols = OLS(y_ex, X_ex, w=w)
slm = SpatialLag(y_ex, X_ex, w)
sem = SpatialError(y_ex, X_ex, w)

coefs = pd.DataFrame({"OLS": ols.beta, "SLM": slm.beta, "SEM": sem.beta}, index=ols.names).drop(index="const")
coefs["|SLM - OLS|"] = (coefs["SLM"] - coefs["OLS"]).abs()
coefs["|SEM - OLS|"] = (coefs["SEM"] - coefs["OLS"]).abs()
coefs.round(4)

	OLS	SLM	SEM	|SLM - OLS|	|SEM - OLS|
Obesity	0.1462	0.1368	0.1478	0.0095	0.0016
PhysInact	0.1520	0.1415	0.1501	0.0105	0.0019
Acc_Exer	-0.0060	-0.0039	-0.0047	0.0021	0.0013
PCP_rate	-0.0003	-0.0005	-0.0001	0.0001	0.0002
FoodEnvIdx	-0.2002	-0.1818	-0.1669	0.0184	0.0333
SVI	1.0327	0.9709	0.9464	0.0618	0.0863
Urban	-0.0144	-0.0159	-0.0657	0.0015	0.0513

Answer: the textbook expectation only holds for part of the coefficient vector. For the two strongest, most significant predictors — Obesity and PhysInact, both with $|t|>10$ in the OLS fit — SEM is dramatically closer to OLS than SLM is (about 6× closer: |SEM-OLS| $\approx 0.002$ vs. |SLM-OLS|$\approx 0.010$ for both). But for FoodEnvIdx, SVI, and especially Urban` (already the least significant covariate in every model), the pattern reverses — SEM moves further from OLS than SLM does. The honest conclusion: SEM's closeness to OLS is a tendency for well-identified coefficients, not a guarantee for every coefficient in the model, and a summary statistic averaged across all seven coefficients (mean $|SLM-OLS|=0.015$ vs. mean $|SEM-OLS|=0.025$) would actually get the overall comparison backwards here. Always check coefficient by coefficient rather than trusting an aggregate summary or a single textbook rule of thumb.

4.2 LM tests vs. AIC: do they agree between SLM and SEM?

Chapter 3.3's LM tests recommended the spatial error model. Does comparing SLM and SEM by AIC alone reach the same conclusion?

lm = ols.lm_tests(w)
print("LM recommendation:", lm.recommendation())
print(f"\nAIC:  SLM={slm.aic:.2f}   SEM={sem.aic:.2f}   (lower is better)")
print(f"BIC:  SLM={slm.bic:.2f}   SEM={sem.bic:.2f}")
print(f"\nAIC prefers: {'SEM' if sem.aic < slm.aic else 'SLM'}")

LM recommendation: Spatial error model (both robust tests significant; larger statistic chosen — consider SARMA/Durbin)

AIC:  SLM=1553.58   SEM=1541.65   (lower is better)
BIC:  SLM=1594.09   SEM=1582.16

AIC prefers: SEM

Answer: yes, they agree. SEM's AIC (1541.65) is lower than SLM's (1553.58), and its BIC is lower too (1582.16 vs. 1594.09) — both criteria independently confirm the LM battery's recommendation. This is the expected, reassuring case (contrast with Chapter 3.5, where AIC and BIC disagreed about SDM specifically — Section 4.5 below explains why that disagreement happens there and not here).

4.3 The impacts identity under a non-row-standardised W

Chapter 3.4 showed the exact identity total impact $=\beta/(1-\rho)$ holds for row-standardised weights. Does it still hold for binary (non-row-standardised) weights?

w_binary = Queen(counties, transform="binary")
slm_binary = SpatialLag(y_ex, X_ex, w_binary)
impacts_binary = slm_binary.impacts(seed=1)

beta_binary = pd.Series(slm_binary.beta, index=slm_binary.names)
check = impacts_binary.assign(
    beta=lambda d: d["variable"].map(beta_binary),
    naive=lambda d: d["variable"].map(beta_binary) / (1 - slm_binary.rho),
)
print(f"rho (binary W) = {slm_binary.rho:.5f}")
check[["variable", "total", "beta", "naive"]].round(5)

rho (binary W) = 0.00463

	variable	total	beta	naive
0	Obesity	0.14581	0.14214	0.14280
1	PhysInact	0.15479	0.15089	0.15159
2	Acc_Exer	-0.00516	-0.00503	-0.00506
3	PCP_rate	-0.00031	-0.00030	-0.00030
4	FoodEnvIdx	-0.20785	-0.20261	-0.20355
5	SVI	1.06058	1.03386	1.03868
6	Urban	-0.01377	-0.01342	-0.01349

Answer: no, the simple scalar identity breaks down. The total column (computed via the general trace-based impacts() formula) no longer matches beta / (1 - rho) once $W$ is not row-standardised — every row shown differs, though only slightly here because $\hat\rho$ happens to be very small (0.0046) for this particular binary-$W$ fit. The reason: the shortcut $\text{total}=\beta/(1-\rho)$ relies on the constant vector $\mathbf{1}$ being an eigenvector of $W$ with eigenvalue 1 — true by construction for row-standardised weights (every row sums to exactly 1), but not for binary weights (row sums equal each county's raw neighbour count, which varies). impacts()'s general formula — built from $\text{tr}(A^{-1})$ and $\mathbf{1}^\top A^{-1} \mathbf{1}$ rather than the scalar shortcut — remains correct regardless; only the simplified identity is transform-specific.

4.4 Refitting with KNN(6): do conclusions change?

Chapter 3.3 and 3.5 already checked this piecemeal. Consolidate: under KNN(6) weights, does the model recommendation, the residual-Moran acceptance test, and the AIC/BIC pattern all still hold?

w_knn = KNN(counties, k=6)
ols_knn = OLS(y_ex, X_ex, w=w_knn)
slm_knn = SpatialLag(y_ex, X_ex, w_knn)
sem_knn = SpatialError(y_ex, X_ex, w_knn)
sdm_knn = SpatialDurbin(y_ex, X_ex, w_knn)

print("LM recommendation (KNN(6)):", ols_knn.lm_tests(w_knn).recommendation())
comp_knn = compare_models({"OLS": ols_knn, "SLM": slm_knn, "SEM": sem_knn, "SDM": sdm_knn}, w=w_knn)
comp_knn[["aic", "bic", "resid_moran_I", "resid_moran_p"]].round(4)

LM recommendation (KNN(6)): Spatial error model (both robust tests significant; larger statistic chosen — consider SARMA/Durbin)

	aic	bic	resid_moran_I	resid_moran_p
model
OLS	1584.4029	1620.4133	0.1751	0.0000
SLM	1547.2836	1587.7952	0.0791	0.0002
SEM	1531.9057	1572.4173	-0.0136	0.5728
SDM	1524.6155	1596.6362	-0.0111	0.6548

Answer: every conclusion is unchanged. LM tests still recommend the spatial error model; AIC still prefers SDM (1524.6) and BIC still prefers SEM (1572.4); SEM and SDM residuals still pass the near-zero acceptance test (SEM: $I=-0.014$, $p=0.57$) while OLS and SLM do not. The specific numbers move somewhat (as they did throughout Part 3 whenever the weights definition changed) but the qualitative story built across Chapters 3.2–3.5 is robust to this particular alternative.

4.5 Why AIC and BIC disagree about the SDM

Chapter 3.5 found AIC prefers SDM while BIC prefers SEM. Explain the mechanism precisely, using the actual parameter counts and log-likelihoods.

sdm = SpatialDurbin(y_ex, X_ex, w)
n = sdm.n
delta_loglik = sdm.loglik - sem.loglik
delta_params = sdm.n_params - sem.n_params
print(f"n = {n}")
print(f"SEM: n_params={sem.n_params}, loglik={sem.loglik:.3f}")
print(f"SDM: n_params={sdm.n_params}, loglik={sdm.loglik:.3f}")
print(f"extra log-likelihood from SDM's {delta_params} extra parameters: {delta_loglik:.3f}")
print(f"\nAIC penalty per extra parameter : 2")
print(f"BIC penalty per extra parameter : ln(n) = {np.log(n):.3f}")
print(f"\nAIC: SDM wins by {2*(-delta_params) + 2*delta_loglik:.2f} (negative = SDM lower/better)" )
print(f"BIC: SEM wins by {np.log(n)*delta_params - 2*delta_loglik:.2f} (positive = SDM's extra BIC cost exceeds its loglik gain)")

n = 666
SEM: n_params=9, loglik=-761.826
SDM: n_params=16, loglik=-753.109
extra log-likelihood from SDM's 7 extra parameters: 8.717

AIC penalty per extra parameter : 2
BIC penalty per extra parameter : ln(n) = 6.501

AIC: SDM wins by 3.43 (negative = SDM lower/better)
BIC: SEM wins by 28.07 (positive = SDM's extra BIC cost exceeds its loglik gain)

Answer: SDM's 7 extra parameters buy a log-likelihood improvement of about 8.7 points ($-753.1$ vs. $-761.8$). AIC charges only 2 points per extra parameter ($2 \times 7 = 14$), so a gain of $2\times8.7\approx17.4$ in likelihood terms comfortably clears that charge, and AIC prefers SDM. BIC charges $\ln(n) = \ln(666) \approx 6.5$ points per extra parameter ($6.5\times7\approx45.6$) — far more than the ~17.4-point likelihood gain can justify, so BIC prefers the more parsimonious SEM. The disagreement is not a contradiction or a sign that either criterion is wrong; it is a direct, computable consequence of how much more heavily BIC penalises additional parameters as $n$ grows ($\ln(n)$ versus AIC's constant 2), and it recurs in exactly the same form any time a richer model's likelihood gain is real but modest relative to its added parameter count.

---

5. Summary and Conclusions

This chapter closed Part 3 with a reusable code skeleton and five exercises that stress-tested specific claims from Chapters 3.2–3.5.

What we did

1. Condensed the full Chapters 3.2–3.5 workflow into a compact, reusable skeleton.
2. Exercise 1 found the "SEM close to OLS" rule of thumb holds strongly for the two best-identified coefficients but reverses for three weaker ones — a genuine complication, not a clean confirmation.
3. Exercise 2 confirmed LM tests and AIC/BIC agree on SLM vs. SEM (unlike the AIC/BIC disagreement about SDM specifically).
4. Exercise 3 showed the impacts identity total$=\beta/(1-\rho)$ is specific to row-standardised weights and breaks down, in general, for binary weights.
5. Exercise 4 consolidated every KNN(6) robustness check from earlier chapters: every qualitative conclusion in Part 3 survives the alternative weights definition.
6. Exercise 5 derived, with actual numbers, exactly why AIC prefers SDM while BIC prefers SEM — a direct consequence of BIC's $\ln(n)$ per-parameter penalty versus AIC's constant 2.

Practical takeaways

• Textbook rules of thumb (like "SEM tracks OLS") are tendencies, not guarantees — verify them coefficient by coefficient on your own data before repeating them as a finding.
• A convenient closed-form identity (impacts under row-standardisation) is often a special case; check whether the general formula is being used before assuming a shortcut applies.
• AIC vs. BIC disagreements are fully explainable from the parameter count and sample size alone — compute the breakeven point rather than treating the disagreement as ambiguous.

Part 3 in review

Across seven chapters, Part 3 took a single diagnosed problem (Part 2's confirmed spatial autocorrelation in diabetes prevalence) through a complete, validated workflow: diagnose the mechanism (Chapter 3.3), fit the recommended model and understand what its parameters mean (Chapter 3.4), confirm the fit actually solved the problem it was chosen for (Chapter 3.5), and stress-test the one assumption every model in the Part still shares (Chapter 3.6). Part 4 (Spatial Cluster Analysis) returns to the where question with a fuller toolkit than Part 2's LISA; Part 6 (Geographically Weighted Models) is the direct continuation of Chapter 3.6's finding that a single global $\beta$ is not always enough.

---

6. Resources

Textbooks

Reference	Notes
Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer Academic.	The reference underlying every method used across Part 3
LeSage, J. & Pace, R. K. (2009). Introduction to Spatial Econometrics. CRC Press.	Impacts, model comparison, and practical estimation guidance

Key papers

Paper	Contribution
Burnham, K. P. & Anderson, D. R. (2004). Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods & Research, 33, 261–304.	The AIC/BIC penalty comparison used in Exercise 4.5

In this library

Object	Role
spatialstats.regression.*	Every object used across Part 3, summarised in the Section 3 skeleton
Part 4	spatialstats.cluster — the fuller "where" toolkit continuing from Part 2's LISA
Part 6	spatialstats.gwmodel — the direct continuation of Chapter 3.6