R cross-check — bayesm::rhierMnlRwMixture on margarine¶
Is the population of tastes non-normal? (Rossi, Allenby & McCulloch, BSM Ch. 5)¶
Companion to hmnl_mixture_python.ipynb (from-scratch). Same model — random-coefficients MNL with a mixture-of-normals heterogeneity distribution — fitted with the textbook routine, comparing ncomp = 1, 2, 3 by log marginal density and inspecting the fitted taste distribution.
In [1]:
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages(library(bayesm))
cp <- read.csv("margarine_choice.csv")
hh <- unique(cp$hhid); N <- length(hh); pricecols <- 3:12
lgt <- vector("list", N)
for (i in 1:N) {
r <- which(cp$hhid == hh[i])
lgt[[i]] <- list(y = cp$choice[r],
X = createX(p=10, na=1, Xa=as.matrix(cp[r, pricecols]),
nd=NULL, Xd=NULL, INT=TRUE, base=1)) # 9 brand intercepts + price
}
cat(sprintf("%d households, 10 products, design k=10 (9 intercepts[base=Parkay stick] + price)\n", N))
516 households, 10 products, design k=10 (9 intercepts[base=Parkay stick] + price)
Model comparison — ncomp = 1, 2, 3 by log marginal density¶
Newton–Raftery harmonic-mean estimate (logMargDenNR — bayesm's own metric; higher = better). Single normal = ncomp=1.
In [2]:
fits <- list()
cat(sprintf("%-8s%14s%10s %s\n", "ncomp", "logMargDenNR", "price", "weights"))
for (K in 1:3) {
set.seed(1)
out <- rhierMnlRwMixture(Data = list(p=10, lgtdata=lgt),
Prior = list(ncomp=K), Mcmc = list(R=8000, keep=2, nprint=0))
fits[[K]] <- out
ll <- out$loglike; ll <- ll[(length(ll)%/%2+1):length(ll)]
bd <- out$betadraw; nd <- dim(bd)[3]; ki <- (nd%/%2+1):nd
pw <- colMeans(out$nmix$probdraw[(nrow(out$nmix$probdraw)%/%2+1):nrow(out$nmix$probdraw), , drop=FALSE])
cat(sprintf("%-8d%14.0f%10.2f %s\n", K, logMargDenNR(ll), mean(bd[,10,ki]),
paste(round(pw,2), collapse=",")))
}
ncomp logMargDenNR price weights
Z not specified
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
MCMC Parms:
s= 0.753 w= 0.1 R= 8000 keep= 2 nprint= 0
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
1 -3857 -9.38 1 Z not specified
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 2 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 5
MCMC Parms:
s= 0.753 w= 0.1 R= 8000 keep= 2 nprint= 0
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
2 -3838 -9.43 0.46,0.54 Z not specified
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 3 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 5 5
MCMC Parms:
s= 0.753 w= 0.1 R= 8000 keep= 2 nprint= 0
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
3 -3851 -9.09 0.21,0.65,0.14
The fitted heterogeneity density¶
The label-invariant population density of a coefficient from the mixture draws,
$\;f(b)=\frac1S\sum_s\sum_k \pi_k^{(s)} N\!\big(b;\ \mu_k^{(s)},\ \Sigma_k^{(s)}\big)$.
Compare the single normal (ncomp=1) with the mixture (ncomp=2) for price and a brand part-worth (Parkay tub).
In [3]:
fitdens <- function(out, coord, g) {
pd <- out$nmix$probdraw; cd <- out$nmix$compdraw; ki <- (nrow(pd)%/%2+1):nrow(pd)
d <- numeric(length(g))
for (r in ki) { pr <- pd[r,]; cp <- cd[[r]]
for (k in seq_along(pr)) {
S <- solve(cp[[k]]$rooti %*% t(cp[[k]]$rooti))
d <- d + pr[k] * dnorm(g, cp[[k]]$mu[coord], sqrt(S[coord, coord])) } }
d / length(ki)
}
mom <- function(out, coord) { # posterior-mean skew & excess kurtosis of the mixture
pd <- out$nmix$probdraw; cd <- out$nmix$compdraw; ki <- (nrow(pd)%/%2+1):nrow(pd); sk <- ek <- numeric(0)
for (r in ki) { pr <- pd[r,]; mu <- sapply(seq_along(pr), function(k) cd[[r]][[k]]$mu[coord])
s2 <- sapply(seq_along(pr), function(k){S<-solve(cd[[r]][[k]]$rooti%*%t(cd[[r]][[k]]$rooti)); S[coord,coord]})
m <- sum(pr*mu); v <- sum(pr*(s2+mu^2))-m^2
m3 <- sum(pr*(mu^3+3*mu*s2))-3*m*v-m^3; m4 <- sum(pr*(mu^4+6*mu^2*s2+3*s2^2))-4*m*m3-6*m^2*v-m^4
sk <- c(sk, m3/v^1.5); ek <- c(ek, m4/v^2-3) }
c(skew=mean(sk), exkurt=mean(ek))
}
par(mfrow=c(1,2))
for (cd in list(c(10,"price sensitivity",-16,-2), c(8,"Parkay-tub part-worth",-7,8))) {
coord <- as.integer(cd[1]); nm <- cd[2]; g <- seq(as.numeric(cd[3]), as.numeric(cd[4]), length=300)
plot(g, fitdens(fits[[1]], coord, g), type="l", col="blue", lwd=2, lty=2,
xlab=nm, ylab="density", main=paste("Heterogeneity:", nm))
lines(g, fitdens(fits[[2]], coord, g), col="red", lwd=2.4)
legend("topright", c("single normal (ncomp=1)","mixture (ncomp=2)"),
col=c("blue","red"), lwd=2, lty=c(2,1), bty="n")
}
par(mfrow=c(1,1))
cat(sprintf("ncomp=2 fitted price: skew=%.2f exkurt=%.2f\n", mom(fits[[2]],10)[1], mom(fits[[2]],10)[2]))
cat(sprintf("ncomp=2 fitted Pk-tub: skew=%.2f exkurt=%.2f\n", mom(fits[[2]],8)[1], mom(fits[[2]],8)[2]))
ncomp=2 fitted price: skew=0.11 exkurt=0.55
ncomp=2 fitted Pk-tub: skew=0.01 exkurt=1.55
Conclusion — bayesm agrees with the from-scratch result¶
| ncomp | logMargDenNR | price | weights |
|---|---|---|---|
| 1 | −3857 | −9.38 | 1 |
| 2 | −3838 | −9.43 | 0.46, 0.54 |
| 3 | −3851 | −9.09 | 0.21, 0.65, 0.14 |
ncomp=2 wins (log marginal density up ≈ +19 over the single normal; ncomp=3 falls back) — exactly the from-scratch Python verdict (+44 on its scale; both favor 2, both reject 3). The price coefficient (≈ −9.4) matches, and the fitted-density shape confirms a non-normal, leptokurtic taste distribution (ncomp=2 fitted excess kurtosis ≈ 0.55 for price, 1.55 for the Parkay-tub part-worth; mild skew) — the mixture density is more peaked / fatter-tailed than the single normal, Rossi et al.'s point.
As in Python, this is the mixture earning its keep as a **flexible heterogeneity *distribution***, not as discrete segments (the two components overlap heavily in price). Two engines, one conclusion.
(Caveat: logMargDenNR is the Newton–Raftery harmonic-mean estimator — known to be noisy and complexity-biased — so the fitted-density shape is the more robust evidence.)