How to fit spatial tobit models

Outcomes censored at a threshold — spending that cannot go below zero, demand recorded only above a reporting limit. SARTobit puts the spatial lag on the latent outcome, SEMTobit puts the spatial structure in the errors, and SDMTobit adds spatially lagged covariates.

All three are NUTS-only: none has a registered Gibbs sampler, so fit() uses NUTS and target_accept applies. Equations and constructor arguments are in Supported Models.

Note

For binary outcomes, use the Pólya–Gamma logit classes SARLogit, SEMLogit and SARLogitStructural are the modern route to spatial binary outcomes, and they have a conjugate Gibbs sampler — see How to fit spatial logit models. The package’s SARProbit is a region-random-effects specification rather than the standard spatial probit, so it is not the class to reach for by default.

import contextlib
import io
import time

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

from neighbayes.dgp import simulate_sar_tobit, simulate_sdm_tobit, simulate_sem_tobit
from neighbayes.models import SARTobit, SDMTobit, SEMTobit

FIT = dict(draws=1000, tune=800, chains=4, random_seed=1)


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

Fit a censored outcome

censoring is the threshold, and it must match your data — pass the value below which observations are recorded at the limit rather than at their true level. The DGP here censors at zero, so about a fifth of the sample piles up there.

data_sar = simulate_sar_tobit(
    n=18, rho=0.4, beta=np.array([1.0, 1.5]), sigma=0.8, censoring=0.0, seed=1
)
y, X, W = data_sar["y"], data_sar["X"], data_sar["W_graph"]

print(f"observations   : {len(y)}   (18 x 18 queen grid)")
print(f"censored at 0  : {(y <= 0).mean():.1%}")
print("true parameters: rho=0.4, beta=[1.0, 1.5], sigma=0.8")

observations   : 324   (18 x 18 queen grid)
censored at 0  : 21.0%
true parameters: rho=0.4, beta=[1.0, 1.5], sigma=0.8

model_sar = SARTobit(y=y, X=X, W=W, censoring=0.0)
idata_sar, secs_sar = fit_timed(model_sar, **FIT)

az.summary(idata_sar, var_names=["rho", "beta", "sigma"]).round(3)
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
rho 0.404 0.044 0.321 0.486 0.001 0.001 2342.0 2731.0 1.0
beta[x0] 1.040 0.080 0.898 1.197 0.002 0.001 2303.0 2329.0 1.0
beta[x1] 1.537 0.057 1.436 1.652 0.001 0.001 2595.0 2614.0 1.0
sigma 0.792 0.037 0.725 0.860 0.001 0.001 2954.0 2966.0 1.0

Important

Get censoring right or the fit is meaningless The threshold is not inferred. Pass censoring=0.0 for data censored at zero, and the actual limit otherwise. Fitting a censored outcome with the wrong threshold — or with an ordinary Gaussian model — biases every coefficient toward zero, because the pile-up at the limit is read as genuine variation in the outcome.

Read marginal effects, not coefficients

This is the part that most often goes wrong. In a linear SAR model \(\beta\) is a marginal effect; under censoring it is not, because a unit change in \(x\) moves the latent outcome and only part of that reaches the observed one — the part that does depends on how close each observation sits to the threshold.

spatial_effects() does the averaging and splits the result into direct, indirect and total components.

model_sar.spatial_effects().round(4)
direct direct_ci_lower direct_ci_upper direct_pvalue indirect indirect_ci_lower indirect_ci_upper indirect_pvalue total total_ci_lower total_ci_upper total_pvalue
variable
x1 1.6126 1.496 1.735 0.0 0.9779 0.6694 1.3501 0.0 2.5906 2.2155 3.0315 0.0

Report those. The gap between direct and the raw coefficient is the combined effect of the censoring and of the spatial multiplier feeding back onto the unit itself; indirect is what reaches other units.

Put the spatial structure in the errors instead

Use SEMTobit when the spatial pattern is a nuisance — clustered unobservables rather than a channel through which one unit’s covariates reach another. The spatial parameter is named lam, and because there are no covariate-mediated spillovers, the effects decomposition collapses to the coefficient itself.

