Choosing a Spatial Regression Model with Lagrange Multiplier Tests

TL;DR: Call spreg.OLS(y, X, w, spat_diag=True, moran=True) and read five numbers: lm_lag, lm_error, rlm_lag, rlm_error and lm_sarma. Neither LM significant means keep OLS. One significant means fit that model. Both significant means follow the more significant robust test. Both robust significant means the specification, not the spatial term, is wrong.

Why This Matters

The choice between a lag specification and an error specification is not a matter of taste, and it is not settled by whichever model reports the higher pseudo R-squared. Fitting a lag model when the process is an error process biases every coefficient, because the endogenous term you introduced does not belong there. Fitting an error model when the process is a lag process leaves an omitted regressor in the equation and understates the total effect of every covariate, since none of the feedback through neighbours is counted. The Lagrange Multiplier sequence set out by Anselin is the standard way of deciding this from the data before either model is estimated, and it costs one extra keyword argument on the OLS call you were going to make anyway.

What makes the sequence worth learning properly, rather than memorising as “pick the smaller p-value”, is that it has an escape hatch. The rule refuses to answer in one specific case, and that refusal is the most informative outcome it produces. The conceptual background to the two specifications is covered in spatial lag vs spatial error model; this page is about running the tests, reading the printed block line by line, and knowing how much of the answer is really a statement about your spatial weight matrix. It sits inside the wider set of spatial regression models guidance.

Written as a table of the two standard tests, the whole rule fits in four cells.

Environment and Version Pinning

The diagnostics live in spreg; the weights come from libpysal; the residual check uses esda. Pin all three, because the diagnostic attribute names on spreg.OLS have moved between major versions.

bash

pip install "geopandas>=1.0.1" "libpysal>=4.12.0" "spreg>=1.8.0" \
            "esda>=2.6.0" "numpy>=1.26.0" "scipy>=1.13.0" "pandas>=2.2.0"

python

import geopandas as gpd
import numpy as np
import pandas as pd
import libpysal
import spreg
from esda.moran import Moran

Step-by-Step Implementation

The worked example uses 211 municipalities in a metropolitan region, with the log of annual mean PM2.5 as the outcome and three covariates: road density in kilometres per square kilometre, the percentage of land zoned industrial, and log population density. Atmospheric transport gives a genuine reason to expect lag dependence, and unmeasured topography gives a genuine reason to expect error dependence, so the test has real work to do.

1. Build and validate the weights

python

gdf = gpd.read_file("data/municipalities.gpkg").to_crs("EPSG:3035")
gdf = gdf.reset_index(drop=True)

w = libpysal.weights.Queen.from_dataframe(gdf, use_index=False)
w.transform = "r"                      # row-standardise before any diagnostic

assert len(w.islands) == 0, f"islands present: {w.islands}"
assert w.n_components == 1, f"{w.n_components} disconnected components"
print(f"n = {w.n}, mean neighbours = {w.mean_neighbors:.2f}")

text

n = 211, mean neighbours = 5.62

An island contributes a row of zeros to $W$ , which enters every LM statistic as a silently dropped observation. Fail the run rather than discovering it in the residuals.

2. Fit OLS with the spatial diagnostics switched on

python

y = gdf["log_pm25"].values.reshape(-1, 1)      # spreg wants a column vector
names = ["road_km_per_km2", "pct_industrial", "log_pop_density"]
X = gdf[names].values                          # no constant; spreg adds it

ols = spreg.OLS(
    y, X, w,
    spat_diag=True,          # computes LM-Lag, LM-Error, robust pair, SARMA
    moran=True,              # adds Moran's I of the OLS residuals
    name_y="log_pm25", name_x=names, name_w="Queen_r", name_ds="municipalities",
)
print(ols.summary)

The header and coefficient block of the report:

text

SUMMARY OF OUTPUT: ORDINARY LEAST SQUARES
-----------------------------------------
Data set            :  municipalities
Weights matrix      :        Queen_r
Dependent Variable  :       log_pm25    Number of Observations:         211
Mean dependent var  :         2.4187    Number of Variables   :           4
S.D. dependent var  :         0.5263    Degrees of Freedom    :         207
R-squared           :         0.6421
Adjusted R-squared  :         0.6369
Sum squared residual:         20.819    F-statistic           :    123.7911
Sigma-square        :          0.101    Prob(F-statistic)     :   3.710e-46
S.E. of regression  :          0.317    Log likelihood        :     -55.1230
Sigma-square ML     :         0.0987    Akaike info criterion :     118.246
S.E of regression ML:         0.3141    Schwarz criterion     :     131.654

