Spatial Lag and Spatial Error Regression Models

This guide is part of the Disease Clustering & Spatial Statistical Modeling section and shows how to move from an ordinary least squares (OLS) model of area-level health rates to a global spatial regression when the OLS residuals are spatially autocorrelated. When neighbouring counties share unexplained variation, OLS violates the independence assumption, its standard errors are wrong, and any coefficient you report to a program office is indefensible. The fix is to model the spatial process explicitly — as a spatial lag (SAR) of the outcome, or as spatially correlated error — and to let the Lagrange Multiplier tests decide which.

Two model families answer two different substantive questions, and confusing them is the most common mistake in applied spatial epidemiology. A spatial lag model says the outcome in one unit is caused by the outcome in its neighbours — a genuine spillover or contagion process. A spatial error model says nothing spills over; instead an omitted, spatially-structured covariate leaves a footprint in the residuals — a nuisance to be absorbed, not a mechanism to be interpreted. The whole decision procedure below exists to keep those two stories separate.

Concept & Epidemiological Alignment

Areal health models — a county-level rate of low birth weight regressed on poverty, insurance coverage, and provider density — routinely produce residuals that are not independent. Two adjacent counties share commuting zones, media markets, hospital referral regions, and unmeasured environmental exposures, so what one model row gets wrong the next row tends to get wrong in the same direction. OLS assumes $\varepsilon_i$ are independent and identically distributed; spatially autocorrelated residuals break that assumption, deflate the estimated standard errors, and inflate the apparent significance of your covariates. A coefficient that looks decisive under OLS may be an artifact of counting the same spatial signal many times.

Global spatial regression repairs this by adding one spatial parameter to the model. The two canonical forms differ in where the spatial dependence enters. The spatial lag, or spatial autoregressive (SAR), model puts a spatially-lagged dependent term on the right-hand side:

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

Here $W$ is a row-standardized spatial weights matrix, $Wy$ is the average outcome among each unit’s neighbours, and $\rho$ measures how strongly a unit’s own outcome responds to its neighbours’ outcomes. Because $y$ appears on both sides, the reduced form is $y = (I - \rho W)^{-1}(X\beta + \varepsilon)$, which means a change in a covariate in one county propagates through the entire connected system. That propagation is the substantive claim of the lag model: the outcome is contagious. Vaccine-hesitancy clustering, communicable-disease incidence, and behaviourally-diffusing risk factors are the cases where a lag term is theoretically justified.

The spatial error, or spatial error model (SEM), instead leaves the mean structure alone and moves the spatial dependence into the disturbance:

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

Now $\lambda$ measures spatial correlation among the errors, not the outcomes. This is the model to reach for when the clustering you see is a symptom of a missing covariate — a regional pollution gradient, an unmeasured cultural boundary, a data-collection artifact that follows administrative geography. Nothing spills over in any causal sense; the spatial term is a nuisance correction that restores valid standard errors on $\beta$ without changing their interpretation. Getting this distinction wrong is not a rounding error: reporting a lag “spillover” when the truth is an omitted variable manufactures a causal story out of a data gap.

Both models are global — a single $\rho$ or $\lambda$ describes the whole study area. When the relationship between covariates and outcome itself varies across space, you want a local method instead; geographically weighted regression estimates a separate coefficient surface per location and answers a different question (“where does this predictor matter more?”) than the global spillover-versus-nuisance question these two models settle.

The Lagrange Multiplier Decision Procedure

You do not choose between lag and error by taste; you choose by the Lagrange Multiplier (LM) test battery that Luc Anselin formalized, run on the OLS residuals before you fit any spatial model. The LM tests ask, one at a time, whether adding a lag term ($\text{LM}_{\rho}$) or an error term ($\text{LM}_{\lambda}$) would significantly improve the fit. Their weakness is that each classic test is sensitive to the other form of misspecification — a true error process will often make the LM-lag test reject too, and vice versa. The robust variants ($\text{RLM}_{\rho}$, $\text{RLM}_{\lambda}$) correct each test for the possible presence of the other, and they are the tie-breaker when both classic tests are significant.

The decision follows a fixed rule so it can be automated and audited rather than argued:

