How to estimate spatial negative binomial models

Overdispersed counts with spatial dependence. SARNegBin is the reduced form and the default choice; SARNegBinStructural adds an explicit latent noise term and a sigma parameter, at a large cost in mixing.

Both sample by Pólya–Gamma Gibbs, which fit() selects for you; SARNegBin also has a NUTS path. Equations and constructor arguments are in Supported Models.

import contextlib
import io
import time

import arviz as az

from neighbayes.dgp import simulate_sar_negbin
from neighbayes.models import SARNegBin

/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

Simulate a known process

simulate_sar_negbin builds a square grid of side n, so n=25 gives 625 observations. As with the spatial logit models, \(\rho\) in a count model is weakly identified at small \(n\) — a few hundred units is the floor for the spatial parameter to be recoverable at all.

RHO_TRUE, ALPHA_TRUE = 0.4, 2.0

data = simulate_sar_negbin(n=25, rho=RHO_TRUE, seed=1)
y, X, W = data["y"], data["X"], data["W_graph"]

print(f"observations   : {len(y)}")
print(f"true parameters: {data['params_true']}")
print(f"mean count     : {y.mean():.2f}")
print(f"variance       : {y.var():.2f}   (>> mean: overdispersed, as intended)")

observations   : 625
true parameters: {'rho': 0.4, 'beta': array([1. , 0.6]), 'alpha': 2.0}
mean count     : 6.52
variance       : 64.00   (>> mean: overdispersed, as intended)

Fit the reduced-form model

def fit_timed(model, **kw):
    start = time.perf_counter()
    with (
        contextlib.redirect_stdout(io.StringIO()),
        contextlib.redirect_stderr(io.StringIO()),
    ):
        idata = model.fit(progressbar=False, **kw)
    return idata, time.perf_counter() - start


COMMON = dict(draws=1500, tune=800, chains=4, random_seed=1)

idata_gibbs, secs_gibbs = fit_timed(SARNegBin(y=y, X=X, W=W), **COMMON)
az.summary(idata_gibbs, var_names=["rho", "beta", "alpha"]).round(3)
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
rho 0.384 0.090 0.212 0.550 0.002 0.001 3533.0 4222.0 1.0
beta[x0] 1.055 0.154 0.768 1.343 0.003 0.002 3556.0 3854.0 1.0
beta[x1] 0.598 0.039 0.523 0.669 0.001 0.000 3117.0 4612.0 1.0
alpha 2.029 0.158 1.738 2.333 0.002 0.002 5764.0 4185.0 1.0 
print(f"rho   : {float(idata_gibbs.posterior['rho'].mean()):.3f}   (true {RHO_TRUE})")
print(
    f"alpha : {float(idata_gibbs.posterior['alpha'].mean()):.3f}   (true {ALPHA_TRUE})"
)

rho   : 0.384   (true 0.4)
alpha : 2.029   (true 2.0)

The dispersion parameter recovers tightly. \(\rho\) recovers to within its posterior standard deviation — read the interval, not the point.

Configure the sampler

gibbs_backend is "auto" (JAX when installed, else NumPy), "jax", or "numpy", as elsewhere in the package. The count families accept these additional options; anything else raises TypeError.

Option

SARNegBin

SARNegBinStructural

Meaning

slice_width

✅

—

initial slice interval for ρ

init_jitter

✅

—

dispersion of the chain starting points around the GLM warm start

n_rho_omega_cycles

✅

—

how many ρ/ω sub-sweeps per iteration

krylov_degree, krylov_dmax

✅

✅

shift-invert Krylov basis for the ρ conditional

krylov_reuse

✅

—

reuse that basis across draws

timeout

✅

—

abort a run exceeding this many seconds

lanczos_deg, n_probes

—

✅

stochastic log-determinant depth and probe count

pg_n_terms

—

✅

truncation of the Pólya–Gamma series

return_eta

—

✅

keep the latent field in the returned InferenceData

 
idata_tuned, _ = fit_timed(
    SARNegBin(y=y, X=X, W=W),
    gibbs_backend="numpy",
    slice_width=0.3,
    init_jitter=0.5,
    krylov_reuse=True,
    draws=800,
    tune=500,
    chains=2,
    random_seed=7,
)
print(
    f"rho = {float(idata_tuned.posterior['rho'].mean()):.3f}, "
    f"ess = {float(az.ess(idata_tuned, var_names=['rho'])['rho']):.0f}"
)

rho = 0.381, ess = 1005

Is the negative binomial the right likelihood?

The dispersion parameter answers this. \(\alpha \to \infty\) is the Poisson, so a posterior for alpha that is large and poorly determined says the data carry no conditional overdispersion for the extra parameter to fit.

alpha = idata_gibbs.posterior["alpha"]
print(f"alpha: mean {float(alpha.mean()):.2f}, sd {float(alpha.std()):.2f}")
print(
    f"94% HDI: {az.hdi(idata_gibbs, var_names=['alpha'], hdi_prob=0.94)['alpha'].values.round(2)}"
)

alpha: mean 2.03, sd 0.16
94% HDI: [1.74 2.33]

Here alpha is tightly determined near its true value of 2, so the negative binomial is earning its parameter.

Caution

Do not approximate a Poisson by fixing alpha large It is tempting to reach for a big fixed alpha when you want Poisson behaviour. That is the one regime where PG augmentation degenerates: the working precision \(E[\omega]\) diverges while the marginal Fisher information stays at \(\mu\), and effective sample size collapses — measured here as ESS on \(\beta\) falling from 57 at \(\alpha = 10\) to 3 at \(\alpha = 10^4\), with \(\hat R = 1.80\). Use a real Poisson sampler instead; the flow case is covered in How to estimate Poisson origin–destination flow models.

See also