Bayesian Spatial Panel Models

This notebook demonstrates the panel-model classes in bayespecon:

  • Fixed effects / pooled panel models: OLSPanelFE, SARPanelFE, SEMPanelFE, SDMPanelFE, SDEMPanelFE

  • Random effects panel models: OLSPanelRE, SARPanelRE, SEMPanelRE

It also demonstrates Bayesian model comparison via WAIC, LOO-CV, and Bayes factors.

Panel Model Equations

Let observations be stacked by time, then unit (panel_g convention).

Fixed-effects / pooled classes:

  • OLS panel: \(y_{it} = x_{it}'\beta + a_i + \tau_t + \epsilon_{it}\)

  • SAR panel: \(y_{it} = \rho (Wy)_{it} + x_{it}'\beta + a_i + \tau_t + \epsilon_{it}\)

  • SEM panel: \(y_{it} = x_{it}'\beta + a_i + \tau_t + u_{it},\; u_{it}=\lambda (Wu)_{it}+\epsilon_{it}\)

  • SDM panel: \(y_{it} = \rho (Wy)_{it} + x_{it}'\beta + (Wx)_{it}'\theta + a_i + \tau_t + \epsilon_{it}\)

  • SDEM panel: \(y_{it} = x_{it}'\beta + (Wx)_{it}'\theta + a_i + \tau_t + u_{it},\; u_{it}=\lambda (Wu)_{it}+\epsilon_{it}\)

Random-effects classes:

  • OLS random effects: \(y_{it} = x_{it}'\beta + \alpha_i + \epsilon_{it},\; \alpha_i \sim N(0, \sigma_\alpha^2)\)

  • SAR random effects: \(y_{it} = \rho (Wy)_{it} + x_{it}'\beta + \alpha_i + \epsilon_{it},\; \alpha_i \sim N(0, \sigma_\alpha^2)\)

  • SEM random effects: \(y_{it} = x_{it}'\beta + \alpha_i + u_{it},\; u_{it}=\lambda (Wu)_{it}+\epsilon_{it},\; \alpha_i \sim N(0, \sigma_\alpha^2)\)

The model argument selects the pooled or fixed-effects specification for the FE classes:

  • 0 pooled

  • 1 unit FE

  • 2 time FE

  • 3 two-way FE

The RE classes are hierarchical models with unit-level random intercepts, so they internally use the pooled design and estimate sigma_alpha directly.

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from libpysal.graph import Graph

from bayespecon import (
    OLSPanelFE,
    OLSPanelRE,
    SARPanelFE,
    SARPanelRE,
    SDEMPanelFE,
    SDMPanelFE,
    SEMPanelFE,
    SEMPanelRE,
    dgp,
)

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

Build a Balanced Panel Dataset

For pedagogy, we generate synthetic panel data from the bayespecon.dgp module on a regular polygon grid and then assemble a formula-ready DataFrame.

# Generate base geometry via cross-sectional DGP, then simulate panel data via panel DGP.
rng = np.random.default_rng(2026)
xcols = ["poverty", "rev_rating", "num_spots", "crowded"]
ycol = "price_pp"

base_gdf = dgp.simulate_sar(n=8, seed=2026, create_gdf=True, geometry_type="polygon")
W = Graph.build_contiguity(base_gdf, rook=False).transform("r")
N = len(base_gdf)
T = 4
beta = np.array([1.5, -0.8, 0.6, 0.4, -0.5], dtype=float)

panel_sim = dgp.simulate_panel_sar_fe(
    N=N,
    T=T,
    rho=0.35,
    beta=beta,
    sigma=1.0,
    sigma_alpha=0.5,
    rng=rng,
    W=W,
)

panel = pd.DataFrame(
    {
        "unit": panel_sim["unit"],
        "time": panel_sim["time"],
        ycol: panel_sim["y"],
    }
)
for j, name in enumerate(xcols, start=1):
    panel[name] = panel_sim["X"][:, j]

Shared Helpers

def fit_panel_model(
    model_cls, formula, data, W, model=3, draws=300, tune=300, chains=2, seed=42
):
    """Fit a panel model and return (model, summary, effects_df).

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


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

Fit: OLSPanelFE

formula = "price_pp ~ poverty + rev_rating + num_spots + crowded"
ols_panel, summary_ols, effects_ols = fit_panel_model(
    OLSPanelFE, formula, panel, W, model=3
)
display(summary_ols)
display(effects_ols)
display(diagnostics_table(ols_panel.inference_data, ["beta", "sigma"]))
show_trace(ols_panel.inference_data, ["sigma"], "OLSPanelFE trace")

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma]
Sampling 2 chains for 300 tune and 300 draw iterations (600 + 600 draws total) took 5 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
/tmp/ipykernel_8142/2596509600.py:42: UserWarning: The figure layout has changed to tight
  plt.tight_layout()
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept 4469.980 937184.658 -1833361.622 1657213.656 30687.678 45792.231 919.547 425.073 1.013
poverty -0.780 0.064 -0.891 -0.653 0.002 0.002 1004.674 420.546 1.002
rev_rating 0.622 0.066 0.491 0.750 0.002 0.003 1222.699 454.771 1.003
num_spots 0.399 0.084 0.238 0.560 0.003 0.003 1016.108 446.547 0.998
crowded -0.484 0.067 -0.601 -0.358 0.002 0.003 1018.383 466.395 0.999
sigma 0.927 0.040 0.858 1.008 0.001 0.002 1034.072 444.583 1.004
direct direct_ci_lower direct_ci_upper direct_pvalue indirect indirect_ci_lower indirect_ci_upper indirect_pvalue total total_ci_lower total_ci_upper total_pvalue
variable
poverty -0.779711 -0.902205 -0.652772 0.0 0.0 0.0 0.0 0.0 -0.779711 -0.902205 -0.652772 0.0
rev_rating 0.621883 0.491038 0.767581 0.0 0.0 0.0 0.0 0.0 0.621883 0.491038 0.767581 0.0
num_spots 0.399134 0.237576 0.572683 0.0 0.0 0.0 0.0 0.0 0.399134 0.237576 0.572683 0.0
crowded -0.483840 -0.615611 -0.358625 0.0 0.0 0.0 0.0 0.0 -0.483840 -0.615611 -0.358625 0.0
mean sd ess_bulk ess_tail r_hat
beta[Intercept] 4469.980 937184.658 919.547 425.073 1.013
beta[poverty] -0.780 0.064 1004.674 420.546 1.002
beta[rev_rating] 0.622 0.066 1222.699 454.771 1.003
beta[num_spots] 0.399 0.084 1016.108 446.547 0.998
beta[crowded] -0.484 0.067 1018.383 466.395 0.999
sigma 0.927 0.040 1034.072 444.583 1.004
../_images/b0a06bd0751397a141bab346de04ef7439285ccdddf5880b13bf8e5a2e8c8cdf.png

Fit: SARPanelFE

sar_panel, summary_sar, effects_sar = fit_panel_model(
    SARPanelFE, formula, panel, W, model=3
)
display(summary_sar)
display(effects_sar)
display(diagnostics_table(sar_panel.inference_data, ["rho", "beta", "sigma"]))
show_trace(sar_panel.inference_data, ["rho", "sigma"], "SARPanelFE trace")

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [rho, beta, sigma]
Sampling 2 chains for 300 tune and 300 draw iterations (600 + 600 draws total) took 6 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
/tmp/ipykernel_8142/2596509600.py:42: UserWarning: The figure layout has changed to tight
  plt.tight_layout()
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept -21947.032 1071440.748 -2177831.912 1837778.029 30117.687 46279.693 1254.124 443.863 0.998
poverty -0.785 0.064 -0.906 -0.674 0.002 0.003 1129.311 420.399 1.003
rev_rating 0.621 0.066 0.505 0.741 0.002 0.003 855.824 431.533 1.012
num_spots 0.374 0.079 0.234 0.524 0.003 0.003 792.083 493.628 0.999
crowded -0.494 0.062 -0.611 -0.374 0.002 0.003 1049.105 419.372 1.001
rho 0.246 0.085 0.079 0.391 0.003 0.004 1049.293 336.693 1.001
sigma 0.908 0.041 0.834 0.985 0.001 0.002 1315.087 488.279 1.002
direct direct_ci_lower direct_ci_upper direct_pvalue indirect indirect_ci_lower indirect_ci_upper indirect_pvalue total total_ci_lower total_ci_upper total_pvalue
variable
poverty -0.794941 -0.921184 -0.672052 0.0 -0.260856 -0.532395 -0.070084 0.006667 -1.055797 -1.393542 -0.796507 0.0
rev_rating 0.629219 0.506198 0.755766 0.0 0.206738 0.053158 0.439318 0.006667 0.835957 0.617998 1.147002 0.0
num_spots 0.378158 0.217951 0.525261 0.0 0.123778 0.028644 0.261822 0.006667 0.501936 0.283112 0.760715 0.0
crowded -0.499755 -0.626853 -0.374754 0.0 -0.164171 -0.343405 -0.043669 0.006667 -0.663925 -0.908107 -0.451911 0.0
mean sd ess_bulk ess_tail r_hat
rho 0.246 0.085 1049.293 336.693 1.001
beta[Intercept] -21947.032 1071440.748 1254.124 443.863 0.998
beta[poverty] -0.785 0.064 1129.311 420.399 1.003
beta[rev_rating] 0.621 0.066 855.824 431.533 1.012
beta[num_spots] 0.374 0.079 792.083 493.628 0.999
beta[crowded] -0.494 0.062 1049.105 419.372 1.001
sigma 0.908 0.041 1315.087 488.279 1.002
../_images/6c2fc98829dee6d0dff921fe51997c64b583a2a8248b44eae5ab12899914bd42.png

Fit: SEMPanelFE

sem_panel, summary_sem, effects_sem = fit_panel_model(
    SEMPanelFE, formula, panel, W, model=3
)
display(summary_sem)
display(effects_sem)
display(diagnostics_table(sem_panel.inference_data, ["lam", "beta", "sigma"]))
show_trace(sem_panel.inference_data, ["lam", "sigma"], "SEMPanelFE trace")

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [lam, beta, sigma]
Sampling 2 chains for 300 tune and 300 draw iterations (600 + 600 draws total) took 6 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
/tmp/ipykernel_8142/2596509600.py:42: UserWarning: The figure layout has changed to tight
  plt.tight_layout()
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept -24246.255 1031031.388 -1974330.009 1939371.889 36558.371 44704.183 770.792 272.065 1.002
poverty -0.772 0.065 -0.898 -0.653 0.002 0.003 1362.188 517.233 0.997
rev_rating 0.600 0.064 0.499 0.728 0.002 0.003 911.255 512.937 1.012
num_spots 0.358 0.076 0.232 0.506 0.003 0.003 847.769 515.956 1.003
crowded -0.493 0.060 -0.607 -0.378 0.002 0.002 1368.274 579.247 1.001
lam 0.237 0.101 0.043 0.423 0.003 0.005 963.662 388.022 1.002
sigma 0.912 0.039 0.846 0.999 0.001 0.002 839.287 415.063 0.999
direct direct_ci_lower direct_ci_upper direct_pvalue indirect indirect_ci_lower indirect_ci_upper indirect_pvalue total total_ci_lower total_ci_upper total_pvalue
variable
poverty -0.771928 -0.898169 -0.644194 0.0 0.0 0.0 0.0 0.0 -0.771928 -0.898169 -0.644194 0.0
rev_rating 0.599947 0.478713 0.715597 0.0 0.0 0.0 0.0 0.0 0.599947 0.478713 0.715597 0.0
num_spots 0.357872 0.207894 0.499055 0.0 0.0 0.0 0.0 0.0 0.357872 0.207894 0.499055 0.0
crowded -0.492793 -0.613448 -0.374854 0.0 0.0 0.0 0.0 0.0 -0.492793 -0.613448 -0.374854 0.0
mean sd ess_bulk ess_tail r_hat
lam 0.237 0.101 963.662 388.022 1.002
beta[Intercept] -24246.255 1031031.388 770.792 272.065 1.002
beta[poverty] -0.772 0.065 1362.188 517.233 0.997
beta[rev_rating] 0.600 0.064 911.255 512.937 1.012
beta[num_spots] 0.358 0.076 847.769 515.956 1.003
beta[crowded] -0.493 0.060 1368.274 579.247 1.001
sigma 0.912 0.039 839.287 415.063 0.999
../_images/5e70645268d53aebdfbe38c1e1a774f69a7da0289ab81b5425ba4d8aead10986.png

Fit: SDMPanelFE

sdm_panel, summary_sdm, effects_sdm = fit_panel_model(
    SDMPanelFE, formula, panel, W, model=3
)
display(summary_sdm)
display(effects_sdm)
display(diagnostics_table(sdm_panel.inference_data, ["rho", "beta", "sigma"]))
show_trace(sdm_panel.inference_data, ["rho", "sigma"], "SDMPanelFE trace")

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [rho, beta, sigma]
Sampling 2 chains for 300 tune and 300 draw iterations (600 + 600 draws total) took 6 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
/tmp/ipykernel_8142/2596509600.py:42: UserWarning: The figure layout has changed to tight
  plt.tight_layout()
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept 34587.630 1019167.264 -1865686.365 1866428.626 38032.659 40070.974 717.266 502.445 1.007
poverty -0.775 0.065 -0.888 -0.643 0.002 0.004 992.792 380.606 1.003
rev_rating 0.645 0.071 0.522 0.783 0.003 0.003 690.804 402.857 1.002
num_spots 0.363 0.075 0.219 0.493 0.002 0.003 1003.892 467.368 1.000
crowded -0.486 0.067 -0.600 -0.345 0.002 0.003 835.617 441.015 1.001
W*poverty 0.040 0.172 -0.285 0.369 0.008 0.006 411.047 419.883 1.003
W*rev_rating 0.254 0.204 -0.104 0.629 0.008 0.007 710.387 504.674 1.003
W*num_spots 0.419 0.195 0.046 0.736 0.007 0.006 827.667 539.682 1.004
W*crowded 0.213 0.165 -0.116 0.478 0.006 0.006 881.069 468.563 0.999
rho 0.206 0.099 0.025 0.388 0.004 0.003 528.233 422.862 0.999
sigma 0.897 0.040 0.827 0.974 0.001 0.002 760.620 392.460 1.002
direct direct_ci_lower direct_ci_upper direct_pvalue indirect indirect_ci_lower indirect_ci_upper indirect_pvalue total total_ci_lower total_ci_upper total_pvalue
variable
poverty -0.779551 -0.918199 -0.648176 0.0 -0.147005 -0.510916 0.232309 0.433333 -0.926555 -1.338995 -0.493122 0.000000
rev_rating 0.658783 0.512025 0.789349 0.0 0.478657 -0.008991 0.979734 0.060000 1.137440 0.583078 1.732045 0.000000
num_spots 0.380620 0.228367 0.522112 0.0 0.610940 0.166183 1.051684 0.003333 0.991560 0.507988 1.504997 0.000000
crowded -0.481824 -0.618741 -0.349318 0.0 0.139059 -0.269562 0.533493 0.496667 -0.342765 -0.790152 0.117110 0.133333
mean sd ess_bulk ess_tail r_hat
rho 0.206 0.099 528.233 422.862 0.999
beta[Intercept] 34587.630 1019167.264 717.266 502.445 1.007
beta[poverty] -0.775 0.065 992.792 380.606 1.003
beta[rev_rating] 0.645 0.071 690.804 402.857 1.002
beta[num_spots] 0.363 0.075 1003.892 467.368 1.000
beta[crowded] -0.486 0.067 835.617 441.015 1.001
beta[W*poverty] 0.040 0.172 411.047 419.883 1.003
beta[W*rev_rating] 0.254 0.204 710.387 504.674 1.003
beta[W*num_spots] 0.419 0.195 827.667 539.682 1.004
beta[W*crowded] 0.213 0.165 881.069 468.563 0.999
sigma 0.897 0.040 760.620 392.460 1.002
../_images/9e3f7645e88dac158c8b6bb3ab07e91490ef098fbfa240d9edec922fa4a6f226.png

Fit: SDEMPanelFE

sdem_panel, summary_sdem, effects_sdem = fit_panel_model(
    SDEMPanelFE, formula, panel, W, model=3
)
display(summary_sdem)
display(effects_sdem)
display(diagnostics_table(sdem_panel.inference_data, ["lam", "beta", "sigma"]))
show_trace(sdem_panel.inference_data, ["lam", "sigma"], "SDEMPanelFE trace")

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [lam, beta, sigma]
Sampling 2 chains for 300 tune and 300 draw iterations (600 + 600 draws total) took 6 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
/tmp/ipykernel_8142/2596509600.py:42: UserWarning: The figure layout has changed to tight
  plt.tight_layout()
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept 14239.715 971945.515 -1694232.195 1990422.215 32678.216 37979.119 895.797 491.463 1.001
poverty -0.782 0.064 -0.900 -0.668 0.002 0.003 959.971 549.410 0.997
rev_rating 0.659 0.075 0.521 0.792 0.003 0.003 611.073 466.179 0.998
num_spots 0.382 0.075 0.246 0.525 0.002 0.003 1077.238 539.092 1.003
crowded -0.489 0.062 -0.611 -0.383 0.002 0.002 850.873 556.673 1.004
W*poverty -0.158 0.165 -0.520 0.123 0.005 0.007 939.327 422.862 1.008
W*rev_rating 0.439 0.221 0.041 0.849 0.008 0.008 817.469 445.688 0.997
W*num_spots 0.590 0.211 0.201 0.975 0.007 0.008 929.699 546.450 1.002
W*crowded 0.116 0.194 -0.217 0.503 0.007 0.007 749.457 463.143 1.003
lam 0.256 0.095 0.069 0.426 0.003 0.004 1094.502 488.178 1.012
sigma 0.894 0.042 0.806 0.969 0.001 0.002 805.225 402.429 1.000
direct direct_ci_lower direct_ci_upper direct_pvalue indirect indirect_ci_lower indirect_ci_upper indirect_pvalue total total_ci_lower total_ci_upper total_pvalue
variable
poverty -0.781545 -0.908611 -0.662342 0.0 -0.158036 -0.504484 0.192397 0.326667 -0.939581 -1.359283 -0.548555 0.0
rev_rating 0.659203 0.514644 0.792674 0.0 0.439298 0.016644 0.854632 0.046667 1.098501 0.575104 1.625128 0.0
num_spots 0.381790 0.241262 0.534002 0.0 0.590301 0.171367 1.004782 0.006667 0.972091 0.522053 1.412150 0.0
crowded -0.488718 -0.608805 -0.368327 0.0 0.116017 -0.255579 0.501576 0.533333 -0.372702 -0.795999 0.070927 0.1
mean sd ess_bulk ess_tail r_hat
lam 0.256 0.095 1094.502 488.178 1.012
beta[Intercept] 14239.715 971945.515 895.797 491.463 1.001
beta[poverty] -0.782 0.064 959.971 549.410 0.997
beta[rev_rating] 0.659 0.075 611.073 466.179 0.998
beta[num_spots] 0.382 0.075 1077.238 539.092 1.003
beta[crowded] -0.489 0.062 850.873 556.673 1.004
beta[W*poverty] -0.158 0.165 939.327 422.862 1.008
beta[W*rev_rating] 0.439 0.221 817.469 445.688 0.997
beta[W*num_spots] 0.590 0.211 929.699 546.450 1.002
beta[W*crowded] 0.116 0.194 749.457 463.143 1.003
sigma 0.894 0.042 805.225 402.429 1.000
../_images/6ab9384debe4c65e672a098a9e6faa0db20e071e9bee2d88da9d1019d54f8c90.png

MCMC adequacy for the spatial parameter

Panel SAR/SDEM models inherit the slow-mixing behaviour of \(\rho\) (or \(\lambda\)) documented by Wolf, Anselin & Arribas-Bel (2018) [Wolf et al., 2018]. The spatial_mcmc_diagnostic helper checks ESS, sampler yield, \(\hat{R}\), and HPDI stability for the spatial scalar.

from bayespecon import spatial_mcmc_diagnostic

spatial_mcmc_diagnostic(sdm_panel, emit_warnings=False).to_frame()
ess_bulk ess_tail r_hat mcse_mean yield_pct hpdi_drift_pct adequate
parameter
rho 528.233427 422.862388 0.999436 0.004291 88.038904 1.505472 False

Random-Effects Panel Models

The random-effects models replace deterministic unit fixed effects with latent unit intercepts \(\alpha_i \sim N(0, \sigma_\alpha^2)\). That lets the notebook estimate the variance of unit heterogeneity directly and compare FE and RE behavior on the same synthetic balanced panel.

Fit: OLSPanelRE

ols_panel_re, summary_ols_re, effects_ols_re = fit_panel_model(
    OLSPanelRE, formula, panel, W, model=0
)
display(summary_ols_re)
display(effects_ols_re)
display(
    diagnostics_table(ols_panel_re.inference_data, ["beta", "sigma", "sigma_alpha"])
)
show_trace(ols_panel_re.inference_data, ["sigma", "sigma_alpha"], "OLSPanelRE trace")

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma, sigma_alpha, alpha]
Sampling 2 chains for 300 tune and 300 draw iterations (600 + 600 draws total) took 6 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
/tmp/ipykernel_8142/2596509600.py:42: UserWarning: The figure layout has changed to tight
  plt.tight_layout()
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept 2.306 0.096 2.134 2.483 0.005 0.003 337.869 493.011 1.003
poverty -0.825 0.067 -0.976 -0.718 0.002 0.003 1179.867 487.587 1.000
rev_rating 0.635 0.078 0.506 0.783 0.002 0.003 1203.406 416.888 1.000
num_spots 0.418 0.084 0.264 0.566 0.002 0.003 1160.332 558.398 0.999
crowded -0.467 0.069 -0.596 -0.342 0.002 0.003 1035.486 378.629 1.004
... ... ... ... ... ... ... ... ... ...
alpha[61] -0.384 0.397 -1.177 0.299 0.011 0.022 1336.922 416.492 1.008
alpha[62] -0.483 0.421 -1.312 0.243 0.015 0.018 737.316 431.464 1.005
alpha[63] -0.061 0.376 -0.792 0.678 0.011 0.020 1266.437 379.581 1.000
sigma 1.076 0.054 0.968 1.170 0.002 0.002 779.873 376.840 1.007
sigma_alpha 0.569 0.101 0.372 0.748 0.007 0.003 197.392 381.333 1.007

71 rows × 9 columns

direct direct_ci_lower direct_ci_upper direct_pvalue indirect indirect_ci_lower indirect_ci_upper indirect_pvalue total total_ci_lower total_ci_upper total_pvalue
variable
poverty -0.824634 -0.962644 -0.691388 0.0 0.0 0.0 0.0 0.0 -0.824634 -0.962644 -0.691388 0.0
rev_rating 0.635404 0.493474 0.781262 0.0 0.0 0.0 0.0 0.0 0.635404 0.493474 0.781262 0.0
num_spots 0.418216 0.252982 0.570942 0.0 0.0 0.0 0.0 0.0 0.418216 0.252982 0.570942 0.0
crowded -0.466801 -0.598624 -0.336878 0.0 0.0 0.0 0.0 0.0 -0.466801 -0.598624 -0.336878 0.0
mean sd ess_bulk ess_tail r_hat
beta[Intercept] 2.306 0.096 337.869 493.011 1.003
beta[poverty] -0.825 0.067 1179.867 487.587 1.000
beta[rev_rating] 0.635 0.078 1203.406 416.888 1.000
beta[num_spots] 0.418 0.084 1160.332 558.398 0.999
beta[crowded] -0.467 0.069 1035.486 378.629 1.004
sigma 1.076 0.054 779.873 376.840 1.007
sigma_alpha 0.569 0.101 197.392 381.333 1.007
../_images/9db76d015ccf60747dc90427ba4d96b04ff4f91b64bb0e20b72ce86a0a5a88c8.png

Fit: SARPanelRE

sar_panel_re, summary_sar_re, effects_sar_re = fit_panel_model(
    SARPanelRE, formula, panel, W, model=0
)
display(summary_sar_re)
display(effects_sar_re)
display(
    diagnostics_table(
        sar_panel_re.inference_data, ["rho", "beta", "sigma", "sigma_alpha"]
    )
)
show_trace(
    sar_panel_re.inference_data, ["rho", "sigma", "sigma_alpha"], "SARPanelRE trace"
)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [rho, beta, sigma, sigma_alpha, alpha]
Sampling 2 chains for 300 tune and 300 draw iterations (600 + 600 draws total) took 6 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
/tmp/ipykernel_8142/2596509600.py:42: UserWarning: The figure layout has changed to tight
  plt.tight_layout()
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept 1.568 0.207 1.253 2.002 0.013 0.009 258.305 265.162 1.002
poverty -0.822 0.065 -0.954 -0.702 0.003 0.002 614.656 515.829 0.999
rev_rating 0.638 0.069 0.506 0.752 0.002 0.002 873.808 437.132 0.998
num_spots 0.390 0.081 0.251 0.538 0.003 0.003 738.920 547.655 1.001
crowded -0.472 0.068 -0.590 -0.343 0.003 0.003 563.675 423.879 1.002
... ... ... ... ... ... ... ... ... ...
alpha[62] -0.470 0.378 -1.227 0.165 0.018 0.014 454.997 354.585 1.007
alpha[63] -0.042 0.351 -0.707 0.603 0.011 0.014 963.985 397.817 0.998
rho 0.311 0.081 0.157 0.446 0.005 0.004 241.374 292.423 1.009
sigma 1.051 0.056 0.953 1.166 0.003 0.002 472.513 439.050 1.002
sigma_alpha 0.501 0.094 0.345 0.708 0.007 0.004 162.027 206.940 1.003

72 rows × 9 columns

direct direct_ci_lower direct_ci_upper direct_pvalue indirect indirect_ci_lower indirect_ci_upper indirect_pvalue total total_ci_lower total_ci_upper total_pvalue
variable
poverty -0.838263 -0.973437 -0.713678 0.0 -0.370082 -0.662114 -0.138628 0.0 -1.208345 -1.543110 -0.906244 0.0
rev_rating 0.650614 0.522126 0.786255 0.0 0.287146 0.107687 0.527340 0.0 0.937760 0.689768 1.273430 0.0
num_spots 0.397257 0.249877 0.554125 0.0 0.174668 0.060997 0.342674 0.0 0.571925 0.336794 0.848725 0.0
crowded -0.481035 -0.610604 -0.344887 0.0 -0.212878 -0.388136 -0.071129 0.0 -0.693913 -0.967211 -0.454681 0.0
mean sd ess_bulk ess_tail r_hat
rho 0.311 0.081 241.374 292.423 1.009
beta[Intercept] 1.568 0.207 258.305 265.162 1.002
beta[poverty] -0.822 0.065 614.656 515.829 0.999
beta[rev_rating] 0.638 0.069 873.808 437.132 0.998
beta[num_spots] 0.390 0.081 738.920 547.655 1.001
beta[crowded] -0.472 0.068 563.675 423.879 1.002
sigma 1.051 0.056 472.513 439.050 1.002
sigma_alpha 0.501 0.094 162.027 206.940 1.003
../_images/b390c3d1643e5d1b5d8226cda27b121d3e2c0bc0076672eaddb7c75bf2ca5c1e.png

Fit: SEMPanelRE

sem_panel_re, summary_sem_re, effects_sem_re = fit_panel_model(
    SEMPanelRE, formula, panel, W, model=0
)
display(summary_sem_re)
display(effects_sem_re)
display(
    diagnostics_table(
        sem_panel_re.inference_data, ["lam", "beta", "sigma", "sigma_alpha"]
    )
)
show_trace(
    sem_panel_re.inference_data, ["lam", "sigma", "sigma_alpha"], "SEMPanelRE trace"
)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [lam, beta, sigma, sigma_alpha, alpha]
Sampling 2 chains for 300 tune and 300 draw iterations (600 + 600 draws total) took 7 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
/tmp/ipykernel_8142/2596509600.py:42: UserWarning: The figure layout has changed to tight
  plt.tight_layout()
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept 2.313 0.116 2.116 2.535 0.004 0.004 678.194 396.660 1.002
poverty -0.816 0.072 -0.950 -0.693 0.003 0.003 741.497 438.429 1.003
rev_rating 0.616 0.068 0.486 0.742 0.003 0.003 668.569 428.846 1.000
num_spots 0.377 0.083 0.231 0.531 0.003 0.004 788.356 340.977 1.018
crowded -0.472 0.063 -0.581 -0.353 0.002 0.002 660.402 378.602 1.001
... ... ... ... ... ... ... ... ... ...
alpha[62] -0.353 0.363 -1.071 0.331 0.019 0.017 354.366 359.628 1.003
alpha[63] -0.007 0.347 -0.591 0.633 0.009 0.020 1666.891 356.491 1.009
lam 0.340 0.105 0.160 0.548 0.004 0.004 574.602 424.964 1.003
sigma 1.074 0.057 0.962 1.180 0.003 0.002 356.763 515.954 1.005
sigma_alpha 0.448 0.122 0.168 0.650 0.020 0.013 41.714 26.368 1.055

72 rows × 9 columns

direct direct_ci_lower direct_ci_upper direct_pvalue indirect indirect_ci_lower indirect_ci_upper indirect_pvalue total total_ci_lower total_ci_upper total_pvalue
variable
poverty -0.816111 -0.953122 -0.679809 0.0 0.0 0.0 0.0 0.0 -0.816111 -0.953122 -0.679809 0.0
rev_rating 0.615899 0.474850 0.739893 0.0 0.0 0.0 0.0 0.0 0.615899 0.474850 0.739893 0.0
num_spots 0.377059 0.219276 0.552455 0.0 0.0 0.0 0.0 0.0 0.377059 0.219276 0.552455 0.0
crowded -0.472267 -0.592100 -0.353982 0.0 0.0 0.0 0.0 0.0 -0.472267 -0.592100 -0.353982 0.0
mean sd ess_bulk ess_tail r_hat
lam 0.340 0.105 574.602 424.964 1.003
beta[Intercept] 2.313 0.116 678.194 396.660 1.002
beta[poverty] -0.816 0.072 741.497 438.429 1.003
beta[rev_rating] 0.616 0.068 668.569 428.846 1.000
beta[num_spots] 0.377 0.083 788.356 340.977 1.018
beta[crowded] -0.472 0.063 660.402 378.602 1.001
sigma 1.074 0.057 356.763 515.954 1.005
sigma_alpha 0.448 0.122 41.714 26.368 1.055
../_images/9098ea18a25fcfd9e7e1f09f9cf89ae5bc918a4936b12c81c28842c6e78c3d4c.png

Model Comparison (WAIC and LOO) Across FE and RE Models

# Collect all fitted panel models for comparison
panel_models = {
    "OLSPanelFE": ols_panel,
    "SARPanelFE": sar_panel,
    "SEMPanelFE": sem_panel,
    "SDMPanelFE": sdm_panel,
    "SDEMPanelFE": sdem_panel,
    "OLSPanelRE": ols_panel_re,
    "SARPanelRE": sar_panel_re,
    "SEMPanelRE": sem_panel_re,
}
idata_dict = {name: m.inference_data for name, m in panel_models.items()}

# WAIC and LOO comparison
for ic in ("waic", "loo"):
    try:
        cmp = az.compare(idata_dict, ic=ic, method="BB-pseudo-BMA")
        print(f"\n{ic.upper()} comparison")
        display(cmp)
    except Exception as e:
        print(f"{ic.upper()} comparison not available: {type(e).__name__}: {e}")

# Per-model IC table
rows = []
for name, model in panel_models.items():
    idata = model.inference_data
    row = {"model": name}
    try:
        waic_res = az.waic(idata)
        row["elpd_waic"] = float(waic_res.elpd_waic)
        row["p_waic"] = float(waic_res.p_waic)
    except Exception:
        row["elpd_waic"] = np.nan
        row["p_waic"] = np.nan
    try:
        loo_res = az.loo(idata)
        row["elpd_loo"] = float(loo_res.elpd_loo)
        row["p_loo"] = float(loo_res.p_loo)
    except Exception:
        row["elpd_loo"] = np.nan
        row["p_loo"] = np.nan
    rows.append(row)

pd.DataFrame(rows).sort_values("model").reset_index(drop=True)

WAIC comparison

LOO comparison

/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.64 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.64 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.64 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: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.64 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: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.64 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: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.64 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(
rank elpd_waic p_waic elpd_diff weight se dse warning scale
SDEMPanelFE 0 -339.893972 9.325117 0.000000 4.501560e-01 11.991857 0.000000 True log
SDMPanelFE 1 -340.533526 9.307454 0.639554 2.344013e-01 11.991790 0.541625 True log
SARPanelFE 2 -341.718297 5.848315 1.824325 1.805254e-01 11.640996 2.958246 True log
SEMPanelFE 3 -343.202940 5.448393 3.308968 1.029303e-01 11.739321 3.635423 False log
OLSPanelFE 4 -345.359753 5.139489 5.465781 3.198706e-02 11.836319 4.162486 False log
SARPanelRE 5 -395.935964 33.745352 56.041992 1.678728e-18 11.962884 7.311718 True log
SEMPanelRE 6 -401.274137 31.303476 61.380165 1.415375e-19 12.082551 8.264962 True log
OLSPanelRE 7 -401.533418 36.094738 61.639446 7.042905e-20 12.189155 7.647900 True log
rank elpd_loo p_loo elpd_diff weight se dse warning scale
SDEMPanelFE 0 -339.972133 9.403278 0.000000 4.568785e-01 11.445295 0.000000 False log
SDMPanelFE 1 -340.626175 9.400104 0.654042 2.362939e-01 11.467946 0.545009 False log
SARPanelFE 2 -341.765300 5.895318 1.793167 1.881267e-01 11.205857 2.959321 False log
SEMPanelFE 3 -343.238728 5.484181 3.266595 8.774769e-02 11.406379 3.636596 False log
OLSPanelFE 4 -345.392961 5.172696 5.420828 3.095314e-02 11.518371 4.163942 False log
SARPanelRE 5 -396.559968 34.369355 56.587835 2.820065e-20 11.329130 7.309362 True log
SEMPanelRE 6 -401.945528 31.974867 61.973394 1.743175e-20 11.761979 8.238647 True log
OLSPanelRE 7 -402.168361 36.729682 62.196228 1.269921e-20 11.893886 7.662275 True log
model elpd_waic p_waic elpd_loo p_loo
0 OLSPanelFE -345.359753 5.139489 -345.392961 5.172696
1 OLSPanelRE -401.533418 36.094738 -402.168361 36.729682
2 SARPanelFE -341.718297 5.848315 -341.765300 5.895318
3 SARPanelRE -395.935964 33.745352 -396.559968 34.369355
4 SDEMPanelFE -339.893972 9.325117 -339.972133 9.403278
5 SDMPanelFE -340.533526 9.307454 -340.626175 9.400104
6 SEMPanelFE -343.202940 5.448393 -343.238728 5.484181
7 SEMPanelRE -401.274137 31.303476 -401.945528 31.974867

Notes

  • This notebook uses modest draw counts for speed. For real analyses, increase draws and tune.

  • If you see divergences or tree-depth warnings, increase target_accept (for example 0.95-0.99).

  • For production inference, validate sensitivity to FE specification (model=0/1/2/3) and priors.

Bayesian Panel LM Diagnostics

The bayespecon.diagnostics module provides a full suite of Bayesian Lagrange Multiplier (LM) tests for panel models, following the framework of Dogan et al. (2021) with panel-specific formulas from Anselin (2008), Elhorst (2014), and Koley & Bera (2024).

These tests help answer the question: which spatial panel specification is most appropriate? They test for omitted spatial lag (\(\rho\)), spatial error (\(\lambda\)), and spatially lagged covariates (\(\gamma\)), both individually and jointly, with robust variants that account for the presence of other spatial effects.

The decision tree mirrors the classical Anselin (1988) / Koley-Bera (2024) approach, but uses the full posterior distribution rather than point estimates — yielding Bayesian p-values and credible intervals for each LM statistic.

from bayespecon.diagnostics.bayesian_lmtests import (
    bayesian_panel_lm_error_test,
    bayesian_panel_lm_lag_test,
    bayesian_panel_lm_sdm_joint_test,
    bayesian_panel_lm_slx_error_joint_test,
    bayesian_panel_lm_wx_test,
    bayesian_panel_robust_lm_error_sdem_test,
    bayesian_panel_robust_lm_error_test,
    bayesian_panel_robust_lm_lag_sdm_test,
    bayesian_panel_robust_lm_lag_test,
    bayesian_panel_robust_lm_wx_test,
)

Standard Panel LM Tests (from OLS null)

The four standard panel LM tests use the OLS panel model as the null. They test for spatial lag, spatial error, and joint SDM/SDEM alternatives.

# Standard panel LM tests (from OLS null)
panel_standard = pd.DataFrame(
    [
        bayesian_panel_lm_lag_test(ols_panel).to_series(),
        bayesian_panel_lm_error_test(ols_panel).to_series(),
        bayesian_panel_lm_sdm_joint_test(ols_panel).to_series(),
        bayesian_panel_lm_slx_error_joint_test(ols_panel).to_series(),
    ],
    index=["LM-Lag", "LM-Error", "LM-SDM Joint", "LM-SLX-Error Joint"],
)

panel_standard

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/bayespecon/diagnostics/bayesian_lmtests.py:1704: RuntimeWarning: X'X (panel LM-lag) is ill-conditioned (cond=2.53e+14); falling back to pseudo-inverse.
  XtX_inv = _safe_inv(X.T @ X, "X'X (panel LM-lag)")
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/bayespecon/diagnostics/bayesian_lmtests.py:1588: RuntimeWarning: X'X (panel info blocks) is ill-conditioned (cond=2.53e+14); falling back to pseudo-inverse.
  XtX_inv = _safe_inv(X.T @ X, "X'X (panel info blocks)")
lm_samples mean median credible_interval bayes_pvalue test_type df n_draws N T k_wx
LM-Lag [7.846839032979865, 10.490111917114326, 7.9766... 9.051813 9.064912 (7.012499679824649, 11.372619935132999) 0.002624 bayesian_panel_lm_lag 1 600 64 4 NaN
LM-Error [4.111440912077958, 6.420019554621757, 4.08610... 5.144037 4.982020 (2.4361821394356338, 8.880377987635839) 0.023326 bayesian_panel_lm_error 1 600 64 4 NaN
LM-SDM Joint [12.272450028161778, 14.847758677815387, 12.46... 13.276470 13.235997 (11.672859758702915, 15.131582797987058) 0.020921 bayesian_panel_lm_sdm_joint 5 600 64 4 4.0
LM-SLX-Error Joint [12.508512249606115, 15.652369984378204, 12.63... 13.758191 13.711836 (11.442624198522244, 16.288082901401626) 0.017221 bayesian_panel_lm_slx_error_joint 5 600 64 4 4.0

Robust Panel LM Tests (from OLS null)

The robust variants (Elhorst, 2014) test one spatial effect while accounting for the possible local presence of the other. For example, the robust LM-Lag tests for \(\rho\) robust to \(\lambda\), and vice versa.

# Robust panel LM tests (from OLS null)
panel_robust = pd.DataFrame(
    [
        bayesian_panel_robust_lm_lag_test(ols_panel).to_series(),
        bayesian_panel_robust_lm_error_test(ols_panel).to_series(),
    ],
    index=["Robust LM-Lag", "Robust LM-Error"],
)

panel_robust

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/bayespecon/diagnostics/bayesian_lmtests.py:1881: RuntimeWarning: X'X (panel robust LM-lag) is ill-conditioned (cond=2.53e+14); falling back to pseudo-inverse.
  XtX_inv = _safe_inv(X.T @ X, "X'X (panel robust LM-lag)")
/home/runner/micromamba/envs/test/lib/python3.14/site-packages/bayespecon/diagnostics/bayesian_lmtests.py:1976: RuntimeWarning: X'X (panel robust LM-error) is ill-conditioned (cond=2.53e+14); falling back to pseudo-inverse.
  XtX_inv = _safe_inv(X.T @ X, "X'X (panel robust LM-error)")
lm_samples mean median credible_interval bayes_pvalue test_type df n_draws N T
Robust LM-Lag [3.1878927333332534, 4.79017592049596, 3.92507... 4.255999 4.078053 (2.274674009417828, 7.073934310539664) 0.039112 bayesian_panel_robust_lm_lag 1 600 64 4
Robust LM-Error [0.044156243612000164, 0.0522672095975571, 0.1... 0.211100 0.099645 (0.0005608562592401218, 0.987332424243747) 0.645907 bayesian_panel_robust_lm_error 1 600 64 4

SDM/SDEM Variant Tests

These tests address the Koley & Bera (2024) decision tree for choosing between SAR/SEM and SDM/SDEM specifications. They require fitting intermediate models (SAR or SLX) as the alternative model.

  • LM-WX (from SAR null): Should we extend SAR → SDM by adding \(\gamma\)?

  • Robust LM-Lag-SDM (from SLX null): Should we extend SLX → SDM by adding \(\rho\), robust to \(\gamma\)?

  • Robust LM-WX (from SAR null): Should we extend SAR → SDM by adding \(\gamma\), robust to \(\rho\)?

  • Robust LM-Error-SDEM (from SLX null): Should we extend SLX → SDEM by adding \(\lambda\), robust to \(\gamma\)?

We need to fit an SLX panel model first to serve as the null for the robust SDM/SDEM tests.

from bayespecon import SLXPanelFE

slx_panel, summary_slx, effects_slx = fit_panel_model(
    SLXPanelFE, formula, panel, W, model=3
)
display(summary_slx)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, sigma]
Sampling 2 chains for 300 tune and 300 draw iterations (600 + 600 draws total) took 6 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
Intercept 23013.637 1119295.980 -2001523.165 2033869.360 33205.778 55061.102 1132.165 288.850 1.005
poverty -0.778 0.065 -0.891 -0.656 0.002 0.003 1050.968 467.836 1.001
rev_rating 0.657 0.068 0.523 0.772 0.003 0.003 697.044 392.270 1.014
num_spots 0.371 0.079 0.231 0.510 0.002 0.004 1060.848 400.256 1.005
crowded -0.479 0.060 -0.583 -0.357 0.002 0.002 1390.022 561.679 0.998
W*poverty -0.128 0.169 -0.481 0.173 0.004 0.008 1639.578 386.282 0.999
W*rev_rating 0.385 0.178 0.066 0.709 0.006 0.008 1027.626 539.242 1.016
W*num_spots 0.538 0.187 0.211 0.885 0.006 0.008 865.475 395.011 1.002
W*crowded 0.114 0.169 -0.218 0.412 0.004 0.006 1444.035 496.294 1.002
sigma 0.905 0.042 0.836 0.988 0.001 0.002 1433.403 466.069 1.002 
# SDM/SDEM variant panel LM tests
panel_sdem = pd.DataFrame(
    [
        bayesian_panel_lm_wx_test(sar_panel).to_series(),
        bayesian_panel_robust_lm_lag_sdm_test(slx_panel).to_series(),
        bayesian_panel_robust_lm_wx_test(sar_panel).to_series(),
        bayesian_panel_robust_lm_error_sdem_test(slx_panel).to_series(),
    ],
    index=["LM-WX", "Robust LM-Lag-SDM", "Robust LM-WX", "Robust LM-Error-SDEM"],
)

panel_sdem

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/bayespecon/diagnostics/bayesian_lmtests.py:1588: RuntimeWarning: X'X (panel info blocks) is ill-conditioned (cond=2.53e+14); falling back to pseudo-inverse.
  XtX_inv = _safe_inv(X.T @ X, "X'X (panel info blocks)")
lm_samples mean median credible_interval bayes_pvalue test_type df n_draws k_wx N T
LM-WX [9.052266749530466, 5.88295092917965, 7.078770... 6.554297 6.426463 (5.494612656398929, 8.285775877938786) 0.161401 bayesian_panel_lm_wx 4 600 4 64 4
Robust LM-Lag-SDM [2.455700587495692, 2.3422889610625552, 2.5800... 2.476347 2.453508 (1.6629588840389784, 3.308331312416007) 0.115570 bayesian_panel_robust_lm_lag_sdm 1 600 4 64 4
Robust LM-WX [8.390384590189568, 6.674678445761224, 6.97317... 6.420691 6.372788 (5.692600783839155, 7.375952689100776) 0.169857 bayesian_panel_robust_lm_wx 4 600 4 64 4
Robust LM-Error-SDEM [4.846749485308717, 3.5910547807901474, 4.5415... 4.355527 4.205693 (3.0815884770940833, 6.574888886195661) 0.036889 bayesian_panel_robust_lm_error_sdem 1 600 4 64 4

Interpreting the Panel LM Decision Tree

The panel LM tests follow the same decision logic as the classical Anselin/Koley-Bera approach, but with Bayesian p-values:

Scenario

Significant?

Interpretation

LM-Lag ✓, LM-Error ✗

Lag only

Choose SAR panel

LM-Lag ✗, LM-Error ✓

Error only

Choose SEM panel

Both significant → check robust tests

Robust LM-Lag ✓, Robust LM-Error ✗

Lag dominates

Choose SAR panel

Robust LM-Lag ✗, Robust LM-Error ✓

Error dominates

Choose SEM panel

Both robust significant

Both effects

Consider SDM or SDEM

Joint LM-SDM ✓

SDM preferred over OLS

Fit SDM panel

Joint LM-SLX-Error ✓

SDEM preferred over OLS

Fit SDEM panel

LM-WX ✓ (from SAR)

WX effects present

Extend SAR → SDM

Robust LM-Lag-SDM ✓ (from SLX)

Lag needed beyond WX

Extend SLX → SDM

Robust LM-WX ✓ (from SAR)

WX needed beyond lag

Extend SAR → SDM

Robust LM-Error-SDEM ✓ (from SLX)

Error needed beyond WX

Extend SLX → SDEM

A Bayesian p-value below 0.05 (or a threshold you choose) indicates evidence against the null hypothesis. The credible_interval field shows the posterior uncertainty in the LM statistic itself.

from scipy import stats as sp_stats

# All panel LM tests in a single table
all_panel = pd.concat([panel_standard, panel_robust, panel_sdem])

# Add chi-squared reference values
all_panel["chi2_ref_95"] = all_panel["df"].apply(lambda df: sp_stats.chi2.ppf(0.95, df))

all_panel
lm_samples mean median credible_interval bayes_pvalue test_type df n_draws N T k_wx chi2_ref_95
LM-Lag [7.846839032979865, 10.490111917114326, 7.9766... 9.051813 9.064912 (7.012499679824649, 11.372619935132999) 0.002624 bayesian_panel_lm_lag 1 600 64 4 NaN 3.841459
LM-Error [4.111440912077958, 6.420019554621757, 4.08610... 5.144037 4.982020 (2.4361821394356338, 8.880377987635839) 0.023326 bayesian_panel_lm_error 1 600 64 4 NaN 3.841459
LM-SDM Joint [12.272450028161778, 14.847758677815387, 12.46... 13.276470 13.235997 (11.672859758702915, 15.131582797987058) 0.020921 bayesian_panel_lm_sdm_joint 5 600 64 4 4.0 11.070498
LM-SLX-Error Joint [12.508512249606115, 15.652369984378204, 12.63... 13.758191 13.711836 (11.442624198522244, 16.288082901401626) 0.017221 bayesian_panel_lm_slx_error_joint 5 600 64 4 4.0 11.070498
Robust LM-Lag [3.1878927333332534, 4.79017592049596, 3.92507... 4.255999 4.078053 (2.274674009417828, 7.073934310539664) 0.039112 bayesian_panel_robust_lm_lag 1 600 64 4 NaN 3.841459
Robust LM-Error [0.044156243612000164, 0.0522672095975571, 0.1... 0.211100 0.099645 (0.0005608562592401218, 0.987332424243747) 0.645907 bayesian_panel_robust_lm_error 1 600 64 4 NaN 3.841459
LM-WX [9.052266749530466, 5.88295092917965, 7.078770... 6.554297 6.426463 (5.494612656398929, 8.285775877938786) 0.161401 bayesian_panel_lm_wx 4 600 64 4 4.0 9.487729
Robust LM-Lag-SDM [2.455700587495692, 2.3422889610625552, 2.5800... 2.476347 2.453508 (1.6629588840389784, 3.308331312416007) 0.115570 bayesian_panel_robust_lm_lag_sdm 1 600 64 4 4.0 3.841459
Robust LM-WX [8.390384590189568, 6.674678445761224, 6.97317... 6.420691 6.372788 (5.692600783839155, 7.375952689100776) 0.169857 bayesian_panel_robust_lm_wx 4 600 64 4 4.0 9.487729
Robust LM-Error-SDEM [4.846749485308717, 3.5910547807901474, 4.5415... 4.355527 4.205693 (3.0815884770940833, 6.574888886195661) 0.036889 bayesian_panel_robust_lm_error_sdem 1 600 64 4 4.0 3.841459

  • LM-Lag is highly significant (p=0.002) — correctly detecting the spatial lag used in the DGP

  • Robust LM-Lag remains significant (p=0.039) while Robust LM-Error is not (p=0.66) — pointing to SAR over SEM

  • Joint LM-SDM and Joint LM-SLX-Error are both significant — suggesting spatial structure beyond OLS

  • LM-WX from SAR null is not significant (p=0.17) — the DGP had no WX effects, so SAR is sufficient

  • Robust LM-Error-SDEM is marginally significant (p=0.036) — slight error signal when WX is present