neighbayes.models.SARZINB

class neighbayes.models.SARZINB(formula=None, data=None, y=None, X=None, Z=None, W=None, W_sel=None, priors=None, logdet_method=None, robust=False, **kwargs)[source]

Bayesian zero-inflated SAR Negative Binomial with PG-Gibbs sampler.

Parameters:
formula : str, optional

Wilkinson-style formula for the count equation, e.g. "y ~ x1 + x2". Requires data.

data : pandas.DataFrame or geopandas.GeoDataFrame, optional

Data source for formula mode.

y : array-like, optional

Non-negative integer counts of shape (n,). Required in matrix mode.

X : array-like, optional

Count covariate matrix of shape (n, k). Required in matrix mode.

Z : array-like, optional

Selection covariate matrix of shape (n, p). If None, defaults to X (same covariates for both equations).

W : libpysal.graph.Graph or scipy.sparse matrix

Spatial weights for the count equation, shape (n, n).

W_sel : libpysal.graph.Graph or scipy.sparse matrix, optional

Spatial weights for the selection equation. If None, uses W (same weights for both equations).

priors : dict, optional

Override default priors. Supported keys:

  • gamma_mu (float, default 0.0): Normal prior mean for γ.

  • gamma_sigma (float, default 1e6): Normal prior std for γ.

  • lam_lower (float, default -0.999): Lower bound for λ.

  • lam_upper (float, default 0.999): Upper bound for λ.

  • beta_mu (float, default 0.0): Normal prior mean for β.

  • beta_sigma (float, default 10.0): Normal prior std for β.

  • rho_lower (float, default -0.999): Lower bound for ρ.

  • rho_upper (float, default 0.999): Upper bound for ρ.

  • alpha_sigma (float, default 2.5): Half-Normal scale for α.

  • alpha_nu (float, default 3.0): Half-Normal ν for α.

logdet_method : str, optional

How to compute log|I − ρW|. None (default) auto-selects.

robust : bool, default False

Not supported. Raises NotImplementedError if True.

Notes

The model is fit with a custom 9-block Gibbs sampler that composes the SAR-logit blocks (ω^sel, η^sel, γ, λ), the zero-allocation block (z), and the reduced-form SAR-NB blocks (ω^cnt, β, ρ, α). No PyMC model is constructed; calling _build_pymc_model raises NotImplementedError.

__init__(formula=None, data=None, y=None, X=None, Z=None, W=None, W_sel=None, priors=None, logdet_method=None, robust=False, **kwargs)[source]

Methods

__init__([formula, data, y, X, Z, W, W_sel, ...])

corridor_probabilities()

Posterior-mean corridor activation probabilities.

fit([draws, tune, chains, target_accept, ...])

Draw samples from the posterior.

fitted_values()

Return fitted values at posterior mean parameters.

residuals()

Return residuals y - fitted_values.

spatial_diagnostics()

Run Bayesian LM specification tests and return a summary table.

spatial_diagnostics_decision([alpha, format])

Return a model-selection decision from Bayesian LM test results.

spatial_effects([equation, ...])

Compute Bayesian inference for direct, indirect, and total impacts.

summary([var_names])

Return posterior summary table.

zero_attribution()

Decompose observed zeros into structural vs sampling zeros.

Attributes

inference_data

Return the ArviZ InferenceData from the most recent fit.

pymc_model

Return the PyMC model object built for the most recent fit.

corridor_probabilities()[source]

Posterior-mean corridor activation probabilities.

Returns π_i = logit⁻¹(η_i^sel) at posterior means of λ and γ.

Returns:

pi – Fitted activation probabilities.

Return type:

ndarray of shape (n,)

fit(draws=2000, tune=1000, chains=4, target_accept=None, random_seed=None, progressbar=True, sampler=None, gibbs_backend='auto', thin=1, n_jobs=-1, idata_kwargs=None, **sample_kwargs)[source]

Draw samples from the posterior.

Dispatches to this model’s Gibbs sampler (sampler="gibbs") or NUTS (sampler="nuts"). When sampler is None (default), Gibbs is used if the model has a registered Gibbs sampler, otherwise NUTS.

Parameters:
draws : int

Post-warmup draws, warmup steps, and number of chains.

tune : int

Post-warmup draws, warmup steps, and number of chains.

