R cross-check — mixture-of-normals error in the PSID earnings regression

bayesm::rnmixGibbs (Bayesian) + flexmix (frequentist BIC)

Companion to mixreg_python.ipynb. We confirm, with off-the-shelf R, that the log-earnings regression error is a mixture of normals, not a single normal.

Approach — isolate the error. Off-the-shelf "mixture regression" packages fit a mixture of regressions (each component its own slopes = a coefficient/latent-class mixture). Our model (Geweke–Keane style) keeps the slope β common and mixes only the error. To match it, we first take OLS residuals (β is stable across components, as the Python fit showed) and then fit a mixture of normals to those residuals — exactly the error distribution.

• flexmix — EM; BIC across k = the model-selection cross-check for our WAIC.
• bayesm::rnmixGibbs — Bayesian mixture of normals (your existing toolchain); the error-shape cross-check.

.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(bayesm); library(flexmix)})

d  <- read.csv("psid_earnings.csv")          # 2,733 positive earners (exported from psid.mat)
le <- log(d$earn)
az <- scale(d$age); ez <- scale(d$educ); a2 <- scale(az^2)
fit <- lm(le ~ az + a2 + ez)                 # common-beta mean function
res <- as.numeric(residuals(fit))
cat(sprintf("n=%d   OLS coefs: const=%.2f age=%.2f age2=%.2f educ=%.2f   resid sd=%.2f\n",
            length(res), coef(fit)[1], coef(fit)[2], coef(fit)[3], coef(fit)[4], sd(res)))
cat(sprintf("residual skew=%.2f  kurtosis=%.2f  (0,0 = normal)\n",
            mean(scale(res)^3), mean(scale(res)^4) - 3))

n=2733   OLS coefs: const=10.40 age=0.09 age2=-0.19 educ=0.30   resid sd=0.97

residual skew=-3.28  kurtosis=22.82  (0,0 = normal)

---

flexmix — how many normal components does BIC want?

stepFlexmix(res ~ 1, k = 1:4) fits a 1..4-component normal mixture of the residuals and compares by BIC. (k>2 collapse back to 2 — empty components are dropped — so BIC is flat past k=2.)

set.seed(1)
sm <- stepFlexmix(res ~ 1, k = 1:4, nrep = 3, verbose = FALSE)
print(round(BIC(sm), 1))
best <- getModel(sm, "BIC")
cat(sprintf("\nBIC selects k = %d components (Delta BIC vs k=1 = %.0f)\n",
            best@k, BIC(sm)[1] - min(BIC(sm))))
pp <- prior(best); pr <- parameters(best)
for (j in 1:best@k)
  cat(sprintf("  comp %d: weight=%.2f  mean=%+.2f  sd=%.2f\n", j, pp[j], pr[1, j], pr[2, j]))

     1      2      3      4 
7577.9 6043.1 6043.1 6043.1 

BIC selects k = 2 components (Delta BIC vs k=1 = 1535)

  comp 1: weight=0.87  mean=+0.17  sd=0.50
  comp 2: weight=0.13  mean=-1.10  sd=1.98

flexmix result matches Python almost exactly. BIC drops ≈1,535 from k=1 to k=2 (cf. the Python WAIC drop ≈1,580) and selects k=2: a tight majority component (≈87%, mean ≈ +0.17, sd ≈ 0.50) plus a wide low-earnings component (≈13%, mean ≈ −1.10, sd ≈ 1.98) — the same fat left tail the Python sampler found (86%/14%, sd 0.49/1.94).

---

bayesm::rnmixGibbs — Bayesian mixture of normals

The Bayesian workhorse (Rossi–Allenby–McCulloch). We fit a 3-component mixture to the residuals and form the posterior-mean error density (averaged over draws, which sidesteps label switching), then compare it to the single normal.

Y <- matrix(res, ncol = 1)
set.seed(1)
out <- rnmixGibbs(Data = list(y = Y), Prior = list(ncomp = 3),
                  Mcmc = list(R = 4000, keep = 2, nprint = 0))
pd <- out$nmix$probdraw; cd <- out$nmix$compdraw
ki <- (nrow(pd) %/% 2 + 1):nrow(pd)                       # drop burn-in half
cat("avg component weights:", round(colMeans(pd[ki, , drop = FALSE]), 2), "\n")

g <- seq(min(res), max(res), length = 300); dens <- numeric(length(g))
for (r in ki) { pr <- pd[r, ]; comp <- cd[[r]]
  for (k in seq_along(pr)) dens <- dens + pr[k] * dnorm(g, comp[[k]]$mu, 1 / abs(comp[[k]]$rooti[1, 1])) }
dens <- dens / length(ki)

hist(res, breaks = 80, freq = FALSE, col = "grey85", border = "grey60",
     main = "Earnings residual: mixture (bayesm) vs single normal",
     xlab = "log-earnings residual")
curve(dnorm(x, 0, sd(res)), add = TRUE, col = "blue", lwd = 2, lty = 2)
lines(g, dens, col = "red", lwd = 2.4)
legend("topleft", c("rnmixGibbs mixture", "single normal"),
       col = c("red", "blue"), lwd = 2, lty = c(1, 2), bty = "n")
cat(sprintf("mixture density peak=%.2f vs single-normal peak=%.2f -> much more peaked + fat-tailed\n",
            max(dens), dnorm(0, 0, sd(res))))

 Starting Gibbs Sampler for Mixture of Normals
  2733  observations on  1  dimensional data
     using  3  mixture components
 
 Prior Parms: 
  mu_j ~ N(mubar,Sigma (x) A^-1)
  mubar = 
     [,1]
[1,]    0
  precision parm for prior variance of mu vectors (A)=  0.01
  Sigma_j ~ IW(nu,V) nu=  4
  V =
     [,1]
[1,]    4
  Dirichlet parameters 
[1] 5 5 5
 
 Mcmc Parms: R=  4000  keep=  2  nprint=  0  LogLike=  FALSE
 
starting value for z
z
  1   2   3 
911 911 911 
 

avg component weights: 0.03 0.25 0.72 

mixture density peak=0.74 vs single-normal peak=0.41 -> much more peaked + fat-tailed

The OLS-residual histogram (grey) with two fits overlaid: a single normal (blue dashed) and the bayesm::rnmixGibbs posterior-mean mixture density (red). The mixture is sharply peaked with a fat left tail, where the single normal is far too flat at the centre and too thin below — the same two-component error (an ≈87% tight component plus an ≈13% wide low-earnings one) that the from-scratch Gibbs and flexmix both found.

---

Conclusion — R confirms the Python result

Evidence	Verdict
flexmix BIC	k=2 beats k=1 by ≈1,535 → reject the single normal
flexmix components	87% (sd 0.50) + 13% low (mean −1.1, sd 2.0) ≈ Python's 86%/14%
rnmixGibbs density	sharply peaked, fat left tail — far from a single normal
Python (mixreg_gibbs)	WAIC drop ≈1,580; same two-component error

Three independent routes — from-scratch Gibbs (Python), EM/BIC (flexmix), and a Bayesian mixture sampler (bayesm) — agree: the earnings regression error is a mixture of normals with a heavy low-earnings tail, exactly the Geweke–Keane "flexible disturbance" story, here on the earnings level.

(No off-the-shelf R package implements the Geweke–Keane mixture-error probit itself — that remains the from-scratch mixprobit_gibbs.py.)