Mixture-of-normals regression — flexible-error modeling of earnings¶
The continuous sibling of Geweke & Keane (1997)¶
References: Geweke & Keane (1997), Mixture of Normals Probit Models (FRB Minneapolis SR 237); Albert & Chib (1993); Geweke (2005), Contemporary Bayesian Econometrics and Statistics (source of the PSID extract).
| Section | Content |
|---|---|
| Model | Regression with a mixture-of-normals error; relation to Geweke–Keane & to coefficient-mixtures |
| Section 1 | Synthetic validation — recover β and the error shape; WAIC favors the mixture |
| Section 2 | PSID earnings — log-earnings on age + education with a flexible mixture error |
From-scratch sampler: mixreg_gibbs.py. (The binary Geweke–Keane probit itself is in mixprobit_gibbs.py.)
Model — regression with a mixture-of-normals error¶
$$y_t = x_t'\beta + \varepsilon_t, \qquad \varepsilon_t \sim \sum_{j=1}^m p_j\, N(\mu_j,\ \sigma_j^2)$$
This is the continuous sibling of Geweke & Keane (1997): GK put a mixture of normals on the disturbance of a binary probit (a latent-utility model); here we observe $y$ directly, so the same flexible error applies to an ordinary regression. A mixture of normals can approximate essentially any density, so this removes the normality assumption on the error — capturing skew and fat tails that a single normal cannot (earnings are a textbook case).
Two different "mixtures of normals" — don't confuse them:
- here / Geweke–Keane: mixture on the error $\varepsilon_t$ → flexible distribution / link (one equation, no hierarchy).
hlm_mixture/rhierMnlRwMixture: mixture on the coefficients $\beta_i$ across units → heterogeneity / segments (needs panel data).
Identification: the error mixture is recentred to mean 0 each sweep (its location is absorbed by the intercept), so $\beta$ is the usual mean function and the mixture describes the disturbance shape.
Gibbs blocks: (1) component indicator $s_t$ — categorical; (2) $\beta$ — weighted Gaussian (GLS, weights $1/\sigma^2_{s_t}$); (3) $(\mu_j,\sigma_j^2)$ — Normal-Inverse-Gamma; (4) $p$ — Dirichlet. Model comparison by WAIC here; the companion mixreg_pymc.ipynb adds Bayes factors / posterior model probabilities via SMC (as GK used).
The Gibbs sampler in detail¶
The mixture likelihood $\sum_j p_j N(\cdot)$ is not conjugate on its own, but data augmentation makes it so. Introduce a latent component label $s_t\in\{1,\dots,m\}$ for each observation — which mixture component generated its disturbance. Conditional on the labels every block has a closed-form full conditional, giving a four-block sweep:
Labels $s_t\mid\beta,\{\mu_j,\sigma_j^2\},p$ — categorical. With residual $r_t=y_t-x_t'\beta$, $$P(s_t=j)\;\propto\;p_j\,N(r_t;\,\mu_j,\sigma_j^2).$$ Each residual is stochastically assigned to the component that best explains it.
Coefficients $\beta\mid s,\{\mu_j,\sigma_j^2\},y$ — Gaussian (weighted least squares). Given the labels, observation $t$ has known variance $\sigma^2_{s_t}$, so this is GLS on the de-meaned response $y_t-\mu_{s_t}$ with weights $1/\sigma^2_{s_t}$: $$\beta\sim N(\bar b,\bar B),\quad \bar B^{-1}=X'\Omega X+B_0^{-1},\ \ \bar b=\bar B\,X'\Omega\,(y-\mu_s),\ \ \Omega=\mathrm{diag}(1/\sigma^2_{s_t}),$$ drawn via its Cholesky factor.
Components $(\mu_j,\sigma_j^2)\mid s,\beta$ — Normal-Inverse-Gamma, per component. The residuals assigned to component $j$ are iid $N(\mu_j,\sigma_j^2)$; with the conjugate NIG prior ($\sigma_j^2\sim\mathrm{IG}(a_0,b_0)$, $\mu_j\mid\sigma_j^2\sim N(0,\sigma_j^2/\kappa_0)$) the posterior is again NIG, so $\sigma_j^2$ then $\mu_j\mid\sigma_j^2$ are drawn in closed form.
Weights $p\mid s$ — Dirichlet$(\alpha+n_1,\dots,\alpha+n_m)$, with $n_j$ the count assigned to component $j$.
Identification (each sweep). Two non-identified directions are fixed: (i) the mixture's overall location is confounded with the intercept, so we recentre $\sum_j p_j\mu_j=0$ and pass the shift into $\beta_0$ — then $\beta$ is the mean function and the mixture is purely the error shape; (ii) the labels are exchangeable (label switching), so we sort components by $\mu_j$.
For model comparison each kept draw also stores the closed-form per-observation log density $\log\sum_j p_j\,N(y_t;\,x_t'\beta+\mu_j,\sigma_j^2)$ (labels integrated out) — the input WAIC needs.
What is WAIC?¶
WAIC (Widely Applicable Information Criterion; Watanabe 2010) measures a model's out-of-sample predictive accuracy directly from the posterior draws — lower is better. With draws $\theta^{(1)},\dots,\theta^{(S)}$ from the posterior,
$$\mathrm{WAIC} = -2\Big(\underbrace{\textstyle\sum_{i} \log \tfrac1S\textstyle\sum_s p(y_i\mid\theta^{(s)})}_{\text{lppd — fit}} \;-\; \underbrace{\textstyle\sum_i \mathrm{Var}_s\big[\log p(y_i\mid\theta^{(s)})\big]}_{p_\mathrm{WAIC}\ \text{— complexity penalty}}\Big).$$
- lppd (log pointwise predictive density): each observation's likelihood averaged over the posterior, then logged and summed — rewards fit.
- $p_\mathrm{WAIC}$ (effective number of parameters): the posterior variance of each observation's log-likelihood, summed — penalizes models that fit only by being flexible.
Unlike AIC/DIC (which plug in a single point estimate), WAIC uses the whole posterior and asymptotically approximates leave-one-out cross-validation. It needs only the per-observation likelihoods — here the closed-form mixture density $p(y_t\mid x_t)=\sum_j p_j\,N(y_t;\,x_t'\beta+\mu_j,\sigma_j^2)$, computed by waic() in mixreg_gibbs.py. It plays the same role as Geweke–Keane's Bayes factors (comparing models on predictive grounds) but is cheaper and more stable. Differences of a few units are meaningful; the gaps we see are in the hundreds–thousands.
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
from scipy import stats
from mixreg_gibbs import mixreg_gibbs, waic, error_density
print('mixreg_gibbs loaded')
mixreg_gibbs loaded
1. Synthetic validation¶
$y = X\beta + \varepsilon$, $n=3000$, true $\beta=(1.0,\,0.5,\,-0.7)$, and a skewed two-component error: 70% $N(-0.3,0.5^2)$ + 30% $N(0.7,1.4^2)$. Expect: $\beta$ recovered by any $m$ (it is the mean function), but WAIC should strongly prefer the mixture, and the fitted error density should match the skewed truth.
rng = np.random.default_rng(0)
n = 3000
X = np.column_stack([np.ones(n), rng.normal(size=n), rng.normal(size=n)])
beta_true = np.array([1.0, 0.5, -0.7])
z = rng.random(n) < 0.7
eps = np.where(z, rng.normal(-0.3, 0.5, n), rng.normal(0.7, 1.4, n)) # skewed mixture
y = X @ beta_true + eps
fits = {m: mixreg_gibbs(y, X, m=m, R=4000, burn=1500, seed=1) for m in (1, 2, 3)}
print(f"{'m':>3}{'WAIC':>10} beta")
for m, o in fits.items():
print(f"{m:>3}{waic(o):>10.1f} {np.round(o['beta'].mean(0), 3)}")
print(f"\ntruth beta = {beta_true}")
m WAIC beta 1 8568.7 [ 1.007 0.486 -0.682] 2 7639.4 [ 1.008 0.485 -0.688]
3 7638.2 [ 1.007 0.484 -0.689] truth beta = [ 1. 0.5 -0.7]
grid = np.linspace(eps.min(), eps.max(), 300)
fig, ax = plt.subplots(figsize=(8, 4.6))
ax.hist(eps, bins=60, density=True, color='lightgray', label='true errors (sample)')
ax.plot(grid, error_density(fits[1], grid), 'b--', lw=2, label='normal (m=1)')
ax.plot(grid, error_density(fits[2], grid), 'r-', lw=2.2, label='mixture (m=2)')
ax.set_xlabel('error'); ax.set_ylabel('density')
ax.set_title('Synthetic: mixture recovers the skewed error a single normal misses')
ax.legend(); plt.tight_layout()
plt.savefig('mixreg_synth.png', dpi=120, bbox_inches='tight'); plt.show()
1. Results — synthetic¶
| m | WAIC | β |
|---|---|---|
| 1 (normal) | 8568.7 | [1.01, 0.49, −0.68] |
| 2 (mixture) | 7639.4 | [1.01, 0.49, −0.69] |
| 3 | 7638.2 | [1.01, 0.48, −0.69] |
β is recovered by every model (it is the conditional mean — a misspecified error does not bias it much here), but WAIC drops by ≈929 from m=1 to m=2 and then flattens — correctly identifying the two-component truth. The density plot shows the single normal forced into a symmetric bell while the mixture captures the skew. This validates the sampler before the earnings application.
2. PSID earnings: a flexible-error wage regression¶
psid.mat (from Geweke 2005, Contemporary Bayesian Econometrics): age, education, earnings. We model log-earnings (the 2,733 positive earners) on age, age², and education with a mixture-of-normals error. Earnings are famously non-normal (here the OLS residual has kurtosis ≈ 23), so this is exactly where the flexible error should pay off — and where Geweke–Keane's machinery shines.
M = loadmat('psid.mat')['psid']
age, educ, earn = M[:, 0], M[:, 1], M[:, 2]
pos = earn > 0
y = np.log(earn[pos]) # log-earnings, positive earners
a, e = age[pos], educ[pos]
az = (a - a.mean()) / a.std()
ez = (e - e.mean()) / e.std()
a2 = (az ** 2 - (az ** 2).mean()) / (az ** 2).std()
Xe = np.column_stack([np.ones(pos.sum()), az, a2, ez]) # [1, age, age^2, educ] (standardized)
resid = y - Xe @ np.linalg.lstsq(Xe, y, rcond=None)[0]
print(f"n={pos.sum()} positive earners mean log-earn={y.mean():.2f} (~${np.exp(y.mean()):,.0f})")
print(f"OLS residual: skew={stats.skew(resid):.2f} kurtosis={stats.kurtosis(resid):.2f} (0,0 = normal)")
n=2733 positive earners mean log-earn=10.40 (~$32,713)
OLS residual: skew=-3.28 kurtosis=22.84 (0,0 = normal)
efits = {m: mixreg_gibbs(y, Xe, m=m, R=5000, burn=2000, seed=1) for m in (1, 2, 3, 4)}
print(f"{'m':>3}{'WAIC':>10} beta=[const, age, age^2, educ]")
for m, o in efits.items():
print(f"{m:>3}{waic(o):>10.1f} {np.round(o['beta'].mean(0), 3)}")
best = min(efits, key=lambda m: waic(efits[m]))
print(f"\nbest by WAIC: m={best} (Delta WAIC vs normal = {waic(efits[1]) - waic(efits[best]):.0f})")
o2 = efits[2]
print("m=2 components (weight, mean, sd):")
for j in range(2):
print(f" {o2['p'].mean(0)[j]:.2f} {o2['mu'].mean(0)[j]:+.2f} {np.sqrt(o2['sig2'].mean(0))[j]:.2f}")
m WAIC beta=[const, age, age^2, educ] 1 7583.7 [10.396 0.086 -0.191 0.297] 2 6001.5 [10.395 0.096 -0.134 0.266]
3 5908.1 [10.395 0.099 -0.123 0.269] 4 5900.7 [10.393 0.1 -0.125 0.268]
best by WAIC: m=4 (Delta WAIC vs normal = 1683) m=2 components (weight, mean, sd): 0.14 -1.11 1.94 0.86 +0.18 0.49
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
resid = y - Xe @ efits[1]['beta'].mean(0)
grid = np.linspace(resid.min(), resid.max(), 400)
ax[0].hist(resid, bins=80, density=True, color='lightgray', label='OLS residuals')
ax[0].plot(grid, error_density(efits[1], grid), 'b--', lw=2, label='normal (m=1)')
ax[0].plot(grid, error_density(efits[3], grid), 'r-', lw=2.2, label='mixture (m=3)')
ax[0].set_xlabel('log-earnings residual'); ax[0].set_ylabel('density')
ax[0].set_title('Earnings error: fat left tail the normal misses'); ax[0].legend()
ax[1].plot(list(efits), [waic(efits[m]) for m in efits], 'o-', color='firebrick')
ax[1].set_xlabel('number of mixture components m'); ax[1].set_ylabel('WAIC (lower = better)')
ax[1].set_xticks(list(efits)); ax[1].set_title('WAIC: the mixture decisively beats the normal')
plt.tight_layout()
plt.savefig('mixreg_earnings.png', dpi=120, bbox_inches='tight'); plt.show()
2. Results — earnings¶
| m | WAIC | β = [const, age, age², educ] |
|---|---|---|
| 1 (normal) | 7583.7 | [10.40, 0.09, −0.19, 0.30] |
| 2 | 6001.5 | [10.40, 0.10, −0.13, 0.27] |
| 3 | 5908.1 | [10.40, 0.10, −0.12, 0.27] |
| 4 | 5900.7 | [10.39, 0.10, −0.13, 0.27] |
The mixture wins decisively: WAIC falls ≈1,580 from the normal to the mixture (m=1→m=2), improving to ≈1,680 by m=4. The conditional mean is stable across m (intercept 10.4 ⇒ mean earnings ≈ $33k; a positive, concave age profile; an education coefficient ≈0.27 per SD — a solid return to schooling). What changes is the error: the m=2 fit is a tight majority component (≈86%, sd ≈ 0.5) plus a wide low-earnings component (≈14%, mean ≈ −1.1, sd ≈ 1.9) — a heavy left tail of people earning far less than age/education predict. A single normal cannot represent this and is strongly rejected.
Why it matters: the normal-error regression would mis-state predictive intervals and the probability of low-earnings outcomes; the mixture gives a calibrated, flexible earnings distribution. This is the continuous counterpart of Geweke & Keane's finding that a mixture of normals dominates the single-normal assumption — here on the earnings level rather than a participation probit (where, on this extract, the 90/10 binary outcome was too uninformative to detect it).
Contrast with the coefficient-mixture projects: this mixture flexes the error distribution in one equation; hlm_mixture / hierarchical-MNL mixtures flex coefficient heterogeneity across units. Different targets, same tool.