How to set priors¶
Every model in neighbayes takes a priors= argument. This guide covers what
the defaults are, how to override them, how the package stops you from
overriding them wrongly, and the one case where leaving them alone actually
changes an answer you might report.
The short version: the defaults are weakly informative and scaled to your data, and you should usually leave them there. Override when you have real prior information, when you need to bound the spatial parameter, or when a referee asks for a sensitivity analysis.
import contextlib
import io
import pandas as pd
from neighbayes.dgp import simulate_sar
from neighbayes.models import OLS, SAR, SLX
from neighbayes.models.priors import SARPriors, SEMPriors
data = simulate_sar(n=20, rho=0.5, seed=3)
y, X, W = data["y"], data["X"], data["W_graph"]
print(f"n = {len(y)}, k = {X.shape[1]}")
n = 400, k = 2
Three ways to pass priors¶
priors= accepts None, a plain dict, or a typed dataclass. All three end up
in the same place — the dict and None forms are coerced to the dataclass at
construction.
m_default = SAR(y=y, X=X, W=W) # None: defaults
m_dict = SAR(y=y, X=X, W=W, priors={"rho_lower": 0.0}) # dict of overrides
m_typed = SAR(y=y, X=X, W=W, priors=SARPriors(rho_lower=0.0)) # typed dataclass
print(m_dict.priors_obj.rho_lower, m_typed.priors_obj.rho_lower)
0.0 0.0
Prefer the typed form in code you intend to keep. It gives you autocomplete on the supported keys, and a type checker catches a bad field before you run anything. The dict form is convenient in a notebook.
Inspect what a model resolved¶
Two public attributes. priors_obj is the dataclass; priors is a dict of the
keys that carry a concrete value.
print(f"class : {type(m_default.priors_obj).__name__}")
print(f"object: {m_default.priors_obj}")
print()
print(f"dict : {m_default.priors}")
class : SARPriors
object: SARPriors(beta_mu=None, beta_sigma=None, sigma2_alpha=2.0, sigma2_beta=None, sigma_sigma=10.0, nu=4.0, rho_lower=-1.0, rho_upper=1.0)
dict : {'sigma2_alpha': 2.0, 'sigma_sigma': 10.0, 'nu': 4.0, 'rho_lower': -1.0, 'rho_upper': 1.0}
type(model.priors_obj) is also how you discover which dataclass a
given model wants — there are 39 of them and the mapping follows the model
hierarchy, so looking it up beats guessing.
What the defaults actually are¶
Three fields come back as None above: beta_mu, beta_sigma, and
sigma2_beta. None here does not mean “flat” — it means resolved from the
data at construction.
Coefficients: a data-scaled weakly-informative normal¶
Following Gelman et al. (2008), each coefficient gets a normal prior scaled to the response and to its own covariate:
Column |
\(\mu\) |
\(\sigma\) |
|---|---|---|
intercept |
\(\overline{y}\) |
\(2.5\,\mathrm{sd}(y)\) |
slope \(j\) |
\(0\) |
\(2.5\,\mathrm{sd}(y) / \mathrm{sd}(x_j)\) |
Dividing by \(\mathrm{sd}(x_j)\) is what makes this behave sensibly regardless of the units your covariates happen to be in — a prior that is weak for income in dollars would be crushing for income in millions.
Noise variance: a scale-aware inverse gamma¶
\(\sigma^2 \sim \mathrm{InverseGamma}(\alpha, \beta)\) with \(\alpha = 2\) (finite mean, weakly informative) and \(\beta\) resolved to \(\mathrm{Var}(y)\), so the prior mean sits at roughly the observed variance. Following LeSage (2009), the conjugate inverse gamma is what keeps the \(\sigma^2\) Gibbs block in closed form — the Gaussian models use it for both NUTS and Gibbs so that the two paths target exactly the same posterior.
Spatial parameter: bounded uniform¶
\(\rho \sim \mathrm{Uniform}(\texttt{rho\_lower}, \texttt{rho\_upper})\), default
\((-1, 1)\); SEMPriors spells the same fields lam_lower / lam_upper.
Overriding¶
Set a scalar and it broadcasts to every coefficient; set a vector of length \(p\) and it applies per coefficient, in design-matrix column order.
m_scalar = SAR(y=y, X=X, W=W, priors={"beta_sigma": 5.0})
m_vector = SAR(y=y, X=X, W=W, priors={"beta_mu": [0.0, 1.0], "beta_sigma": [10.0, 0.5]})
print("scalar :", m_scalar.priors_obj.beta_sigma)
print("vector :", m_vector.priors_obj.beta_mu, m_vector.priors_obj.beta_sigma)
scalar : 5.0
vector : [0.0, 1.0] [10.0, 0.5]
Bounding the spatial parameter¶
The common substantive override. If theory rules out negative spatial dependence — diffusion processes, contagion — say so, rather than letting the sampler spend draws in a region you would not believe anyway.
m_positive = SAR(y=y, X=X, W=W, priors={"rho_lower": 0.0})
idata_positive = m_positive.fit(
draws=1000, tune=500, chains=4, random_seed=1, progressbar=False
)
print(f"rho posterior minimum: {float(idata_positive.posterior['rho'].min()):.4f}")
print("(the prior bound is a hard constraint, not a nudge)")
rho posterior minimum: 0.3603
(the prior bound is a hard constraint, not a nudge)
/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)
Mistakes fail loudly¶
A misspelled key is silently ignored in most libraries. Here it raises at construction, with the allowed keys in the message.
try:
SAR(y=y, X=X, W=W, priors={"rho_lo": 0.0})
except TypeError as err:
print(f"TypeError: {err}")
TypeError: Unknown prior key(s) for SARPriors: ['rho_lo']. Allowed keys: ['beta_mu', 'beta_sigma', 'nu', 'rho_lower', 'rho_upper', 'sigma2_alpha', 'sigma2_beta', 'sigma_sigma'].
try:
SAR(y=y, X=X, W=W, priors=SEMPriors())
except TypeError as err:
print(f"TypeError: {err}")
TypeError: priors must be None, a dict, or a SARPriors instance; got SEMPriors.
The second case is why the typed form is worth using: SEMPriors and
SARPriors differ only in whether the spatial bounds are named lam_* or
rho_*, which is exactly the kind of thing that is easy to get wrong and hard
to notice.
The case where the default matters¶
Warning
Wide priors are not “uninformative” for model comparison
Setting beta_sigma to a large value to be “objective” reproduces the legacy
near-improper prior, and it will change your Bayes factors by many nats. This
is the Bartlett–Lindley paradox: the marginal likelihood integrates the
likelihood against the prior, so spreading the prior over a wider region
dilutes the evidence for the model carrying more coefficients — regardless of
what the data say.
Compare a non-spatial OLS against SLX, which adds spatially lagged
covariates, under the default prior and under a deliberately diffuse one:
FIT_KW = dict(draws=1000, tune=500, chains=4, random_seed=1, progressbar=False)
def log_bf_slx_vs_ols(priors):
from neighbayes.diagnostics.bayesfactor import bayes_factor_compare_models
with (
contextlib.redirect_stdout(io.StringIO()),
contextlib.redirect_stderr(io.StringIO()),
):
ols = OLS(y=y, X=X, priors=priors)
ols.fit(**FIT_KW)
slx = SLX(y=y, X=X, W=W, priors=priors)
slx.fit(**FIT_KW)
table = bayes_factor_compare_models(
[ols, slx], model_labels=["OLS", "SLX"], log=True
)
return float(table.loc["SLX", "OLS"])
rows = {
"default (data-scaled)": log_bf_slx_vs_ols(None),
"beta_sigma = 1e6": log_bf_slx_vs_ols({"beta_sigma": 1e6}),
}
pd.DataFrame({"log BF (SLX vs OLS)": rows}).round(3)
| log BF (SLX vs OLS) | |
|---|---|
| default (data-scaled) | 29.196 |
| beta_sigma = 1e6 | 17.803 |
Same data, same models, same sampler. The only thing that changed is how much prior mass sits in regions the data rule out — and the evidence swings by more than ten log units.
The direction is the tell. SLX carries the extra spatially-lagged covariates,
so it is the model that pays when the prior is spread thin: widening
beta_sigma dilutes the likelihood mass that its additional coefficients can
claim, and its Bayes factor falls. Push the prior wide enough and the richer
model loses on prior width alone, whatever the data say. That is the paradox,
and it does not go away by making the prior wider still — it gets worse.
Neither number is “wrong”; they answer different questions. But the diffuse one is not the neutral choice it looks like. If you compare models with Bayes factors, keep the default and say so, or report the comparison across a range of prior widths.
Sensitivity as routine practice¶
Refitting across a grid of prior widths costs little and answers the question a reader will ask. Do it on the smallest dataset you plan to report, because that is where a prior can still move the answer — at \(n = 400\) above, nothing reasonable would.
small = simulate_sar(n=7, rho=0.5, seed=3) # 7 x 7 grid = 49 observations
ys, Xs, Ws = small["y"], small["X"], small["W_graph"]
rows = []
for width in (0.25, 1.0, 5.0, None):
label = "default" if width is None else f"beta_sigma={width:g}"
priors = None if width is None else {"beta_sigma": width}
idata = SAR(y=ys, X=Xs, W=Ws, priors=priors).fit(**FIT_KW)
rows.append(
pd.Series(
{
"rho mean": float(idata.posterior["rho"].mean()),
"beta[0] mean": float(idata.posterior["beta"][..., 0].mean()),
"beta[1] mean": float(idata.posterior["beta"][..., 1].mean()),
},
name=label,
)
)
pd.DataFrame(rows).round(3)
/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)
/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)
/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)
/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)
| rho mean | beta[0] mean | beta[1] mean | |
|---|---|---|---|
| beta_sigma=0.25 | 0.459 | 1.229 | 1.315 |
| beta_sigma=1 | 0.459 | 0.937 | 2.058 |
| beta_sigma=5 | 0.459 | 0.917 | 2.106 |
| default | 0.459 | 0.917 | 2.107 |
Two things to read off this table.
The coefficients shrink when the prior is tight and then stop. At
beta_sigma=0.25 the slope is pulled hard toward zero — the prior is narrow
relative to the scale of the data. By beta_sigma=5 the posterior has stopped
caring, and the default sits on that plateau. Reaching the plateau is what tells
you an estimate is data-driven rather than an artefact of the prior; it is the
thing worth reporting, and it is cheap to check.
rho does not move at all. beta_sigma is the prior on the coefficients;
the spatial parameter has its own, and shrinking \(\beta\) does not touch it. If
you want to probe sensitivity in \(\rho\), vary rho_lower / rho_upper —
varying the wrong hyperparameter and finding nothing is not a sensitivity
analysis.
See also¶
How to choose and configure the Gibbs sampler — the σ² prior is what keeps that block conjugate
How to run Bayesian LM specification tests — model choice without Bayes factors, and so without the prior-width sensitivity
Supported Models — the full catalogue and each family’s parameters