Bayesian Spatial Models: A Pedagogical Walkthrough¶

This notebook demonstrates all five models currently implemented in bayespecon:

1. SLX

2. SAR

3. SEM

4. SDM

5. SDEM

Each section explains the model idea, fits the model, and inspects posterior summaries and spatial effects.

Conceptual Roadmap¶

Let \(W\) denote the spatial weights matrix and \(WX\) the spatial lag of covariates.

• SLX: \(y = X\beta + WX\theta + \epsilon\)

• SAR: \(y = \rho Wy + X\beta + \epsilon\)

• SEM: \(y = X\beta + u\), \(u = \lambda Wu + \epsilon\)

• SDM: \(y = \rho Wy + X\beta + WX\theta + \epsilon\)

• SDEM: \(y = X\beta + WX\theta + u\), \(u = \lambda Wu + \epsilon\)

In bayespecon, all models accept either:

• formula mode: model_class(formula=..., data=..., W=...), or

• matrix mode: model_class(y=..., X=..., W=...).

import arviz as az
import geodatasets
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
from libpysal.graph import Graph

from bayespecon import OLS, SAR, SDEM, SDM, SEM, SLX
from bayespecon.diagnostics.bayesfactor import bayes_factor_compare_models

az.style.use("arviz-white")

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

# Load sample spatial dataset
gdf = gpd.read_file(geodatasets.data.geoda.airbnb.url)
xcols = ["poverty", "rev_rating", "num_spots", "crowded"]
ycol = "price_pp"
gdf = gdf.dropna(subset=xcols + [ycol]).copy()

# Build contiguity graph (queen)
W = Graph.build_contiguity(gdf, rook=False).transform("r")

print(f"Observations: {len(gdf)}")
print(f"Predictors: {xcols}")

Observations: 67
Predictors: ['poverty', 'rev_rating', 'num_spots', 'crowded']

Helper Functions

To keep each model section focused, we use utility functions to:

• fit with consistent MCMC settings,

• print posterior summaries,

• tabulate direct/indirect/total effects.

For a quick pedagogical run, we use small draw counts. Increase these for real inference.

def fit_and_report(model_cls, formula, data, W, draws=400, tune=400, chains=4, seed=42):
    """Fit a spatial model and return (model, summary, effects_df).

    Uses minimal MCMC settings for pedagogical demonstration.
    Increase draws/tune for real analyses.
    """
    model = model_cls(formula=formula, data=data, W=W, logdet_method="eigenvalue")
    model.fit(
        draws=draws,
        tune=tune,
        chains=chains,
        target_accept=0.9,
        random_seed=seed,
        progressbar=False,
        idata_kwargs={"log_likelihood": True},
    )
    summary = model.summary(round_to=3)
    effects_df = pd.DataFrame(model.spatial_effects())
    return model, summary, effects_df

Convergence Diagnostics Helper

For each fitted model, we will inspect:

• r_hat (target close to 1.00)

• effective sample sizes (ess_bulk, ess_tail)

• trace plots for key parameters.

def diagnostics_table(idata, var_names):
    """Show key MCMC diagnostics for the given parameters."""
    cols = ["mean", "sd", "ess_bulk", "ess_tail", "r_hat"]
    return az.summary(idata, var_names=var_names, round_to=3)[cols]


def show_trace(idata, var_names, title):
    """Plot trace plots for the given parameters."""
    az.plot_trace(idata, var_names=var_names)
    plt.suptitle(title, y=1.02)
    plt.tight_layout()
    plt.show()

Models¶

0) OLS: What could go wrong?¶

Model:

ols = OLS(
    formula="price_pp ~ poverty + rev_rating + num_spots + crowded",
    data=gdf,
    W=W,
)
olsfit = ols.fit(
    draws=400,
    tune=400,
    chains=4,
    target_accept=0.9,
    random_seed=42,
    progressbar=False,
    idata_kwargs={"log_likelihood": True},
)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 2 jobs)
NUTS: [beta, sigma2]
Sampling 4 chains for 400 tune and 400 draw iterations (1_600 + 1_600 draws total) took 14 seconds.