This table is the operational heart of the section, and the two guides below drill into the halves that trip people up most: diagnosing spatial autocorrelation in regression residuals covers reading the LM battery and Moran’s I of the residuals correctly, and choosing between spatial lag and error models covers applying the robust-LM tie-breaker and reconciling it with the substantive question. The residual autocorrelation itself is quantified with Moran’s I on the OLS residuals, which is the smoke that tells you to run the LM tests in the first place.

Spatial Data Prerequisites

Every number in this section is downstream of the weights matrix, so it must be built, versioned, and justified before a single model runs. The weights object $W$ encodes who neighbours whom; change it from queen contiguity to a distance band and $\rho$, $\lambda$, and every LM statistic change with it. Fix the following before modelling:

• A single projected CRS. Contiguity is topological and safe on any CRS, but any distance- or k-nearest-based weights must be computed on a projected, distance-preserving CRS — EPSG:5070 for CONUS-wide work — not on WGS84 degrees.
• Row-standardized weights. Standardize $W$ so each row sums to one; then $Wy$ is a neighbour average and $\rho$ is bounded in an interpretable range. Record the transform (transform="r") in the run metadata.
• Valid, gap-free topology. A single fractured shared boundary drops a genuine neighbour link and biases the autocorrelation downward. Assert geometry validity and handle island polygons explicitly rather than letting them become disconnected rows.
• A modelled outcome that matches the geometry. Rates per administrative unit pair with contiguity weights; the response and every covariate must be aligned to the same stable unit ID before the design matrix is assembled.

The construction choices — contiguity order, k, distance band, island handling — are the subject of a dedicated guide; the constraint here is only that the same weights object feeds the OLS diagnostics and the spatial model, because comparing an LM test built on one $W$ against a model fit on another is a silent invalidation.

Production Implementation

The workflow is fixed: fit OLS in spreg with spat_diag=True and moran=True, read the LM battery, then fit exactly the model the battery indicates. The code below loads areal data, builds queen-contiguity weights, runs the diagnostics, applies the decision rule programmatically, and reports $\rho$ or $\lambda$ with a scalar impact summary. Every step logs a structured line for the audit trail.

# Pinned: spreg==1.4.0, libpysal==4.12.1, geopandas==0.14.4, numpy==1.26.4
import logging
import json
import numpy as np
import geopandas as gpd
from libpysal.weights import Queen
from spreg import OLS, GM_Lag, GM_Error_Het

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("spreg.select")

ALPHA = 0.05  # significance threshold; justify in the run record, do not tune post hoc

def load_areal(path, y_col, x_cols, crs="EPSG:5070"):
    """Load areal health data, project, and enforce a stable row order."""
    gdf = gpd.read_file(path).to_crs(crs)
    if not gdf.is_valid.all():
        raise ValueError("Invalid geometries; heal topology before building weights.")
    gdf = gdf.sort_values("geoid").reset_index(drop=True)  # stable => reproducible W
    y = gdf[[y_col]].to_numpy(float)
    X = gdf[x_cols].to_numpy(float)
    return gdf, y, X

def build_weights(gdf):
    """Row-standardized queen contiguity; islands must be resolved upstream."""
    w = Queen.from_dataframe(gdf, use_index=True)
    if w.islands:
        raise ValueError(f"{len(w.islands)} island polygons; resolve before modelling.")
    w.transform = "r"  # row-standardize: Wy is a neighbour average
    return w

def ols_diagnostics(y, X, w, x_names):
    """OLS with the full LM spatial-diagnostic battery on the residuals."""
    m = OLS(y, X, w=w, spat_diag=True, moran=True,
            name_y="rate", name_x=x_names)
    diag = {
        "moran_res_I": float(m.moran_res[0]), "moran_res_p": float(m.moran_res[2]),
        "lm_lag": float(m.lm_lag[0]),   "lm_lag_p": float(m.lm_lag[1]),
        "lm_error": float(m.lm_error[0]), "lm_error_p": float(m.lm_error[1]),
        "rlm_lag": float(m.rlm_lag[0]),   "rlm_lag_p": float(m.rlm_lag[1]),
        "rlm_error": float(m.rlm_error[0]), "rlm_error_p": float(m.rlm_error[1]),
    }
    log.info("ols_diagnostics=%s", json.dumps(diag))
    return m, diag

