Spatial Regression Models: Lag, Error, and GWR Workflows in Python

Spatial regression models extend ordinary least squares by explicitly parameterising the dependence structure that arises when geographic neighbours share unobserved factors or directly influence each other’s outcomes. They are the standard tool in spatial econometrics, environmental epidemiology, urban analytics, and real estate valuation — any domain where residual maps show pattern rather than noise. This page covers the full workflow within the broader context of Python Workflows for Spatial Modeling & Regression: from OLS diagnostics through GMM and ML estimation to production deployment.

Prerequisites #

• Python 3.9+
• numpy>=1.22, scipy>=1.9, geopandas>=0.13, libpysal>=4.8, spreg>=1.4, esda>=2.4, mgwr>=2.2
• A GeoDataFrame with valid, non-overlapping geometries reprojected to a metric CRS (e.g., EPSG:32633 for Europe, EPSG:26918 for eastern North America). Geographic coordinates (degrees) distort distance calculations and invalidate neighbourhood construction.
• A validated spatial weight matrix with no isolated observations (W.islands == []) and a single connected component (W.n_components == 1). Disconnected graphs prevent the matrix inversion $(I - \rho W)^{-1}$ required by lag and error estimators.
• Response vector y and design matrix X with row order matching the GeoDataFrame index after reset_index(drop=True). Row misalignment is the most common source of silent coefficient corruption.

Mathematical Core #

Why OLS Fails: Two Spatial Dependence Processes #

Tobler’s First Law — nearby observations resemble each other more than distant ones — manifests in regression residuals when spatial structure is unmodelled. Two distinct data-generating processes require different corrections.

Spatial Lag Model (SLM). The outcome at location $i$ is directly influenced by neighbours’ outcomes through a feedback process:

$$\mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \quad \boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I})$$

where $\rho$ is the spatial autoregressive coefficient and $\mathbf{W}$ is the row-standardised weight matrix. Ignoring $\rho \mathbf{W}\mathbf{y}$ creates an omitted-variable that is correlated with $\mathbf{X}$ whenever $\mathbf{X}$ is spatially clustered — producing biased, inconsistent OLS estimates.

Spatial Error Model (SEM). Spatial dependence enters through an unobserved, spatially structured disturbance:

$$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{u}, \quad \boldsymbol{u} = \lambda \mathbf{W}\boldsymbol{u} + \boldsymbol{\varepsilon}, \quad \boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I})$$

where $\lambda$ is the spatial error autocorrelation coefficient. OLS remains unbiased here but is inefficient — conventional standard errors are incorrect, leading to invalid $t$-statistics and confidence intervals.

Spatial Durbin Model (SDM). Combines spatially lagged outcomes and spatially lagged covariates:

$$\mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \mathbf{W}\mathbf{X}\boldsymbol{\theta} + \boldsymbol{\varepsilon}$$

SDM nests both SLM and SEM as special cases and is appropriate when theory motivates spillover effects in both the outcome and predictors.

Reduced Form and Effect Decomposition #

For the SLM, solving for $\mathbf{y}$ gives the reduced form:

$$\mathbf{y} = (I - \rho W)^{-1} \mathbf{X}\boldsymbol{\beta} + (I - \rho W)^{-1} \boldsymbol{\varepsilon}$$

The matrix $(I - \rho W)^{-1}$ encodes feedback paths of all orders through the network. A one-unit change in covariate $k$ at location $i$ generates a direct effect (impact on $y_i$) and an indirect spillover (impact on all $y_j, j \neq i$) whose sum is the total effect. These three quantities replace the scalar $\beta_k$ from OLS interpretation.

Annotated Implementation #

Step 1 — Build a Production-Ready Weights Matrix #

import geopandas as gpd
import libpysal
import numpy as np

# Load, repair, and project geometry
gdf = gpd.read_file("data/parcels.gpkg").to_crs("EPSG:32633")
if not gdf.geometry.is_valid.all():
    gdf["geometry"] = gdf.geometry.make_valid()

# libpysal requires a 0-based sequential integer index
gdf = gdf.reset_index(drop=True)

# K-nearest neighbours (k=5): avoids islands in sparse or irregular point clouds
w = libpysal.weights.KNN.from_dataframe(gdf, k=5)

# Row-standardise: each row sums to 1.0, Wy becomes a neighbourhood average
w.transform = "r"

# Structural checks — must pass before any model fitting
assert len(w.islands) == 0, f"Islands found: {w.islands}"
assert w.n_components == 1, f"Disconnected graph: {w.n_components} components"

Queen or Rook contiguity is appropriate for polygon areal data; KNN is the safer default for point observations because it guarantees every observation has exactly $k$ neighbours. Full topology selection guidance is in the spatial weight matrices reference.

