Bayesian binary logit via bayesm (R) — independence Metropolis¶

Companion to binary_logit_python.ipynb¶

Reference: Rossi, Allenby & McCulloch (2005), Bayesian Statistics and Marketing (bayesm).

Section Content
Model Logit as a 2-alternative MNL; why bayesm uses Metropolis (not Gibbs); rmnlIndepMetrop; refs
Section 1 Synthetic data (shared binary_logit_synth.csv) — posterior vs. truth, MLE, and the Python samplers
Section 2 Real data — Mroz (1987) women’s labor-force participation

Fits the same synthetic data as the Python notebook (binary_logit_synth.csv) → a direct cross-language check. bayesm estimates a binary logit as a multinomial logit with 2 alternatives by independence Metropolis.

Model — binary logit as MNL(2), independence Metropolis¶

$$y_i \sim \text{Bernoulli}(\sigma(x_i'\beta))$$

bayesm has no binary-logit Gibbs sampler — because, unlike probit's Albert–Chib, logit has no conjugate data augmentation built in — so it fits logit as a multinomial logit with 2 alternatives via independence Metropolis (rmnlIndepMetrop):

  • alt 1 = "success" ($y=1$), alt 2 = base ($y=0$); utilities $u_1 = x'\beta$, $u_2 = 0 \Rightarrow \Pr(y=1)=\sigma(x'\beta)$.
  • createX(p=2, nd=3, Xd=[x1,x2,x3], INT=TRUE, base=2) builds the stacked $(2n)\times k$ design.
  • Independence Metropolis: proposes $\beta$ from a multivariate-$t$ centred at the MLE with the MLE covariance — so proposals already resemble the posterior and acceptance is high (~0.8), unlike the random-walk Metropolis in the Python notebook (acceptance ~0.3).

The symmetry in bayesm: probit → rbprobitGibbs (Albert–Chib Gibbs); logit → rmnlIndepMetrop (Metropolis). The Pólya–Gamma trick (Section 2 of the Python notebook) is exactly what restores Gibbs for logit.

References: Rossi, Allenby & McCulloch (2005); McCulloch & Rossi (1994).

In [5]:
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages(library(bayesm))

d  <- read.csv("C:/Users/user/project/binary_logit_synth.csv")
cat(sprintf("n=%d  mean(y)=%.3f\n", nrow(d), mean(d$y)))

Xd   <- as.matrix(d[, c("x1", "x2", "x3")])
ymnl <- ifelse(d$y == 1, 1L, 2L)            # alt 1 = success, alt 2 = base
X    <- createX(p = 2, na = NULL, nd = 3, Xa = NULL, Xd = Xd, INT = TRUE, base = 2)
cat("createX -> design dim:", dim(X), " (expect", nrow(d) * 2, "x 4)\n")

n=2000  mean(y)=0.569
createX -> design dim: 4000 4  (expect 4000 x 4)
In [6]:
set.seed(0)
out <- rmnlIndepMetrop(Data = list(p = 2, y = ymnl, X = X),
                       Mcmc = list(R = 20000, keep = 5, nprint = 0))
bd   <- out$betadraw
burn <- nrow(bd) %/% 4
bd   <- bd[(burn + 1):nrow(bd), ]

truth <- c(0.5, -1.0, 1.5, -0.75)
mle   <- c(0.445, -1.005, 1.535, -0.619)    # Python notebook, identical data
nm    <- c("const", "x1", "x2", "x3")
pm <- colMeans(bd); ps <- apply(bd, 2, sd)

cat(sprintf("\nindependence-Metropolis acceptance = %.3f\n\n", mean(out$acceptr)))
cat(sprintf("%-7s %7s %9s %12s %10s\n", "param", "truth", "MLE", "bayesm mean", "bayesm sd"))
for (j in 1:4)
  cat(sprintf("%-7s %7.2f %9.3f %12.3f %10.3f\n", nm[j], truth[j], mle[j], pm[j], ps[j]))

 table of y values
y
   1    2 
1138  862 
 
Starting Independence Metropolis Sampler for Multinomial Logit Model
   2000  obs with  2  alternatives
 
Table of y Values
y
   1    2 
1138  862 
Prior Parms: 
betabar
[1] 0 0 0 0
A
     [,1] [,2] [,3] [,4]
[1,] 0.01 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00
[3,] 0.00 0.00 0.01 0.00
[4,] 0.00 0.00 0.00 0.01
 
MCMC parms: 
R=  20000  keep=  5  nprint=  0  nu (df for st candidates) =  6
 

independence-Metropolis acceptance = 0.818

param     truth       MLE  bayesm mean  bayesm sd
const      0.50     0.445        0.444      0.058
x1        -1.00    -1.005       -1.008      0.066
x2         1.50     1.535        1.538      0.080
x3        -0.75    -0.619       -0.619      0.060
In [7]:
options(repr.plot.width = 12, repr.plot.height = 3.4)
truth <- c(0.5, -1.0, 1.5, -0.75); mle <- c(0.445, -1.005, 1.535, -0.619)
nm <- c("const", "x1", "x2", "x3")
par(mfrow = c(1, 4), mar = c(4, 3, 2, 1))
for (j in 1:4) {
  hist(bd[, j], breaks = 40, freq = FALSE, col = "lightsteelblue", border = "white",
       main = nm[j], xlab = "beta")
  abline(v = truth[j], col = "black", lty = 2, lwd = 2)
  abline(v = mle[j],   col = "firebrick", lwd = 2)
}
legend("topright", c("truth", "MLE"), col = c("black", "firebrick"),
       lty = c(2, 1), lwd = 2, bty = "n", cex = 0.8)
No description has been provided for this image

Results — bayesm logit vs. the Python samplers¶

param truth MLE bayesm (R) Python RWM Python PG bambi PyMC
const 0.50 0.445 0.444 0.446 0.449 0.444 0.444
x1 −1.00 −1.005 −1.008 −1.008 −1.016 −1.000 −1.005
x2 1.50 1.535 1.538 1.542 1.552 1.527 1.539
x3 −0.75 −0.619 −0.619 −0.623 −0.625 −0.616 −0.620

(bayesm posterior sd ≈ [0.058, 0.066, 0.080, 0.060], matching the MLE se.)

Six methods agree across two languages — the from-scratch RW-Metropolis and Pólya–Gamma Gibbs, bambi, PyMC (Python), the bayesm independence Metropolis (R), and the MLE baseline. Independence-Metropolis acceptance is ~0.82, far above the random walk's ~0.30, because the proposal is the MLE-centred multivariate-$t$ — proposals already look like the posterior.

The probit/logit symmetry, made concrete: bayesm does probit by Gibbs (rbprobitGibbs, Albert–Chib) but logit by Metropolis (rmnlIndepMetrop), precisely because logit lacks a built-in conjugate augmentation. The Pólya–Gamma sampler in Section 2 of the Python notebook is the missing piece that restores Gibbs for logit.

This completes the binary-logit cross-validation. Next: choose a real-data application.

2. Real data: Mroz (1987) women's labor-force participation¶

The same Mroz LFP model as binary_logit_python.ipynb Section 4 (mroz_lfp.csv, $n=753$):

$$\text{inlf} \sim \text{nwifeinc} + \text{educ} + \text{exper} + \text{exper}^2 + \text{age} + \text{kidslt6} + \text{kidsge6}$$

fit by bayesm's independence Metropolis — a cross-language check against the Python samplers and Wooldridge's published estimates.

Variables (Mroz 1987; married white women aged 30–60, $n=753$; shared mroz_lfp.csv):

Variable Definition Units mean range exp. sign
inlf (outcome) =1 if in the labor force in 1975 0/1 0.568 {0,1} —
nwifeinc Non-wife income = (family income − wife's earnings) $1000s 20.13 −0.03…96 − (income effect)
educ Years of schooling years 12.29 5…17 +
exper Actual labor-market experience years 10.63 0…45 +
expersq exper² years² 178.0 0…2025 − (concave)
age Age years 42.54 30…60 −
kidslt6 Number of children < 6 count 0.238 0…3 − (childcare)
kidsge6 Number of children 6–18 count 1.353 0…8 ≈ 0

In the MNL(2) coding, alt 1 = "in labor force" with utility $x'\beta$ and alt 2 (out) as base, so $\Pr(\text{inlf})=\Lambda(x'\beta)$ and the coefficients are read exactly as a binary logit.

In [8]:
dm <- read.csv("C:/Users/user/project/mroz_lfp.csv")
cols <- c("nwifeinc", "educ", "exper", "expersq", "age", "kidslt6", "kidsge6")
Xdm  <- as.matrix(dm[, cols])
ymnl2 <- ifelse(dm$inlf == 1, 1L, 2L)
Xm <- createX(p = 2, na = NULL, nd = length(cols), Xa = NULL, Xd = Xdm, INT = TRUE, base = 2)

set.seed(0)
outm <- rmnlIndepMetrop(Data = list(p = 2, y = ymnl2, X = Xm),
                        Mcmc = list(R = 20000, keep = 5, nprint = 0))
bdm <- outm$betadraw; bdm <- bdm[(nrow(bdm) %/% 4 + 1):nrow(bdm), ]
nm  <- c("const", cols)
wool <- c(0.425, -0.021, 0.221, 0.206, -0.003, -0.088, -1.443, 0.060)
pm <- colMeans(bdm); ps <- apply(bdm, 2, sd)

cat(sprintf("independence-Metropolis acceptance = %.3f\n\n", mean(outm$acceptr)))
cat(sprintf("%-9s %11s %12s %10s\n", "param", "Wooldridge", "bayesm mean", "bayesm sd"))
for (j in 1:8) cat(sprintf("%-9s %11.3f %12.3f %10.3f\n", nm[j], wool[j], pm[j], ps[j]))

 table of y values
y
  1   2 
428 325 
 
Starting Independence Metropolis Sampler for Multinomial Logit Model
   753  obs with  2  alternatives
 
Table of y Values
y
  1   2 
428 325 
Prior Parms: 
betabar
[1] 0 0 0 0 0 0 0 0
A
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00
[7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00
[8,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01
 
MCMC parms: 
R=  20000  keep=  5  nprint=  0  nu (df for st candidates) =  6
 
independence-Metropolis acceptance = 0.492

param      Wooldridge  bayesm mean  bayesm sd
const           0.425        0.437      0.855
nwifeinc       -0.021       -0.022      0.008
educ            0.221        0.224      0.043
exper           0.206        0.208      0.033
expersq        -0.003       -0.003      0.001
age            -0.088       -0.089      0.015
kidslt6        -1.443       -1.467      0.204
kidsge6         0.060        0.062      0.075

2. Results — Mroz LFP (R bayesm vs. Python & Wooldridge)¶

param Wooldridge bayesm (R) Python PG odds ratio
const 0.425 0.437 0.441 —
nwifeinc −0.021 −0.022 −0.022 0.98 / $1k
educ 0.221 0.224 0.225 1.25 / yr
exper 0.206 0.208 0.208 1.23
expersq −0.003 −0.003 −0.003 concave
age −0.088 −0.089 −0.090 0.91 / yr
kidslt6 −1.443 −1.467 −1.467 0.23 / child < 6
kidsge6 0.060 0.062 0.059 1.06 (n.s.)

bayesm's independence Metropolis (acceptance ≈ 0.49) reproduces Wooldridge's logit estimates and the Python from-scratch / bambi / PyMC results — the same economic story (a young child is the dominant drag on participation; education and experience raise it). This closes the binary-logit project: two languages, six samplers, synthetic + the canonical Mroz dataset, all in agreement.

In [9]:
# --- Analysis of findings: odds ratios, predicted probabilities, average partial effects ---
B  <- bdm
lo <- apply(B, 2, quantile, 0.025); hi <- apply(B, 2, quantile, 0.975); mn <- colMeans(B)
cat(sprintf("%-9s %8s  %-20s %4s %8s\n", "coef", "mean", "95% CI", "sig", "OR"))
for (j in 1:8)
  cat(sprintf("%-9s %8.3f  [%7.3f, %7.3f] %4s %8.3f\n", nm[j], mn[j], lo[j], hi[j],
              ifelse(lo[j] > 0 | hi[j] < 0, "yes", "no"), exp(mn[j])))

Xf   <- cbind(1, Xdm); xbar <- colMeans(Xf)
Pavg <- function(x) mean(plogis(as.numeric(B %*% x)))
prof <- function(k, v) { x <- xbar; x[which(nm == k)] <- v; x }
cat(sprintf("\nP(inlf) at sample means = %.3f\n", Pavg(xbar)))
cat(sprintf("  educ 12 -> 16:  %.3f -> %.3f\n", Pavg(prof("educ", 12)), Pavg(prof("educ", 16))))
cat(sprintf("  kids<6 0 -> 1:  %.3f -> %.3f\n", Pavg(prof("kidslt6", 0)), Pavg(prof("kidslt6", 1))))

base <- plogis(Xf %*% t(B))
ame  <- function(k, d) { X2 <- Xf; X2[, which(nm == k)] <- X2[, which(nm == k)] + d
  e <- colMeans(plogis(X2 %*% t(B)) - base); c(mean(e), quantile(e, 0.025), quantile(e, 0.975)) }
ak <- ame("kidslt6", 1); ae <- ame("educ", 1)
cat(sprintf("  AME +1 child<6 = %+.3f [%+.3f, %+.3f]\n", ak[1], ak[2], ak[3]))
cat(sprintf("  AME +1 yr educ = %+.3f [%+.3f, %+.3f]\n", ae[1], ae[2], ae[3]))
cat(sprintf("\nin-sample accuracy (0.5 cutoff) = %.3f  (base rate %.3f)\n",
            mean((plogis(Xf %*% mn) > 0.5) == dm$inlf), max(mean(dm$inlf), 1 - mean(dm$inlf))))

coef          mean  95% CI                sig       OR
const        0.437  [ -1.212,   2.106]   no    1.548
nwifeinc    -0.022  [ -0.038,  -0.005]  yes    0.979
educ         0.224  [  0.142,   0.309]  yes    1.251
exper        0.208  [  0.144,   0.271]  yes    1.232
expersq     -0.003  [ -0.005,  -0.001]  yes    0.997
age         -0.089  [ -0.118,  -0.062]  yes    0.915
kidslt6     -1.467  [ -1.876,  -1.082]  yes    0.231
kidsge6      0.062  [ -0.082,   0.208]   no    1.063

P(inlf) at sample means = 0.583
  educ 12 -> 16:  0.567 -> 0.761
  kids<6 0 -> 1:  0.664 -> 0.315
  AME +1 child<6 = -0.258 [-0.314, -0.199]
  AME +1 yr educ = +0.039 [+0.025, +0.052]

in-sample accuracy (0.5 cutoff) = 0.736  (base rate 0.568)

2. Analysis of findings — Mroz LFP (bayesm posterior)¶

Computed from the bayesm independence-Metropolis draws — and matching the Python Section 4 results (same data, same model).

  • Significance: six of seven predictors are credibly nonzero (nwifeinc, educ, exper, expersq, age, kidslt6); only kidsge6 is not (95% CI ≈ [−0.08, 0.21]) — school-age children don't measurably affect participation once the rest is controlled.
  • Odds ratios: kidslt6 0.23 (each child under 6 multiplies the odds of working by ~0.23 — the dominant effect), educ 1.25/yr, exper 1.23 (concave via expersq), age 0.91/yr, nwifeinc 0.98/$1k.
  • Predicted probabilities (others at means): ≈ 0.58 at the average profile; HS→college 0.57 → 0.76; first child under 6 0.66 → 0.32.
  • Average partial effects: +1 child < 6 ≈ −0.26 probability, +1 year educ ≈ +0.039, each with a credible interval.
  • Fit: ~73.6% in-sample accuracy vs. a 56.8% base rate.

Cross-language agreement: these match the Python from-scratch / bambi / PyMC findings to ~0.01 — the bayesm posterior, the two from-scratch Python samplers, and the production tools all tell the identical economic story. The Bayesian framing delivers every quantity (odds ratios, predicted probabilities, partial effects) with a credible interval directly from the posterior draws.