def choose_model(diag, alpha=ALPHA):
    """Apply the Anselin LM decision rule; returns 'ols' | 'lag' | 'error' | 'both'."""
    lag_sig = diag["lm_lag_p"] < alpha
    err_sig = diag["lm_error_p"] < alpha
    if not lag_sig and not err_sig:
        return "ols"
    if lag_sig and not err_sig:
        return "lag"
    if err_sig and not lag_sig:
        return "error"
    # both classic tests reject -> defer to the robust variants
    rlag_sig = diag["rlm_lag_p"] < alpha
    rerr_sig = diag["rlm_error_p"] < alpha
    if rlag_sig and not rerr_sig:
        return "lag"
    if rerr_sig and not rlag_sig:
        return "error"
    return "both"  # ambiguous: compare magnitudes / consider SARAR

def total_impact_lag(model):
    """Average total impact of each X in a lag model: beta / (1 - rho).

    In a lag model the marginal effect of X on y is NOT beta; a change
    propagates through (I - rho W)^-1, so the scalar summary divides by (1 - rho).
    """
    rho = float(model.rho)
    betas = model.betas[1:-1].ravel()  # drop intercept and the trailing rho
    return rho, betas / (1.0 - rho)

def run(path, y_col, x_cols):
    gdf, y, X = load_areal(path, y_col, x_cols)
    w = build_weights(gdf)
    ols, diag = ols_diagnostics(y, X, w, x_names=list(x_cols))
    choice = choose_model(diag)
    log.info("indicated_model=%s (residual Moran's I p=%.4g)", choice, diag["moran_res_p"])

    if choice == "ols":
        return {"model": "ols", "result": ols}
    if choice in ("lag", "both"):
        m = GM_Lag(y, X, w=w, name_y="rate", name_x=list(x_cols))
        rho, impacts = total_impact_lag(m)
        log.info("GM_Lag rho=%.4f total_impacts=%s", rho, np.round(impacts, 5).tolist())
        return {"model": "lag", "result": m, "rho": rho, "impacts": impacts}
    # error path — use the heteroskedasticity-robust GMM error estimator by default
    m = GM_Error_Het(y, X, w=w, name_y="rate", name_x=list(x_cols))
    lam = float(m.betas[-1][0])
    log.info("GM_Error_Het lambda=%.4f", lam)
    return {"model": "error", "result": m, "lambda": lam}

Three implementation choices in that code are deliberate and defensible. First, the decision rule is coded, not eyeballed — the same diag dictionary that drives the branch is the one written to the audit log, so a reviewer sees exactly why a lag model was chosen. Second, total_impact_lag never reports the raw $\beta$ as a marginal effect; in a lag model the effect of a covariate is diffused through the inverse $(I - \rho W)^{-1}$, and quoting $\beta$ alone understates the true effect and misleads the reader. Third, the error branch defaults to GM_Error_Het rather than the homoskedastic estimator, because areal health rates over populations of wildly different size are almost always heteroskedastic, and ignoring that produces error variances — and therefore inference — that are simply wrong.

For maximum-likelihood estimation of the lag model when you need the log-likelihood for an AIC comparison, swap GM_Lag for ML_Lag (also in spreg); the GMM estimators above are faster and robust to non-normal errors, while ML gives you the likelihood-based fit statistics the comparison guide uses.

Parameter Selection & Tuning

The parameters that move the results are not hyperparameters in the machine-learning sense; they are modelling decisions that each need a written justification.

• The weights specification. This is the single most consequential choice. Queen contiguity is the sensible default for polygon health data; k-nearest-neighbours guarantees a fixed neighbour count but can link across natural barriers; a distance band respects geography but strands sparse rural units with empty rows. Run the diagnostics under at least two plausible $W$ definitions and confirm the indicated model is stable — if the choice between lag and error flips with the weights, neither is well identified.
• Row standardization. Keep it on for interpretability. Binary weights make $\rho$ and $\lambda$ scale with neighbour count and are hard to compare across units of differing connectivity.
• The significance threshold $\alpha$. Fix it once, in the config, at 0.05 or 0.01, and never lower it after seeing the p-values. Post-hoc threshold tuning to force a preferred model is the spatial-regression equivalent of p-hacking.
• Estimator family. GMM (GM_*) for speed and robustness to non-normality; ML (ML_*) when you need a likelihood for AIC. For the error model prefer the _Het variants unless a formal test clears homoskedasticity.

Edge Cases & Failure Modes

