Spatial Epi Public Health GIS

Diagnosing Spatial Autocorrelation in Regression Residuals

A clean-looking OLS on county health rates — high R2R^2, tidy coefficients, small p-values — routinely hides residuals that cluster in space, and when they do its standard errors are understated and its inference is wrong. This guide, part of Spatial Lag and Spatial Error Regression Models, shows how to test the OLS residuals for that hidden structure with spreg and read the diagnostics that tell you which spatial model to fit next.

Problem Context & Constraints

The naive workflow fits OLS, checks R2R^2 and the coefficient table, and stops. On areal epidemiological data that is dangerous, because the OLS estimator can be unbiased in its point estimates while its variance estimates are badly wrong. The independence assumption — that residual εi\varepsilon_i carries no information about residual εj\varepsilon_j — fails whenever adjacent units share unmeasured drivers: a regional pollution gradient, a common labor market, a data-collection boundary. The residuals then cluster, the effective sample size is smaller than nn, and the reported standard errors are too small. Coefficients that look decisive are not.

Three constraints make residual diagnosis distinct from generic regression checking:

  • The signal is in the residuals, not the outcome. The raw outcome is supposed to vary across space; that is what the covariates explain. What matters is whether the part the model fails to explain still clusters. So every test here runs on residuals, never on yy.
  • Every diagnostic depends on the weights matrix. Moran’s I of the residuals and the whole Lagrange Multiplier (LM) battery are computed against a spatial weights matrix. A different definition of “neighbour” can flip a result, so the weights specification is part of the diagnosis, not a detail.
  • A significant result must point somewhere. Detecting autocorrelation is only half the job; the diagnostics have to discriminate between a lag process and an error process so the next model is chosen for a reason, not by habit.

The diagnostic path from a fitted OLS to a specific indicated model runs through two gates — the residual Moran’s I, then the LM battery:

Residual Autocorrelation Diagnostic Flow A top-down diagnostic flow. An OLS model is fit, then Moran's I is computed on the OLS residuals. A decision checks whether the residual Moran's I is significant. If not, OLS inference is valid and the process stops. If it is significant, the full Lagrange Multiplier battery of lag, error, and robust tests is run, and the outcome routes to one of two diagnoses: if only one LM test is significant, fit that indicated model, or if both are significant, apply the robust LM tie-break. OLS model fit areal health rates Moran's I on residuals permutation inference residual Moran's I significant? no OLS inference valid stop — no spatial term yes Run LM battery lag · error · robust variants One LM significant fit that indicated model Both LM significant robust LM tie-break A significant residual Moran's I sends the LM battery to a specific diagnosis

Prerequisites

Pin the stack so the diagnostics reproduce exactly across reviewers:

  • python 3.11
  • geopandas 0.14.4 (pulls shapely 2.0.x, pyproj 3.6.x)
  • libpysal 4.12.1 — weights construction
  • spreg 1.4.0 — OLS with the spatial diagnostic battery
  • numpy 1.26.4

Input state assumed by the code:

  • Areal polygons (counties, tracts) with a modelled rate and its covariates already joined on a stable geoid.
  • CRS: a single projected, distance-preserving CRS (EPSG:5070 for CONUS) if any distance- or k-nearest-based weights are used; contiguity weights are CRS-agnostic but the layer should still be projected for consistency.
  • Weights: a row-standardized weights object with island polygons already resolved — an unresolved island is a zero row that quietly biases every test.
  • A stable row order: sort by geoid before building weights so the neighbour graph, and therefore every diagnostic, is byte-reproducible.

Step-by-Step Solution

Fit OLS once with spat_diag=True and moran=True. That single call returns Moran’s I of the residuals and the full LM battery; the interpretation logic then reads those tuples and reports the indicated direction.

# python 3.11
# geopandas==0.14.4  libpysal==4.12.1  spreg==1.4.0  numpy==1.26.4
import logging
import json
import geopandas as gpd
from libpysal.weights import Queen
from spreg import OLS

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

ALPHA = 0.05  # fixed in advance; never tuned after seeing p-values

def build_weights(gdf):
    """Row-standardized queen contiguity on a stably-ordered frame."""
    gdf = gdf.sort_values("geoid").reset_index(drop=True)
    w = Queen.from_dataframe(gdf, use_index=True)
    if w.islands:
        raise ValueError(f"{len(w.islands)} island rows; resolve before diagnosing.")
    w.transform = "r"
    return gdf, w

def diagnose(gdf, y_col, x_cols):
    """Fit OLS, extract residual Moran's I and the LM battery, and interpret."""
    gdf, w = build_weights(gdf)
    y = gdf[[y_col]].to_numpy(float)
    X = gdf[x_cols].to_numpy(float)

    m = OLS(y, X, w=w, spat_diag=True, moran=True,
            name_y=y_col, name_x=list(x_cols))

    # spreg returns each diagnostic as a tuple; unpack statistic and p-value
    diag = {
        "moran_res_I": float(m.moran_res[0]),  "moran_res_p": float(m.moran_res[2]),
        "lm_lag_p":    float(m.lm_lag[1]),      "lm_error_p":  float(m.lm_error[1]),
        "rlm_lag_p":   float(m.rlm_lag[1]),     "rlm_error_p": float(m.rlm_error[1]),
    }
    log.info("residual_diagnostics=%s", json.dumps(diag))
    return m, diag

