Bayesian LM tests against classical `spreg`¶

Cross-implementation check for the Lagrange-Multiplier specification tests. neighbayes evaluates the LM statistic at every posterior draw and reports the posterior mean; spreg computes the classical point statistic from an OLS fit. On data simulated under \(H_0\) with weakly-informative priors the two should tell the same story, and where they diverge the reason should be understood rather than assumed.

This is evidence about the implementation, not a guide to using it — for that see How to run Bayesian LM specification tests.

import warnings

import libpysal
import numpy as np
import pandas as pd
from spreg import OLS as SpregOLS
from spreg import LMtests as SpregLMtests

from neighbayes.diagnostics.lmtests import (
    bayesian_lm_error_test,
    bayesian_lm_lag_test,
    bayesian_lm_sdm_joint_test,
    bayesian_lm_slx_error_joint_test,
    bayesian_lm_wx_sem_test,
    bayesian_lm_wx_test,
    bayesian_robust_lm_error_sdem_test,
    bayesian_robust_lm_error_test,
    bayesian_robust_lm_lag_sdm_test,
    bayesian_robust_lm_lag_test,
    bayesian_robust_lm_wx_test,
)
from neighbayes.models import OLS, SAR, SEM, SLX

warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning)

/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

import geopandas as gpd
from libpysal.graph import Graph
from libpysal.weights import Rook

# Generate data under H0: no spatial effects
np.random.seed(42)

# Load Columbus dataset for spatial weights
columbus_path = libpysal.examples.get_path("columbus.shp")
gdf = gpd.read_file(columbus_path)

# Create a Graph (modern libpysal API) for neighbayes models
# Row-standardize the graph so spatial models work correctly
g = Graph.build_contiguity(gdf, rook=True).transform("r")
n = g.n

# Legacy W for spreg comparison
w_spreg = Rook.from_shapefile(columbus_path)
w_spreg.transform = "r"

# Get sparse and dense W matrices
W_sparse = g.sparse.tocsr().astype(np.float64)
W_dense = np.array(W_sparse.todense())

# Design matrix
k = 3
X = np.column_stack([np.ones(n), np.random.normal(size=(n, k - 1))])
beta_true = np.array([1.0, 2.0, -1.5])
y = X @ beta_true + np.random.normal(scale=1.0, size=n)

print(f"W shape: {W_sparse.shape}, nnz: {W_sparse.nnz}")

W shape: (49, 49), nnz: 200

# Fit OLS model (null for joint tests)
ols_model = OLS(y=y, X=X, W=g)
ols_model.fit(draws=5000, tune=5000, chains=4, random_seed=42)

# Fit SAR model (null for LM-WX and robust LM-WX)
sar_model = SAR(y=y, X=X, W=g)
sar_model.fit(draws=5000, tune=5000, chains=4, random_seed=42)

# Fit SLX model (null for robust LM-Lag-SDM and robust LM-Error-SDEM)
slx_model = SLX(y=y, X=X, W=g)
slx_model.fit(draws=5000, tune=5000, chains=4, random_seed=42)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 2 jobs)
NUTS: [beta, sigma2]
Sampling 4 chains for 5_000 tune and 5_000 draw iterations (20_000 + 20_000 draws total) took 17 seconds.
/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)
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 2 jobs)
NUTS: [beta, sigma2]
Sampling 4 chains for 5_000 tune and 5_000 draw iterations (20_000 + 20_000 draws total) took 18 seconds.

Gibbs sampling (sar): 4 chains for 5,000 tune and 5,000 draw iterations (4 x 10,000 = 40,000 draws total)

Sampling took 4s (9,732 draws/s)

arviz.InferenceData

Bayesian LM statistics¶

Evaluated at the null model each test requires.

# Fit the SEM null required by LM-WX-SEM (Bayesian analogue of Koley–Bera 2024).
sem_model = SEM(y=y, X=X, W=g)
sem_model.fit(draws=2500, tune=2500, chains=2, random_seed=42)

# Non-robust Bayesian LM tests evaluated at the appropriate null model.
non_robust_results = pd.DataFrame(
    [
        bayesian_lm_lag_test(ols_model).to_series(),
        bayesian_lm_error_test(ols_model).to_series(),
        bayesian_lm_wx_test(sar_model).to_series(),
        bayesian_lm_wx_sem_test(sem_model).to_series(),
        bayesian_lm_sdm_joint_test(ols_model).to_series(),
        bayesian_lm_slx_error_joint_test(ols_model).to_series(),
    ],
    index=[
        "LM-Lag",
        "LM-Error",
        "LM-WX",
        "LM-WX-SEM",
        "LM-SDM Joint",
        "LM-SLX-Error Joint",
    ],
)

non_robust_results

/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 (sem): 2 chains for 2,500 tune and 2,500 draw iterations (2 x 5,000 = 10,000 draws total)

Sampling took 4s (2,638 draws/s)

	lm_samples	mean	median	credible_interval	bayes_pvalue	test_type	df	n_draws	k_wx