chains : int

Post-warmup draws, warmup steps, and number of chains.

target_accept : float, optional

Target acceptance rate for NUTS. NUTS-only: passing it with the Gibbs sampler raises TypeError. Defaults to 0.9 for NUTS.

random_seed : int, optional

Seed for reproducibility.

progressbar : bool, default True

Show progress bar(s) during sampling.

sampler : {"gibbs", "nuts", None}, default None

Sampling method. None auto-selects Gibbs when this model has one, else NUTS.

gibbs_backend : {"auto", "jax", "numpy"}, default "auto"

Execution backend for the Gibbs sampler. "auto" uses JAX when installed and supported by the family, else NumPy. Ignored for NUTS.

thin : int, default 1

Keep every thin-th post-warmup Gibbs draw (Gibbs only).

n_jobs : int, default -1

Parallel workers for the NumPy Gibbs path (Gibbs only).

idata_kwargs : dict, optional

{"log_likelihood": True} stores the complete Jacobian-corrected pointwise log-likelihood that az.loo / az.waic / az.compare need, for Gibbs and NUTS alike. Off by default, as in PyMC: it holds one value per draw, chain, and observation (16 GB at n = 250,000 with 4 × 2,000 draws). For NUTS the dict is also passed to pm.sample.

**sample_kwargs

For NUTS, forwarded to pm.sample (nuts_sampler=...); for Gibbs, the family’s declared options (an unsupported key raises).

Return type:

arviz.InferenceData

fitted_values()[source]

Return fitted values at posterior mean parameters.

Returns:

Posterior-mean fitted values (on the model’s native scale; fixed-effects-transformed for panel models).

Return type:

np.ndarray

property inference_data : arviz.data.inference_data.InferenceData | None[source]

Return the ArviZ InferenceData from the most recent fit.

Returns:

The inference data object, or None if the model has not been fit yet.

Return type:

arviz.InferenceData or None

property pymc_model : pymc.model.core.Model | None[source]

Return the PyMC model object built for the most recent fit.

For Gibbs-fitted models the PyMC model is not constructed during sampling; it is built lazily on first access so that downstream consumers (e.g. bridge sampling for marginal likelihoods) can evaluate logp and the prior under the same model definition used by the NUTS path.

Returns:

The model object used by fit(), or None if the instance has not been fit yet.

Return type:

pymc.Model or None

residuals()[source]

Return residuals y - fitted_values.

Returns:

Residual vector y - fitted_values on the same scale as fitted_values().

Return type:

np.ndarray

spatial_diagnostics()[source]

Run Bayesian LM specification tests and return a summary table.

Looks up the diagnostic suite registered for this model class and calls each test function on this fitted model, collecting the results into a tidy DataFrame. The set of tests depends on the model type — for example, an OLS model runs LM-Lag, LM-Error, LM-SDM-Joint, and LM-SLX-Error-Joint, while an SAR model runs LM-Error, LM-WX, and Robust-LM-WX. Panel models run the Panel--prefixed analogues (e.g. Panel-LM-Lag).

Requires the model to have been fit (.fit() called) and a spatial weights matrix W to have been supplied at construction time.

Returns:

DataFrame indexed by test name with columns:

Column

Description

statistic

Posterior mean of the LM statistic

median

Posterior median of the LM statistic

df

Degrees of freedom for the \(\chi^2\) reference

p_value

Bayesian p-value: 1 - chi2.cdf(mean, df)

ci_lower

Lower bound of 95% credible interval (2.5%)

ci_upper

Upper bound of 95% credible interval (97.5%)

The DataFrame has attrs["model_type"] (class name) and attrs["n_draws"] (total posterior draws) metadata.

Return type:

pandas.DataFrame

Raises:
  • RuntimeError – If the model has not been fit yet.

  • ValueError – If no spatial weights matrix W was supplied.

See also

spatial_diagnostics_decision

Model-selection decision based on the test results.

spatial_effects

Posterior inference for direct/indirect/total impacts.

Examples

>>> ols = OLS(formula="price ~ income + crime", data=df, W=w)
>>> ols.fit()
>>> ols.spatial_diagnostics()
                 statistic  median  df  p_value  ci_lower  ci_upper