------------------------------------------------------------------------------------
            Variable     Coefficient       Std.Error     t-Statistic     Probability
------------------------------------------------------------------------------------
            CONSTANT         1.86420         0.11740        15.87904         0.00000
     road_km_per_km2         0.07210         0.01190         6.05882         0.00000
      pct_industrial         0.01840         0.00430         4.27907         0.00003
     log_pop_density         0.10930         0.02060         5.30583         0.00000
------------------------------------------------------------------------------------

3. Read the spatial dependence block

The block you actually make the decision from is printed near the foot of the same report.

text

DIAGNOSTICS FOR SPATIAL DEPENDENCE
TEST                           MI/DF       VALUE           PROB
Moran's I (error)             0.2841         7.9312          0.0000
Lagrange Multiplier (lag)          1        48.9271          0.0000
Robust LM (lag)                    1        12.4413          0.0004
Lagrange Multiplier (error)        1        38.6641          0.0000
Robust LM (error)                  1         2.1783          0.1400
Lagrange Multiplier (SARMA)        2        51.1054          0.0000

Each of these is also available as an attribute, which is what you use in a pipeline. Every one is a (statistic, p-value) tuple except moran_res, which is (I, z, p).

python

print("Moran's I (error) :", tuple(round(v, 4) for v in ols.moran_res))
print("LM-Lag            :", tuple(round(v, 4) for v in ols.lm_lag))
print("Robust LM-Lag     :", tuple(round(v, 4) for v in ols.rlm_lag))
print("LM-Error          :", tuple(round(v, 4) for v in ols.lm_error))
print("Robust LM-Error   :", tuple(round(v, 4) for v in ols.rlm_error))
print("SARMA             :", tuple(round(v, 4) for v in ols.lm_sarma))

text

Moran's I (error) : (0.2841, 7.9312, 0.0)
LM-Lag            : (48.9271, 0.0)
Robust LM-Lag     : (12.4413, 0.0004)
LM-Error          : (38.6641, 0.0)
Robust LM-Error   : (2.1783, 0.14)
SARMA             : (51.1054, 0.0)

4. Encode the decision rule

Never apply the rule by eye. Write it once, so the choice is reproducible and so the “respecify” outcome cannot be quietly skipped.

python

def anselin_rule(fitted_ols, alpha=0.05):
    """Return 'ols', 'lag', 'error', 'respecify' or 'ambiguous'."""
    lag_sig = fitted_ols.lm_lag[1] < alpha
    err_sig = fitted_ols.lm_error[1] < 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 standard tests fire: each is contaminated by the other process,
    # so the robust variants are the only informative statistics left.
    rlag_sig = fitted_ols.rlm_lag[1] < alpha
    rerr_sig = fitted_ols.rlm_error[1] < alpha
    if rlag_sig and not rerr_sig:
        return "lag"
    if rerr_sig and not rlag_sig:
        return "error"
    if rlag_sig and rerr_sig:
        return "respecify"     # both survive: the model, not the term, is wrong
    return "ambiguous"         # neither survives: no spatial term is identified


choice = anselin_rule(ols)
print("selected specification:", choice)

text

selected specification: lag

5. Fit the selected model and recheck the residuals

With 211 observations, exact maximum likelihood is comfortably affordable, so ML_Lag is preferable to the instrumental-variables GM_Lag. The full mechanics of that estimator are in implementing spatial lag models in Python.

python

if choice == "lag":
    model = spreg.ML_Lag(y, X, w, name_y="log_pm25", name_x=names)
elif choice == "error":
    model = spreg.ML_Error(y, X, w, name_y="log_pm25", name_x=names)
else:
    raise RuntimeError(f"LM sequence returned '{choice}' — respecify before fitting")

rho = float(np.ravel(model.rho)[0])
print(f"rho = {rho:.4f}   AIC = {model.aic:.3f}   (OLS AIC = {ols.aic:.3f})")

text