LM-Lag	[0.2142964999714562, 1.0784033337488124, 1.280...	1.063234	0.656830	(0.00204675682175708, 4.2687677969034805)	0.302479	bayesian_lm_lag	1	20000	NaN
LM-Error	[0.7079187163674764, 0.7756572487282847, 0.744...	0.669052	0.683859	(0.014393878917328799, 1.3853911845443372)	0.413382	bayesian_lm_error	1	20000	NaN
LM-WX	[1.4814399069845414, 3.360987014281438, 0.8809...	2.112357	1.413838	(0.09593622926595723, 8.194730280544038)	0.347782	bayesian_lm_wx	2	20000	2.0
LM-WX-SEM	[3.2315485732635807, 1.0743489513363125, 0.196...	1.351305	0.898337	(0.05536302605253931, 5.110249054287485)	0.508824	bayesian_lm_wx_sem	2	5000	2.0
LM-SDM Joint	[0.7240693263411846, 2.7662057621909844, 2.499...	3.926935	2.786218	(0.42552310256782655, 14.128384100207635)	0.269463	bayesian_lm_sdm_joint	3	20000	2.0
LM-SLX-Error Joint	[1.3998088300832858, 3.5034345838334366, 1.597...	2.113251	1.633809	(0.7011929236378855, 6.1492522670831296)	0.549237	bayesian_lm_slx_error_joint	3	20000	2.0

# Robust Bayesian LM tests (Neyman orthogonal score)
robust_results = pd.DataFrame(
    [
        bayesian_robust_lm_lag_test(ols_model).to_series(),
        bayesian_robust_lm_error_test(ols_model).to_series(),
        bayesian_robust_lm_lag_sdm_test(slx_model).to_series(),
        bayesian_robust_lm_wx_test(sar_model).to_series(),
        bayesian_robust_lm_error_sdem_test(slx_model).to_series(),
    ],
    index=[
        "Robust LM-Lag",
        "Robust LM-Error",
        "Robust LM-Lag-SDM",
        "Robust LM-WX",
        "Robust LM-Error-SDEM",
    ],
)

robust_results

	lm_samples	mean	median	credible_interval	bayes_pvalue	test_type	df	n_draws	k_wx	tr_MzWMzW_pair
Robust LM-Lag	[0.04470652999465163, 0.06039321001634612, 0.0...	0.052556	0.049121	(0.021080445691699274, 0.10259303831611782)	0.818673	bayesian_robust_lm_lag	1	20000	NaN	NaN
Robust LM-Error	[0.2592088833186119, 0.35016040229998335, 0.23...	0.304723	0.284806	(0.12222462329905084, 0.5948354055071119)	0.580937	bayesian_robust_lm_error	1	20000	NaN	NaN
Robust LM-Lag-SDM	[1.497849382284626, 1.7175227063115825, 1.8123...	2.003181	1.974842	(1.2853662789897897, 2.877608504387878)	0.156969	bayesian_robust_lm_lag_sdm	1	20000	2.0	19.966602
Robust LM-WX	[2.548722904936673, 2.5345872759229016, 0.1998...	3.448256	1.918213	(0.08724517028584271, 15.464011910652697)	0.178328	bayesian_robust_lm_wx	2	20000	2.0	NaN
Robust LM-Error-SDEM	[1.4121519102807134, 1.6051095505017356, 1.687...	1.844711	1.826631	(1.222114877277616, 2.565960814166178)	0.174400	bayesian_robust_lm_error_sdem	1	20000	2.0	19.966602

Classical `spreg` statistics on the same data¶

# Classical spreg LM tests
ols_spreg = SpregOLS(y, X)
lm_spreg = SpregLMtests(ols_spreg, w_spreg)

spreg_results = {
    "LM-Lag": lm_spreg.lml,
    "LM-Error": lm_spreg.lme,
    "LM-WX": lm_spreg.lmwx,
    "LM-SDM Joint": lm_spreg.lmspdurbin,
    "LM-SLX-Error Joint": lm_spreg.lmslxerr,
    "Robust LM-Lag": lm_spreg.rlml,
    "Robust LM-Error": lm_spreg.rlme,
    "Robust LM-Lag-SDM": lm_spreg.rlmdurlag,
    "Robust LM-WX": lm_spreg.rlmwx,
}

all_results = pd.concat([non_robust_results, robust_results])
comparison_rows = []
for bname in all_results.index:
    row = all_results.loc[bname]
    if bname in spreg_results:
        s_stat, s_pval = spreg_results[bname]
    else:
        s_stat, s_pval = np.nan, np.nan
    comparison_rows.append(
        {
            "bayes_mean": row["mean"],
            "spreg_stat": s_stat,
            "bayes_pvalue": row["bayes_pvalue"],
            "spreg_pvalue": s_pval,
        }
    )

comparison_df = pd.DataFrame(comparison_rows, index=all_results.index)
comparison_df

/home/runner/micromamba/envs/test/lib/python3.14/site-packages/spreg/diagnostics.py:620: ComplexWarning: Casting complex values to real discards the imaginary part
  ci_result = sqrt(max_eigval / min_eigval)

	bayes_mean	spreg_stat	bayes_pvalue	spreg_pvalue
LM-Lag	1.063234	0.886074	0.302479	0.346543
LM-Error	0.669052	1.341625	0.413382	0.246748
LM-WX	2.112357	0.616118	0.347782	0.734872
LM-WX-SEM	1.351305	NaN	0.508824	NaN
LM-SDM Joint	3.926935	1.957743	0.269463	0.581224
LM-SLX-Error Joint	2.113251	1.957743	0.549237	0.581224
Robust LM-Lag	0.052556	0.057811	0.818673	0.809990
Robust LM-Error	0.304723	0.513362	0.580937	0.473687
Robust LM-Lag-SDM	2.003181	1.341625	0.156969	0.246748
Robust LM-WX	3.448256	1.071669	0.178328	0.585181
Robust LM-Error-SDEM	1.844711	NaN	0.174400	NaN

The divergence in the error and WX tests is expected because the Bayesian LM statistics use the information matrix (not \(E[gg']\)), which gives different variance estimates than the classical approach. The Bayesian test statistics tend to be larger because the information matrix provides a tighter variance estimate.

Bayesian LM tests against classical spreg¶

Bayesian LM statistics¶

Classical spreg statistics on the same data¶

Bayesian LM tests against classical `spreg`¶

Classical `spreg` statistics on the same data¶