Getting Started with neighbayes¶

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 neighbayes.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 neighbayes.models 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 5 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

/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

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.0000, 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 neighbayes.models 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/neighbayes/_logdet/_jax.py:188: ComplexWarning: Casting complex values to real discards the imaginary part
  W_arr = np.asarray(W, dtype=np.float64)

Gibbs sampling (sar): 4 chains for 100 tune and 2,000 draw iterations (4 x 2,100 = 8,400 draws total)

/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 took 4s (1,994 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