How to fit spatial tobit panel models¶

Censored outcomes observed for the same units over several periods. SARPanelTobit puts the spatial lag on the latent outcome; SEMPanelTobit puts the spatial structure in the disturbances.

Both are NUTS-only, and both need N, T and censoring stated explicitly — the panel shape and the censoring threshold are not inferred. Equations and constructor arguments are in Supported Models.

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

from neighbayes import dgp
from neighbayes.models import SARPanelTobit, SEMPanelTobit

N, T, SIDE = 100, 5, 10  # 10x10 rook grid, 5 periods
RHO, LAM, SIGMA = 0.35, 0.35, 0.8
BETA = np.array([1.0, 1.4])
FIT = dict(draws=1000, tune=1000, chains=4, random_seed=42, progressbar=False)
rng = np.random.default_rng(42)

Fit a censored panel with a spatial lag¶

y and X are stacked long — \(N \cdot T\) rows, unit-major — and you tell the model the shape with N and T.

sar_data = dgp.simulate_panel_sar_tobit_fe(
    N=N,
    T=T,
    n=SIDE,
    contiguity="rook",
    rho=RHO,
    beta=BETA,
    sigma=SIGMA,
    censoring=0.0,
    rng=rng,
)
y_sar, X_sar, W = sar_data["y"], sar_data["X"], sar_data["W_graph"]

print(f"rows          : {len(y_sar)}   ({N} units x {T} periods)")
print(f"censored at 0 : {(y_sar == 0).mean():.1%}")

sar_model = SARPanelTobit(y=y_sar, X=X_sar, W=W, N=N, T=T, censoring=0.0)
sar_idata = sar_model.fit(**FIT)

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

rows          : 500   (100 units x 5 periods)
censored at 0 : 20.6%

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 2 jobs)
NUTS: [rho, beta, sigma, y_cens_gap]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 17 seconds.

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
rho	0.307	0.039	0.234	0.380	0.001	0.001	3155.0	2972.0	1.0
beta[x0]	0.983	0.072	0.853	1.114	0.001	0.001	3061.0	3014.0	1.0
beta[x1]	1.435	0.047	1.349	1.521	0.001	0.001	3295.0	2566.0	1.0
sigma	0.889	0.032	0.826	0.946	0.000	0.000	4685.0	3245.0	1.0

Important

censoring must match your data The threshold is not inferred. Pass the value below which observations are recorded at the limit rather than at their true level — 0.0 here. Get it wrong, or fit an uncensored model instead, and every coefficient is biased toward zero because the pile-up at the limit reads as genuine variation.

Interpreting the coefficients¶

Warning

The effects decomposition is not available for panel tobit spatial_effects() raises NotImplementedError on SARPanelTobit and SEMPanelTobit. That matters for interpretation, because \(\beta\) is not a marginal effect once the outcome is censored — a unit change in \(x\) moves the latent outcome, and only the uncensored part of that reaches the observed one.

Read the coefficients as latent-scale parameters until the decomposition lands. Where the direct/indirect split is central to your argument, the cross-sectional SARTobit does provide it.

try:
    sar_model.spatial_effects()
except NotImplementedError as err:
    print(f"NotImplementedError: {err}")

NotImplementedError: Spatial effects not yet implemented for panel Tobit models.

Put the spatial structure in the errors instead¶

SEMPanelTobit treats spatial correlation as a nuisance in the disturbances rather than a channel between units. The spatial parameter is lam, and because there are no covariate-mediated spillovers the effects decomposition collapses to the coefficient.

sem_data = dgp.simulate_panel_sem_tobit_fe(
    N=N,
    T=T,
    n=SIDE,
    contiguity="rook",
    lam=LAM,
    beta=BETA,
    sigma=SIGMA,
    censoring=0.0,
    rng=rng,
)
sem_model = SEMPanelTobit(
    y=sem_data["y"],
    X=sem_data["X"],
    W=sem_data["W_graph"],
    N=N,
    T=T,
    censoring=0.0,
)
sem_idata = sem_model.fit(**FIT)

print(f"censored at 0 : {(sem_data['y'] == 0).mean():.1%}")
az.summary(sem_idata, var_names=["lam", "beta", "sigma"]).round(3)

censored at 0 : 31.8%

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 2 jobs)
NUTS: [lam, beta, sigma, y_cens_gap]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 18 seconds.

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
lam	0.324	0.065	0.195	0.440	0.001	0.001	4443.0	3143.0	1.0
beta[x0]	0.852	0.069	0.720	0.977	0.001	0.001	6055.0	2920.0	1.0
beta[x1]	1.351	0.055	1.253	1.460	0.001	0.001	3233.0	3119.0	1.0
sigma	0.952	0.038	0.883	1.024	0.001	0.001	5269.0	3444.0	1.0

What to check before trusting the output¶

The spatial parameter mixes slowest, and these are NUTS models, so divergences matter. Any at all are worth chasing down before you read the posterior.

pd.DataFrame(
    [
        pd.Series(
            {
                "spatial parameter": p,
                "posterior mean": float(i.posterior[p].mean()),
                "truth": t,
                "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, p, t, i in (
            ("SARPanelTobit", "rho", RHO, sar_idata),
            ("SEMPanelTobit", "lam", LAM, sem_idata),
        )
    ]
).round(3)

	spatial parameter	posterior mean	truth	divergences	min ess_bulk	max rhat
SARPanelTobit	rho	0.307	0.35	0	1739.0	1.0
SEMPanelTobit	lam	0.324	0.35	0	2188.0	1.0

Short chains are used here so the docs build quickly. If ess_bulk on the spatial parameter is in the low hundreds on your own data, raise draws and tune before concluding anything from the interval.