data_sem = simulate_sem_tobit(
    n=18, lam=0.4, beta=np.array([1.0, 1.5]), sigma=0.8, censoring=0.0, seed=1
)
model_sem = SEMTobit(
    y=data_sem["y"], X=data_sem["X"], W=data_sem["W_graph"], censoring=0.0
)
idata_sem, secs_sem = fit_timed(model_sem, **FIT)

az.summary(idata_sem, var_names=["lam", "beta", "sigma"]).round(3)
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
lam 0.439 0.071 0.307 0.571 0.001 0.001 3489.0 3178.0 1.0
beta[x0] 1.100 0.079 0.953 1.250 0.001 0.001 4434.0 2802.0 1.0
beta[x1] 1.505 0.058 1.402 1.617 0.001 0.001 2449.0 2800.0 1.0
sigma 0.757 0.037 0.691 0.825 0.001 0.001 4311.0 3180.0 1.0

Add spatially lagged covariates

SDMTobit includes \(WX\) alongside \(X\), so a neighbour’s covariates enter your outcome directly rather than only through the lag. It is the right choice when you think the characteristics of neighbours matter, not just their outcomes — and it is the specification the LM tests point to when both the lag and the error channel register.

data_sdm = simulate_sdm_tobit(
    n=18, rho=0.4, beta1=np.array([1.0, 1.5]), sigma=0.8, censoring=0.0, seed=1
)
model_sdm = SDMTobit(
    y=data_sdm["y"], X=data_sdm["X"], W=data_sdm["W_graph"], censoring=0.0
)
idata_sdm, secs_sdm = fit_timed(model_sdm, **FIT)

print(f"coefficients: {idata_sdm.posterior['beta'].sizes['coefficient']} (X and WX)")
az.summary(idata_sdm, var_names=["rho", "sigma"]).round(3)

coefficients: 3 (X and WX)
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
rho 0.385 0.062 0.273 0.507 0.001 0.001 1747.0 2351.0 1.0
sigma 0.778 0.036 0.715 0.846 0.001 0.001 3882.0 2955.0 1.0 
pd.DataFrame(
    [
        pd.Series(
            {
                "spatial parameter": p,
                "posterior mean": float(i.posterior[p].mean()),
                "posterior sd": float(i.posterior[p].std()),
                "seconds": s,
                "ess": float(az.ess(i, var_names=[p])[p]),
                "rhat": float(az.rhat(i, var_names=[p])[p]),
            },
            name=n,
        )
        for n, p, i, s in (
            ("SARTobit", "rho", idata_sar, secs_sar),
            ("SEMTobit", "lam", idata_sem, secs_sem),
            ("SDMTobit", "rho", idata_sdm, secs_sdm),
        )
    ]
).round(3)
spatial parameter posterior mean posterior sd seconds ess rhat
SARTobit rho 0.404 0.044 27.408 2342.069 1.003
SEMTobit lam 0.439 0.071 20.066 3488.875 1.001
SDMTobit rho 0.385 0.062 21.722 1747.178 1.000

What to check before trusting the output

ess_bulk and r_hat on the spatial parameter first — it is the slowest-mixing quantity in all four models. Then check that the censored share in your data matches what you told the model, and look at spatial_effects() rather than the coefficient table.

Because these are NUTS models, sample_stats also carries divergence counts; any divergences at all are worth investigating before reading the posterior.

pd.DataFrame(
    [
        pd.Series(
            {
                "divergences": int(i.sample_stats["diverging"].sum()),
                "min ess_bulk": float(az.summary(i)["ess_bulk"].min()),
                "max rhat": float(az.summary(i)["r_hat"].max()),
            },
            name=n,
        )
        for n, i in (
            ("SARTobit", idata_sar),
            ("SEMTobit", idata_sem),
            ("SDMTobit", idata_sdm),
        )
    ]
)
divergences min ess_bulk max rhat
SARTobit 0.0 1366.0 1.00
SEMTobit 0.0 1831.0 1.01
SDMTobit 0.0 1624.0 1.00

See also