Step 2 — OLS Baseline with Spatial Diagnostics #

import spreg

y = gdf["log_price"].values.reshape(-1, 1)
X = gdf[["sqft", "age", "dist_cbd", "crime_rate"]].values
feature_names = ["sqft", "age", "dist_cbd", "crime_rate"]

# OLS with Lagrange Multiplier tests and Moran's I on residuals
ols = spreg.OLS(
    y, X, w,
    spat_diag=True,         # computes LM-Lag, LM-Error, Robust variants
    moran=True,             # Moran's I on OLS residuals
    name_y="log_price",
    name_x=feature_names
)
print(ols.summary)

# Extract decision-critical statistics
lm_lag_stat,  lm_lag_p  = ols.lm_lag
lm_err_stat,  lm_err_p  = ols.lm_error
rlm_lag_stat, rlm_lag_p = ols.lm_lag_robust
rlm_err_stat, rlm_err_p = ols.lm_error_robust

print(f"LM-Lag:        {lm_lag_stat:.3f}  p={lm_lag_p:.4f}")
print(f"LM-Error:      {lm_err_stat:.3f}  p={lm_err_p:.4f}")
print(f"Robust LM-Lag: {rlm_lag_stat:.3f}  p={rlm_lag_p:.4f}")
print(f"Robust LM-Err: {rlm_err_stat:.3f}  p={rlm_err_p:.4f}")

The Robust LM tests are the decision-relevant statistics. Standard LM-Lag and LM-Error are biased whenever both processes coexist — the robust variants partial out the contamination from the alternative process. Never build a model selection decision on the non-robust statistics alone.

Step 3 — Fit the Spatial Lag Model (GMM) #

# Generalized Method of Moments: robust to heteroskedasticity, scales to large n
slm = spreg.GM_Lag(
    y, X, w,
    robust="hac",           # HAC standard errors for spatially correlated residuals
    name_y="log_price",
    name_x=feature_names
)
print(slm.summary)

# rho: spatial autoregressive coefficient (range: -1 to 1)
print(f"rho = {slm.rho:.4f}  (p={slm.z_stat[-1][1]:.4f})")

GM_Lag instruments the endogenous spatial lag $\mathbf{W}\mathbf{y}$ using spatially lagged exogenous variables. Setting robust="hac" adjusts standard errors for residual spatial autocorrelation — recommended when ols.moran_res is significant. For samples under approximately 5,000 observations with approximately normal residuals, spreg.ML_Lag provides asymptotically efficient estimates.

Step 4 — Fit the Spatial Error Model (GMM) #

# Spatial Error Model: bias-corrects OLS when spatial dependence is in the disturbance
sem = spreg.GM_Error(
    y, X, w,
    name_y="log_price",
    name_x=feature_names
)
print(sem.summary)

# lambda: spatial error autocovariance (range: -1 to 1)
print(f"lambda = {sem.e_lam:.4f}")

# Model comparison via AIC (lower is better)
print(f"SLM AIC: {slm.aic:.2f}  |  SEM AIC: {sem.aic:.2f}")

When $\lambda$ is large and significant, it confirms that important spatially structured predictors are missing from X. The SEM coefficient vector $\boldsymbol{\beta}$ remains unbiased and interpretable in the standard OLS sense — unlike SLM, there are no spillover corrections needed.

Step 5 — Fit the Spatial Durbin Model #

# SDM: spatially lagged y AND spatially lagged X covariates
# w_lags=2 includes both Wy and WX (the Durbin terms)
sdm = spreg.GM_Lag(
    y, X, w,
    w_lags=2,
    name_y="log_price",
    name_x=feature_names
)
print(sdm.summary)

Use the SDM when theory motivates spatial spillovers in the covariates themselves — for instance, when a neighbourhood’s average income (W @ income) independently affects local property values beyond a unit’s own income.

Diagnostic Configuration and Model Selection Rationale #

The Anselin-Florax decision sequence, implemented in the diagram above, is the standard empirical guide:

Beyond the LM sequence, always inspect:

• Breusch-Pagan statistic (ols.bp_stat): significant heteroskedasticity justifies robust="hac" or robust="white" in GMM models.
• Moran’s I on model residuals (esda.Moran(slm.u, w).I): if residual autocorrelation persists after fitting SLM or SEM, the model is misspecified — consider additional covariates or SDM.
• Condition number (ols.std_err and np.linalg.cond(X)): near-collinear predictors inflate coefficient variance and obscure spatial parameter estimates.

Output Interpretation #

Reading Spatial Lag Model Results #

For the SLM, $\hat{\boldsymbol{\beta}}$ coefficients alone do not represent marginal effects. The correct quantities are:

$$\text{Direct} = \frac{1}{n}\text{tr}\left[(I - \hat{\rho}W)^{-1}\right] \hat{\beta}_k$$

$$\text{Total} = \frac{1}{n}\mathbf{1}^\top (I - \hat{\rho}W)^{-1} \mathbf{1} \hat{\beta}_k$$

$$\text{Indirect} = \text{Total} - \text{Direct}$$

import numpy as np
from scipy.sparse.linalg import inv as sparse_inv
from scipy.sparse import eye

# Effect decomposition for the spatial lag model
rho_hat = slm.rho
S = sparse_inv(eye(len(y), format="csc") - rho_hat * w.sparse.tocsc())

# For covariate k (0-indexed), beta_k is slm.betas[k, 0]
# (betas order: intercept, then feature_names order, then rho last)
for k, name in enumerate(feature_names):
    beta_k = slm.betas[k + 1, 0]   # +1 skips intercept
    S_dense_diag = S.diagonal()
    direct_k   = float(S_dense_diag.mean() * beta_k)
    total_k    = float(S.sum() / len(y) * beta_k)
    indirect_k = total_k - direct_k
    print(f"{name:15s}  direct={direct_k:+.4f}  indirect={indirect_k:+.4f}  total={total_k:+.4f}")

A positive indirect effect means that an increase in covariate $k$ anywhere in the network raises $y$ at all other locations — the classic spatial spillover. This is the mechanism that housing price changes in one neighbourhood diffuse into adjacent neighbourhoods.

Reading Spatial Error Model Results #

SEM coefficients are standard regression slopes in the filtered residual space. Interpretation is identical to OLS: $\hat{\beta}_k$ is the expected change in $y$ per unit increase in $x_k$, holding other covariates fixed, with correct standard errors. The $\lambda$ parameter is a nuisance quantity that cleans the error structure; it does not affect the substantive interpretation of $\boldsymbol{\beta}$.

Residual Map as a Diagnostic #

After fitting, always visualise the model residuals on a map:

import matplotlib.pyplot as plt

gdf["residuals"] = slm.u.flatten()

fig, ax = plt.subplots(figsize=(10, 8))
gdf.plot(
    column="residuals",
    cmap="RdBu_r",
    legend=True,
    legend_kwds={"label": "Residual (log-price units)"},
    ax=ax
)
ax.set_title("Spatial Lag Model Residuals")
plt.tight_layout()
plt.savefig("outputs/slm_residuals.png", dpi=150)

A residual map with remaining spatial clusters — contiguous patches of dark red or dark blue — signals that the model has not captured all spatial structure. This is a warning sign to: add spatially structured covariates, switch to SDM, or investigate geographically varying coefficients via the GWR extension below.

Geographically Weighted Regression for Non-Stationary Coefficients #

When diagnostic maps reveal that the relationship between $X$ and $y$ varies systematically across space — a different $\beta_k$ in the north versus the south — a global SLM or SEM applies a single parameter vector that masks this heterogeneity. Geographically Weighted Regression (GWR) relaxes this by estimating a local coefficient surface $\hat{\boldsymbol{\beta}}(u_i, v_i)$ at each location $(u_i, v_i)$:

from mgwr.gwr import GWR
from mgwr.sel_bw import Sel_BW
import numpy as np

# mgwr requires coordinate arrays, not GeoDataFrame
coords = list(zip(gdf.geometry.centroid.x, gdf.geometry.centroid.y))
y_gwr = gdf["log_price"].values.reshape(-1, 1)
X_gwr = gdf[["sqft", "age", "dist_cbd", "crime_rate"]].values

# Bandwidth selection via cross-validation (AICc criterion)
selector = Sel_BW(coords, y_gwr, X_gwr)
bw = selector.search(bw_min=20, bw_max=500)
print(f"Optimal bandwidth: {bw:.0f} neighbours")

# Fit GWR
gwr_model = GWR(coords, y_gwr, X_gwr, bw=bw)
gwr_results = gwr_model.fit()

# Local coefficient surface: shape (n, n_covariates + 1)
local_betas = gwr_results.params
print(f"Local beta range for sqft: [{local_betas[:, 1].min():.4f}, {local_betas[:, 1].max():.4f}]")

The bandwidth controls the degree of spatial smoothing: a smaller bandwidth produces highly localised estimates that capture fine-grained heterogeneity but may overfit; a larger bandwidth approaches the global OLS solution. AICc-based selection balances fit against model complexity.

GWR produces a $t$-statistic surface for each coefficient — map the pseudo $t$-values to identify regions where a covariate is locally statistically significant. Stationarity analysis using mgwr.utils.compare_surfaces tests whether local variance is genuine heterogeneity or sampling noise; see stationarity and trend analysis for the conceptual framing.