ols.summary()

1) SLX: Spatially Lagged Covariates Only¶

Model:

Interpretation:

• beta captures local covariate effects.

• theta captures neighbor-covariate spillovers.

• No spatial lag on \(y\), so no autoregressive propagation through outcomes.

slx, summary_slx, effects_slx = fit_and_report(
    SLX,
    formula="price_pp ~ poverty + rev_rating + num_spots + crowded",
    data=gdf,
    W=W,
)
display(summary_slx)
display(effects_slx)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 2 jobs)
NUTS: [beta, sigma2]
Sampling 4 chains for 400 tune and 400 draw iterations (1_600 + 1_600 draws total) took 17 seconds.

az.plot_forest(slx.inference_data)
display(diagnostics_table(slx.inference_data, ["beta", "sigma"]))
show_trace(slx.inference_data, ["sigma"], "SLX Trace: sigma")

	mean	sd	ess_bulk	ess_tail	r_hat
beta[Intercept]	90.613	77.431	732.657	918.725	1.005
beta[poverty]	-0.659	0.472	951.327	892.080	1.003
beta[rev_rating]	0.592	0.773	1069.421	1194.587	1.000
beta[num_spots]	0.021	0.032	1017.610	1030.384	1.000
beta[crowded]	-0.953	1.188	993.532	933.014	1.005
beta[W*poverty]	1.347	0.634	1023.534	876.226	1.005
beta[W*rev_rating]	-0.971	0.967	944.270	951.677	1.004
beta[W*num_spots]	0.196	0.047	1226.908	870.559	1.001
beta[W*crowded]	-1.797	1.628	1007.256	1074.373	1.004
sigma	24.443	2.155	1197.693	1052.758	1.002

/tmp/ipykernel_7302/3581438262.py:11: UserWarning: The figure layout has changed to tight
  plt.tight_layout()

2) SAR: Spatial Lag of Outcome¶

Model:

Interpretation:

• \(\rho\) controls feedback from neighboring outcomes.

• Effects are amplified through \((I - \rho W)^{-1}\), so direct and indirect impacts differ from raw \(\beta\).

sar, summary_sar, effects_sar = fit_and_report(
    SAR,
    formula="price_pp ~ poverty + rev_rating + num_spots + crowded",
    data=gdf,
    W=W,
)
summary_sar

effects_sar

az.plot_forest(sar.inference_data)
display(diagnostics_table(sar.inference_data, ["rho", "beta", "sigma"]))
show_trace(sar.inference_data, ["rho", "sigma"], "SAR Trace: rho, sigma")

	mean	sd	ess_bulk	ess_tail	r_hat
rho	0.421	0.145	1590.247	1184.079	1.000
beta[Intercept]	14.767	61.025	1414.628	1400.482	1.003
beta[poverty]	0.023	0.316	1532.303	1457.932	0.999
beta[rev_rating]	0.310	0.646	1363.433	1486.916	1.003
beta[num_spots]	0.081	0.026	1652.805	1258.488	0.999
beta[crowded]	-1.699	0.902	1536.060	1394.467	1.000
sigma	25.386	2.359	1464.873	1571.805	0.999

/tmp/ipykernel_7302/3581438262.py:11: UserWarning: The figure layout has changed to tight
  plt.tight_layout()

3) SEM: Spatial Correlation in Errors¶

Model:

Interpretation:

• Spatial dependence is moved to latent shocks rather than outcome feedback.

• Useful when omitted spatial factors induce correlated residuals.

sem, summary_sem, effects_sem = fit_and_report(
    SEM,
    formula="price_pp ~ poverty + rev_rating + num_spots + crowded",
    data=gdf,
    W=W,
)
display(summary_sem)
display(effects_sem)

az.plot_forest(sem.inference_data)
display(diagnostics_table(sem.inference_data, ["lam", "beta", "sigma"]))
show_trace(sem.inference_data, ["lam", "sigma"], "SEM Trace: lam, sigma")

	mean	sd	ess_bulk	ess_tail	r_hat
