neighbayes.models.SEM¶

class neighbayes.models.SEM(formula=None, data=None, y=None, X=None, W=None, priors=None, logdet_method=None, robust=False, w_vars=None, logdet_refit=True, logdet_refit_pad_sd=10.0, logdet_aaa_check=True, logdet_probe_check=True)[source]¶

Bayesian Spatial Error Model.

Spatial dependence enters through the disturbance via the autoregressive parameter \(\lambda\):

\[y = X\beta + u, \quad u = \lambda Wu + \varepsilon, \quad \varepsilon \sim N(0, \sigma^2 I).\]

The likelihood includes the spatial Jacobian \(\log|I - \lambda W|\) so that posterior inference on \(\lambda\) is exact.

Parameters:¶

formula : str, optional¶

Wilkinson-style formula, e.g. "y ~ x1 + x2". Requires data. Intercept is included by default; suppress with "y ~ x - 1".

data : pandas.DataFrame or geopandas.GeoDataFrame, optional¶

Data source for formula mode.

y : array-like, optional¶

Dependent variable of shape (n,). Required in matrix mode.

X : array-like or pandas.DataFrame, optional¶

Design matrix. Required in matrix mode. DataFrame columns are preserved as feature names.

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

Spatial weights of shape (n, n). Accepts a libpysal.graph.Graph or any scipy.sparse matrix. The legacy libpysal.weights.W object is not accepted; pass w.sparse or libpysal.graph.Graph.from_W(w). Should be row-standardized; a UserWarning is raised otherwise.

priors : dict, optional¶

Override default priors. Supported keys:

lam_lower (float, default -1.0): Lower bound of the Uniform prior on \(\lambda\).
lam_upper (float, default 1.0): Upper bound of the Uniform prior on \(\lambda\).
beta_mu (float, default 0.0): Normal prior mean for \(\beta\).
beta_sigma (float, default 1e6): Normal prior std for \(\beta\).
sigma2_alpha (float, default 2.0): Shape of the InverseGamma prior on \(\sigma^2\).
sigma2_beta (float, default Var(y)): Scale of the InverseGamma prior on \(\sigma^2\).
nu (float, default 4.0): Fixed Student-t degrees of freedom (only used when robust=True).

logdet_method : str, optional¶

How to compute \(\log|I - \lambda W|\). None (default) auto-selects by size: "eigenvalue" for n <= 500; for 500 < n <= 60000, "cheb_cholesky" (exact, sparse Cholesky at Chebyshev nodes) when W is symmetric else "aaa" (AAA rational approximation); "cheb_stochastic" for n > 60000. Explicit opt-ins: "chebyshev" (Barry-Pace) and "slq" (stochastic Lanczos quadrature).

logdet_refit : bool, default True¶

Rebuild the log-determinant interpolant halfway through warmup, on the range the chains have found rather than the interval implied by the prior. A warmup posterior is typically one to two orders of magnitude narrower than the prior, which needs far fewer interpolation nodes and drives the approximation error over the posterior’s support down to the factorization’s roundoff floor. Applies to "cheb_cholesky", "lu_cheb", "aaa" and "chol_aaa"; ignored otherwise.

The interpolant is only valid on its interval, so the refit window becomes the sampler’s support. The window is padded by logdet_refit_pad_sd warmup standard deviations, recorded in idata.attrs["logdet_refit_window"], and a warning is raised if the retained draws ever reach an edge the refit introduced.

logdet_refit_pad_sd : float, default 10.0¶

Padding for the refit window, in warmup posterior standard deviations. At the default the truncated tail is ~1e-23 under normality, and the padding costs a node or two at most.

logdet_aaa_check : bool, default True¶

For the AAA methods ("aaa", "chol_aaa"), set the number of exact factorizations from where the posterior lies. Warmup starts on a 14-node fit; halfway through, nodes are added only if the warmup posterior lies closer to a singularity of the Jacobian than four times the distance at which the fit’s poles resolve it. The count used is recorded in idata.attrs["logdet_aaa_nodes"]. Off, or on a path without a warmup midpoint, the count is fixed by the prior interval: 14 nodes within |ρ| ≤ 0.9, 18 otherwise.

logdet_probe_check : bool, default True¶

For "cheb_stochastic", set the number of Hutchinson probes from where the posterior lies. Warmup starts on 50 probes; halfway through, the probes’ own spread prices the bias they leave in the posterior mean of the spatial parameter, and the pool grows, to at most 200, until that bias is expected to stay under 0.0225 posterior sd. The count used is recorded in idata.attrs["logdet_probes"], with a warning if the cap is reached first. Probes cost setup time only; the cost of each draw does not depend on how many there are.

robust : bool, default False¶

If True, replace the Normal disturbance with Student-t. See Robust regression below.

Notes

Because spatial dependence enters only through the disturbance, direct effects equal \(\beta\) and indirect effects are zero.

Robust regression

When robust=True, the spatially-filtered innovation is Student-t:

\[\varepsilon = (I - \lambda W)(y - X\beta) \sim t_\nu(0, \sigma^2 I)\]

where \(\nu\) is a fixed hyperparameter set by priors={"nu": value} (default 4, LeSage’s rval); larger values approach the Normal. Values must exceed 2 so the variance exists.

__init__(formula=None, data=None, y=None, X=None, W=None, priors=None, logdet_method=None, robust=False, w_vars=None, logdet_refit=True, logdet_refit_pad_sd=10.0, logdet_aaa_check=True, logdet_probe_check=True)[source]¶

Methods

`__init__`([formula, data, y, X, W, priors, ...])
`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`([return_posterior_samples])	Compute Bayesian inference for direct, indirect, and total impacts.
`summary`([var_names])	Return posterior summary table.

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.

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.