Implementing Spatial Lag Models in Python

TL;DR #

Use spreg.ML_Lag(y, X, w) for datasets up to ~2,500 observations, or spreg.GM_Lag(y, X, w) for larger ones. Always call w.transform = 'r' first to row-standardize the weights matrix. Run spreg.OLS(..., spat_diag=True) beforehand and confirm ols.lm_lag[1] < 0.05 before fitting — that one diagnostic check prevents silent OLS mis-specification.

Why this matters #

A spatial lag model (also called a simultaneous autoregressive, or SAR, model) adds a spatially lagged copy of the dependent variable — $\rho W y$ — to the right-hand side of the regression equation. This term captures the direct feedback between neighboring observations: pollution at one monitoring station influences readings at adjacent stations; house prices in one block spill over into the next. Ignoring this process causes OLS to absorb the spillover into the residuals, producing omitted-variable bias across the entire coefficient vector.

This page covers the exact implementation steps within the spatial regression models workflow. It sits inside the broader Python workflows for spatial modeling and regression pipeline, where the same GeoDataFrame and weights matrix you prepare here feed directly into cross-validation and model comparison tasks.

Environment and version pinning #

These exact versions are tested against the code below. Pin them in your requirements.txt or environment.yml to avoid silent API changes in spreg.

pip install "geopandas>=0.14.0" "libpysal>=4.8.0" "spreg>=1.5.0" \
            "numpy>=1.24.0" "scipy>=1.10.0" "pandas>=2.0.0"

import geopandas as gpd        # >= 0.14.0
import numpy as np             # >= 1.24.0
import pandas as pd            # >= 2.0.0
import libpysal                # >= 4.8.0
import spreg                   # >= 1.5.0

CRS requirement: all geometry must be in a projected CRS (e.g. EPSG:3857 or a local UTM zone) before building distance-based weights. Queen and Rook contiguity weights are topology-based and CRS-agnostic, but any neighbor definition that uses distance thresholds requires metric coordinates.

Import path: spreg dropped the legacy pysal.spreg namespace in v2.0. Always import directly from spreg.

Step-by-step implementation #

Step 1 — Load data and build the spatial weights matrix #

# Load a projected GeoDataFrame
gdf = gpd.read_file("your_spatial_data.gpkg")
assert gdf.crs.is_projected, "Reproject to a metric CRS before continuing"

# Drop rows with nulls in the model variables
model_vars = ["y", "X1", "X2"]
df = gdf[model_vars].dropna().copy()

# Queen contiguity: standard for polygon lattice data
# For point data use: libpysal.weights.KNN.from_dataframe(gdf, k=8)
w = libpysal.weights.Queen.from_dataframe(gdf)

# Row-standardize — mandatory before ML_Lag
w.transform = "r"

Row-standardization rescales each row of $W$ so its entries sum to 1. This bounds the spatial autoregressive coefficient $\rho$ to $(-1, 1)$ and stabilizes the Jacobian determinant computation that maximum-likelihood estimation requires. See spatial weight matrices for the full treatment of row-standardization and its effect on eigenvalue bounds.

Step 2 — Extract arrays #

y = df["y"].values.reshape(-1, 1)   # spreg expects a column vector
X = df[["X1", "X2"]].values         # no constant needed; spreg adds it

Step 3 — OLS baseline with Lagrange Multiplier diagnostics #

Run OLS first. The Lagrange Multiplier tests embedded in spreg.OLS tell you whether the lag specification, the error specification, or both are needed — without them you are guessing.

ols = spreg.OLS(
    y, X, w,
    spat_diag=True,                      # activates LM-Lag and LM-Error
    name_y="y",
    name_x=["X1", "X2"],
    name_w="Queen_w"
)
print(ols.summary)

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    # robust version
rlm_err_stat, rlm_err_p = ols.lm_error_robust

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

Decision rule (Anselin–Florax procedure):

• If only LM-Lag is significant → use spatial lag model.
• If only LM-Error is significant → use spatial error model.
• If both are significant, compare robust variants: whichever robust test is more significant indicates the dominant process.

Step 4 — Fit the spatial lag model #

try:
    model = spreg.ML_Lag(
        y, X, w,
        name_y="y",
        name_x=["X1", "X2"],
        name_w="Queen_w"
    )
    estimator = "ML"
except Exception as exc:
    print(f"ML failed ({exc}). Falling back to GMM.")
    model = spreg.GM_Lag(
        y, X, w,
        name_y="y",
        name_x=["X1", "X2"],
        name_w="Queen_w"
    )
    estimator = "GMM"

# Normalize rho: ML_Lag returns a scalar; GM_Lag returns a 1-D array
rho = float(model.rho) if np.isscalar(model.rho) else float(model.rho[0])

print(f"\n--- {estimator} Spatial Lag Results ---")
print(f"rho  : {rho:.4f}")
print(f"R²   : {model.r2:.4f}")
print(f"AIC  : {model.aic:.2f}")
print(model.summary)

Step 5 — Compute total marginal effects #

Coefficients from model.betas are direct effects only. The spatial multiplier propagates each unit change in $X$ through the entire neighborhood network:

$$\text{Total effect} = (I - \rho W)^{-1} \beta$$

