Spatial Epi Public Health GIS

Choosing Between Spatial Lag and Spatial Error Models

Analysts too often pick a spatial lag or spatial error model by habit — whichever they fit last time — when the choice should follow the robust Lagrange Multiplier (LM) tests and the substantive question the analysis is asking. This guide, part of Spatial Lag and Spatial Error Regression Models, applies the Anselin-Florax decision rule programmatically from the robust LM p-values, fits the selected model, and confirms it with AIC and spatial pseudo R2R^2.

Problem Context & Constraints

Once the residuals of an OLS model on areal health data test positive for spatial autocorrelation, you have to commit to a spatial specification. The two candidates encode incompatible stories. A spatial lag (SAR) model, y=ρWy+Xβ+εy = \rho W y + X\beta + \varepsilon, says the outcome is contagious — a county’s rate responds to its neighbours’ rates, a genuine spillover. A spatial error (SEM) model, y=Xβ+uy = X\beta + u with u=λWu+εu = \lambda W u + \varepsilon, says nothing spills over; the clustering is the footprint of an omitted, spatially-structured covariate, and the spatial term is a nuisance correction that restores valid standard errors. Report a lag “spillover” when the truth is a missing variable and you have invented a causal mechanism from a data gap.

Three constraints shape the decision:

  • The classic LM tests cannot separate the two. Each classic test has power against the other form of misspecification, so on real data both usually reject. The classic pair tells you a spatial term is needed; it does not tell you which. Only the robust variants discriminate.
  • The statistics narrow, the mechanism decides. The robust LM rule is a strong prior, not a proof. When it is ambiguous — or when interpretability matters more than a marginal fit gain — the substantive question about whether the outcome plausibly diffuses has to break the tie.
  • The choice must be reproducible. A model selected by a coded rule from logged p-values can be audited; a model selected by inspection cannot.

The Anselin-Florax rule turns the robust LM comparison into a fixed routing from the ambiguous “both classic tests reject” state to a single specification:

Robust LM Decision Rule for Lag vs Error A decision flow starting from the robust LM tests. A central decision asks which robust LM test rejects. If only the robust LM-lag rejects, fit a spatial lag SAR model, interpreted as contagion or spillover. If only the robust LM-error rejects, fit a spatial error SEM model, interpreted as an omitted covariate. If both robust tests reject, the process is higher-order and a combined SARAR model or a revisited covariate set is warranted. Robust LM tests both classic rejected Which robust LM rejects? lag only error only both Spatial Lag (SAR) contagion / spillover Spatial Error (SEM) omitted covariate SARAR / revisit covariates higher-order process The robust LM tests route the ambiguous case to one specification

Prerequisites

  • python 3.11
  • geopandas 0.14.4, libpysal 4.12.1, spreg 1.4.0, numpy 1.26.4
  • Input: areal polygons with the modelled rate and covariates joined on a stable geoid, projected to a single distance-preserving CRS (EPSG:5070 for CONUS).
  • Weights: one row-standardized weights object with islands resolved, built on a stable row order so the neighbour graph — and every LM statistic — is reproducible.
  • A completed OLS diagnostic run whose robust LM p-values you have already logged; the diagnosis itself is covered in Diagnosing Spatial Autocorrelation in Regression Residuals.

Step-by-Step Solution

Apply the robust-LM rule to the logged p-values, fit exactly the indicated model with maximum likelihood so an AIC is available, and confirm the choice with the likelihood-based fit statistics. Using ML rather than GMM here is deliberate: the AIC comparison needs a log-likelihood, which the ML_* estimators expose and the GM_* estimators do not.

# 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, ML_Lag, ML_Error

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("lag.vs.error")

ALPHA = 0.05  # fixed before results are seen

def prep(gdf):
    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)} islands; resolve before selection.")
    w.transform = "r"
    return gdf, w

def robust_lm_choice(diag, alpha=ALPHA):
    """Anselin-Florax rule applied to the ROBUST LM p-values only."""
    rlag = diag["rlm_lag_p"] < alpha
    rerr = diag["rlm_error_p"] < alpha
    if rlag and not rerr:
        return "lag"
    if rerr and not rlag:
        return "error"
    if rlag and rerr:
        return "both"       # higher-order: SARAR or missing covariate
    return "neither"        # robust tests disagree with classic; re-diagnose

def select_and_fit(gdf, y_col, x_cols):
    gdf, w = prep(gdf)
    y = gdf[[y_col]].to_numpy(float)
    X = gdf[x_cols].to_numpy(float)

    ols = OLS(y, X, w=w, spat_diag=True, moran=True,
              name_y=y_col, name_x=list(x_cols))
    diag = {"rlm_lag_p": float(ols.rlm_lag[1]),
            "rlm_error_p": float(ols.rlm_error[1]),
            "moran_res_p": float(ols.moran_res[2])}
    choice = robust_lm_choice(diag)
    log.info("robust_lm=%s -> choice=%s", json.dumps(diag), choice)

    if choice == "lag":
        m = ML_Lag(y, X, w=w, name_y=y_col, name_x=list(x_cols))
        log.info("ML_Lag rho=%.4f aic=%.2f pr2=%.4f", float(m.rho), m.aic, m.pr2)
    elif choice == "error":
        m = ML_Error(y, X, w=w, name_y=y_col, name_x=list(x_cols))
        log.info("ML_Error lambda=%.4f aic=%.2f pr2=%.4f", float(m.lam), m.aic, m.pr2)
    else:
        m = None  # ambiguous: escalate to SARAR / covariate review, do not force one
        log.warning("ambiguous robust LM (%s); no single model fit", choice)
    return ols, m, choice