LM-Lag                3.21    2.98   1    0.073      0.12      8.54
LM-Error              5.67    5.34   1    0.017      0.34     12.10
LM-SDM-Joint          7.89    7.12   4    0.096      1.23     18.32
LM-SLX-Error-Joint    6.45    5.98   4    0.168      0.89     15.67
spatial_diagnostics_decision(alpha=0.05, format='graphviz')[source]

Return a model-selection decision from Bayesian LM test results.

Implements the decision tree from Koley and Bera [2024] (the Bayesian analogue of the classical stge_kb procedure in Anselin et al. [1996]). Panel models use the Panel--prefixed test analogues and the panel decision specs, following Elhorst [2014]. The decision logic depends on the current model type and the pattern of significant tests:

From OLS (6-test decision tree):

  1. If only LM-Lag is significant → SAR.

  2. If only LM-Error is significant → SEM.

  3. If both are significant → use the Anselin–Florax / Koley–Bera robust pair: Robust-LM-Lag → SAR, Robust-LM-Error → SEM, both → SARAR. If neither robust test is significant, fall back to the lower raw p-value.

  4. If neither naive test is significant → OLS.

From SAR (3-test decision tree):

  • LM-Error significant → SARAR; LM-WX significant → SDM; Robust-LM-WX significant → SDM.

From SEM (2-test decision tree):

  • LM-Lag significant → SARAR; LM-WX significant → SDEM.

From SLX (4-test decision tree):

  • Robust-LM-Lag-SDM significant → SDM; Robust-LM-Error-SDEM significant → SDEM; both → MANSAR; neither → SLX.

From SDM: LM-Error-SDM significant → MANSAR; else SDM.

From SDEM: LM-Lag-SDEM significant → MANSAR; else SDEM.

Parameters:
alpha : float, default 0.05

Significance level for the Bayesian p-values.

format : {"graphviz", "ascii", "model"}, default "graphviz"

Output format. "model" returns the recommended-model name string. "ascii" returns an indented box-drawing rendering of the full decision tree with the chosen path highlighted. "graphviz" returns a graphviz.Digraph object that renders inline in Jupyter; if the optional graphviz package is not installed a UserWarning is issued and the ASCII rendering is returned instead.

Returns:

Recommended model name when format="model", an ASCII tree string when format="ascii", or a graphviz.Digraph when format="graphviz" (with ASCII fallback on missing dep).

Return type:

str or graphviz.Digraph

See also

spatial_diagnostics

Compute the Bayesian LM test statistics.

References

Koley and Bera [2024], Anselin et al. [1996], Elhorst [2014]

spatial_effects(equation='count', return_posterior_samples=False)[source]

Compute Bayesian inference for direct, indirect, and total impacts.

The ZINB-SAR model has two spatial lag equations — a count equation with parameter ρ and a selection (logit) equation with parameter λ — each with its own LeSage–Pace impact decomposition. Use the equation parameter to select which equation’s impacts to report.

Parameters:
equation : {"count", "selection"}, default "count"

Which equation to compute impacts for.

  • "count": impacts of X on E[y | d=1] through (I − ρW)⁻¹ Xβ.

  • "selection": impacts of Z on P(d=1) on the log-odds scale through (I − λW_sel)⁻¹ Zγ.

return_posterior_samples : bool, default False

If True, return a (DataFrame, dict) tuple where the dict contains the full posterior draws.

Returns:

Impact summary table, optionally with posterior draws.

Return type:

pd.DataFrame or tuple of (pd.DataFrame, dict)

summary(var_names=None, **kwargs)[source]

Return posterior summary table.

Parameters:
var_names : list, optional

Variable names to include in the summary.

**kwargs

Additional arguments passed to arviz.summary().

Returns:

Posterior summary statistics.

Return type:

pandas.DataFrame

zero_attribution()[source]

Decompose observed zeros into structural vs sampling zeros.

For each observation with y_i = 0, computes the posterior probability that the zero is structural (d_i = 0) versus sampling (d_i = 1 but NB draw was zero).

Returns:

structural_probndarray of shape (n_zero,)

P(d_i = 0 | y_i = 0) for each zero observation.

sampling_probndarray of shape (n_zero,)

P(d_i = 1, NB zero | y_i = 0) for each zero observation.

zero_indicesndarray of shape (n_zero,)

Indices of zero observations.

Return type:

dict with keys