3.5 Comparison, Reporting and Diagnostics¶
Part: 3 — Spatial Regression (Chapter 3.5 of 7)
Dataset: diabetes_atlantic.shp (666 counties)
Focus: Side-by-side model comparison, likelihood-ratio tests, the residual-Moran acceptance criterion, and reporting standards
Library: spatialstats.regression (PySpatialStats)
1. Overview¶
Chapter 3.4 fitted five models on the same data. This chapter puts them side by side, formally, and closes the loop back to Chapter 3.3's recommendation: did fitting the spatial error model actually solve the problem it was recommended for?
| Step | What we do |
|---|---|
| Comparison table | compare_models — log-likelihood, AIC, BIC, pseudo-R², spatial parameter, residual Moran, all in one table |
| LR tests | Formal nested-model comparisons: OLS vs. SLM, OLS vs. SEM, SLM vs. SDM |
| Pseudo-R² caveat | Why it is not comparable to OLS's ordinary R², demonstrated directly |
| Acceptance test | The standard the whole exercise is judged against: residual Moran's I should collapse toward zero |
| Sensitivity | Repeat the comparison under a different W |
| Reporting | What a complete, reproducible spatial regression report states |
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, SLX, 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")
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)
ols = OLS(y, X, w=w)
slx = SLX(y, X, w=w)
slm = SpatialLag(y, X, w)
sem = SpatialError(y, X, w)
sdm = SpatialDurbin(y, X, w)
print(f"spatialstats : version {spatialstats.__version__}")
print("re-fit OLS, SLX, SLM, SEM, SDM from Chapters 3.2-3.4")
spatialstats : version 0.2.0 re-fit OLS, SLX, SLM, SEM, SDM from Chapters 3.2-3.4
3. Comparing Every Model Side by Side¶
compare_models assembles every fitted model's log-likelihood, information criteria, pseudo-R², spatial
parameter, and residual Moran's I into a single table — likelihood-based columns are only meaningfully
comparable because every model here shares the same $y$, $X$, and $W$ (Chapter 3.2, Section 4).
comparison = compare_models({"OLS": ols, "SLX": slx, "SLM": slm, "SEM": sem, "SDM": sdm}, w=w)
comparison.round(4)
| n_params | loglik | aic | bic | pseudo_r2 | rho | lambda | resid_moran_I | resid_moran_p | |
|---|---|---|---|---|---|---|---|---|---|
| model | |||||||||
| OLS | 8 | -784.2015 | 1584.4029 | 1620.4133 | 0.7470 | NaN | NaN | 0.1812 | 0.0000 |
| SLX | 15 | -772.8100 | 1575.6200 | 1643.1394 | 0.7555 | NaN | NaN | 0.1699 | 0.0000 |
| SLM | 9 | -767.7914 | 1553.5828 | 1594.0944 | 0.7607 | 0.1784 | NaN | 0.0753 | 0.0018 |
| SEM | 9 | -761.8263 | 1541.6527 | 1582.1643 | 0.7454 | NaN | 0.3521 | -0.0216 | 0.4140 |
| SDM | 16 | -753.1089 | 1538.2178 | 1610.2385 | 0.7744 | 0.3208 | NaN | -0.0208 | 0.4324 |
Every spatial model beats OLS on log-likelihood and AIC. AIC prefers SDM (1538.2, the lowest) but BIC prefers SEM (1582.2, versus SDM's 1610.2) — BIC's heavier penalty per parameter outweighs SDM's extra log-likelihood once its 7 additional $WX$ terms are counted. This is not a contradiction to resolve; it is the two criteria trading off fit against parsimony differently, and it is worth reporting both rather than picking whichever supports a preferred conclusion.
4. Likelihood-Ratio Tests for Nested Pairs¶
lr_test(restricted, full) formally compares two nested models — the simpler one recoverable from the more
complex one by fixing a parameter to zero.
print("OLS vs SLM (restricting rho=0):", lr_test(ols, slm))
print("OLS vs SEM (restricting lambda=0):", lr_test(ols, sem))
print("SLM vs SDM (restricting theta=0):", lr_test(slm, sdm))
OLS vs SLM (restricting rho=0): {'statistic': 32.82010968464101, 'df': 1, 'p_value': 1.0109360594916218e-08}
OLS vs SEM (restricting lambda=0): {'statistic': 44.750264218103894, 'df': 1, 'p_value': 2.2383673975788455e-11}
SLM vs SDM (restricting theta=0): {'statistic': 29.364994482358952, 'df': 7, 'p_value': 0.00012410813694179224}
All three reject their restriction decisively. lr_test only accepts pairs where the second model has more
parameters than the first (ValueError otherwise) — SLM and SEM are not nested in each other (neither's
parameter set is a subset of the other's), so lr_test(slm, sem) is not a meaningful comparison; AIC/BIC
(Section 3) are the right tool for comparing non-nested models like these two.
5. The Pseudo-R² Caveat¶
pseudo_r2 (squared correlation of $y$ with the model's prediction) is defined identically for every model, so
it is comparable across model types — but it is not the same quantity as OLS's ordinary $R^2$, even
though the two coincide exactly for OLS itself (Chapter 3.2).
print(f"OLS pseudo_r2 = {ols.pseudo_r2:.4f} (== ols.r2 = {ols.r2:.4f})")
print(f"SEM pseudo_r2 = {sem.pseudo_r2:.4f}")
print(f"SLM pseudo_r2 = {slm.pseudo_r2:.4f}")
OLS pseudo_r2 = 0.7470 (== ols.r2 = 0.7470) SEM pseudo_r2 = 0.7454 SLM pseudo_r2 = 0.7607
SEM's pseudo-R² (0.7454) is actually lower than OLS's (0.7470) despite SEM having a substantially better log-likelihood — the two metrics answer different questions. Log-likelihood (and AIC/BIC) measure fit to the full probability model, including the spatial error structure; pseudo-R² measures only how well the point predictions $\hat y$ track $y$, which for SEM uses $X\beta$ alone (not the fitted spatial error term). Do not use pseudo-R² to argue a spatial model "fits better" than OLS — use the likelihood-based criteria for that, and treat pseudo-R² only as a rough, comparable-across-model-types summary of point-prediction accuracy.
6. The Acceptance Test: Residual Moran's I Near Zero¶
The entire exercise, from Chapter 3.2's diagnosis through Chapter 3.4's fits, is judged against one concrete question: did the chosen model's residuals stop being spatially autocorrelated?
acceptance = comparison[["resid_moran_I", "resid_moran_p"]].round(4)
acceptance
| resid_moran_I | resid_moran_p | |
|---|---|---|
| model | ||
| OLS | 0.1812 | 0.0000 |
| SLX | 0.1699 | 0.0000 |
| SLM | 0.0753 | 0.0018 |
| SEM | -0.0216 | 0.4140 |
| SDM | -0.0208 | 0.4324 |
OLS (0.181) and SLX (0.170) both leave clear, significant residual autocorrelation — neither actually addresses the problem, despite SLX adding seven extra covariates. SLM reduces it substantially but not to zero (0.075, still significant at $p=0.002$) — consistent with Chapter 3.3's diagnosis that error, not lag, was the dominant mechanism here. SEM (−0.022, $p=0.414$) and SDM (−0.021, $p=0.432$) both pass: residual autocorrelation is statistically indistinguishable from zero. This is the formal confirmation that Chapter 3.3's recommendation (spatial error) was correct, and that SDM — which nests a lag term alongside the $WX$ covariates — also succeeds by capturing enough of the same structure a different way.
7. Sensitivity to the Weights Definition¶
Repeating the full comparison under KNN(6) weights (Chapter 3.3's robustness check, extended to every model):
w_knn = KNN(counties, k=6)
ols_k = OLS(y, X, w=w_knn)
slm_k = SpatialLag(y, X, w_knn)
sem_k = SpatialError(y, X, w_knn)
sdm_k = SpatialDurbin(y, X, w_knn)
comparison_knn = compare_models({"OLS": ols_k, "SLM": slm_k, "SEM": sem_k, "SDM": sdm_k}, w=w_knn)
comparison_knn[["aic", "bic", "resid_moran_I", "resid_moran_p"]].round(4)
| 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 |
The same qualitative pattern holds: AIC still prefers SDM, BIC still prefers SEM, and SEM/SDM residuals are still the ones that pass the near-zero acceptance test — the conclusion is not an artefact of the Queen weights choice specifically.
8. Reporting Standards¶
A reproducible spatial regression report states, at minimum:
- The weights definition (type, transform, $k$/threshold if applicable) and whether it changes the conclusion under a reasonable alternative (Section 7).
- The estimator and its assumptions (ML here; note the Gaussian error assumption).
- Full diagnostics leading to the model choice: the LM battery and its recommendation (Chapter 3.3), not just the final model in isolation.
- Impacts, not raw coefficients, for any lag-containing model (Chapter 3.4) — reporting SLM/SDM $\beta$ alone without direct/indirect/total effects misstates the substantive finding.
- The acceptance check (Section 6): did the chosen model's residual Moran's I actually reach non-significance?
- AIC and BIC together, and an explicit note when they disagree (Section 3), rather than silently reporting only the criterion that favours a preferred model.
A summary table combining Sections 3 and 6 is usually the single most useful artefact to include:
report_table = comparison[["loglik", "aic", "bic", "pseudo_r2", "rho", "lambda", "resid_moran_I"]].round(3)
report_table
| loglik | aic | bic | pseudo_r2 | rho | lambda | resid_moran_I | |
|---|---|---|---|---|---|---|---|
| model | |||||||
| OLS | -784.201 | 1584.403 | 1620.413 | 0.747 | NaN | NaN | 0.181 |
| SLX | -772.810 | 1575.620 | 1643.139 | 0.756 | NaN | NaN | 0.170 |
| SLM | -767.791 | 1553.583 | 1594.094 | 0.761 | 0.178 | NaN | 0.075 |
| SEM | -761.826 | 1541.653 | 1582.164 | 0.745 | NaN | 0.352 | -0.022 |
| SDM | -753.109 | 1538.218 | 1610.238 | 0.774 | 0.321 | NaN | -0.021 |
9. A Note on Prediction¶
Every model in this chapter is evaluated in-sample only: predy (SLM/SEM/SDM) uses the observed $y$ of
neighbouring counties (through $Wy$ or the filtered residual structure) to help predict each county's own
value. This library does not currently support prediction at a genuinely new, out-of-sample location — doing so
correctly requires deciding how to obtain $Wy$ for a location with no observed neighbours in the original
sample, which is a nontrivial extension beyond the scope of what is implemented here. Every number in this
Part, including every model comparison, should be read as a statement about fit to the observed 666
counties, not as a validated predictive model for a 667th.
10. Summary and Conclusions¶
This chapter closed the loop from Chapter 3.3's diagnosis to a validated, comparable set of fitted models.
What we did¶
- Built the full model comparison table: AIC prefers SDM (1538.2), BIC prefers SEM (1582.2 vs. SDM's 1610.2).
- Ran likelihood-ratio tests for the three genuinely nested pairs (OLS vs. SLM, OLS vs. SEM, SLM vs. SDM), all strongly significant; flagged that SLM and SEM are not nested and so not LR-comparable.
- Demonstrated the pseudo-R² caveat directly: SEM's pseudo-R² (0.745) is lower than OLS's (0.747) despite a far better log-likelihood — the two metrics are not interchangeable.
- Applied the acceptance test: OLS and SLX both fail (residual Moran's I still significant); SLM improves but does not pass; SEM and SDM both pass (I ≈ −0.02, not significant) — the formal confirmation of Chapter 3.3's recommendation.
- Confirmed the whole pattern is robust to a KNN(6) alternative weights definition.
- Laid out six concrete reporting standards for a reproducible spatial regression analysis.
- Flagged that every result in this Part is in-sample only; out-of-sample prediction is not currently supported.
Practical takeaways¶
- Report AIC and BIC, and say explicitly when they disagree — do not cherry-pick the one that favours a preferred model.
- The residual-Moran acceptance test is the real finish line, not the LM recommendation from Chapter 3.3 — a model can be "recommended" and still fail this check if the recommendation was marginal.
- Never compare pseudo-R² across model types as if it were ordinary $R^2$; use the likelihood-based criteria for fit comparisons instead.
Next steps¶
Chapter 3.6 (Limits and Connections) steps back from model selection to ask what OLS, SLM, and SEM all still share: a single, region-wide $\beta$ — and previews Part 6's geographically weighted models as the response.
11. Resources¶
Textbooks¶
| Reference | Notes |
|---|---|
| Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer Academic. | Model comparison and diagnostic standards for spatial regression |
| Burnham, K. P. & Anderson, D. R. (2002). Model Selection and Multimodel Inference (2nd ed.). Springer. | General AIC/BIC theory and the trade-off explored in Section 3 |
Key papers¶
| Paper | Contribution |
|---|---|
| Schwarz, G. (1978). Estimating the dimension of a model. Annals of Statistics, 6, 461–464. | Origin of BIC, whose heavier parameter penalty explains Section 3's AIC/BIC disagreement |
| Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19, 716–723. | Origin of AIC |
In this library¶
| Object | Role |
|---|---|
spatialstats.regression.compare_models |
Side-by-side comparison table used throughout this chapter |
spatialstats.regression.lr_test |
Formal nested-model comparison |
| Chapter 3.4 | The five fitted models compared here |
| Chapter 3.6 | What every model in this chapter still cannot do: vary its coefficients over space |