from scipy.sparse import eye, csr_matrix
from scipy.sparse.linalg import inv as sparse_inv

n = w.n
W_sparse = csr_matrix(w.full()[0])           # dense → sparse conversion
multiplier = sparse_inv(eye(n) - rho * W_sparse)

# Direct effect of X1 (index 1 in betas, 0 is the constant)
beta_x1 = float(model.betas[1])
total_effect_x1 = multiplier.mean(axis=1).A1.mean() * beta_x1
print(f"Average total effect of X1: {total_effect_x1:.4f}")

For large $N$, full matrix inversion is expensive. Use the average of the diagonal of the multiplier (direct effect) and the row-sum average (total effect) as a computationally tractable approximation.

Interpreting the output #

Key output fields:

A $\rho$ of 0.35, for example, means that 35% of the dependent variable’s variation at each location can be attributed to a spatially weighted average of its neighbors, after conditioning on $X$. This is a substantively large spillover in most policy contexts.

Critical best practices #

Row-standardize before maximum-likelihood estimation #

w.transform = 'r' is not optional for ML_Lag. Without it the eigenvalues of $W$ are not bounded, the Jacobian determinant overflows, and $\rho$ will be estimated outside $(-1, 1)$ — producing nonsensical results with no warning. GMM estimators are more tolerant but still converge faster with standardized weights.

Choose the estimator by dataset size #

ML_Lag performs exact eigenvalue decomposition on the full $n \times n$ weight matrix. The operation is $O(n^2)$ in memory and $O(n^3)$ in compute time, making it impractical beyond roughly 2,500 observations on typical hardware. GM_Lag uses a two-stage least-squares approach with spatial lags of $X$ as instruments — it runs in sparse matrix operations and comfortably handles tens of thousands of observations.

Handle spatial islands before fitting #

An island (a polygon with no neighbors in the weight graph) causes w.transform = 'r' to divide by zero. Detect them with w.islands and resolve before fitting:

if w.islands:
    # Option A: fall back to k=1 nearest-neighbor for islands only
    w_knn = libpysal.weights.KNN.from_dataframe(gdf, k=1)
    w = libpysal.weights.w_union(w, w_knn)
    w.transform = "r"
    # Option B: drop island observations
    # keep_idx = [i for i in range(len(gdf)) if i not in w.islands]
    # gdf = gdf.iloc[keep_idx].reset_index(drop=True)

Understanding how neighbor definitions propagate through a model is covered in depth on the building custom spatial weights matrices page.

Gate on diagnostics programmatically #

Embed the LM check as an assertion or guard clause in pipelines:

if lm_lag_p > 0.05:
    raise ValueError(
        f"LM-Lag not significant (p={lm_lag_p:.3f}). "
        "Spatial lag specification is not warranted. Check for spatial error or OLS."
    )

This is the same diagnostic-gating pattern used throughout the cross-validation strategies workflow to prevent silent mis-specification before model scoring.

Interpret coefficients as direct effects only #

Unlike OLS, where each $\beta_k$ is the full marginal effect of $X_k$, the spatial lag model splits effects into direct (own-observation) and indirect (spillover-through-neighbors) components. Reporting only model.betas without the multiplier understates the true policy impact. For robust inference on marginal effects, use Monte Carlo simulation over draws from model.vm to construct confidence intervals for the total effect.

Troubleshooting #

FAQ #

When should I use ML_Lag versus GM_Lag?

Use ML_Lag for datasets with fewer than about 2,500 observations where exact maximum-likelihood is feasible. Switch to GM_Lag for larger datasets; it uses instrumental variables and sparse algebra and scales to tens of thousands of observations.

Why must spatial weights be row-standardized for ML_Lag?

Row-standardization bounds the eigenvalue spectrum of $W$ between $-1$ and $1$, which keeps $\rho$ interpretable and prevents numerical overflow in the Jacobian determinant that ML_Lag must compute.

What does a significant rho actually mean?

A significant positive $\rho$ means that high values in the dependent variable cluster next to other high values after controlling for the covariates — a direct spatial spillover effect. The OLS coefficients in a spatial lag model represent only direct effects; total effects require multiplying by the spatial multiplier $(I - \rho W)^{-1}$.

What happens if my dataset has spatial islands (polygons with no neighbors)?

Islands make row-standardization undefined (division by zero). Either assign $k=1$ nearest-neighbor fallback weights or remove the island observations before fitting.

Next steps #

For a comparison of spatial lag versus spatial error specifications and the decision criteria for choosing between them, see the spatial regression models overview. For validating any fitted spatial model against held-out data, the spatial k-fold cross-validation setup guide covers geographically blocked folds that respect spatial autocorrelation when scoring.

← Back to Spatial Regression Models

Related

• Spatial Regression Models — parent overview covering SAR, SEM, and GWR model selection
• Building Custom Spatial Weights Matrices — construct Queen, Rook, k-NN, and distance-band weights
• Spatial Weight Matrices — theory and standardization options
• Spatial K-Fold Cross-Validation Setup — validate spatial regression models with geographic blocking
• How to Calculate Moran’s I in PySAL — detect and quantify spatial autocorrelation before model fitting