Production Considerations #

Computational scaling. GMM estimators (GM_Lag, GM_Error) scale roughly $O(n^2)$ in matrix operations and handle datasets up to approximately 50,000 observations in RAM. ML estimators require an $O(n^2)$ or $O(n^3)$ log-determinant computation — use the ML_Lag Ord approximation only below approximately 5,000 observations. GWR is $O(n^2)$ per bandwidth evaluation; for $n > 10{,}000$, use mgwr.gwr.MGWR with parallelism or tile the study area.

Memory layout. Always pass W.sparse (CSR format) rather than the dense W.full() array. At $n = 20{,}000$, the dense matrix is 3.2 GB; the sparse equivalent with average $k=5$ neighbours is 1.6 MB.

Serialising the pipeline. The spatial weight matrix must be reconstructed identically during inference — it is not a model parameter but a data-derived artefact. Store the construction metadata alongside the model:

import joblib

artifact = {
    "model": slm,
    "feature_names": feature_names,
    "crs": "EPSG:32633",
    "weights_config": {"type": "knn", "k": 5, "transform": "r"}
}
joblib.dump(artifact, "models/slm_v1.joblib")

Spatial drift monitoring. At inference time, compute Moran’s I on incoming residuals. A sudden rise in spatial autocorrelation — tracked as a control chart metric — signals covariate or neighbourhood structure drift that warrants retraining.

Spatially honest validation. Standard random $k$-fold cross-validation leaks spatial information because training and test folds share contiguous neighbours. Always use spatially blocked folds or buffered holdouts. Methodologies for these resampling strategies are covered in detail in cross-validation strategies.

Troubleshooting #

FAQ #

When should I use Spatial Lag vs Spatial Error? Choose the Spatial Lag Model when the feedback between neighbouring outcomes is theoretically meaningful and estimating the spillover magnitude matters — housing prices, disease transmission, economic multipliers. Choose the Spatial Error Model when spatial dependence is a nuisance from unmeasured spatially structured factors and your goal is unbiased coefficient estimation. The Anselin-Florax Robust LM sequence provides the empirical guide; theory provides the confirmation.

Why can I not read spatial lag beta coefficients like OLS betas? The SLM embeds a network feedback: changing $x_k$ at one location propagates through $\mathbf{W}$, altering neighbours’ $y$, which feeds back to alter $y$ everywhere. The scalar $\hat{\beta}_k$ captures only the local impulse; the indirect and total effects computed via $(I - \hat{\rho}W)^{-1}$ capture the full network response. Ignoring this understates total impact whenever $\hat{\rho} > 0$.

GMM or ML — which to use? GMM (GM_Lag, GM_Error) is the production default: distribution-free, heteroskedasticity-robust, and computationally tractable for large $n$. ML (ML_Lag, ML_Error) is asymptotically more efficient under correct specification with normal errors, but breaks down for $n > 5{,}000$ due to the log-determinant cost and is sensitive to distributional misspecification.

How do I predict out-of-sample with a Spatial Lag Model? True out-of-sample prediction requires observed or imputed neighbour $y$ values to compute $\mathbf{W}\mathbf{y}$ for new locations. Use the reduced-form solution $(I - \hat{\rho}W)^{-1} \mathbf{X}\hat{\boldsymbol{\beta}}$ when the full prediction-set $W$ is available, or the truncated power-series approximation $\sum_{k=0}^{K} \hat{\rho}^k W^k X\hat{\beta}$ for large $n$. Reconstruct $W$ using the exact parameters stored in the deployment artifact.

Next Steps #

For a hands-on walkthrough of direct, indirect, and total effect calculation and prediction in production, see Implementing Spatial Lag Models in Python. Spatially honest performance evaluation — blocking strategies that prevent neighbourhood leakage — is detailed in cross-validation strategies. When datasets grow beyond memory, the chunked and parallelised patterns in memory-efficient processing apply directly to weight construction and GMM fitting.

Related

• Implementing Spatial Lag Models in Python — direct, indirect, and total effect decomposition with full worked examples
• Cross-Validation Strategies — spatial blocking and buffered holdout methods for unbiased model evaluation
• Spatial Weight Matrices — construction, row-standardisation, and diagnostic validation of the W matrix that all spatial regression estimators require
• Spatial Autocorrelation Metrics — Moran’s I, Geary’s C, and LISA statistics for quantifying and locating spatial dependence before and after model fitting
• Stationarity & Trend Analysis — testing whether global model assumptions hold or whether GWR-style local estimation is warranted

← Back to Python Workflows for Spatial Modeling & Regression