How to fit spatial logit models¶

Binary outcomes with spatial dependence. SARLogit is the reduced form and the one to reach for; SARLogitStructural is the same model parameterised without inverting \((I - \rho W)\); SEMLogit puts the spatial structure in the errors instead, as a nuisance rather than a channel for spillovers.

All three sample by Pólya–Gamma Gibbs and have no NUTS path — fit(sampler="nuts") raises. Equations and constructor arguments are in Supported Models; if you have not yet decided between a lag and an error specification, start with the Bayesian LM specification tests.

import arviz as az
import numpy as np
import pandas as pd

from neighbayes.dgp import simulate_sar_logit, simulate_sem_logit
from neighbayes.models import SARLogit, SARLogitStructural, SEMLogit

/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_logit builds a square grid of side n, so n=35 gives \(35^2 = 1225\) observations. That matters: binary outcomes carry far less information per observation than continuous ones, and \(\rho\) in a spatial logit is weakly identified below roughly a thousand units. Fitting one of these to \(n = 100\) and concluding the sampler is broken is a common mistake.

RHO_TRUE, BETA_TRUE = 0.5, np.array([-0.3, 0.5])

data_sar = simulate_sar_logit(n=35, rho=RHO_TRUE, beta=BETA_TRUE, seed=42)
y_sar, X_sar, W_sar = data_sar["y"], data_sar["X"], data_sar["W_graph"]

print(f"observations : {len(y_sar)}  (35 x 35 queen grid)")
print(f"covariates   : {X_sar.shape[1]}")
print(f"outcome      : {int(y_sar.sum())} ones, {int(len(y_sar) - y_sar.sum())} zeros")

observations : 1225  (35 x 35 queen grid)
covariates   : 2
outcome      : 463 ones, 762 zeros

Fit the reduced-form model¶

The sweep has four blocks, only one of which is non-conjugate:

1. \(\omega \mid \eta\) — Pólya–Gamma augmentation

2. \(\eta \mid \omega, \rho, \beta\) — spatial-normal draw

3. \(\beta \mid \eta, \rho\) — conjugate normal

4. \(\rho \mid \beta, \omega, y\) — collapsed one-dimensional slice, with \(\eta\) integrated out

Collapsing block 4 is what makes this mix. Conditioning \(\rho\) on the current \(\eta\) draw instead couples the two so tightly that the chain crawls.

model_sar = SARLogit(y=y_sar, X=X_sar, W=W_sar)
idata_sar = model_sar.fit(
    draws=1500, tune=800, chains=4, random_seed=42, progressbar=False
)

az.summary(idata_sar, var_names=["rho", "beta"]).round(3)

print(f"rho  : {float(idata_sar.posterior['rho'].mean()):.3f}   (true {RHO_TRUE})")
print(
    f"beta : {idata_sar.posterior['beta'].mean(dim=['chain', 'draw']).values.round(3)}"
    f"   (true {BETA_TRUE})"
)

rho  : 0.432   (true 0.5)
beta : [-0.287  0.437]   (true [-0.3  0.5])

\(\rho\) comes back with a wide credible interval even at \(n = 1225\). That is the model, not the sampler — check ess_bulk and r_hat above to confirm the chain is healthy, then read the width as what the data support.

Read the spillovers¶

In a linear SAR model \(\beta\) is a marginal effect. Here it is not. The link is nonlinear, so a change in \(x_i\) moves \(y_j\) through both the spatial multiplier \((I - \rho W)^{-1}\) and the logistic derivative, which varies with \(\eta_i\). spatial_effects() does the averaging and returns the LeSage–Pace decomposition — direct, indirect, total — on the probability scale.

model_sar.spatial_effects().round(4)

Report these, not the raw coefficients, whenever the question is “how much does \(x\) change the outcome”. The gap between direct and the coefficient is the part of the effect that the spatial multiplier feeds back onto the unit itself.

Fit the error model¶

SEMLogit swaps the lag for spatially correlated disturbances. The collapsed density for \(\lambda\) carries an extra term the SAR version does not: the correction \(-\tfrac{1}{2}(X\beta)^\top A_\lambda^\top A_\lambda (X\beta)\) depends on \(\lambda\), whereas the SAR analogue is constant in \(\rho\) and drops out.