lam	0.490	0.184	1508.653	918.977	1.002
beta[Intercept]	19.475	58.940	1451.461	1156.616	1.000
beta[poverty]	-0.165	0.440	1521.688	1706.191	1.002
beta[rev_rating]	0.630	0.628	1472.748	1257.599	1.001
beta[num_spots]	0.073	0.036	1375.208	1377.750	1.001
beta[crowded]	-1.795	1.146	1560.680	1598.345	0.999
sigma	25.856	2.362	1479.718	1570.485	1.001

/tmp/ipykernel_7302/3581438262.py:11: UserWarning: The figure layout has changed to tight
  plt.tight_layout()

4) SDM: SAR + SLX¶

Model:

Interpretation:

• Includes both outcome feedback and neighbor-covariate channels.

• Often used as a flexible nesting model for SAR/SLX.

sdm, summary_sdm, effects_sdm = fit_and_report(
    SDM,
    formula="price_pp ~ poverty + rev_rating + num_spots + crowded",
    data=gdf,
    W=W,
)
display(summary_sdm)
display(effects_sdm)

az.plot_forest(sdm.inference_data)
display(diagnostics_table(sdm.inference_data, ["rho", "beta", "sigma"]))
show_trace(sdm.inference_data, ["rho", "sigma"], "SDM Trace: rho, sigma")

	mean	sd	ess_bulk	ess_tail	r_hat
rho	0.223	0.171	1721.619	1299.748	1.000
beta[Intercept]	87.065	78.479	1647.108	1361.040	0.999
beta[poverty]	-0.643	0.469	1777.358	1611.313	1.000
beta[rev_rating]	0.688	0.770	1474.167	1515.338	1.000
beta[num_spots]	0.019	0.031	1849.856	1659.838	1.000
beta[crowded]	-0.998	1.134	1640.344	1335.392	1.003
beta[W*poverty]	1.173	0.632	1609.890	1571.707	1.000
beta[W*rev_rating]	-1.168	0.976	1791.192	1595.063	1.000
beta[W*num_spots]	0.162	0.052	1706.260	1534.498	0.999
beta[W*crowded]	-1.239	1.630	1733.541	1568.763	1.002
sigma	24.095	2.253	1229.013	1530.625	1.001

/tmp/ipykernel_7302/3581438262.py:11: UserWarning: The figure layout has changed to tight
  plt.tight_layout()

5) SDEM: SLX + Spatial Error¶

Model:

Interpretation:

• Neighbor covariates matter directly (through \(WX\)),

• and unmodeled shocks are spatially correlated (through \(u\)).

sdem, summary_sdem, effects_sdem = fit_and_report(
    SDEM,
    formula="price_pp ~ poverty + rev_rating + num_spots + crowded",
    data=gdf,
    W=W,
)
display(summary_sdem)
display(effects_sdem)

az.plot_forest(sdem.inference_data)
display(diagnostics_table(sdem.inference_data, ["lam", "beta", "sigma"]))
show_trace(sdem.inference_data, ["lam", "sigma"], "SDEM Trace: lam, sigma")

	mean	sd	ess_bulk	ess_tail	r_hat
lam	0.332	0.189	1539.489	1235.231	1.004
beta[Intercept]	89.945	79.819	1581.511	1350.080	1.001
beta[poverty]	-0.559	0.437	1496.614	1618.919	1.001
beta[rev_rating]	0.759	0.707	1559.760	1497.604	0.998
beta[num_spots]	0.029	0.031	1543.692	1616.363	1.000
beta[crowded]	-1.035	1.123	1448.302	1425.393	1.001
beta[W*poverty]	1.137	0.643	1641.260	1425.588	1.000
beta[W*rev_rating]	-1.094	0.889	1612.120	1261.345	1.007
beta[W*num_spots]	0.180	0.052	1547.216	1573.058	1.002
beta[W*crowded]	-1.905	1.865	1585.542	1518.769	1.002
sigma	24.203	2.263	1080.165	1484.121	1.001

/tmp/ipykernel_7302/3581438262.py:11: UserWarning: The figure layout has changed to tight
  plt.tight_layout()

