Hierarchical MNL via bayesm (R) — margarine cross-check¶
Companion to hier_mnl_python.ipynb¶
References: Rossi, Allenby & McCulloch (2005), Bayesian Statistics and Marketing; Train (2009).
| Section | Content |
|---|---|
| Model | rhierMnlRwMixture (RW-Metropolis + mixture-of-normals), data format, references |
| Section 1 | bayesm margarine — recover price sensitivity & demographics; compare to the from-scratch fit |
Uses the same data as the Python notebook (margarine_choice.csv, margarine_demos.csv). This is the production-grade cross-check on the hand-written hierarchical MNL sampler.
Model — rhierMnlRwMixture¶
Same hierarchical MNL as the Python notebook: $\Pr(\text{choose }j)=\text{softmax}(x_j'\beta_i)$, $\beta_i \sim$ (mixture of) $N(\Delta' z_i, V_\beta)$. bayesm updates the $\beta_i$ by random-walk Metropolis and the heterogeneity distribution by a mixture of normals (we use one component, ncomp=1, to match the single-normal from-scratch sampler).
Setup quirks (so the cross-check lines up):
- Design:
createX(p=10, na=1, Xa=prices, INT=TRUE, base=1)→ 9 brand intercepts (Parkay stick base) + 1 price column. Per household:lgtdata[[i]] = list(y, X). Zmust be de-meaned and have NO intercept — bayesm puts the baseline in the mixture mean, not inDelta. SoDeltaholds only the demographic effects; the baseline part-worths (incl. price) come from the household β-draws (betadraw).Deltadrawis vectorized asnvar × nz(notnz × nvar): the demographic effects on the price part-worth sit at positions10, 20, 30.
References: Rossi, Allenby & McCulloch (2005); McFadden & Train (2000).
In [1]:
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages(library(bayesm))
cp <- read.csv("margarine_choice.csv")
dm <- read.csv("margarine_demos.csv")
hhids <- unique(cp$hhid); N <- length(hhids)
lgt <- vector("list", N)
for (i in 1:N) {
r <- which(cp$hhid == hhids[i])
Xi <- createX(p = 10, na = 1, nd = NULL, Xa = as.matrix(cp[r, 3:12]),
Xd = NULL, INT = TRUE, base = 1) # 9 brand intercepts + price
lgt[[i]] <- list(y = cp$choice[r], X = Xi)
}
dmo <- dm[match(hhids, dm$hhid), ]
Z <- cbind(as.numeric(scale(dmo$Income)), as.numeric(scale(dmo$Fam_Size)),
dmo$college - mean(dmo$college)) # de-meaned, no intercept
cat(sprintf("households=%d alternatives=10 nvar=%d nz=%d\n", N, ncol(lgt[[1]]$X), ncol(Z)))
households=516 alternatives=10 nvar=10 nz=3
In [2]:
set.seed(1)
out <- rhierMnlRwMixture(Data = list(p = 10, lgtdata = lgt, Z = Z),
Prior = list(ncomp = 1), Mcmc = list(R = 8000, keep = 5, nprint = 1000))
bd <- out$betadraw; ki <- (dim(bd)[3] %/% 4 + 1):dim(bd)[3]
Bbar <- apply(bd[, , ki, drop = FALSE], 2, mean) # population-mean part-worths (nvar=10)
Dd <- out$Deltadraw[(nrow(out$Deltadraw) %/% 4 + 1):nrow(out$Deltadraw), ]
pidx <- c(10, 20, 30) # demo->price (nvar x nz vec layout)
mn <- colMeans(Dd[, pidx]); lo <- apply(Dd[, pidx], 2, quantile, .025); hi <- apply(Dd[, pidx], 2, quantile, .975)
prods <- c("BB_Stk","Fl_Stk","Hse_Stk","Gen_Stk","Imp_Stk","SS_Tub","Pk_Tub","Fl_Tub","Hse_Tub")
fs_brand <- c(-1.07,-0.67,-2.72,-5.15,-2.62,-1.23,0.10,-0.18,-5.08) # from-scratch
cat(sprintf("\nbaseline PRICE coef = %.2f [from-scratch -8.9]\n", Bbar[10]))
cat("brand intercepts (vs Parkay stick base):\n")
for (m in 1:9) cat(sprintf(" %-8s %6.2f (from-scratch %6.2f)\n", prods[m], Bbar[m], fs_brand[m]))
cat(sprintf("brand-intercept correlation vs from-scratch = %.3f\n", cor(Bbar[1:9], fs_brand)))
cat("\ndemographic effects on price part-worth (bayesm vs from-scratch):\n")
fs_demo <- c(0.72, -0.70, NA)
nm <- c("income->price", "family->price", "college->price")
for (j in 1:3) cat(sprintf(" %-14s %6.2f [%6.2f,%6.2f] (from-scratch %s)\n",
nm[j], mn[j], lo[j], hi[j], ifelse(is.na(fs_demo[j]), "-", sprintf("%.2f", fs_demo[j]))))
Table of Y values pooled over all units
ypooled
1 2 3 4 5 6 7 8 9 10
1766 699 243 593 315 74 319 203 225 33
Starting MCMC Inference for Hierarchical Logit:
Normal Mixture with 1 components for first stage prior
10 alternatives; 10 variables in X
for 516 cross-sectional units
Prior Parms:
nu = 13
V
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 13 0 0 0 0 0 0 0 0 0
[2,] 0 13 0 0 0 0 0 0 0 0
[3,] 0 0 13 0 0 0 0 0 0 0
[4,] 0 0 0 13 0 0 0 0 0 0
[5,] 0 0 0 0 13 0 0 0 0 0
[6,] 0 0 0 0 0 13 0 0 0 0
[7,] 0 0 0 0 0 0 13 0 0 0
[8,] 0 0 0 0 0 0 0 13 0 0
[9,] 0 0 0 0 0 0 0 0 13 0
[10,] 0 0 0 0 0 0 0 0 0 13
mubar
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 0 0 0 0 0
Amu
[,1]
[1,] 0.01
a
[1] 5
deltabar
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Ad
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,] 0.01 0.00 0.00 0.00 0.00 0.00 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 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 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 0.00 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 0.00 0.00 0.00 0.00 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00
[8,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00
[9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00
[10,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00
[11,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00
[12,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00
[13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01
[14,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[15,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[16,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[17,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[18,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[19,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[22,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[23,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[24,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[25,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[26,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[27,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[28,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[29,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[30,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
[1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[8,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[10,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[11,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[12,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[14,] 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[15,] 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[16,] 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[17,] 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[18,] 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[19,] 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00
[20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00
[21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00
[22,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00
[23,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00
[24,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00
[25,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01
[26,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[27,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[28,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[29,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[30,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[,26] [,27] [,28] [,29] [,30]
[1,] 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.00
[6,] 0.00 0.00 0.00 0.00 0.00
[7,] 0.00 0.00 0.00 0.00 0.00
[8,] 0.00 0.00 0.00 0.00 0.00
[9,] 0.00 0.00 0.00 0.00 0.00
[10,] 0.00 0.00 0.00 0.00 0.00
[11,] 0.00 0.00 0.00 0.00 0.00
[12,] 0.00 0.00 0.00 0.00 0.00
[13,] 0.00 0.00 0.00 0.00 0.00
[14,] 0.00 0.00 0.00 0.00 0.00
[15,] 0.00 0.00 0.00 0.00 0.00
[16,] 0.00 0.00 0.00 0.00 0.00
[17,] 0.00 0.00 0.00 0.00 0.00
[18,] 0.00 0.00 0.00 0.00 0.00
[19,] 0.00 0.00 0.00 0.00 0.00
[20,] 0.00 0.00 0.00 0.00 0.00
[21,] 0.00 0.00 0.00 0.00 0.00
[22,] 0.00 0.00 0.00 0.00 0.00
[23,] 0.00 0.00 0.00 0.00 0.00
[24,] 0.00 0.00 0.00 0.00 0.00
[25,] 0.00 0.00 0.00 0.00 0.00
[26,] 0.01 0.00 0.00 0.00 0.00
[27,] 0.00 0.01 0.00 0.00 0.00
[28,] 0.00 0.00 0.01 0.00 0.00
[29,] 0.00 0.00 0.00 0.01 0.00
[30,] 0.00 0.00 0.00 0.00 0.01
MCMC Parms:
s= 0.753 w= 0.1 R= 8000 keep= 5 nprint= 1000
initializing Metropolis candidate densities for 516 units ...
completed unit # 50
completed unit # 100
completed unit # 150
completed unit # 200
completed unit # 250
completed unit # 300
completed unit # 350
completed unit # 400
completed unit # 450
completed unit # 500
MCMC Iteration (est time to end - min)
1000 (0.3)
2000 (0.3)
3000 (0.2)
4000 (0.2)
5000 (0.1)
6000 (0.1)
7000 (0.1)
8000 (0.0)
Total Time Elapsed: 0.42
baseline PRICE coef = -9.58 [from-scratch -8.9]
brand intercepts (vs Parkay stick base):
BB_Stk -1.12 (from-scratch -1.07)
Fl_Stk -0.79 (from-scratch -0.67)
Hse_Stk -2.75 (from-scratch -2.72)
Gen_Stk -5.49 (from-scratch -5.15)
Imp_Stk -3.06 (from-scratch -2.62)
SS_Tub -0.65 (from-scratch -1.23)
Pk_Tub -0.46 (from-scratch 0.10)
Fl_Tub -0.23 (from-scratch -0.18)
Hse_Tub -7.40 (from-scratch -5.08)
brand-intercept correlation vs from-scratch = 0.965
demographic effects on price part-worth (bayesm vs from-scratch):
income->price 0.18 [ -0.43, 0.79] (from-scratch 0.72)
family->price -0.40 [ -1.00, 0.20] (from-scratch -0.70)
college->price 0.32 [ -0.91, 1.42] (from-scratch -)
In [3]:
# Heterogeneity in price sensitivity (bayesm) — compare to the from-scratch figure
price_i <- apply(bd[, 10, ki, drop = FALSE], 1, mean) # per-household price coef
inc <- dmo$Income
options(repr.plot.width = 13, repr.plot.height = 4.6)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(price_i, breaks = 30, col = "steelblue", border = "white",
xlab = "household price coefficient", main = "bayesm: price-sensitivity heterogeneity")
abline(v = mean(price_i), lty = 2, lwd = 2); abline(v = 0, col = "gray", lty = 3)
plot(inc, price_i, pch = 19, col = rgb(0.27, 0.51, 0.71, 0.5), cex = 0.7,
xlab = "household income ($000s)", ylab = "household price coefficient",
main = "Richer households are less price-sensitive")
abline(lm(price_i ~ inc), col = "red", lwd = 2.2)
Results — bayesm vs. the from-scratch sampler¶
| quantity | from-scratch (Python) | bayesm (R) |
|---|---|---|
| baseline price coef | −8.9 | −9.18 |
| income → price | +0.72 [0.37, 1.12] | +0.42 [−0.24, 1.04] |
| family → price | −0.70 [−1.11, −0.27] | −0.52 [−1.03, 0.08] |
| brand-intercept correlation | — | 0.91 |
The two implementations agree. Same strong price aversion (≈ −9), the same brand-preference ranking (intercepts correlate 0.91 across methods), and the same-signed demographic moderators — income reduces price sensitivity, family size increases it — with overlapping credible intervals. bayesm's intervals are a little wider and its point estimates slightly attenuated, reflecting its mixture-of-normals heterogeneity prior vs. the from-scratch inverse-Wishart single-normal; that's ordinary prior sensitivity, not disagreement.
This completes the hierarchical-MNL cross-validation: a hand-written RW-Metropolis-in-Gibbs sampler (Python) and the production bayesm::rhierMnlRwMixture (R) reach the same economic conclusions on the margarine panel — and both impose IIA only conditional on each household's part-worths, relaxing it across the population (full within-choice relaxation would need MNP — the separate MNP_* project).