data_sem = simulate_sem_logit(n=35, lam=0.5, beta=BETA_TRUE, seed=42)
y_sem, X_sem, W_sem = data_sem["y"], data_sem["X"], data_sem["W_graph"]

model_sem = SEMLogit(y=y_sem, X=X_sem, W=W_sem)
idata_sem = model_sem.fit(
    draws=1500, tune=800, chains=4, random_seed=42, progressbar=False
)

az.summary(idata_sem, var_names=["lam", "beta"]).round(3)

Reduced form or structural form¶

SARLogitStructural samples the same SAR-logit model without inverting \((I - \rho W)\) into the mean. Fit both on the same data and they should agree; what differs is the sampler’s internals.

model_struct = SARLogitStructural(y=y_sar, X=X_sar, W=W_sar)
idata_struct = model_struct.fit(
    draws=1500, tune=800, chains=4, random_seed=42, progressbar=False
)

pd.DataFrame(
    {
        "reduced (SARLogit)": [
            float(idata_sar.posterior["rho"].mean()),
            float(idata_sar.posterior["rho"].std()),
            float(az.ess(idata_sar, var_names=["rho"])["rho"]),
        ],
        "structural": [
            float(idata_struct.posterior["rho"].mean()),
            float(idata_struct.posterior["rho"].std()),
            float(az.ess(idata_struct, var_names=["rho"])["rho"]),
        ],
    },
    index=["rho mean", "rho sd", "rho ess"],
).round(3)

The two parameterisations are algebraically the same model — multiply the structural form through by \((I - \rho W)^{-1}\) and you have the reduced one — so the check is whether the posteriors are consistent, not whether the means match to three decimals. Compare the gap between the means against the posterior standard deviations: a difference well inside one standard deviation is two independent implementations agreeing, and that is worth more than either fit alone.

They are not interchangeable in practice. The structural sampler conditions on \(\eta\) rather than integrating through the inverse, which changes both the cost per sweep and the width of the \(\rho\) posterior it reports. Prefer SARLogit unless you have a reason to want the latent field itself.

Configure the sampler¶

gibbs_backend behaves as it does everywhere in the package: "auto" (the default — JAX when installed, else NumPy), "jax", or "numpy".

Beyond the shared fit() arguments, the two families accept different options. Anything else raises TypeError rather than being silently ignored.

idata_tuned = SARLogit(y=y_sar, X=X_sar, W=W_sar).fit(
    draws=800,
    tune=500,
    chains=2,
    random_seed=7,
    progressbar=False,
    gibbs_backend="numpy",
    slice_width=0.4,
    krylov_reuse=True,
)
print(
    f"rho = {float(idata_tuned.posterior['rho'].mean()):.3f}, "
    f"ess = {float(az.ess(idata_tuned, var_names=['rho'])['rho']):.0f}"
)

rho = 0.433, ess = 1277

Fitted probabilities¶

fitted_probabilities() returns \(P(y_i = 1)\) at the posterior mean. The two model types differ in what the spatial parameter does to the mean:

• SARLogit — \(\hat\eta = (I - \hat\rho W)^{-1} X\hat\beta\); the multiplier is in the mean

• SEMLogit — \(\hat\eta = X\hat\beta\); the spatial parameter shapes the variance, leaving the mean alone

probs_sar = model_sar.fitted_probabilities()
probs_sem = model_sem.fitted_probabilities()

pd.DataFrame(
    {
        "SARLogit": [probs_sar.mean(), probs_sar.min(), probs_sar.max()],
        "SEMLogit": [probs_sem.mean(), probs_sem.min(), probs_sem.max()],
    },
    index=["mean", "min", "max"],
).round(3)

What to check before trusting the output¶

print("groups:", idata_sar.groups())
az.summary(idata_sar, var_names=["rho"]).round(3)

groups: ['posterior', 'observed_data']

See also¶

• How to run Bayesian LM specification tests — deciding between lag and error before you fit

• How to set priors — SARLogitPriors and SEMLogitPriors

• How to estimate spatial negative binomial models — the same Pólya–Gamma machinery for counts