def confirm_with_aic(gdf, y_col, x_cols):
    """Fit BOTH models and compare AIC / spatial pseudo R^2 as a confirmation."""
    gdf, w = prep(gdf)
    y = gdf[[y_col]].to_numpy(float)
    X = gdf[x_cols].to_numpy(float)
    lag = ML_Lag(y, X, w=w)
    err = ML_Error(y, X, w=w)
    table = {
        "lag":   {"aic": round(lag.aic, 2), "pr2": round(lag.pr2, 4)},
        "error": {"aic": round(err.aic, 2), "pr2": round(err.pr2, 4)},
    }
    log.info("aic_comparison=%s", json.dumps(table))
    # lower AIC is better; it should AGREE with the robust LM verdict
    return "lag" if lag.aic < err.aic else "error", table

if __name__ == "__main__":
    gdf = gpd.read_file("county_health.gpkg").to_crs("EPSG:5070")
    ols, model, choice = select_and_fit(gdf, "lbw_rate",
                                        ["poverty", "uninsured", "pcp_rate"])
    aic_pick, table = confirm_with_aic(gdf, "lbw_rate",
                                       ["poverty", "uninsured", "pcp_rate"])
    log.info("robust_lm_choice=%s aic_choice=%s agree=%s",
             choice, aic_pick, choice == aic_pick)

The design encodes a discipline. robust_lm_choice reads only the robust p-values, because the classic pair has already done its job (telling you a spatial term is needed) and cannot separate the two forms. confirm_with_aic fits both models and compares AIC — where lower is better — and spatial pseudo R2R^2 as an independent check; the AIC winner should agree with the robust-LM verdict. When they disagree, that is a signal to investigate, not a licence to pick the answer you like. The ambiguous branches return no model at all rather than defaulting to one, so an unresolved decision cannot silently become a fitted result.

Two subtleties keep the confirmation honest. First, the spatial pseudo R2R^2 that spreg reports is not the ordinary coefficient of determination and is not directly comparable across the lag and error forms — the lag model’s pseudo R2R^2 is computed on the predicted values including the spatial multiplier, while the error model’s is computed on the mean structure alone. Use AIC, which is likelihood-based and comparable across the two specifications, as the primary confirmation and treat pseudo R2R^2 as a secondary sanity check. Second, an AIC difference of one or two points is not a decisive win for either model; when the gap is that small the two specifications fit the data almost equally well, and the choice should fall back on the robust LM verdict and the substantive mechanism rather than on chasing the marginally lower AIC. A rule of thumb worth logging: treat an AIC gap under about two as a tie, in which case interpretability — discussed below — governs.

Validation & Edge Cases

1. Residual Moran’s I that survives the fit. After fitting the selected model, re-test its residuals for spatial autocorrelation. A significant remaining residual Moran’s I is the single most important validation signal: it means the specification you chose did not capture the process. If a lag model leaves autocorrelated residuals, the mechanism is probably an error process (or both), and vice versa. The autocorrelation statistic itself is the same Moran’s I you ran on the OLS residuals, now applied to the spatial model’s residuals.

INFO robust_lm={"rlm_lag_p": 0.021, "rlm_error_p": 0.44, "moran_res_p": 0.001} -> choice=lag
INFO ML_Lag rho=0.31 aic=612.4 pr2=0.58
WARNING post_fit_moran={"moran_res_I": 0.12, "moran_res_p": 0.03}  # residual structure remains

A remaining significant Moran’s I after a lag fit, as above, means a single lag term was insufficient — escalate to a combined SARAR model or revisit the covariate set rather than reporting the lag result as final.

2. When to prefer SEM for interpretability. The AIC gap between the two models is often small. When the robust LM rule is genuinely ambiguous and the fit statistics are close, prefer the spatial error model unless you have a concrete spillover mechanism to defend. The reason is interpretive honesty: in the SEM the β\beta coefficients retain their ordinary marginal-effect meaning, so “a one-point rise in poverty is associated with β\beta more low-birth-weight cases” reads exactly as it would in OLS. In the lag model that same statement is wrong — the true effect is diffused through (I−ρW)−1(I - \rho W)^{-1} — so an audience that misreads β\beta as a marginal effect will overstate or understate the relationship. Choosing SEM when the mechanism is unclear avoids asserting a contagion story you cannot support and keeps the coefficients directly interpretable.

3. Robust tests disagree with the classic pair. Occasionally neither robust test rejects even though both classic tests did. This neither branch usually means the autocorrelation is weak or the weights are misspecified; do not force a lag or error model onto it. Re-examine the weights specification and covariate set before committing.

Compliance Notes

  • Log the rule, the inputs, and the confirmation. Record the robust LM p-values, the fixed α\alpha, the model the rule selected, the AIC and pseudo R2R^2 of both candidates, and whether the AIC check agreed. The selection must be re-derivable from the logged statistics, not from an analyst’s recollection.
  • Record the post-fit residual test. Store the residual Moran’s I of the fitted spatial model; it is the evidence that the chosen specification actually absorbed the spatial structure.
  • Deterministic inputs. Sort by stable geoid, pin every library version, and version-control the weights specification so a reviewer reproduces the same choice and the same coefficients.