rho = 0.4127   AIC = 79.612   (OLS AIC = 118.246)

Then confirm that the dependence has actually gone:

python

mi_ols = Moran(ols.u.flatten(), w, permutations=999)
mi_lag = Moran(model.u.flatten(), w, permutations=999)

for label, mi in (("OLS", mi_ols), ("ML_Lag", mi_lag)):
    print(f"{label:>7} residuals: I = {mi.I: .4f}  E[I] = {mi.EI: .4f}  "
          f"z = {mi.z_sim: .2f}  p_sim = {mi.p_sim:.3f}")

text

    OLS residuals: I =  0.2841  E[I] = -0.0048  z =  7.93  p_sim = 0.001
 ML_Lag residuals: I =  0.0271  E[I] = -0.0048  z =  0.76  p_sim = 0.221

6. Re-run the rule under alternative weights

Every LM statistic is a quadratic form in $W$ . Change the neighbour definition and all five numbers change, so a decision reported from a single weights matrix is a decision reported from a single untested assumption.

python

cent = gpd.GeoDataFrame(geometry=gdf.geometry.centroid, crs=gdf.crs)

specs = {
    "Queen":          libpysal.weights.Queen.from_dataframe(gdf, use_index=False),
    "Rook":           libpysal.weights.Rook.from_dataframe(gdf, use_index=False),
    "KNN k=6":        libpysal.weights.KNN.from_dataframe(cent, k=6),
    "Distance 25 km": libpysal.weights.DistanceBand.from_dataframe(
                          cent, threshold=25_000, binary=True, silence_warnings=True),
}

rows = []
for label, wi in specs.items():
    wi.transform = "r"
    m = spreg.OLS(y, X, wi, spat_diag=True)
    rows.append({
        "weights":  label,
        "LM-Lag":   f"{m.lm_lag[0]:6.2f} ({m.lm_lag[1]:.4f})",
        "RLM-Lag":  f"{m.rlm_lag[0]:6.2f} ({m.rlm_lag[1]:.4f})",
        "LM-Err":   f"{m.lm_error[0]:6.2f} ({m.lm_error[1]:.4f})",
        "RLM-Err":  f"{m.rlm_error[0]:6.2f} ({m.rlm_error[1]:.4f})",
        "choice":   anselin_rule(m),
    })

print(pd.DataFrame(rows).to_string(index=False))

text

       weights          LM-Lag         RLM-Lag          LM-Err         RLM-Err choice
         Queen  48.93 (0.0000)  12.44 (0.0004)  38.66 (0.0000)   2.18 (0.1400)    lag
          Rook  46.28 (0.0000)  13.05 (0.0003)  34.79 (0.0000)   1.56 (0.2116)    lag
       KNN k=6  44.17 (0.0000)   9.86 (0.0017)  36.02 (0.0000)   1.71 (0.1911)    lag
Distance 25 km  31.55 (0.0000)   1.02 (0.3125)  34.61 (0.0000)   4.08 (0.0434)  error

Three of the four specifications select the lag model; the 25 km distance band selects the error model. That is not a contradiction to be averaged away. A 25 km band merges each municipality with a wide, mostly non-adjacent set of neighbours, which smears any genuine adjacency feedback into what then looks like a broad, low-frequency disturbance — exactly the signature of an error process. The honest reading is that the lag conclusion holds for contiguity-scale interaction, and that a defensible choice of weights has to be argued for on substantive grounds rather than chosen after seeing which answer it gives.

Interpreting the Output

Take the six rows in order. Moran’s I (error) of 0.2841 with $z = 7.9312$ says the OLS residuals are strongly and positively autocorrelated, which is a reason to keep reading and nothing more; it does not distinguish a lag process from an error process, and it never will. The mechanics of that statistic are in how to calculate Moran’s I in PySAL.

LM-Lag of 48.9271 and LM-Error of 38.6641 are both far past the one-degree-of-freedom critical value of 3.841, so the standard pair is uninformative on its own. This is the ordinary situation, not the exception: each standard LM test is derived under the assumption that the other form of dependence is absent, so a real lag process inflates LM-Error and a real error process inflates LM-Lag. Reading the raw magnitudes and picking the larger is the single most common way this rule is misapplied.