MCMC Sampling Adequacy¶

Bayesian estimation of spatial models has a well-known efficiency pitfall: the spatial dependence parameter \(\rho\) (or \(\lambda\)) often mixes slowly, so a chain that looks converged in its point estimate can still produce posterior credible intervals that are 10–12 % too narrow because the sampler has not visited the tails enough times [Wolf et al., 2018]. The spatial_mcmc_diagnostic helper checks effective sample size, sampler yield, \(\hat{R}\), and HPDI stability for the spatial parameter and warns when any threshold is violated.

from bayespecon import spatial_mcmc_diagnostic

# Run on the SDM (most parameters; both rho and impact summaries depend on it)
report = spatial_mcmc_diagnostic(sdm, emit_warnings=False)
report.to_frame()

sdm.spatial_diagnostics()

ols.spatial_diagnostics()

sar.spatial_diagnostics_decision()

sem.spatial_diagnostics()

report

SpatialMCMCReport(parameters=['rho'], ess_bulk={'rho': 1721.618975952847}, ess_tail={'rho': 1299.7476399296434}, r_hat={'rho': 1.0003077944547376}, mcse_mean={'rho': 0.004136803152868457}, nominal_size=1600, yield_pct={'rho': 107.60118599705292}, hpdi_drift_pct={'rho': 0.07560087735308905}, warnings_triggered=[], adequate=True, adequate_by_param={'rho': True})

Model Comparison¶

# Collect all fitted models for comparison
model_dict = {
    "OLS": ols,
    "SLX": slx,
    "SAR": sar,
    "SEM": sem,
    "SDM": sdm,
    "SDEM": sdem,
}
idata_dict = {name: m.inference_data for name, m in model_dict.items()}

Bayes Factor Model Comparison¶

Bayes factors provide an alternative to information criteria for comparing competing models. While WAIC and LOO assess predictive performance, Bayes factors compare marginal likelihoods — the probability of the data under each model, integrated over all parameter values weighted by the prior:

This makes Bayes factors sensitive to the prior: models with diffuse priors on unnecessary parameters are penalized more heavily, because the marginal likelihood averages over all possible parameter values weighted by the prior. In spatial models, this means that models with many WX coefficients under wide priors (e.g., SLX, SDM, SDEM) can receive substantially lower marginal likelihoods than more parsimonious alternatives (e.g., OLS, SAR, SEM) — even if their log-likelihoods are similar.

Interpreting Bayes factors (Kass & Raftery, 1995):

Method. We use bridge sampling (Meng & Wong, 1996) to estimate each model’s log marginal likelihood, following the R bridgesampling package (Gronau et al., 2020). The implementation uses ESS weighting, two-phase convergence, and MCSE diagnostics. For reliable estimates, 40,000+ posterior samples are recommended.

Important caveat. Bayes factors can be very sensitive to prior specification. Models with diffuse priors on many parameters will tend to have lower marginal likelihoods. This is by design — it is the Bayesian Occam’s razor at work — but it means that Bayes factors and information criteria may disagree when priors are uninformative.

bayes_factor_compare_models(model_dict, method="bridge", log=True).round(3)

/tmp/ipykernel_7302/497883469.py:1: UserWarning: Bridge sampling with 1600 posterior samples for 'OLS' may yield imprecise marginal-likelihood estimates. A conservative rule of thumb is 40,000+ samples (Gronau, Singmann, & Wagenmakers, 2017).
  bayes_factor_compare_models(model_dict, method="bridge", log=True).round(3)

bayes_factor_compare_models(model_dict, method="bic", log=True).round(3)

Bridge Sampling vs. BIC¶

• BIC: SAR > SLX > OLS > SEM (SAR wins, SLX is moderate, SEM is worst)

• Bridge: SAR > SEM > OLS (SAR wins, SEM is moderate)

The two methods can produce qualitatively different model rankings. This is expected and reflects a fundamental difference in what they assume about priors:

• BIC approximates \(\log(ML) \approx \hat\ell_{\max} - \frac{k}{2}\log(n)\), which assumes unit-information priors (priors containing as much information as a single observation). The penalty per parameter is fixed at \(\frac{1}{2}\log(n) \approx 2.1\) for \(n = 77\).

• Bridge sampling integrates over the actual priors in the model. When priors are wide (e.g., Normal(0, 100) on WX coefficients), the marginal likelihood penalizes each such parameter by roughly \(\log(\sigma_{\text{prior}} / \sigma_{\text{post}})\), which can be 5–10× larger than the BIC penalty.

This is why models with many WX terms (SLX, SDM, SDEM) may look reasonable under BIC but receive extreme Bayes factors under bridge sampling: the wide priors on the WX coefficients are “wasted” — they spread probability mass over implausible parameter values, reducing the marginal likelihood. This is Bayesian Occam’s razor at work, and bridge sampling is generally more trustworthy because it accounts for the actual prior specification.

Model Comparison: WAIC and LOO¶

We compare fitted models using information criteria from ArviZ.

# WAIC and LOO comparison
for ic in ("waic", "loo"):
    cmp = az.compare(idata_dict, ic=ic)
    az.plot_compare(cmp)

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:1652: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. 
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:1652: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. 
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:1652: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. 
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:1652: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. 
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:1652: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. 
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:1652: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. 
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.69 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.69 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.69 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.69 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/arviz/stats/stats.py:782: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.69 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

Bayesian LM Specification Tests¶

The bayespecon.diagnostics module provides Bayesian Lagrange-multiplier tests that operate on posterior draws rather than point estimates. These replace the frequentist LM tests that were previously available.

Below we run the Bayesian LM tests from the fitted OLS model:

• bayesian_lm_lag_test — tests \(H_0: \rho = 0\) (no spatial lag)

• bayesian_lm_error_test — tests \(H_0: \lambda = 0\) (no spatial error)

• bayesian_lm_wx_test — tests \(H_0: \gamma = 0\) (no WX, from SAR null)

• bayesian_lm_sdm_joint_test — joint test for SDM (\(H_0: \rho = 0\) and \(\gamma = 0\))

• bayesian_lm_slx_error_joint_test — joint test for SDEM (\(H_0: \lambda = 0\) and \(\gamma = 0\))

• bayesian_robust_lm_lag_test / bayesian_robust_lm_error_test — Anselin–Florax robust pair (SAR vs SEM)

• bayesian_robust_lm_lag_sdm_test / bayesian_robust_lm_wx_test / bayesian_robust_lm_error_sdem_test — Neyman-orthogonal robust tests

• bayesian_lm_error_sdm_test / bayesian_lm_lag_sdem_test — SDM/SDEM-aware tests using correctly filtered residuals

In practice, prefer the method API: model.spatial_diagnostics() runs the tests wired to that model class, and model.spatial_diagnostics_decision() returns a recommended specification.

Each returns a BayesianLMTestResult with posterior summary statistics and a Bayesian p-value.

# `spatial_diagnostics()` runs all tests wired to a model class
# and returns a tidy DataFrame.
ols.spatial_diagnostics()

Compare Spatial Parameters Across Models¶

This cell compares posterior means/intervals for the spatial scalar where present:

• rho in SAR/SDM

• lam in SEM/SDEM

# Compare spatial parameters across models
spatial_rows = []
for name, model, var in [
    ("SAR", sar, "rho"),
    ("SEM", sem, "lam"),
    ("SDM", sdm, "rho"),
    ("SDEM", sdem, "lam"),
]:
    if var in model.inference_data.posterior:
        summary = az.summary(model.inference_data, var_names=[var], round_to=3)
        summary.insert(0, "model", name)
        spatial_rows.append(summary)

pd.concat(spatial_rows)

Matrix-Mode API Example (Optional)¶

The package also supports direct matrix inputs. This is useful if your design matrix is already engineered.

y = gdf[ycol]
X = gdf[xcols]

sar_matrix = SAR(y=y, X=X, W=W)
sar_matrix.fit(draws=200, tune=200, chains=2, random_seed=7, progressbar=False)
display(sar_matrix.summary(round_to=3))