• Weights misspecification. A $W$ that is too dense (large distance band, high k) smears real local structure and can wash out a true lag signal; one that is too sparse fractures it. There is no autocorrelation “truth” independent of $W$, so the weights choice must be argued from geography and mobility, then held fixed.
• Endogeneity in the lag term. $Wy$ is correlated with the error by construction, so OLS on a lag model is biased and inconsistent — this is exactly why GM_Lag uses spatial two-stage least squares with $WX$ as instruments rather than plain OLS. Never fit a lag model by adding $Wy$ as an ordinary regressor.
• Heteroskedasticity. Population-weighted rates have variance that scales inversely with the denominator. Standard SEM inference assumes homoskedastic $\varepsilon$; use GM_Error_Het to get a heteroskedasticity-consistent $\lambda$ and standard errors, or the reported significance is unreliable.
• Both robust tests significant. When $\text{RLM}_{\rho}$ and $\text{RLM}_{\lambda}$ both reject, the process is likely higher-order; a combined SARAR (lag + error) model or a re-examination of the omitted-covariate set is warranted rather than forcing a single term.
• Residual autocorrelation that survives the fit. After fitting a lag model, re-test the residuals; a lag term that leaves significant residual Moran’s I means the spillover specification did not capture the process and an error component may still be needed.

Compliance & Audit Controls

A spatial regression that informs resource allocation must be reconstructable on demand. Bind every run to an immutable record: a SHA-256 hash of the input geometry, the exact weights specification (type, order/k/band, transform, island policy), the full LM battery with p-values, the significance threshold, the chosen model and its $\rho$ or $\lambda$, and the pinned library versions. Because the model choice is deterministic given those inputs, a reviewer re-running the pipeline must reproduce the same branch and the same coefficients.

• Input geometry hashed (SHA-256) and stored with the run record
• Single projected CRS enforced before any distance-based weights
• Weights specification version-controlled; islands resolved, not silently dropped
• OLS fit with spat_diag=True and moran=True; full LM battery logged
• Model chosen by the coded LM decision rule, not by inspection
• α fixed in config before results are seen
• Lag models estimated with spatial 2SLS/ML, never OLS on Wy
• Error models use the heteroskedasticity-robust estimator unless homoskedasticity is tested and cleared
• Impacts reported as diffused effects, not raw β, for lag models
• Residuals re-tested for autocorrelation after fitting
• Library versions pinned and recorded alongside results

Frequently Asked Questions

When is a spatial lag model the right choice rather than a spatial error model? Choose the lag model when the outcome itself plausibly spills across boundaries — a contagious disease, a diffusing behaviour, a shared-market effect — and the robust LM-lag test is significant while robust LM-error is not. Choose the error model when the clustering is a footprint of an omitted, spatially-structured covariate rather than a real spillover, indicated by a significant robust LM-error. The statistics narrow it down; the substantive mechanism confirms it.

Why can’t I just add the neighbour average of the outcome as a regressor in OLS? Because the spatially-lagged outcome $Wy$ is endogenous — it is correlated with the model error by construction — so OLS estimates of $\rho$ and $\beta$ are biased and inconsistent. Spatial lag models are estimated with spatial two-stage least squares (using $WX$ as instruments) or by maximum likelihood, which is exactly what spreg’s GM_Lag and ML_Lag do.

My OLS looks great — good $R^2$, significant covariates. Why check for spatial autocorrelation at all? Because a strong-looking OLS with spatially autocorrelated residuals has understated standard errors, so its significance is inflated and its inference is invalid. A good $R^2$ does not rescue you; run Moran’s I on the residuals and the LM battery before trusting any coefficient.

What does it mean if both classic LM tests are significant? Each classic test reacts to the other form of misspecification, so both rejecting is common and uninformative on its own. Defer to the robust LM tests, which correct each for the presence of the other; if both robust tests also reject, the process is likely higher-order and a combined SARAR model or a better covariate set is warranted.

Related Topics

• Diagnosing Spatial Autocorrelation in Regression Residuals — read the LM battery and residual Moran’s I that trigger this whole procedure.
• Choosing Between Spatial Lag and Spatial Error Models — apply the robust-LM tie-breaker and reconcile it with the substantive question.
• Spatial Weights Matrix Construction — the weights object every LM test and spatial coefficient depends on.
• Global & Local Moran’s I Implementation — the residual autocorrelation test that signals a spatial model is needed.
• Geographically Weighted Regression — the local alternative when the covariate relationships themselves vary across space.