Robust LM-Lag of 12.4413 at $p = 0.0004$ and Robust LM-Error of 2.1783 at $p = 0.1400$ resolve it. The robust variants correct each test for the local presence of the other, so a robust statistic that stays significant is evidence that form of dependence exists in its own right. Here lag survives and error does not, so the decision is a lag model. SARMA of 51.1054 on two degrees of freedom is the joint test; note that it equals 48.9271 + 2.1783 and equally 38.6641 + 12.4413. That identity holds in every run and is a useful sanity check that you have read the right rows off the report.

After fitting, $\rho = 0.4127$ with the AIC falling from 118.246 to 79.612, and residual Moran’s I falls from 0.2841 to 0.0271 with a pseudo p-value of 0.221. That is what success looks like: the dependence has been absorbed rather than merely modelled around. Warning signs are the mirror image — a residual Moran’s I that barely moves, or a $\rho$ estimate close to the boundary of $(-1, 1)$ , both of which usually mean the weights matrix is at the wrong scale.

Critical Best Practices

The attribute names are `rlm_lag` and `rlm_error`

spreg.OLS exposes the robust statistics as rlm_lag and rlm_error, not lm_lag_robust or robust_lm_lag. The joint test is lm_sarma and the residual Moran’s I is moran_res, a three-element tuple rather than a pair. Guessing these names produces an AttributeError at the worst moment, in a pipeline that has already spent minutes building weights.

Never read the robust tests unless both standard tests fired

The robust statistics are conditional corrections. When only one standard test is significant there is no contamination to correct for, and the robust variant of the other test can be significant for purely mechanical reasons. Applying the robust comparison unconditionally will occasionally invert a decision the standard pair had already made correctly, which is why anselin_rule above reaches them only inside the both-significant branch.

Treat “both robust significant” as a stop, not a coin toss

If rlm_lag and rlm_error are both significant, the data are telling you that neither spatial term explains the other away. In practice that almost always means a missing spatially patterned covariate, a wrong functional form, or a weights matrix built at the wrong range. spreg.GM_Combo will happily estimate a model with both $\rho$ and $\lambda$ , but the two parameters are weakly identified from the same neighbour structure and the fit will be unstable. Fix the specification first.

$\rho$ is not a marginal effect

Once the lag model is fitted, the coefficient on road_km_per_km2 of 0.0508 is not the total effect of road density. The feedback through $W$ multiplies it, and the total impact is $0.0508 / (1 - 0.4127) = 0.0865$ — close to the OLS estimate of 0.0721, which is exactly why the OLS number looked reasonable while resting on a wrong specification. Ask spreg.ML_Lag for spat_impacts="simple" and report direct, indirect and total effects rather than raw betas.

Residual Moran’s I after a spatial fit is indicative, not exact

esda.Moran assumes the values it is given are not residuals from a model that already contains $W$ , so the pseudo p-value on model.u is approximate. It remains the fastest smell test. For a statement you can defend, fit the same specification with spreg.GM_Lag(y, X, w, spat_diag=True) and read ak_test, the Anselin–Kelejian statistic for remaining dependence in an instrumented lag model.

python

gm = spreg.GM_Lag(y, X, w, spat_diag=True, name_y="log_pm25", name_x=names)
print(f"Anselin-Kelejian: stat = {gm.ak_test[0]:.4f}   p = {gm.ak_test[1]:.4f}")

text

Anselin-Kelejian: stat = 1.9024   p = 0.1678

Troubleshooting

Symptom	Likely cause	Fix
`AttributeError: 'OLS' object has no attribute 'lm_lag'`	`spat_diag=True` omitted, or `w` not passed to `spreg.OLS`	Pass both: `spreg.OLS(y, X, w, spat_diag=True)`
`AttributeError` on `lm_lag_robust`	Wrong attribute name	Use `rlm_lag` and `rlm_error`; the joint test is `lm_sarma`
`moran_res` missing while the LM tests are present	`moran=True` not set	Add `moran=True`; it requires `spat_diag=True` as well
All five statistics look implausibly large	Weights not row-standardised, or duplicated geometries	Set `w.transform = "r"`; drop exact duplicate rows before building `w`
Decision flips between Queen and KNN weights	Interaction range differs from the contiguity range	Report the full weights table; choose the range from theory or a variogram, not from the outcome
Both robust tests significant on every weights specification	Omitted spatially patterned covariate	Add the covariate or a trend surface, then re-run; do not fit a spatial model to compensate
Residual Moran’s I still significant after the lag fit	Weights too coarse, or an error process on top of the lag process	Rebuild `w` at a shorter range; check `ak_test` before considering `spreg.GM_Combo`

