Bayesian hierarchical multinomial logit — from scratch¶
Random-coefficients MNL; the rich-state idea applied to choice (bayesm margarine)¶
References: McFadden (1974); McFadden & Train (2000); Rossi, Allenby & McCulloch (2005), Bayesian Statistics and Marketing; Train (2009), Discrete Choice Methods with Simulation.
| Section | Content |
|---|---|
| Model | Random-coefficients MNL, the IIA question (vs. MNP), RW-Metropolis-in-Gibbs sampler, references |
| Section 1 | Synthetic validation — recover known heterogeneity (Δ, V_β, per-respondent β_i) |
| Section 2 | bayesm margarine panel (516 households, 10 products) — price-sensitivity heterogeneity |
| R companion | hier_mnl_bayesm.ipynb — bayesm::rhierMnlRwMixture cross-check on the same data |
From-scratch sampler: hier_mnl_gibbs.py — per-respondent RW-Metropolis (the binary-logit move) inside the conjugate Gibbs hierarchy (the hier-logit Δ / V_β blocks).
Model — hierarchical (random-coefficients) multinomial logit¶
Respondent $i$ chooses among $J$ alternatives in each of $T_i$ tasks; alternative $j$ has attribute vector $x_{j}$:
$$\Pr(\text{choose } j) = \frac{\exp(x_j'\beta_i)}{\sum_{l=1}^J \exp(x_l'\beta_i)}, \qquad \beta_i \sim N(\Delta' z_i,\ V_\beta)$$
with $z_i$ the respondent-level covariates (e.g. $[1,\text{demographics}]$), $\Delta$ the group-level regression, $V_\beta$ the heterogeneity covariance. The $\beta_i$ are individual part-worths / preferences; $V_\beta$ captures how much they vary across people.
Does IIA still hold?¶
Standard MNL has IIA (Independence of Irrelevant Alternatives): $\Pr(j)/\Pr(k)$ doesn't depend on the other alternatives — the red-bus/blue-bus problem — a consequence of i.i.d. Gumbel errors.
In the hierarchical / random-coefficients model:
- Conditional on $\beta_i$: each respondent is still a plain logit, so IIA holds at the individual level.
- Marginally (integrating over the population distribution of $\beta_i$): IIA is broken — heterogeneity induces realistic, flexible substitution. McFadden & Train (2000): mixed logit can approximate any random-utility model.
- Within a single respondent facing genuinely correlated alternatives, conditional IIA remains a limitation. Fully relaxing that needs correlated errors → multinomial probit.
So the two mechanisms are complementary: hierarchical MNL relaxes IIA across people (coefficient heterogeneity); MNP relaxes it within a choice (error covariance $\Sigma$). The MNP side is built out in the separate MNP project (MNP_gibbs.ipynb, MNP_pymc.ipynb, MNP_bayesm.ipynb).
Sampler — RW-Metropolis within Gibbs¶
MNL has no clean Pólya–Gamma/Gibbs step (unlike binary logit), so $\beta_i$ is updated by Metropolis:
| Block | Full conditional / update |
|---|---|
| $\beta_i$ | random-walk Metropolis, proposal scaled by each respondent's MNL Hessian |
| $\Delta$ | Matrix-normal (group-level regression of $\beta_i$ on $z_i$) |
| $V_\beta$ | Inverse-Wishart |
This is the structure of bayesm::rhierMnlRwMixture (we use a single-normal heterogeneity prior; bayesm uses a mixture of normals).
References¶
- McFadden (1974) — conditional logit / IIA
- McFadden & Train (2000) — mixed logit approximates any RUM
- Rossi, Allenby & McCulloch (2005) —
rhierMnlRwMixture - Train (2009) — Discrete Choice Methods with Simulation
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from hier_mnl_gibbs import simulate_hier_mnl, hier_mnl_gibbs, summary
print('hier_mnl_gibbs loaded')
hier_mnl_gibbs loaded
1. Synthetic validation¶
300 respondents × 20 choice tasks, $J=3$ alternatives, $k=3$ attributes, one respondent-level covariate $u$ (so $z_i=[1,u_i]$). Truth:
- group-level regression $\Delta$ (rows $[1,u]$, cols = 3 attributes): baseline part-worths $[1.0,-1.0,0.5]$; demographic effects $[0.4,\,0.0,-0.3]$
- heterogeneity covariance $V_\beta = 0.5\,I_3$
Expect: recover $\Delta$, $V_\beta$, and the per-respondent $\beta_i$ (which are only weakly identified from 20 tasks each, so they shrink toward the group mean — recovery improves with more tasks).
sim = simulate_hier_mnl(N=300, T=20, J=3, k=3, seed=0)
data, Z = sim['data'], sim['Z']
print(f"respondents={len(data)} tasks each={data[0][0].shape[0]} "
f"alternatives={data[0][0].shape[1]} attributes={data[0][0].shape[2]}")
print("true Delta (rows [intercept, u]; cols = 3 attributes):")
print(np.round(sim['Delta_true'], 2))
print("true V_beta diagonal:", np.round(np.diag(sim['Vb_true']), 2))
respondents=300 tasks each=20 alternatives=3 attributes=3 true Delta (rows [intercept, u]; cols = 3 attributes): [[ 1. -1. 0.5] [ 0.4 0. -0.3]] true V_beta diagonal: [0.5 0.5 0.5]
out = hier_mnl_gibbs(data, Z, R=3000, burn=1000, seed=1)
s = summary(out)
print(f"RW-Metropolis acceptance = {out['accept']:.3f}\n")
lab = [('baseline a1 (D00)', 0, 0, 1.0), ('baseline a2 (D01)', 0, 1, -1.0), ('baseline a3 (D02)', 0, 2, 0.5),
('demo->a1 (D10)', 1, 0, 0.4), ('demo->a2 (D11)', 1, 1, 0.0), ('demo->a3 (D12)', 1, 2, -0.3)]
print(f"{'parameter':>20}{'truth':>7}{'post mean':>11}{'95% CI':>20}")
print('-' * 58)
for nm, a, b, tr in lab:
d = out['Delta'][:, a, b]
print(f"{nm:>20}{tr:7.2f}{d.mean():11.3f} [{np.percentile(d,2.5):6.3f},{np.percentile(d,97.5):6.3f}]")
print(f"\nV_beta diagonal: truth [0.5, 0.5, 0.5] est {np.round(np.diag(s['Vb']), 3)}")
Bt, Bh = sim['B_true'], s['B']
print("per-respondent beta_i recovery corr:",
[round(np.corrcoef(Bh[:, j], Bt[:, j])[0, 1], 3) for j in range(3)])
RW-Metropolis acceptance = 0.379
parameter truth post mean 95% CI
----------------------------------------------------------
baseline a1 (D00) 1.00 0.885 [ 0.792, 0.976]
baseline a2 (D01) -1.00 -0.998 [-1.100,-0.903]
baseline a3 (D02) 0.50 0.528 [ 0.435, 0.622]
demo->a1 (D10) 0.40 0.431 [ 0.343, 0.521]
demo->a2 (D11) 0.00 0.100 [-0.002, 0.195]
demo->a3 (D12) -0.30 -0.320 [-0.408,-0.226]
V_beta diagonal: truth [0.5, 0.5, 0.5] est [0.422 0.456 0.464]
per-respondent beta_i recovery corr: [np.float64(0.88), np.float64(0.801), np.float64(0.862)]
# Per-respondent coefficient recovery: posterior-mean beta_i vs. truth
Bt, Bh = sim['B_true'], s['B']
fig, ax = plt.subplots(1, 3, figsize=(13, 4.2))
for j in range(3):
ax[j].scatter(Bt[:, j], Bh[:, j], s=14, alpha=0.55, color='steelblue', edgecolor='none')
lim = [min(Bt[:, j].min(), Bh[:, j].min()), max(Bt[:, j].max(), Bh[:, j].max())]
ax[j].plot(lim, lim, 'k--', lw=1)
ax[j].set_title(f'attribute {j+1} (corr = {np.corrcoef(Bt[:,j], Bh[:,j])[0,1]:.2f})')
ax[j].set_xlabel('true $\\beta_i$')
ax[0].set_ylabel('posterior mean $\\beta_i$')
fig.suptitle('Hierarchical MNL — per-respondent part-worth recovery (shrunk toward the group mean)',
fontsize=12)
plt.tight_layout()
plt.savefig(r'C:\Users\user\project\hier_mnl_synth.png', dpi=120, bbox_inches='tight')
plt.show()
1. Results — synthetic recovery¶
| parameter | truth | post. mean | 95% CI |
|---|---|---|---|
| baseline a1 (Δ₀₀) | 1.00 | 0.885 | [0.792, 0.976] |
| baseline a2 (Δ₀₁) | −1.00 | −0.998 | [−1.100, −0.903] |
| baseline a3 (Δ₀₂) | 0.50 | 0.528 | [0.435, 0.622] |
| demo→a1 (Δ₁₀) | 0.40 | 0.431 | [0.343, 0.521] |
| demo→a2 (Δ₁₁) | 0.00 | 0.100 | [−0.002, 0.195] |
| demo→a3 (Δ₁₂) | −0.30 | −0.320 | [−0.408, −0.226] |
$V_\beta$ diagonal: truth [0.5, 0.5, 0.5] → est [0.42, 0.46, 0.46]. Per-respondent recovery corr: [0.88, 0.80, 0.86]. RW-Metropolis acceptance ≈ 0.38.
The sampler recovers the heterogeneity structure. Both the baseline part-worths and the demographic effects land at (or near) their true values with credible intervals that cover them; the heterogeneity covariance $V_\beta$ is slightly underestimated and the per-respondent $\beta_i$ are pulled toward the group mean — the usual shrinkage when each respondent contributes only ~20 choices. With more tasks per respondent this tightens: at 100 tasks, $V_\beta$ → ~[0.43, 0.46, 0.47] and per-respondent recovery → ~0.97, confirming the from-scratch sampler is unbiased.
The recovery scatter shows it: posterior-mean $\beta_i$ track the truth along the 45° line, with the mild flattening that shrinkage produces.
Next (Section 2): the bayesm margarine household scanner panel — 516 households choosing among 10 margarine products, with demographics as the upper-level $z_i$.
2. Real data: bayesm margarine household scanner panel¶
516 households, 4,470 purchases, each a choice among 10 margarine products (6 stick + 4 tub brands) with the 10 shelf prices observed at each occasion. Household demographics: income, family size, college.
- Design: 9 brand intercepts (Parkay stick = base) + a generic price coefficient → $k=10$ part-worths per household.
- Hierarchy: $\beta_i \sim N(\Delta' z_i, V_\beta)$ with $z_i = [1,\ \text{income},\ \text{family size},\ \text{college}]$.
- Only ~8.7 purchases per household, so the $\beta_i$ are heavily pooled through the hierarchy — exactly what the model is for.
The key questions: how price-sensitive are households, how much does that sensitivity vary, and do demographics explain the variation?
# bayesm margarine: build per-household (X, y) and household covariates Z
cp = pd.read_csv('margarine_choice.csv')
dm = pd.read_csv('margarine_demos.csv')
pricecols = ['PPk_Stk', 'PBB_Stk', 'PFl_Stk', 'PHse_Stk', 'PGen_Stk',
'PImp_Stk', 'PSS_Tub', 'PPk_Tub', 'PFl_Tub', 'PHse_Tub']
prods = ['Pk_Stk', 'BB_Stk', 'Fl_Stk', 'Hse_Stk', 'Gen_Stk',
'Imp_Stk', 'SS_Tub', 'Pk_Tub', 'Fl_Tub', 'Hse_Tub']
J, k = 10, 10
D = np.zeros((J, J - 1)) # brand dummies, product 1 (Parkay stick) = base
for j in range(1, J):
D[j, j - 1] = 1
data, hh = [], []
for h, grp in cp.groupby('hhid'):
P = grp[pricecols].values; y = grp['choice'].values - 1; T = len(grp)
X = np.zeros((T, J, k))
for t in range(T):
X[t, :, :9] = D
X[t, :, 9] = P[t] # price of each alternative
data.append((X, y.astype(int))); hh.append(h)
dm2 = dm.set_index('hhid').loc[hh]
income = dm2['Income'].values
Iz = (income - income.mean()) / income.std()
Fz = (dm2['Fam_Size'].values - dm2['Fam_Size'].mean()) / dm2['Fam_Size'].std()
Z = np.column_stack([np.ones(len(hh)), Iz, Fz, dm2['college'].values])
print(f"households={len(data)} occasions={sum(len(y) for _, y in data)} "
f"avg purchases/hh={np.mean([len(y) for _, y in data]):.1f} k={k} q={Z.shape[1]}")
households=516 occasions=4470 avg purchases/hh=8.7 k=10 q=4
out_m = hier_mnl_gibbs(data, Z, R=4000, burn=1000, seed=1) # ~45s
sm = summary(out_m); Dm = sm['Delta']; Bm = sm['B']
print(f"RW-Metropolis acceptance = {out_m['accept']:.3f}\n")
print("baseline brand intercepts (relative to Parkay stick):")
print(f" {prods[0]+'(base)':16s} 0.00")
for m in range(9):
print(f" {prods[m+1]:16s} {Dm[0, m]:6.2f}")
pr = out_m['Delta'][:, :, 9] # group-level coefficients on the price part-worth
print(f"\nprice coefficient (baseline) : {pr[:,0].mean():6.2f} "
f"[{np.percentile(pr[:,0],2.5):.2f}, {np.percentile(pr[:,0],97.5):.2f}]")
print(f" income -> price : {pr[:,1].mean():+.2f} "
f"[{np.percentile(pr[:,1],2.5):+.2f}, {np.percentile(pr[:,1],97.5):+.2f}] (>0: richer = LESS price-sensitive)")
print(f" family -> price : {pr[:,2].mean():+.2f} "
f"[{np.percentile(pr[:,2],2.5):+.2f}, {np.percentile(pr[:,2],97.5):+.2f}] (<0: bigger family = MORE price-sensitive)")
print(f"\nprice-sensitivity heterogeneity sd = {np.sqrt(out_m['Vb'][:, 9, 9].mean()):.2f}")
print(f"all {len(Bm)} household price coefficients negative: {(Bm[:,9] < 0).mean()*100:.0f}%")
RW-Metropolis acceptance = 0.310 baseline brand intercepts (relative to Parkay stick): Pk_Stk(base) 0.00 BB_Stk -1.07 Fl_Stk -0.67 Hse_Stk -2.72 Gen_Stk -5.15 Imp_Stk -2.62 SS_Tub -1.23 Pk_Tub 0.10 Fl_Tub -0.18 Hse_Tub -5.08 price coefficient (baseline) : -8.90 [-9.45, -8.43] income -> price : +0.72 [+0.37, +1.12] (>0: richer = LESS price-sensitive) family -> price : -0.70 [-1.11, -0.27] (<0: bigger family = MORE price-sensitive) price-sensitivity heterogeneity sd = 2.75 all 516 household price coefficients negative: 100%
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
# (a) heterogeneity in price sensitivity across households
ax[0].hist(Bm[:, 9], bins=30, color='steelblue', alpha=0.75, edgecolor='white')
ax[0].axvline(Bm[:, 9].mean(), color='k', ls='--', lw=1.5, label=f'mean = {Bm[:,9].mean():.1f}')
ax[0].axvline(0, color='gray', lw=0.8, ls=':')
ax[0].set_xlabel('household price coefficient $\\beta_{i,\\mathrm{price}}$')
ax[0].set_ylabel('# households')
ax[0].set_title('Every household is price-averse — by widely varying degrees')
ax[0].legend(fontsize=8)
# (b) demographic gradient: price sensitivity vs income
ax[1].scatter(income, Bm[:, 9], s=14, alpha=0.5, color='steelblue', edgecolor='none')
b1, b0 = np.polyfit(income, Bm[:, 9], 1)
xs = np.linspace(income.min(), income.max(), 50)
ax[1].plot(xs, b0 + b1 * xs, 'r-', lw=2.2, label=f'slope = {b1:+.3f}')
ax[1].set_xlabel('household income ($000s)')
ax[1].set_ylabel('household price coefficient (closer to 0 = less sensitive)')
ax[1].set_title('Richer households are LESS price-sensitive')
ax[1].legend(fontsize=8)
fig.suptitle('Hierarchical MNL on margarine — price-sensitivity heterogeneity & its demographics',
fontsize=12)
plt.tight_layout()
plt.savefig(r'C:\Users\user\project\hier_mnl_margarine.png', dpi=120, bbox_inches='tight')
plt.show()
2. Results — margarine¶
Price sensitivity (the headline):
- Baseline price coefficient ≈ −8.9 [−9.4, −8.4] — strong aversion; raising a product's price sharply lowers its choice probability.
- Income → price = +0.72 [0.37, 1.12] — credibly positive: higher-income households are less price-sensitive (their price coefficient is closer to zero).
- Family size → price = −0.70 [−1.11, −0.27] — credibly negative: bigger families are more price-sensitive.
- Heterogeneity sd ≈ 2.75; every one of the 516 households has a negative price coefficient (mean −9.0, ranging roughly −2 to −15).
Brand preferences (intercepts vs. Parkay stick): Parkay (stick & tub) and Fleischmann are the most-preferred; generics and the cheap house tub sit far below (−5) — households need a big price discount to switch to them.
What the hierarchy buys us: with only ~8.7 purchases per household, no household's preferences are estimable alone; the multilevel model pools them — shrinking noisy individual estimates toward demographically-predicted means while still letting price sensitivity vary widely across households (left panel). The right panel shows the demographic gradient: the cloud of household price coefficients trends toward zero as income rises.
On IIA: the model still imposes IIA conditional on each household's $\beta_i$, but the population mixes a wide range of price sensitivities and brand tastes, so aggregate substitution between products does not obey IIA. To also relax IIA within a household's choice (e.g., correlated stick-vs-tub errors) you'd move to multinomial probit — see the MNP project (MNP_* notebooks).
Next (Section 3): cross-check against R bayesm::rhierMnlRwMixture on the same data.
# Part-worths in dollars: willingness to pay = brand part-worth / |price part-worth|
# (scale-free, so directly comparable across logit and probit models)
zbar = Z.mean(0) # average household covariate profile
popb = np.einsum('gqk,q->gk', out_m['Delta'], zbar) # population-mean part-worths, per draw
wtp = popb[:, :9] / np.abs(popb[:, 9:10]) # $ premium vs Parkay stick, per brand
print("Willingness to pay vs Parkay stick (dollars; + = preferred over Parkay):\n")
print(f" {'Parkay_Stk':<10} 0.00 (base)")
for m in range(9):
print(f" {prods[m+1]:<10} {wtp[:,m].mean():+.2f} [{np.percentile(wtp[:,m],2.5):+.2f}, {np.percentile(wtp[:,m],97.5):+.2f}]")
print(f"\n(price part-worth = {popb[:,9].mean():.2f} utils per dollar)")
Willingness to pay vs Parkay stick (dollars; + = preferred over Parkay): Parkay_Stk 0.00 (base) BB_Stk -0.12 [-0.15, -0.10] Fl_Stk -0.07 [-0.12, -0.02] Hse_Stk -0.31 [-0.35, -0.28] Gen_Stk -0.59 [-0.64, -0.54] Imp_Stk -0.28 [-0.32, -0.21] SS_Tub -0.12 [-0.17, -0.07] Pk_Tub +0.01 [-0.04, +0.06] Fl_Tub -0.01 [-0.08, +0.07] Hse_Tub -0.56 [-0.62, -0.50] (price part-worth = -9.00 utils per dollar)
Part-worths in dollars (willingness to pay)¶
Dividing each brand's part-worth by the price part-worth converts utils into dollars — the premium an average household would pay for that brand over Parkay stick:
| brand | WTP vs Parkay stick |
|---|---|
| Parkay tub | ≈ $0.00 (interchangeable with Parkay stick) |
| Fleischmann (stick/tub) | −$0.01 to −$0.07 |
| Blue Bonnet, Shedd's | −$0.12 |
| Imperial, House stick | −$0.28 to −$0.31 |
| Generic, House tub | −$0.56 to −$0.59 |
So brand equity is worth up to ~60¢: a generic or house-brand tub needs a ~$0.57 discount before the average household will switch from Parkay. WTP is scale-free, which makes it the right quantity for comparing across model families.
Comparison with the MNP project (multinomial probit, same margarine data)¶
The MNP_* notebooks fit multinomial probit to the same data (pooled, top-4 then all-10, with an estimated error covariance $\Sigma$). Three points:
1. Brand valuations agree (scale-free WTP).
| brand | hier-MNL (logit) | MNP M3 (probit) |
|---|---|---|
| House stick | −$0.31 | −$0.29 |
| Blue Bonnet | −$0.12 | −$0.15 |
| Shedd's tub | −$0.12 | −$0.07 |
| Parkay tub | +$0.01 | −$0.09 |
| House tub | −$0.56 | −$0.75 |
Despite totally different error structures, the two model families put essentially the same dollar value on brands.
2. Price scale differs as expected. Logit price part-worth ≈ −9.0 vs. probit ≈ −4.1 — a ratio of ~2.2, the usual logit-vs-probit scaling (logit's larger error variance inflates raw coefficients). Scale-free quantities (WTP, choice probabilities) are unaffected, which is why the WTP table lines up.
3. IIA — the real difference. This is the payoff of the comparison. MNP estimates a full error covariance $\Sigma$ and finds non-proportional substitution: cutting Blue Bonnet's price pulls volume mostly from the other challengers (House, Shedd), with Parkay losing only ~17% of the diverted volume — versus the ~53% a plain logit (IIA) would wrongly predict.
The hierarchical MNL relaxes IIA a different way: random coefficients ($V_\beta$) make challenger-buyers and Parkay-loyalists distinct segments, so aggregate substitution is non-proportional too — but conditional on a household it is still IIA, and it cannot represent residual format correlation (stick vs. tub) the way MNP's $\Sigma$ does directly.
Two complementary routes past IIA: heterogeneous tastes (hierarchical MNL, $V_\beta$) vs. correlated errors (MNP, $\Sigma$). Same brand valuations; different — and complementary — substitution machinery.