Censored Gaussian regression — the Tobit model

From-scratch data-augmentation Gibbs; the PyMC "Censored Data" example

Module: tobit_gibbs.py. When an outcome is only observed within a known region — a sensor that saturates, a value below a detection limit, or hours worked that pile up at 0 — ordinary regression is biased. The censored-Gaussian (Tobit) model fixes it. This is the PyMC "Censored Data" example's domain and the classic econometrics Tobit model (Tobin 1958). The new ingredient over the survival projects is left and interval censoring (not just right).

Censoring types and the likelihood

A latent $y^*_i\sim N(x_i'\beta,\sigma^2)$ is observed only within a known region:

type	what we know	likelihood contribution
exact	$y^*=y$	density $\phi((y-x'\beta)/\sigma)/\sigma$
left	$y^*\le c$ (e.g. Tobit at 0)	$\Phi((c-x'\beta)/\sigma)$
right	$y^*\ge c$	$1-\Phi((c-x'\beta)/\sigma)$
interval	$a<y^*\le b$	$\Phi((b-x'\beta)/\sigma)-\Phi((a-x'\beta)/\sigma)$

(The data-augmentation Gibbs algorithm is detailed in the next section.)

The data-augmentation Gibbs sampler in detail

The censored likelihood is awkward because it mixes two kinds of terms — a Gaussian density for the exactly-observed points and a normal-CDF tail probability for each censored point — so it is not conjugate and cannot be sampled in one clean step. The trick (Tanner–Wong / Chib 1992) is to treat the unobserved latent values as missing data: if we knew the true $y^_i$ behind every censored observation, the completed dataset would be an ordinary Gaussian linear regression* — fully conjugate. So we alternate imputing the latents and updating the parameters.

Let $C$ index the censored observations (region $(\text{lo}_i,\text{hi}_i]$: $(-\infty,c]$ for left, $[c,\infty)$ for right, $(a,b]$ for interval), and keep $y^*_i=y_i$ for exact points.

Block 1 — impute the latents $\;y^*_i\mid\beta,\sigma^2\ \sim\ N(x_i'\beta,\sigma^2)$ truncated to $(\text{lo}_i,\text{hi}_i]$ for each $i\in C$. Each draw is the model's guess of where the censored value actually fell, made consistent with both the regression and the censoring constraint. Sampled by the inverse-CDF method: $$u\sim\mathrm{Unif}\big(\Phi(\alpha_i),\Phi(\beta_i)\big),\quad y^*_i=x_i'\beta+\sigma\,\Phi^{-1}(u),\qquad \alpha_i=\tfrac{\text{lo}_i-x_i'\beta}{\sigma},\ \beta_i=\tfrac{\text{hi}_i-x_i'\beta}{\sigma}.$$

Block 2 — coefficients $\;\beta\mid y^*,\sigma^2$. With the latents filled in this is Bayesian linear regression; under the conjugate prior $\beta\sim N(0,B_0)$ the posterior is Gaussian, $$\beta\sim N(\bar m,\bar V),\qquad \bar V=\big(X'X/\sigma^2+B_0^{-1}\big)^{-1},\quad \bar m=\bar V\,X'y^*/\sigma^2,$$ drawn via its Cholesky factor.

Block 3 — variance $\;\sigma^2\mid y^*,\beta$. Conjugate Inverse-Gamma on the completed residuals, $$\sigma^2\sim\mathrm{IG}\!\Big(a_0+\tfrac n2,\ \ b_0+\tfrac12\textstyle\sum_i (y^*_i-x_i'\beta)^2\Big).$$

Why it works. Cycling these three conditionals is a valid Gibbs sweep on the joint posterior $p(y^*_C,\beta,\sigma^2\mid y_{\text{obs}})$; dropping the imputed $y^*_C$ leaves draws from the target $p(\beta,\sigma^2\mid y_{\text{obs}})$. The intractable censored likelihood has been replaced by a loop of standard, closed-form draws — the same data-augmentation idea as the ZIP latent indicator, except here the latent is the unobserved value itself. Averaging the imputed $y^*_i$ across sweeps also gives the model's estimate of where each censored point "really" sat (the latent-variable figure below).

Mixing caveat. The imputed latents and $(\beta,\sigma^2)$ are correlated, so the autocorrelation of the chain grows with the censored fraction — for Mroz (≈43% censored at 0) this is non-trivial. A marginalised sampler that integrates the latents out (PyMC's pm.Censored + NUTS) avoids it; we contrast the two below.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from tobit_gibbs import simulate_tobit, tobit_gibbs, tobit_mle, summary

# ---- Section 1: synthetic recovery, and why naive OLS is biased ----
y, X, lo, hi, exact, ystar = simulate_tobit(n=3000, beta=(1.0, 2.0, -1.5), sigma=2.0, L=0.0, seed=1)
print('synthetic: %d obs, %.0f%% LEFT-censored at 0' % (len(y), 100*(~exact).mean()))
bm, sm = tobit_mle(y, X, lo, hi, exact)
o = tobit_gibbs(y, X, lo, hi, exact, seed=1)
b_ols = np.linalg.lstsq(X, y, rcond=None)[0]
print('  truth        beta = [1.0, 2.0, -1.5]  sigma = 2.0')
print('  Tobit MLE    beta = %s  sigma %.2f' % (np.round(bm,3), sm))
print('  Tobit Bayes  beta = %s  sigma %.2f' % (np.round(o['beta'].mean(0),3), o['sigma'].mean()))
print('  naive OLS    beta = %s   <-- BIASED toward 0 (slopes attenuated)' % np.round(b_ols,3))

# the SAME sampler handles INTERVAL censoring: observe only which unit-width bin y* falls in
binw = 1.0; a_iv = np.floor(ystar/binw)*binw; b_iv = a_iv + binw
oi = tobit_gibbs(ystar*0.0, X, a_iv, b_iv, np.zeros(len(y), bool), seed=1)
print('\n  interval-censored (only the bin known): Tobit Bayes beta = %s  sigma %.2f  (truth recovered)'
      % (np.round(oi['beta'].mean(0),3), oi['sigma'].mean()))

synthetic: 3000 obs, 37% LEFT-censored at 0

  truth        beta = [1.0, 2.0, -1.5]  sigma = 2.0
  Tobit MLE    beta = [ 0.99   2.016 -1.485]  sigma 1.98
  Tobit Bayes  beta = [ 0.988  2.017 -1.486]  sigma 1.98
  naive OLS    beta = [ 1.825  1.253 -0.914]   <-- BIASED toward 0 (slopes attenuated)

  interval-censored (only the bin known): Tobit Bayes beta = [ 0.999  1.98  -1.496]  sigma 2.00  (truth recovered)

Section 2 — Mroz: married women's hours worked (the classic Tobit)

753 married women (Mroz 1987); the outcome is annual hours worked, which is left-censored at 0 — 325 of 753 (43%) did not work, so their hours pile up at zero while their desired hours are unobserved. Regressing hours on covariates by OLS is biased toward zero; Tobit models the latent desired hours.

field	meaning
hours	annual hours worked (0–4950); 0 for 43% (left-censored)
nwifeinc	non-wife household income ($000)
educ	years of schooling
exper, expersq	labor-market experience and its square
age	age in years
kidslt6, kidsge6	number of children <6 and 6–18

d = pd.read_csv('mroz_hours.csv')
Xm = np.column_stack([np.ones(len(d)), d.nwifeinc, d.educ, d.exper, d.expersq, d.age, d.kidslt6, d.kidsge6])
ym = d.hours.to_numpy(float); ex = ym > 0
loM = np.where(ex, ym, -np.inf); hiM = np.where(ex, ym, 0.0)
om = tobit_gibbs(ym, Xm, loM, hiM, ex, R=15000, burn=4000, seed=1)
b_ols_m = np.linalg.lstsq(Xm, ym, rcond=None)[0]
nm = ['intercept','nwifeinc','educ','exper','expersq','age','kidslt6','kidsge6']
print('%-10s %12s %12s' % ('', 'Tobit (Bayes)', 'naive OLS'))
for j, r in enumerate(summary(om['beta'], nm)):
    print('%-10s %9.2f (%.0f) %11.2f' % (r[0], r[1], r[2], b_ols_m[j]))
print('sigma (Tobit) = %.0f' % om['sigma'].mean())
print('\n(AER tobit target: educ 80.6, exper 131.6, age -54.4, kidslt6 -894, sigma 1122)')

           Tobit (Bayes)    naive OLS
intercept     953.33 (454)     1330.48
nwifeinc       -8.77 (5)       -3.45
educ           81.55 (22)       28.76
exper         132.90 (17)       65.67
expersq        -1.89 (1)       -0.70
age           -54.81 (8)      -30.51
kidslt6      -899.21 (115)     -442.09
kidsge6       -15.22 (39)      -32.78
sigma (Tobit) = 1137

(AER tobit target: educ 80.6, exper 131.6, age -54.4, kidslt6 -894, sigma 1122)

fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
# (left) synthetic: censoring pile-up + truth/Tobit/OLS regression lines (vs x1, others at 0)
x1 = X[:,1]; ax[0].scatter(x1[exact], y[exact], s=8, color='steelblue', alpha=.4, label='observed')
ax[0].scatter(x1[~exact], y[~exact], s=10, color='firebrick', alpha=.5, label='censored (at 0)')
gx = np.linspace(x1.min(), x1.max(), 50)
ax[0].plot(gx, 1.0 + 2.0*gx, 'k-', lw=1.5, label='truth')
ax[0].plot(gx, o['beta'].mean(0)[0] + o['beta'].mean(0)[1]*gx, '--', color='seagreen', lw=1.5, label='Tobit')
ax[0].plot(gx, b_ols[0] + b_ols[1]*gx, ':', color='darkorange', lw=1.8, label='naive OLS')
ax[0].set_xlabel('$x_1$'); ax[0].set_ylabel('y'); ax[0].set_title('Censoring biases OLS; Tobit recovers truth'); ax[0].legend(fontsize=7)
# (right) Mroz: OLS coefficients vs Tobit coefficients (OLS attenuated toward 0)
sl = slice(1, None); tb = om['beta'].mean(0)[sl]; ol = b_ols_m[sl]
ax[1].scatter(tb, ol, s=45, color='purple', zorder=3, edgecolor='k', lw=.3)
lim = [min(tb.min(),ol.min())*1.1, max(tb.max(),ol.max())*1.1]; ax[1].plot(lim, lim, 'k:', lw=1, label='OLS = Tobit')
for j, lab in enumerate(nm[1:]): ax[1].annotate(lab, (tb[j], ol[j]), fontsize=7, xytext=(3,3), textcoords='offset points')
ax[1].axhline(0,color='grey',lw=.6); ax[1].axvline(0,color='grey',lw=.6)
ax[1].set_xlabel('Tobit coefficient'); ax[1].set_ylabel('naive OLS coefficient')
ax[1].set_title('Mroz: OLS attenuated toward 0'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('tobit.png', dpi=120, bbox_inches='tight'); plt.show()

Showing the latent variable

The Tobit's engine is a latent $y^*\sim N(X\beta,\sigma)$ that we see only when positive; the Gibbs sampler imputes the censored $y^*$, so we can look at it directly:

• left (Mroz) — observed hours pile up at 0, but the fitted latent Gaussian extends well below 0: that hidden left tail is what censoring collapses onto the boundary.
• middle (Mroz) — the model's reconstructed "desired hours" for the 325 women at 0 (all negative) — the quantity naive OLS wrongly treats as a genuine 0.
• right (synthetic, truth known) — the fitted latent density matches the true latent across the whole range, including the censored region below 0 — the Tobit recovers the unseen tail from the observed part alone.

# --- showing the LATENT variable y* (what the augmentation reconstructs) ---
from scipy.stats import norm
fig, ax = plt.subplots(1, 3, figsize=(15, 4.4))
# (A) Mroz: observed hours (spike at 0) vs the fitted latent Gaussian (which extends below 0)
bmh = om['beta'].mean(0); smh = om['sigma'].mean(); grid = np.linspace(-3000, 5500, 400)
ax[0].hist(ym, bins=40, density=True, color='lightgray', edgecolor='k', lw=.3, label='observed hours')
ax[0].plot(grid, np.mean([norm.pdf(grid, mu, smh) for mu in Xm @ bmh], axis=0), 'firebrick', lw=2, label='fitted latent $N(X\\beta,\\sigma)$')
ax[0].axvline(0, color='k', ls=':', lw=1); ax[0].set_xlabel('hours'); ax[0].set_yticks([])
ax[0].set_title('Observed (spike at 0) vs latent desired hours'); ax[0].legend(fontsize=7)
# (B) Mroz: imputed latent "desired hours" for the 325 censored women (all <= 0)
ax[1].hist(om['ystar'][~ex], bins=30, color='steelblue', edgecolor='k', lw=.3)
ax[1].axvline(0, color='firebrick', ls='--', lw=1, label='censoring point (0)')
ax[1].set_xlabel('imputed latent "desired hours"'); ax[1].set_ylabel('count')
ax[1].set_title('Reconstructed latent for the\n325 censored women'); ax[1].legend(fontsize=7)
# (C) synthetic (truth known): the fitted latent recovers the FULL distribution, incl. the censored region
bms = o['beta'].mean(0); sms = o['sigma'].mean(); gs = np.linspace(ystar.min()-1, ystar.max()+1, 300)
ax[2].hist(ystar, bins=40, density=True, color='lightgray', edgecolor='k', lw=.3, label='true latent $y^*$')
ax[2].plot(gs, np.mean([norm.pdf(gs, mu, sms) for mu in X @ bms], axis=0), 'seagreen', lw=2, label='fitted latent')
ax[2].axvline(0, color='k', ls=':', lw=1, label='censoring point')
ax[2].set_xlabel('$y^*$'); ax[2].set_yticks([]); ax[2].set_title('Synthetic: latent recovered\nincl. censored region (<0)'); ax[2].legend(fontsize=7)
plt.tight_layout(); plt.savefig('tobit_latent.png', dpi=120, bbox_inches='tight'); plt.show()

Results & interpretation

• The Tobit recovers the truth; OLS does not. On the synthetic data the slopes are recovered (2.0, −1.5) while naive OLS on the censored outcome attenuates them toward 0 (≈1.25, −0.91) — the classic censoring bias. The same tobit_gibbs recovers the truth under interval censoring too (only the bin known).
• Mroz: the Tobit estimates match the AER::tobit benchmark (educ ≈ +81 hrs/yr, experience strongly positive but concave, a child under 6 ≈ −900 hours, age negative). The right panel shows OLS coefficients lying inside the diagonal — every significant effect shrunk toward zero by ignoring the 43% censoring. (The lone exception is kidsge6: small and statistically insignificant in both fits, its OLS value is if anything a touch larger — the point sitting just outside the diagonal — a reminder that censoring attenuation is a strong tendency for the real effects, not a guarantee for every noisy coefficient.)
• σ ≈ 1130 is the latent-hours residual SD.

Takeaways

• Censoring is a general idea, not a survival one. Here it's a Gaussian outcome (Tobit / detection limits); the likelihood is just density-for-exact, CDF-mass-for-censored — covering left, right, and interval in one framework.
• Augmentation makes it trivial. Imputing the censored values from truncated normals turns the awkward censored likelihood into an ordinary conjugate regression — the same data-augmentation trick used for the ZIP indicator.
• Cross-checks: PyMC (pm.Censored) and R (AER::tobit / survreg) follow.