Next Steps

With the specification selected, move on to estimating and interpreting it properly — implementing spatial lag models in Python covers estimator choice, $\rho$ interpretation and impacts, while spatial lag vs spatial error model supplies the substantive reasoning that should confirm, or overrule, what the tests decided.

Frequently Asked Questions

What do I do when both robust LM tests are significant?

Treat it as evidence that the problem is not the spatial term. Both robust statistics firing means each form of dependence survives after partialling out the other, which usually points to a missing spatially patterned covariate, a wrong functional form, or a weights matrix that does not match the interaction range. Fit no spatial model on that basis. Add the missing covariate, transform the outcome, or rebuild the weights, then re-run the diagnostics. A combined model is estimable but weakly identified and will not repair a specification error.

Does the Lagrange Multiplier decision depend on the spatial weights matrix?

Entirely. Every LM statistic is a quadratic form in $W$ , so changing the neighbour definition changes all five numbers and can reverse the decision. In the worked example the Queen, Rook and k-nearest results all select the lag model, while a 25 km distance band selects the error model instead. Run the rule under two or three defensible weights specifications, report the whole table rather than one row, and treat a decision that flips as a finding about the weights rather than a result to be quietly picked from.

Why is the SARMA statistic not part of the decision rule?

SARMA is a joint two-degree-of-freedom test of lag and error dependence together, and it is arithmetically redundant: it always equals LM-Lag plus Robust LM-Error, and equally LM-Error plus Robust LM-Lag. It therefore tells you that some spatial dependence exists without indicating which form. Use it as a gate, in that a non-significant SARMA means neither individual test should be trusted, and as a consistency check that you have read the right rows, but never as the statistic that chooses the model.

Should I check residual Moran’s I after fitting the selected model?

Yes, but read it as a diagnostic rather than an exact test. The null distribution of Moran’s I on residuals from a fitted spatial model is not the one esda assumes, so the pseudo p-value is approximate. It is still the fastest way to see whether the dependence was absorbed: a drop from about 0.28 to about 0.03 with a non-significant p is what success looks like. For a formal statement use the Anselin–Kelejian test exposed as ak_test by spreg.GM_Lag when spat_diag=True.

Related

Spatial Lag vs Spatial Error Model — the two specifications the tests are choosing between
Implementing Spatial Lag Models in Python — estimating and interpreting the model once selected
Spatial Weight Matrices — the W that every LM statistic is a quadratic form in

← Back to Spatial Regression Models

Choosing a Spatial Regression Model with Lagrange Multiplier Tests

Why This Matters #

Environment and Version Pinning #

Step-by-Step Implementation #

1. Build and validate the weights #

2. Fit OLS with the spatial diagnostics switched on #

3. Read the spatial dependence block #

4. Encode the decision rule #

5. Fit the selected model and recheck the residuals #

6. Re-run the rule under alternative weights #

Interpreting the Output #

Critical Best Practices #

The attribute names are rlm_lag and rlm_error #

Never read the robust tests unless both standard tests fired #

Treat “both robust significant” as a stop, not a coin toss #

ρ\rhoρ is not a marginal effect #

Residual Moran’s I after a spatial fit is indicative, not exact #

Troubleshooting #

Next Steps #

Frequently Asked Questions #

Related

Why This Matters

Environment and Version Pinning

Step-by-Step Implementation

1. Build and validate the weights

2. Fit OLS with the spatial diagnostics switched on

3. Read the spatial dependence block

4. Encode the decision rule

5. Fit the selected model and recheck the residuals

6. Re-run the rule under alternative weights

Interpreting the Output

Critical Best Practices

The attribute names are `rlm_lag` and `rlm_error`

Never read the robust tests unless both standard tests fired

Treat “both robust significant” as a stop, not a coin toss

$\rho$ is not a marginal effect

Residual Moran’s I after a spatial fit is indicative, not exact

Troubleshooting

Next Steps

Frequently Asked Questions