Chapter 13: Spatial Regression¶
Part II — Geographic Data Science
Learning Objectives¶
By the end of this chapter you will be able to:
- Explain why ordinary least squares (OLS) can be the wrong tool for geographic data, and diagnose the problem using the same Moran's I test from Chapter 12
- Fit and interpret a spatial lag model, where an outcome depends directly on its neighbors' outcomes, reading its effects as direct, indirect and total impacts rather than as a raw coefficient
- Fit and interpret a spatial error model, where unmeasured spatially clustered factors — not the outcome itself — are what's contagious
- Use Lagrange Multiplier diagnostics to choose between the two rather than guessing — and recognize when neither is enough and a Spatial Durbin model is the candidate
Why Regression Needs to Know About Space¶
Every regression course starts with the same assumption, usually stated once and then never revisited: errors are independent of each other. For geographic data, that assumption is rarely true. A municipality's unexplained income, its unexplained disease rate, its unexplained anything tends to resemble its neighbors' unexplained anything — not because of a coding mistake, but because space itself carries information: shared infrastructure, shared climate, shared history, spillover effects that no explanatory variable in the model captured.
Chapter 12 built the tool for detecting exactly this: Moran's I. This chapter puts it to work on regression residuals specifically, then introduces two ways of fixing what it finds.
# Standard imports
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import box
import statsmodels.api as sm
from libpysal.weights import Queen
from esda.moran import Moran
import spreg
np.random.seed(42)
A Synthetic Spatial Process¶
Chapter 12 used a synthetic 15x15 grid instead of the Maranhão municipalities for its checkerboard and clustered patterns, for a specific reason: proving a statistic works requires knowing the ground truth, which real data never hands you. The same reasoning applies here, more strongly — to check that a spatial regression model recovers the right answer, the data has to be built from a known answer in the first place. This chapter reuses the same grid convention (and the same Queen weights machinery from Chapter 10) before the final section confirms the exact same code runs unmodified on maranhao.
# Build a 15x15 regular grid, same convention as Chapter 12
n_side = 15
cells_geom = [box(j, i, j + 1, i + 1) for i in range(n_side) for j in range(n_side)]
grid = gpd.GeoDataFrame({"geometry": cells_geom})
# A single explanatory variable, uncorrelated with position
grid["x"] = np.random.uniform(0, 10, len(grid))
w = Queen.from_dataframe(grid, use_index=False)
w.transform = "r" # row-standardized, exactly as in Chapter 10
n = len(grid)
I_n = np.eye(n)
W_full = w.full()[0]
A Quick OLS Refresher¶
Before diagnosing what's wrong with OLS on spatial data, a one-paragraph reminder of what OLS does: given an outcome y and one or more explanatory variables X, it finds the coefficients that minimize the sum of squared residuals — the gaps between predicted and observed y. statsmodels reports each coefficient, its standard error, and R², the share of y's variance the model explains.
To generate a y with a known spatial lag structure, this section builds it deliberately: $y = \rho W y + \beta_0 + \beta_1 x + \epsilon$. Solving that equation for y (rather than sampling it directly) is what guarantees the dependency is really there — rho_true = 0.6 means every cell's value is genuinely pulled toward its neighbors' values by that amount, not just correlated with them by coincidence.
rho_true = 0.6
beta0_true, beta1_true = 2.0, 1.5
eps = np.random.normal(0, 1, n)
rhs = beta0_true + beta1_true * grid["x"].values + eps
# Solve (I - rho*W) y = rhs for y, so the lag dependency is exact by construction
grid["y_lag"] = np.linalg.solve(I_n - rho_true * W_full, rhs)
X_const = sm.add_constant(grid["x"].values)
ols_lag = sm.OLS(grid["y_lag"].values, X_const).fit()
print(ols_lag.summary().tables[1])
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 15.9968 0.215 74.491 0.000 15.574 16.420
x1 1.5366 0.038 40.153 0.000 1.461 1.612
==============================================================================
Look at both coefficients. The intercept is the loudest symptom: about 16.0 against a true 2.0, because OLS has nowhere else to put the neighbor multiplier. The slope, about 1.54, looks reassuringly close to the true β = 1.5 — but that closeness is misleading, for a reason the spatial lag model makes explicit below: in a lag process, β is not the effect of x on y at all. That's the real danger: a model with a plausible-looking R² and significant coefficients that is nonetheless mis-specified.
Are the Residuals Spatially Autocorrelated?¶
Chapter 12's exact tool answers this directly: run Moran.I on the OLS residuals themselves, using the same w built above.
mi_resid = Moran(ols_lag.resid, w)
print(f"Moran's I on OLS residuals: {mi_resid.I:.3f} (p_sim = {mi_resid.p_sim:.4f})")
Moran's I on OLS residuals: 0.390 (p_sim = 0.0010)
A Moran's I this far from zero, with p_sim far below any conventional threshold, is the same signal Chapter 12 used to detect clustering in a plain variable — except here it's a diagnosis, not just a description: the model is missing something structural, and that something has a spatial shape.
The Spatial Lag Model¶
The spatial lag model builds the missing term directly into the regression: $y = \rho W y + X\beta + \epsilon$. Because y appears on both sides, ordinary OLS can't fit it — spreg.ML_Lag uses maximum likelihood instead, which is why the model needs the weights object w passed in explicitly rather than just an extra column of data.
y_arr = grid[["y_lag"]].values
X_arr = grid[["x"]].values
ml_lag = spreg.ML_Lag(y_arr, X_arr, w=w, name_y="y", name_x=["x"])
print(f"True rho: {rho_true}")
print(f"Estimated rho: {ml_lag.rho:.3f}")
print(f"True beta1: {beta1_true} OLS beta1: {ols_lag.params[1]:.3f} ML_Lag beta1: {ml_lag.betas.flatten()[1]:.3f}")
# Impacts: in a lag model, x moves y through the multiplier (I - rho W)^-1, not through beta alone
M = np.linalg.inv(I_n - ml_lag.rho * W_full)
b1 = ml_lag.betas.flatten()[1]
direct = b1 * np.trace(M) / n # average own-cell effect, feedback included
total = b1 / (1 - ml_lag.rho) # row sums of M are all 1/(1 - rho) for row-standardized W
print(f"Impacts of x — direct: {direct:.3f} indirect: {total - direct:.3f} total: {total:.3f}")
print(f"True total impact: {beta1_true / (1 - rho_true):.3f} OLS slope: {ols_lag.params[1]:.3f}")
ML_Lag True rho: 0.6 Estimated rho: 0.603 True beta1: 1.5 OLS beta1: 1.537 ML_Lag beta1: 1.453 Impacts of x — direct: 1.569 indirect: 2.095 total: 3.664 True total impact: 3.750 OLS slope: 1.537
ml_lag.rho lands close to the true 0.6 — the model recovered a parameter that never appeared anywhere in the input data, only in how y was constructed.
The impacts are the second payoff, and the one that matters for interpretation. In a lag model, raising x in one cell raises that cell's y, which raises its neighbors' y, which feeds back again: the full effect is β multiplied by (I − ρW)⁻¹, not β alone. Following LeSage and Pace (2009), the direct impact averages the diagonal of that multiplier (a cell's own x on its own y, feedback included), the total impact is β/(1 − ρ) for row-standardized weights, and the indirect impact — the spillover onto neighbors — is the difference.
Here the true total effect of x is 1.5 / (1 − 0.6) = 3.75, and ML_Lag's total recovers it closely. OLS's slope of about 1.54, which looked so reassuring above, is estimating something close to the direct effect only: it misses more than half of what x actually does to y. That is the rule to carry into the real example below — compare an OLS slope with the lag model's impacts, never with its raw beta1.
The Spatial Error Model¶
A spatial lag process says the outcome itself is contagious — a municipality's income genuinely depends on its neighbors' income. A spatial error process tells a different story: x fully explains y on average, but the leftover noise — everything x doesn't capture — happens to cluster in space anyway, because of some omitted variable (soil quality, an unmapped historical boundary, anything) that itself varies smoothly across the map. The model is $y = X\beta + u$, with $u = \lambda W u + \epsilon$.
The distinction matters practically: a lag process means neighbors' y belongs in the model; an error process means the coefficients on x were fine all along, but the standard errors around them were not.
lambda_true = 0.6
eps2 = np.random.normal(0, 1, n)
u = np.linalg.solve(I_n - lambda_true * W_full, eps2)
grid["y_err"] = beta0_true + beta1_true * grid["x"].values + u
ols_err = sm.OLS(grid["y_err"].values, X_const).fit()
mi_resid_err = Moran(ols_err.resid, w)
print(f"OLS beta1: {ols_err.params[1]:.3f} (true: {beta1_true})")
print(f"Moran's I on OLS residuals: {mi_resid_err.I:.3f} (p_sim = {mi_resid_err.p_sim:.4f})")
OLS beta1: 1.525 (true: 1.5) Moran's I on OLS residuals: 0.372 (p_sim = 0.0010)
ml_err = spreg.ML_Error(grid[["y_err"]].values, X_arr, w=w, name_y="y", name_x=["x"])
print(f"True lambda: {lambda_true}")
print(f"Estimated lambda: {ml_err.lam:.3f}")
print(f"OLS beta1: {ols_err.params[1]:.3f} ML_Error beta1: {ml_err.betas.flatten()[1]:.3f}")
ML_Error True lambda: 0.6 Estimated lambda: 0.681 OLS beta1: 1.525 ML_Error beta1: 1.512
The pattern confirms the distinction from above: ml_err's beta1 barely moves from OLS's — because the error model was never trying to fix x's coefficient, only to account for the fact that the leftover noise isn't independent across neighboring cells. What actually changes (not shown by a point estimate) is the standard errors around that coefficient, which ML_Error computes correctly and plain OLS does not.
Choosing Between Lag and Error: The LM Tests¶
Both models above were fit knowing, by construction, which process generated the data — a luxury real analysis never has. spreg.OLS with spat_diag=True reports a set of Lagrange Multiplier (LM) tests that make the same call from data alone, by checking each model's diagnostic in a way that's robust to the other kind of dependence also being present.
ols_diag = spreg.OLS(
grid[["y_lag"]].values, X_arr, w=w,
spat_diag=True, moran=True,
name_y="y", name_x=["x"],
)
print(ols_diag.summary[ols_diag.summary.find("DIAGNOSTICS FOR SPATIAL DEPENDENCE"):])
DIAGNOSTICS FOR SPATIAL DEPENDENCE - SARERR - TEST MI/DF VALUE PROB Moran's I (error) 0.3904 11.227 0.0000 Lagrange Multiplier (lag) 1 216.128 0.0000 Robust LM (lag) 1 113.000 0.0000 Lagrange Multiplier (error) 1 119.517 0.0000 Robust LM (error) 1 16.389 0.0001 Lagrange Multiplier (SARMA) 2 232.517 0.0000 - Spatial Durbin - TEST DF VALUE PROB LM test for WX 1 113.000 0.0000 Robust LM WX test 1 16.389 0.0001 Lagrange Multiplier (lag) 1 216.128 0.0000 Robust LM Lag - SDM 1 119.517 0.0000 Joint test for SDM 2 232.517 0.0000 ================================ END OF REPORT =====================================
On data generated by a genuine lag process, Robust LM (lag) comes out far larger than Robust LM (error) — the robust versions matter here specifically because the plain (non-robust) LM tests tend to both look significant whenever any spatial dependence is present, lag or error alike. The robust pair is what actually discriminates between the two. Run the same diagnostic on y_err as an exercise below, and the ranking flips.
The full decision rule (Anselin et al., 1996) has four branches. If neither plain LM test is significant, keep OLS. If only one is, fit that model. If both are, look at the robust pair and fit the model whose robust test stays significant. And if both robust tests stay significant, neither the lag nor the error model alone is enough: a model with both terms (SARMA), or a Spatial Durbin model that also adds the neighbors' explanatory variables WX, is the candidate. spreg prints an LM test for WX in its Spatial Durbin block for exactly that last case — and it is the branch the real data below lands on.
From Grid to Real Geography — and a Bridge Ahead¶
Nothing above depends on the grid being synthetic. Swap grid[["x"]].values for any column of maranhao and w for the Queen.from_dataframe(maranhao) weights already built in Chapters 10–12, and spreg.ML_Lag / spreg.ML_Error run exactly as shown — the synthetic version exists only so the true rho and lambda are known well enough to confirm the models are recovering them correctly, not because the mechanics themselves are any different on real municipalities.
One extension is deliberately left for later: everything in this chapter assumes a continuous outcome. Chapter 26's land-use change model instead asks a yes/no question — does a given cell convert from forest to another use? — which calls for a logistic variant of these same ideas rather than the linear regression used here. The spatial reasoning carries over unchanged; only the outcome's shape does.
A Real Example: GDP per Capita in Maranhão¶
A teaching example, not an economic study
This section uses one year of data and one or two explanatory variables, with no economic theory behind the specification. Its job is to show how the diagnostics and models above behave on real, messy data — and how easily a spatial model's story changes when the specification changes. Read the numbers as illustrations of the method, not as findings about Maranhão's economy.
The grid above proves the mechanism works on data with a known answer. Real Maranhão economic data has no such guarantee — but it lets the same three tools (Moran's I on residuals, ML_Lag, ML_Error) be run on a question with a real geographic flavor: is a municipality's GDP per capita related to its neighbors', or only to its own economic structure?
The explanatory variable is agro_share — agriculture's share of each municipality's GDP, computed from IBGE's PIB dos Municípios value-added breakdown (the sectoral columns of the same table that publishes the headline PIB figures). The outcome is GDP per capita for 2021, the last year with a complete sectoral breakdown. One municipality, Tasso Fragoso, has a GDP per capita over R$260,000 — nearly 8 times São Luís's, on a population of under 9,000 driven almost entirely by soy agribusiness — so, as with any strongly right-skewed economic variable, the regression below uses log(pib_percapita) rather than the raw value.
url = (
"https://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/"
"malhas_municipais/municipio_2025/UFs/MA/MA_Municipios_2025.zip"
)
maranhao = gpd.read_file(url)
pib = pd.read_csv("pib_municipios_ma.csv")
pop = pd.read_csv("populacao_municipios_ma.csv")
pib_2021 = pib[pib["ano"] == 2021].merge(pop[pop["ano"] == 2021], on=["id_municipio", "ano"])
pib_2021["pib_percapita"] = pib_2021["pib"] / pib_2021["populacao"]
pib_2021["agro_share"] = pib_2021["va_agropecuaria"] / pib_2021["pib"]
pib_2021["ind_share"] = pib_2021["va_industria"] / pib_2021["pib"] # used at the end of the chapter
maranhao["CD_MUN"] = maranhao["CD_MUN"].astype(int)
maranhao = maranhao.merge(
pib_2021[["id_municipio", "pib_percapita", "agro_share", "ind_share"]],
left_on="CD_MUN", right_on="id_municipio", how="left",
)
maranhao["log_pib_pc"] = np.log(maranhao["pib_percapita"])
fig, ax = plt.subplots(figsize=(6, 8))
maranhao.plot(column="log_pib_pc", cmap="viridis", legend=True, ax=ax)
ax.set_title("Maranhão: log(GDP per capita), 2021")
ax.set_axis_off()
No missing values and no islands here — every one of the 217 municipalities has at least one Queen neighbor, unlike Chapter 12's São Paulo example, which had to drop one island municipality before fitting Queen.from_dataframe.
w_ma = Queen.from_dataframe(maranhao, use_index=False)
print("Islands:", w_ma.islands)
w_ma.transform = "r"
X_ma = sm.add_constant(maranhao["agro_share"].values)
ols_ma = sm.OLS(maranhao["log_pib_pc"].values, X_ma).fit()
print(ols_ma.summary().tables[1])
Islands: []
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 9.0727 0.062 146.016 0.000 8.950 9.195
x1 1.3991 0.269 5.203 0.000 0.869 1.929
==============================================================================
agro_share's coefficient comes out positive and significant: in this simple model, municipalities with a larger farming share of GDP tend to have higher GDP per capita — the pattern of a few low-population, agribusiness-heavy municipalities in the south, like Tasso Fragoso, posting per-capita figures far above the coastal capital. Keep in mind how little this model contains: one explanatory variable explains about 11% of the variation (R² ≈ 0.11). The last section of the chapter comes back to what that leaves out.
mi_ma = Moran(ols_ma.resid, w_ma)
print(f"Moran's I on OLS residuals: {mi_ma.I:.3f} (p_sim = {mi_ma.p_sim:.4f})")
Moran's I on OLS residuals: 0.269 (p_sim = 0.0010)
The residuals are far from independent — an I of this size, this significant, says agro_share alone is leaving a genuinely spatial signal on the table, exactly the diagnosis that opened this chapter.
diag_ma = spreg.OLS(
maranhao[["log_pib_pc"]].values, maranhao[["agro_share"]].values, w=w_ma,
spat_diag=True, moran=True,
name_y="log_pib_pc", name_x=["agro_share"],
)
print(diag_ma.summary[diag_ma.summary.find("DIAGNOSTICS FOR SPATIAL DEPENDENCE"):])
DIAGNOSTICS FOR SPATIAL DEPENDENCE - SARERR - TEST MI/DF VALUE PROB Moran's I (error) 0.2694 6.483 0.0000 Lagrange Multiplier (lag) 1 57.483 0.0000 Robust LM (lag) 1 37.728 0.0000 Lagrange Multiplier (error) 1 38.319 0.0000 Robust LM (error) 1 18.564 0.0000 Lagrange Multiplier (SARMA) 2 76.047 0.0000 - Spatial Durbin - TEST DF VALUE PROB LM test for WX 1 37.728 0.0000 Robust LM WX test 1 18.564 0.0000 Lagrange Multiplier (lag) 1 57.483 0.0000 Robust LM Lag - SDM 1 38.319 0.0000 Joint test for SDM 2 76.047 0.0000 ================================ END OF REPORT =====================================
Both robust tests stay significant — Robust LM (lag) ≈ 37.7 and Robust LM (error) ≈ 18.6, both with p < 0.0001 — which is the fourth branch of the decision rule above. The lag term is the stronger of the two, so ML_Lag is the natural first model to fit. But the diagnostics already warn that it will not be the last: the LM test for WX in the Spatial Durbin block is significant too, meaning the neighbors' agro_share probably belongs in the model as well.
ml_lag_ma = spreg.ML_Lag(
maranhao[["log_pib_pc"]].values, maranhao[["agro_share"]].values, w=w_ma,
name_y="log_pib_pc", name_x=["agro_share"],
)
n_ma = len(maranhao)
M_ma = np.linalg.inv(np.eye(n_ma) - ml_lag_ma.rho * w_ma.full()[0])
b_ma = ml_lag_ma.betas.flatten()[1]
direct_ma = b_ma * np.trace(M_ma) / n_ma
total_ma = b_ma / (1 - ml_lag_ma.rho)
print(f"Estimated rho: {ml_lag_ma.rho:.3f}")
print(f"OLS slope (agro_share): {ols_ma.params[1]:.3f}")
print(f"ML_Lag impacts of agro_share — direct: {direct_ma:.3f} indirect: {total_ma - direct_ma:.3f} total: {total_ma:.3f}")
ML_Lag Estimated rho: 0.613 OLS slope (agro_share): 1.399 ML_Lag impacts of agro_share — direct: 0.369 indirect: 0.490 total: 0.859
With ρ ≈ 0.61, the lag model's raw beta1 (about 0.33) is not comparable with OLS's 1.40 — its impacts are. The total impact of agro_share on log GDP per capita is about 0.86: smaller than OLS's 1.40, so — within this one-variable specification — OLS overstates the association. And more than half of that total (about 0.49 of 0.86) is indirect — in this model, a municipality's agricultural share goes with higher GDP per capita in its neighbors, propagated through the lag term, more than in itself.
Fitting ML_Error on the same data, as a check, gives a negative beta1 (about −0.51), but with z ≈ −1.54 (p ≈ 0.12) it cannot be told apart from zero. That is a sign that the error model is the wrong specification here, not evidence of a real negative effect. The diagnostics pointed one step further, to the Spatial Durbin model:
sdm_ma = spreg.ML_Lag(
maranhao[["log_pib_pc"]].values, maranhao[["agro_share"]].values, w=w_ma,
slx_lags=1, spat_impacts="full", # slx_lags=1 adds W * agro_share: a Spatial Durbin model
name_y="log_pib_pc", name_x=["agro_share"],
)
for name, b, (z, p) in zip(sdm_ma.name_x, sdm_ma.betas.flatten(), sdm_ma.z_stat):
print(f"{name:14s} {b:7.3f} z = {z:6.2f} p = {p:.4f}")
lr = 2 * (sdm_ma.logll - ml_lag_ma.logll)
print(f"\nLR test, Spatial Durbin vs. lag model: {lr:.1f} on 1 df")
print(sdm_ma.summary[sdm_ma.summary.find("SPATIAL DURBIN MODEL IMPACTS"):])
ML_Lag
CONSTANT 4.229 z = 6.47 p = 0.0000
agro_share -0.911 z = -2.63 p = 0.0085
W_agro_share 2.333 z = 5.05 p = 0.0000
W_log_pib_pc 0.520 z = 7.11 p = 0.0000
LR test, Spatial Durbin vs. lag model: 24.5 on 1 df
SPATIAL DURBIN MODEL IMPACTS
Impacts computed using the 'full' method.
Variable Direct Indirect Total
agro_share -0.9755 3.9352 2.9597
================================ END OF REPORT =====================================
With the neighbors' agricultural share (W_agro_share) in the model, a municipality's own share turns negative (about −0.91) while its neighbors' share is strongly positive (about +2.33), and the likelihood-ratio test against the lag model (about 24.5 on one degree of freedom) favors the Durbin specification. It is tempting to read this geographically — being surrounded by agribusiness matters, being agricultural does not. Before believing a story like that, though, it is worth checking two things every applied analyst checks: whether a few observations are driving the result, and whether an obvious variable is missing. The next section does both, and the story does not survive.
Checking the Story: Influential Municipalities and a Missing Variable¶
Cook's distance measures how much the fitted OLS model would change if one observation were removed; a common rule of thumb flags values above 4/n (about 0.018 here). statsmodels computes it directly from the model already fitted:
influence = ols_ma.get_influence()
cooks_d = influence.cooks_distance[0]
flagged = (
maranhao.assign(cooks_d=cooks_d)
.sort_values("cooks_d", ascending=False)
[["NM_MUN", "cooks_d", "agro_share", "ind_share", "pib_percapita"]]
.head(5)
)
print(f"Rule-of-thumb threshold 4/n = {4 / len(maranhao):.3f}")
flagged.round(3)
Rule-of-thumb threshold 4/n = 0.018
| NM_MUN | cooks_d | agro_share | ind_share | pib_percapita | |
|---|---|---|---|---|---|
| 114 | Tasso Fragoso | 0.997 | 0.792 | 0.029 | 263760.179 |
| 96 | Godofredo Viana | 0.241 | 0.004 | 0.817 | 254030.957 |
| 174 | Sambaíba | 0.223 | 0.726 | 0.030 | 83033.584 |
| 170 | Santo Antônio dos Lopes | 0.210 | 0.006 | 0.900 | 207732.341 |
| 213 | Alto Parnaíba | 0.075 | 0.626 | 0.059 | 50964.877 |
Two kinds of municipality top the list. Tasso Fragoso, Sambaíba and Alto Parnaíba are the agribusiness pattern the model was built around. Godofredo Viana and Santo Antônio dos Lopes have almost no agriculture and still very high GDP per capita: industry accounts for 82% and 90% of their GDP (mining and energy projects show up this way in municipal GDP accounts). A model with only agro_share has no way to explain them, and a spatial term is exactly where that unexplained, spatially clustered signal ends up.
ind_share, industry's share of GDP, comes from the same IBGE table and was loaded with the data at the start of this section. Adding it changes the picture:
for cols in (["agro_share"], ["agro_share", "ind_share"]):
ols_check = spreg.OLS(maranhao[["log_pib_pc"]].values, maranhao[cols].values, name_x=cols)
sdm_check = spreg.ML_Lag(
maranhao[["log_pib_pc"]].values, maranhao[cols].values, w=w_ma,
slx_lags=1, name_y="log_pib_pc", name_x=cols,
)
print(f"Explanatory variables: {cols} OLS R² = {ols_check.r2:.2f}")
for name, b, (z, p) in zip(sdm_check.name_x, sdm_check.betas.flatten(), sdm_check.z_stat):
print(f" {name:16s} {b:7.3f} p = {p:.3f}")
ML_Lag Explanatory variables: ['agro_share'] OLS R² = 0.11 CONSTANT 4.229 p = 0.000 agro_share -0.911 p = 0.008 W_agro_share 2.333 p = 0.000 W_log_pib_pc 0.520 p = 0.000
ML_Lag Explanatory variables: ['agro_share', 'ind_share'] OLS R² = 0.69 CONSTANT 5.374 p = 0.000 agro_share 0.876 p = 0.000 ind_share 3.231 p = 0.000 W_agro_share 1.281 p = 0.001 W_ind_share -0.784 p = 0.089 W_log_pib_pc 0.362 p = 0.000
With industry in the model, the variance explained by OLS jumps from about 11% to about 69%, the municipality's own agricultural share turns positive again (about +0.88), and the neighbors' share shrinks to roughly half its earlier value. The striking "own negative, neighbors positive" pattern of the previous section was, to a large extent, the missing industry variable showing up where the model could put it.
That is the real lesson of this section, and it applies to every spatial model in this book: a lag or Durbin term absorbs anything spatially clustered that the explanatory variables leave out. The diagnostics told us, correctly, that the one-variable model was missing something with a spatial shape; they could not tell us what. Answering that took a look at the influential observations and some knowledge of the place. Even the two-variable model is still a teaching example — one year, a cross-section, no theory of regional growth — so its coefficients are best read as "what this specification says", not as effects.
Exercises¶
- Confirm the flip. Run
spreg.OLS(..., spat_diag=True)ony_errinstead ofy_lag. DoesRobust LM (error)now outrankRobust LM (lag), as the error-process construction would predict? - A weaker process. Rebuild
y_lagwithrho_true = 0.2instead of0.6. DoesML_Lagstill recover a value close to the truerho? Does the Moran's I on the OLS residuals stay as clearly significant as it was atrho_true = 0.6? - No spatial process at all. Set
rho_true = 0(equivalently, just userhsdirectly asy, skipping thenp.linalg.solvestep). Confirm that Moran's I on the OLS residuals is no longer significant, and thatML_Lag's estimatedrholands close to zero — a spatial regression model shouldn't invent a spatial effect that was never there. - Rook instead of Queen. Refit
ML_Lagony_lagusingRookweights instead ofQueen(both fromlibpysal.weights, following Chapter 10's pattern). Does the estimatedrhochange much? What does that suggest about how sensitive spatial regression is to the exact weights specification, compared to Moran's I in Chapter 12? - Impacts, by hand. Using
ml_lagfrom the synthetic example, build the full matrixM = np.linalg.inv(I_n - ml_lag.rho * W_full) * b1. Check that the mean of its diagonal equalsdirectand the mean of its row sums equalstotal. Then rebuildy_lagwithrho_true = 0.2: what share of the total effect is indirect now? - Drop the influential municipalities. Refit the one-variable OLS and Durbin models without the five municipalities in
flagged. How much do the coefficients move? Then do the same withind_shareincluded. Which specification is more stable, and why?
Summary¶
Key concepts introduced¶
- Why OLS's independence-of-errors assumption is routinely violated by geographic data, and why that matters beyond a technicality
- Reusing Chapter 12's Moran's I as a diagnostic on regression residuals, not just on a raw variable
- The spatial lag model (
spreg.ML_Lag): the outcome itself depends on neighbors' outcomes, so a coefficient is read through direct, indirect and total impacts (LeSage and Pace), never as a raw slope - The spatial error model (
spreg.ML_Error): unmeasured spatially clustered factors, leaving coefficients roughly intact but correcting their standard errors - Robust Lagrange Multiplier tests as a principled way to choose between the two from data alone — including the branch where both stay significant and a Spatial Durbin model (
slx_lags=1) is the candidate - A teaching example with 2021 GDP per capita in Maranhão: with only
agro_share, the diagnostics pointed to a Spatial Durbin model and a striking "neighbors matter, own share does not" pattern — which largely disappeared once influential municipalities were inspected and industry's share was added. Spatial terms absorb whatever clustered signal the explanatory variables leave out; checking influence and specification comes before any geographic story
Chapter 14 stays with regression's toolkit a little longer, but turns from explaining a continuous outcome to grouping observations by similarity — clustering and regionalization. Chapter 15 then turns from statistics back to communication — the mapping and visualization techniques for presenting exactly the kind of pattern this chapter and Chapter 12 just learned to detect, measure, and model.
Further Reading¶
- Sergio J. Rey, Dani Arribas-Bel, and Levi J. Wolf, Geographic Data Science with Python (CRC Press) — freely available at https://geographicdata.science/book/; see their chapter on spatial regression for the same lag/error distinction applied to a real Airbnb dataset
- Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer Academic — the foundational text behind the lag/error distinction used throughout this chapter
- LeSage, J., & Pace, R. K. (2009). Introduction to Spatial Econometrics. CRC Press — the source of the direct/indirect/total impacts used to interpret lag and Durbin models
- Anselin, L., Bera, A. K., Florax, R., & Yoon, M. J. (1996). "Simple diagnostic tests for spatial dependence." Regional Science and Urban Economics, 26(1), 77-104 — the robust LM tests and the decision rule above
- spreg documentation, Spatial Regression: https://pysal.org/spreg/
- PySAL, Spatial Regression notebooks: https://pysal.org/spreg/generated/spreg.ML_Lag.html
- IBGE, Produto Interno Bruto dos Municípios — the GDP and sectoral value-added source behind this chapter's Maranhão example: https://www.ibge.gov.br/estatisticas/economicas/contas-nacionais/9088-produto-interno-bruto-dos-municipios.html