Bayesian spatial models: a walkthrough¶

By the end of this lesson you will have fitted six models to one dataset, watched OLS fail a diagnostic that the spatial models pass, and used the package’s own tests to arrive at the specification the data support.

You do not need to know the spatial econometrics beforehand — each model is introduced when the previous one runs out of road. Work through the cells in order; every step produces something you can look at.

Every model here accepts a formula and a weights graph:

SAR(formula="y ~ x1 + x2", data=gdf, W=W)

The equations behind each one are in Supported Models if you want them; you will not need them to follow along.

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 neighbayes.diagnostics.bayesfactor import bayes_factor_compare_models
from neighbayes.models import OLS, SAR, SDEM, SDM, SEM, SLX

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,
        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()

The models, one at a time¶

Start with OLS, and see what it misses¶

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,
    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 13 seconds.

ols.summary()

Add the neighbours’ covariates (SLX)¶

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 16 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]	89.774	76.776	739.464	731.272	1.007
beta[poverty]	-0.609	0.456	1035.774	737.892	1.001
beta[rev_rating]	0.608	0.752	1131.814	1101.190	1.001
beta[num_spots]	0.023	0.033	1038.612	1045.020	1.002
beta[crowded]	-1.029	1.214	894.497	824.635	1.004
beta[W*poverty]	1.287	0.629	1145.391	820.396	1.000
beta[W*rev_rating]	-0.975	0.941	904.702	883.591	1.009
beta[W*num_spots]	0.194	0.049	814.179	847.357	1.008
beta[W*crowded]	-1.684	1.615	1015.970	1055.031	1.004
sigma	24.433	2.242	1382.910	992.499	1.002

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

Let the outcome itself spill over (SAR)¶

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

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/neighbayes/_logdet/_jax.py:188: ComplexWarning: Casting complex values to real discards the imaginary part
  W_arr = np.asarray(W, dtype=np.float64)

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.430	0.135	1551.439	1317.324	1.001
beta[Intercept]	16.073	60.068	1500.723	1618.080	1.001
beta[poverty]	0.029	0.316	1614.876	1574.109	1.001
beta[rev_rating]	0.288	0.641	1499.921	1475.097	1.001
beta[num_spots]	0.079	0.025	1544.044	1486.345	1.002
beta[crowded]	-1.712	0.902	1598.993	1516.523	1.000
sigma	25.429	2.211	1449.005	1390.197	1.001

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

Treat the clustering as nuisance instead (SEM)¶

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)

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/neighbayes/_logdet/_jax.py:188: ComplexWarning: Casting complex values to real discards the imaginary part
  W_arr = np.asarray(W, dtype=np.float64)

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.489	0.176	1525.126	1343.762	1.001
beta[Intercept]	20.974	58.338	1522.358	1611.833	1.001
beta[poverty]	-0.159	0.441	1532.719	1380.956	1.000
beta[rev_rating]	0.611	0.630	1515.520	1589.968	1.002
beta[num_spots]	0.073	0.034	1502.598	1488.721	1.001
beta[crowded]	-1.831	1.152	1606.060	1251.815	1.000
sigma	26.087	2.296	1452.437	1358.397	1.000

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

Combine both channels (SDM)¶

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)

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/neighbayes/_logdet/_jax.py:188: ComplexWarning: Casting complex values to real discards the imaginary part
  W_arr = np.asarray(W, dtype=np.float64)

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.229	0.163	1585.399	1249.928	1.001
beta[Intercept]	88.353	77.536	1392.868	1421.182	1.002
beta[poverty]	-0.634	0.440	1611.936	1598.968	1.001
beta[rev_rating]	0.653	0.758	1576.420	1659.903	1.000
beta[num_spots]	0.016	0.031	1666.324	1565.305	1.001
beta[crowded]	-0.968	1.183	1466.797	1521.188	1.002
beta[W*poverty]	1.167	0.613	1504.523	1502.548	1.002
beta[W*rev_rating]	-1.149	0.953	1384.228	1450.976	1.001
beta[W*num_spots]	0.163	0.050	1657.056	1647.218	1.003
beta[W*crowded]	-1.292	1.679	1556.008	1467.845	1.002
sigma	24.098	2.138	1367.823	1534.073	1.001

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

The error-side combination (SDEM)¶

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)

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/neighbayes/_logdet/_jax.py:188: ComplexWarning: Casting complex values to real discards the imaginary part
  W_arr = np.asarray(W, dtype=np.float64)

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.325	0.187	1451.228	1162.854	1.003
beta[Intercept]	89.304	78.324	1385.310	1457.200	1.003
beta[poverty]	-0.551	0.424	1674.855	1615.808	1.000
beta[rev_rating]	0.728	0.688	1569.599	1572.163	1.001
beta[num_spots]	0.028	0.030	1563.386	1533.317	1.000
beta[crowded]	-1.052	1.149	1544.621	1461.216	1.001
beta[W*poverty]	1.140	0.658	1504.776	1502.020	1.002
beta[W*rev_rating]	-1.057	0.859	1430.407	1447.225	1.003
beta[W*num_spots]	0.180	0.049	1681.046	1574.454	1.000
beta[W*crowded]	-1.947	1.881	1507.133	1479.839	1.002
sigma	24.159	2.142	1359.933	1459.839	1.001

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

Check the chains before comparing anything¶

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 neighbayes.diagnostics 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': 1585.3993209654486}, ess_tail={'rho': 1249.92785492861}, r_hat={'rho': 1.000685177748173}, mcse_mean={'rho': 0.0041244539862753124}, nominal_size=1600, yield_pct={'rho': 99.08745756034054}, hpdi_drift_pct={'rho': 7.980914100963144}, warnings_triggered=["95% HPDI width for 'rho' drifts by 8.0% between the last third and the full chain (> 5%); the credible interval has not yet stabilized. Consider doubling `draws` and/or `tune`, or re-running with more chains. See Wolf, Anselin & Arribas-Bel (2018), Geographical Analysis 50:97-119."], adequate=False, adequate_by_param={'rho': False})

Compare what you have fitted¶

# 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_11205/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.

Let the tests pick the specification¶

The neighbayes.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()

See how the spatial parameter moves 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)

Where to go next¶

• How to run Bayesian LM specification tests — the tests used above, in full

• How to run spatial block cross-validation — an out-of-sample check on the specification you chose

• How to set priors

• Supported Models — every model and its equations