def interpret(diag, alpha=ALPHA):
    """Turn the diagnostic tuples into a plain-language direction."""
    if diag["moran_res_p"] >= alpha:
        return "OLS inference valid: residual Moran's I not significant."

    lag, err = diag["lm_lag_p"] < alpha, diag["lm_error_p"] < alpha
    if lag and not err:
        return "Residuals autocorrelated; LM points to a spatial LAG model."
    if err and not lag:
        return "Residuals autocorrelated; LM points to a spatial ERROR model."
    if lag and err:
        # both classic tests reject -> the robust variants discriminate
        rlag, rerr = diag["rlm_lag_p"] < alpha, diag["rlm_error_p"] < alpha
        if rlag and not rerr:
            return "Both LM significant; robust LM favours the LAG model."
        if rerr and not rlag:
            return "Both LM significant; robust LM favours the ERROR model."
        return ("Both classic and both robust LM significant; likely higher-order "
                "process — consider a combined SARAR model or a missing covariate.")
    return ("Residual Moran's I significant but neither LM test rejects; re-examine "
            "the weights specification and covariate set.")

if __name__ == "__main__":
    gdf = gpd.read_file("county_health.gpkg").to_crs("EPSG:5070")
    model, d = diagnose(gdf, y_col="lbw_rate", x_cols=["poverty", "uninsured", "pcp_rate"])
    log.info("interpretation: %s", interpret(d))
    print(model.summary)  # full OLS + diagnostic report for the audit record

Read the output in a fixed order. First, moran_res_p: if the residual Moran’s I is not significant, the residuals show no spatial structure and OLS inference is valid — stop, and do not add a spatial term you cannot justify. Second, if it is significant, read the LM p-values: a single significant classic test names the model directly. Third, if both classic tests reject — which is common, because each reacts to the other’s misspecification — the robust LM p-values discriminate. The print(model.summary) line dumps the complete spreg report (coefficients, diagnostics, and the LM battery) into the record so the diagnosis is fully traceable.

The sign of the residual Moran’s I matters as much as its p-value. A significant positive value — the usual case for health rates — means the model under-predicts in some regions and over-predicts in others in spatially coherent blocks, which is what deflates the OLS standard errors: neighbouring residuals reinforce rather than cancel, so the effective sample size is smaller than nn and the naive variance is too small. A significant negative residual Moran’s I is rarer and usually signals a modelling artifact — a covariate measured at the wrong scale, or a checkerboard pattern induced by an aggregation boundary — rather than a spillover, and it should send you back to the covariate set before any spatial model. The LM tests are score tests evaluated at the OLS fit, so they are cheap by design: you get the full battery from the single OLS call above without estimating a spatial model, which is exactly why the diagnosis precedes the fit rather than following it.

Validation & Edge Cases

1. Both LM-lag and LM-error are significant. This is the default confusing case, not an anomaly. The classic LM-lag test has power against a true error process and the classic LM-error test has power against a true lag process, so on genuinely autocorrelated data both usually reject. Do not pick the one with the smaller p-value from the classic pair; that is exactly the mistake the robust variants exist to prevent. Read rlm_lag_p and rlm_error_p instead — whichever robust test remains significant when the other does not is the indicated model. If both robust tests still reject, the process is higher-order and a single spatial term will leave structure behind.

INFO residual_diagnostics={"moran_res_I": 0.28, "moran_res_p": 0.001,
  "lm_lag_p": 0.004, "lm_error_p": 0.002, "rlm_lag_p": 0.312, "rlm_error_p": 0.019}
INFO interpretation: Both LM significant; robust LM favours the ERROR model.

Here both classic tests reject, but the robust lag test is far from significant (0.312) while the robust error test holds (0.019) — an omitted spatially-structured covariate, not a spillover.

2. Weights-sensitivity of the diagnosis. Because every statistic is computed against WW, a diagnosis that flips when the weights change is not trustworthy. Re-run diagnose under at least one alternative specification — queen versus k-nearest-neighbours, or a different distance band — and confirm the indicated direction is stable.

from libpysal.weights import KNN

def weights_sensitivity(gdf, y_col, x_cols, k=6):
    gdf_s = gdf.sort_values("geoid").reset_index(drop=True)
    w_knn = KNN.from_dataframe(gdf_s, k=k); w_knn.transform = "r"
    y = gdf_s[[y_col]].to_numpy(float); X = gdf_s[x_cols].to_numpy(float)
    m = OLS(y, X, w=w_knn, spat_diag=True, moran=True)
    d = {"moran_res_p": float(m.moran_res[2]),
         "lm_lag_p": float(m.lm_lag[1]), "lm_error_p": float(m.lm_error[1]),
         "rlm_lag_p": float(m.rlm_lag[1]), "rlm_error_p": float(m.rlm_error[1])}
    log.info("knn_weights_diagnostics=%s", json.dumps(d))
    return interpret(d)

If queen contiguity points to an error model and k-nearest-neighbours points to a lag model, neither is well identified; report the instability rather than picking the answer you prefer.

3. Significant Moran’s I but no LM rejection. Occasionally the residual Moran’s I rejects while neither LM test does. This usually signals that the autocorrelation is being driven by a handful of high-leverage units or a covariate that is itself spatially misaligned. Map the residuals, check for a small cluster of outliers, and revisit the covariate set before forcing a spatial model onto what may be a data problem.

Compliance Notes

Three artifacts make a residual diagnosis defensible in review:

  • The weights specification, logged in full. Record the weights type, contiguity order or k, the row-standardization transform, and the island policy. Because every diagnostic is conditional on WW, a diagnosis without its weights definition is unreproducible.
  • The complete diagnostic battery, not just the verdict. Store the residual Moran’s I with its p-value and all four LM p-values, plus the fixed α\alpha. A reviewer must be able to re-derive the indicated model from the raw statistics.
  • Deterministic inputs. Sort by stable geoid before building weights and pin every library version, so a re-run reproduces the same neighbour graph, the same p-values, and the same interpretation.