Getting Started with bayespecon¶

This notebook walks through a complete spatial econometric workflow:

1. Simulate spatial data from a known DGP

2. Fit a baseline OLS model

3. Run Bayesian LM diagnostics to test for spatial dependence

4. Use the decision tree to select a better model

5. Fit the recommended SAR model

6. Compute direct, indirect, and total spatial effects

7. Peek at the underlying PyMC model for customization

import libpysal

from bayespecon.dgp import simulate_sar

gdf = simulate_sar(n=20, beta=[1, 0.4, 2.5], rho=0.6, create_gdf=True)

gdf

gdf.plot("X_1").set_title("X1")

Text(0.5, 1.0, 'X1')

gdf.plot("y").set_title("y")

Text(0.5, 1.0, 'y')

1. Create a Spatial Weights Matrix¶

Spatial econometric models require a spatial weights matrix \(W\) that encodes the spatial relationships between observations. Here we use Rook contiguity (regions that share a boundary or vertex are neighbors).

G = libpysal.graph.Graph.build_contiguity(gdf, rook=False).transform("r")

2. Fit a Baseline OLS Model¶

We start with a non-spatial OLS regression to establish a baseline.

# remove the intercept since we have one (-1)
form = "y ~ -1 + X_0 + X_1 + X_2"

from bayespecon import OLS

ols = OLS(formula=form, W=G, data=gdf)
ols.fit(draws=1000, tune=500, chains=2, random_seed=42)
ols.summary()

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma2]
Sampling 2 chains for 500 tune and 1_000 draw iterations (1_000 + 2_000 draws total) took 6 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" 
for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Notice all the estimated beta parameters are biased upward.

Use arviz to inspect fit¶

the fit method attaches an inference_data object to each model that can be used directly with arviz functions

import arviz as az

az.plot_forest(ols.inference_data)

array([<Axes: title={'center': '94.0% HDI'}>], dtype=object)

az.plot_trace(ols.inference_data)

array([[<Axes: title={'center': 'beta'}>,
        <Axes: title={'center': 'beta'}>],
       [<Axes: title={'center': 'sigma2'}>,
        <Axes: title={'center': 'sigma2'}>],
       [<Axes: title={'center': 'sigma'}>,
        <Axes: title={'center': 'sigma'}>]], dtype=object)

3. Run Bayesian LM Diagnostics¶

The spatial_diagnostics() method runs a battery of Bayesian LM tests (Doğan, Taşpınar & Bera 2021) that check whether spatial dependence is present in the residuals.

• LM-Lag: Tests \(H_0: \rho = 0\) (no spatial lag on \(y\))

• LM-Error: Tests \(H_0: \lambda = 0\) (no spatial error correlation)

• LM-SDM-Joint: Joint test for \(\rho = 0\) and \(\gamma = 0\)

• Robust-LM-Lag / Robust-LM-Error: Neyman-orthogonal robust versions

diag = ols.spatial_diagnostics()
diag

4. Model Selection Decision Tree¶

The spatial_diagnostics_decision() method walks a Koley & Bera decision tree using the Bayesian p-values above and recommends the next model to fit.

decision = ols.spatial_diagnostics_decision(alpha=0.05, format="ascii")
print(decision)

LM-Lag *  (p=0.0000, alpha=0.05)
├── <sig> LM-Error *  (p=0.0000, alpha=0.05)
│   ├── <sig> Robust-LM-Lag *  (p=0.0000, alpha=0.05)
│   │   ├── <sig> Robust-LM-Error *  (p=0.0001, alpha=0.05)
│   │   │   ├── <sig> Robust-LM-Lag p <= Robust-LM-Error p *
│   │   │   │   ├── [SAR] * ← SELECTED
│   │   │   │   └── [SEM]
│   │   │   └── [SAR]
│   │   └── <not sig> Robust-LM-Error
│   │       ├── [SEM]
│   │       └── <not sig> LM-Lag p <= LM-Error p
│   │           ├── [SAR]
│   │           └── [SEM]
│   └── [SAR]
└── <not sig> LM-Error
    ├── [SEM]
    └── [OLS]

ols.spatial_diagnostics_decision(alpha=0.05)

5. Fit the Recommended SAR Model¶

If the diagnostics suggest spatial dependence, we fit a Spatial Autoregressive (SAR) model:

The likelihood includes the Jacobian \(\log|I - \rho W|\) so that posterior inference on \(\rho\) is exact.

from bayespecon import SAR

sar = SAR(formula=form, W=G, data=gdf)
sar.fit(draws=2000, tune=100, chains=4, random_seed=42)
sar.summary()

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/rich/live.py:260: UserWarning: install "ipywidgets" 
for Jupyter support
  warnings.warn('install "ipywidgets" for Jupyter support')

Sampling 4 chains for 100 tune and 2000 draw iterations, 4 x 2,100 draws total took 5s (1758 draws/s)

az.plot_trace(sar.inference_data)

array([[<Axes: title={'center': 'rho'}>, <Axes: title={'center': 'rho'}>],
       [<Axes: title={'center': 'sigma'}>,
        <Axes: title={'center': 'sigma'}>],
       [<Axes: title={'center': 'sigma2'}>,
        <Axes: title={'center': 'sigma2'}>],
       [<Axes: title={'center': 'beta'}>,
        <Axes: title={'center': 'beta'}>]], dtype=object)

6. Compute Spatial Effects¶

For the SAR model, a change in \(X\) propagates through the spatial multiplier \((I - \rho W)^{-1}\). The effects decompose into:

• Direct: Effect of \(X_i\) on \(y_i\) (own-region)

• Indirect: Effect of \(X_j\) on \(y_i\) for \(j \neq i\) (spillover)

• Total: Direct + Indirect

effects = sar.spatial_effects()
effects