The math behind choicer: hierarchical Bayesian multinomial logit

This article gives a self-contained description of the hierarchical Bayesian multinomial logit (HMNL) implemented by run_hmnlogit() and src/hmnlogit.cpp. Shared hierarchical conditionals are in src/hb_internal.h; sampling and RNG primitives are in src/bayes_samplers.h and src/rng.h. Equations below describe the implementation, including its identification convention, adaptive Metropolis targets, retained draws, and numerical approximations.

1. Indices, data, and dimensions

Respondents are $i=1,\ldots,N$. Respondent $i$ faces tasks $t=1,\ldots,T_i$. Let $C_{it}\subseteq\{1,\ldots,J\}$ denote the observed set of inside alternatives and let $o$ denote an implicit outside option. Choice sets may be unbalanced. For an inside row $(i,t,j)$,

• $x_{ijt}\in\mathbb R^K$ is the structural design vector; no ASC dummies are included;
• $z_j\in\mathbb R^P$ is the time-invariant alternative-level design, whose first element is an intercept;
• $y_{it}\in C_{it}\cup\{o\}$ is the observed choice.

person_col defines respondents. If it is absent, every task is assigned its own respondent and $T_i=1$. The data engine stores inside rows only, plus task sizes, respondent task counts, the chosen within-task position (zero for outside), and the global alternative index of every row.

2. Utility, heterogeneity, and likelihood

For coordinate $k$, define the map from the Gaussian chain scale to the utility scale

$$g_k(\beta_{ik})= \begin{cases} \beta_{ik},&\texttt{rc\_dist}_k=0,\\ \exp(\beta_{ik}),&\texttt{rc\_dist}_k=1. \end{cases}$$

Write $\gamma_i=g(\beta_i)$ elementwise. Utilities are

$$U_{ijt}=x_{ijt}'\gamma_i+\delta_j+\varepsilon_{ijt},\qquad U_{iot}=\varepsilon_{iot},$$

where every $\varepsilon$, including the outside shock, is iid type-I extreme value with unit scale and variance $\pi^2/6$. Let

$$v_{ijt}=x_{ijt}'\gamma_i+\delta_j,\qquad D_{it}=1+\sum_{k\in C_{it}}\exp(v_{ikt}).$$

Then

$$p_{ijt}=\Pr(y_{it}=j\mid\beta_i,\delta)=\frac{\exp(v_{ijt})}{D_{it}}, \qquad p_{iot}=\frac1{D_{it}}.$$

With $c_{it}$ the chosen inside row, or absent when outside is chosen, respondent $i$ contributes

$$\ell_i(\beta_i,\delta) =\sum_{t=1}^{T_i}\left[ \mathbf 1\{y_{it}\ne o\}v_{ic_{it}t}-\log D_{it} \right].$$

There are two deliberately distinct random-effect levels:

$$\beta_i\mid b,W\sim N_K(b,W),$$

and

$$\delta_j=z_j'\theta+\xi_j,\qquad \xi_j\mid\sigma_d^2\stackrel{ind}{\sim}N(0,\sigma_d^2).$$

$W$ describes cross-respondent heterogeneity in structural tastes. It never contains the $J$ alternative effects. The scalar $\sigma_d^2$ describes dispersion of alternative quality around the maintained mean function $z_j'\theta$. This separation avoids a dense $(K+J)$-dimensional respondent covariance and supplies partial pooling and a predictive distribution for a new alternative. It is a maintained exchangeability model, not a guarantee that sparse-alternative estimates are consistent. Entry results can be sensitive to $Z$ and to the prior on $\sigma_d$.

3. Identification and empirical content

The EV1 scale is fixed, so utility scale is identified by assumption. The systematic utility of the outside option is fixed to zero. Consequently, shifting every inside $\delta_j$ by a common constant changes inside-versus-outside probabilities and is not an invariance. The package therefore estimates all $J$ inside effects without a base-alternative or sum-to-zero restriction. $\delta_j$ is quality relative to the outside good; the intercept in $\theta$ locates average inside quality relative to outside.

Without the outside option, adding a constant to all available alternatives would leave the softmax unchanged. For this reason v0.2.0 requires include_outside_option = TRUE. Even with the formal anchor, little observed outside choice can make the common level weakly informed while leaving inside-to-inside contrasts much more precise.

Panel repetition reveals persistent respondent heterogeneity because the same $\beta_i$ governs all $T_i$ tasks. With $T_i=1$, $W$ is learned through the mixture shape, covariate and choice-set variation, and the parametric Gaussian assumption. It is therefore generally less empirically disciplined than in a rich panel. Similarly, $\theta$ and $\sigma_d^2$ require enough alternatives and variation in $Z$; many observations of a few products do not substitute for cross-product variation in the second-level regression.

If a price-like regressor is correlated with $\xi_j$, the likelihood is exogenous only conditional on $Z$. cf_residual_col appends a user-supplied first-stage residual to $x_{ijt}$ as a normal-coordinate control function. The package does not estimate that first stage and does not propagate its estimation uncertainty.

4. Priors and posterior kernel

The hierarchy uses

$$b\sim N_K(\bar b,A^{-1}),\qquad W\sim IW_K(\nu,V),\qquad \theta\sim N_P(\bar\theta,A_\theta^{-1}).$$

Here $A$ and $A_\theta$ are precisions. The inverse-Wishart convention is

$$p(W)\propto |W|^{-(\nu+K+1)/2} \exp\{-\tfrac12\operatorname{tr}(VW^{-1})\},$$

which is proper for the engine’s allowed $\nu\ge K$ and has mean $V/(\nu-K-1)$ only when $\nu>K+1$. Defaults are $\bar b=0$, $A=0.01I_K$, $\nu=K+3$, $V=\nu I_K$, $\bar\theta=0$, and $A_\theta=0.01I_P$. Thus the default raw prior mean of $W$ is $(K+3)I_K/2$, not $I_K$.

By default $\sigma_d\sim C^+(0,s_d)$. The implementation uses the Makalic–Schmidt representation, in the inverse-gamma shape–scale convention

$$\sigma_d^2\mid a_d\sim IG(1/2,1/a_d),\qquad a_d\sim IG(1/2,1/s_d^2),$$

where $X\sim IG(a,d)$ means $1/X\sim\operatorname{Gamma}(a,\text{rate}=d)$. Setting sd_prior$half_cauchy = FALSE instead gives $\sigma_d^2\sim IG(c_0,d_0)$.

Up to constants, the posterior density is

$$\begin{aligned} p(&\beta_{1:N},b,W,\delta,\theta,\sigma_d^2,a_d\mid y)\propto \exp\left\{\sum_i\ell_i(\beta_i,\delta)\right\} \prod_i\phi_K(\beta_i;b,W)\\ &\times\prod_j\phi(\delta_j;z_j'\theta,\sigma_d^2) p(b)p(W)p(\theta)p(\sigma_d^2\mid a_d)p(a_d), \end{aligned}$$

with the last two factors replaced by the plain inverse-gamma prior when requested. Prior-predictive checks should cover taste magnitudes, implied shares, $W$, the quality ladder, and entrant quality—not merely individual hyperparameters in isolation.

5. One Metropolis-within-Gibbs iteration

The chain is a systematic scan. It rebuilds $v_{ijt}$, $\exp(v_{ijt})$, $D_{it}$, task log likelihoods, and respondent log likelihoods at the start of every iteration.

5.1 Respondent coefficients

For respondent $i$, the exact conditional target is

$$\log\pi_i(\beta_i)=\ell_i(\beta_i,\delta) -\tfrac12(\beta_i-b)'W^{-1}(\beta_i-b)+\text{constant}.$$

The proposal is symmetric:

$$\beta_i^*=\beta_i+s_iL_i^{-T}z,\qquad z\sim N_K(0,I_K),$$

where $L_iL_i'=H_i+W^{-1}$. Thus its covariance is $s_i^2(H_i+W^{-1})^{-1}$. $H_i$ is fixed once at startup. At the pooled initialization, for each task,

$$H_{it}=\sum_{j\in C_{it}}p_{ijt}x_{ijt}x_{ijt}' -\bar x_{it}\bar x_{it}',\qquad \bar x_{it}=\sum_{j\in C_{it}}p_{ijt}x_{ijt},$$

and $H_i=\sum_tH_{it}$. The omitted outside vector is zero but its probability remains in the denominator. For log-normal coordinate $k$, row and column $k$ of $H_i$ are multiplied by $\exp(\beta_{k,\mathrm{pooled}})$, implementing the local Jacobian transformation from utility to chain scale.

Because the proposal is symmetric,

$$\log\alpha_i=\ell_i(\beta_i^*,\delta)-\ell_i(\beta_i,\delta) +\tfrac12\left[q_i-q_i^*\right],$$

where $q_i=(\beta_i-b)'W^{-1}(\beta_i-b)$. Non-finite candidate likelihoods are rejected. During burn-in only,

$$\log s_i^{(r+1)}=\operatorname{clip}_{[\log .01,\log 10]} \left(\log s_i^{(r)}+(r+1)^{-0.6}\{\min(1,e^{\log\alpha_i})-a_\beta\}\right),$$

where $a_\beta$ defaults to 0.234. Adaptation stops after burn-in.

5.2 Alternative effects

The $\delta_j$ conditionals are not independent because alternatives in the same task share $D_{it}$. They are therefore updated in a strictly serial random-walk sweep. Propose $\delta_j^*=\delta_j+s_{\delta j}z$. Let $e_{ijt}=\exp(v_{ijt})$, $d=\delta_j^*-\delta_j$, and

$$D_{it}^*=D_{it}-e_{ijt}+e_{ijt}e^d.$$

Only tasks containing $j$ change, giving

$$\Delta\ell_j=\sum_{(i,t):j\in C_{it}} \left[d\,\mathbf1\{y_{it}=j\}+\log D_{it}-\log D_{it}^*\right].$$

With $m_j=z_j'\theta$,

$$\log\alpha_{\delta j}=\Delta\ell_j+ \frac{(\delta_j-m_j)^2-(\delta_j^*-m_j)^2}{2\sigma_d^2}.$$

On acceptance the utility and denominator caches are updated immediately, so the next alternative conditions on the new state. $\log s_{\delta j}$ follows the same burn-in Robbins–Monro rule with fixed target 0.44 and the same $[0.01,10]$ scale clamp.

5.3 Conjugate hierarchy

Let $P_b=A+NW^{-1}$ and $r_b=A\bar b+W^{-1}\sum_i\beta_i$. Then

$$b\mid\cdot\sim N_K(P_b^{-1}r_b,P_b^{-1}),$$

$$W\mid\cdot\sim IW_K\left(\nu+N,V+\sum_i(\beta_i-b)(\beta_i-b)'\right).$$

For the alternative mean function, let

$$P_\theta=A_\theta+Z'Z/\sigma_d^2,\qquad r_\theta=A_\theta\bar\theta+Z'\delta/\sigma_d^2.$$

Then $\theta\mid\cdot\sim N_P(P_\theta^{-1}r_\theta,P_\theta^{-1})$. Writing $\xi=\delta-Z\theta$, the half-Cauchy mixture gives

$$\sigma_d^2\mid\cdot\sim IG\left(\frac{J+1}{2},\frac1{a_d}+\frac12\xi'\xi\right), \qquad a_d\mid\cdot\sim IG\left(1,\frac1{s_d^2}+\frac1{\sigma_d^2}\right).$$

Under the fallback prior, $\sigma_d^2\mid\cdot\sim IG(c_0+J/2,d_0+\xi'\xi/2)$.

6. Initialization

The R wrapper first estimates a pooled MNL over the structural columns, with the alternative effects excluded from that pooled optimization. Normal coordinates start at that estimate. A log-normal coordinate is transformed to $\log\max(\hat\gamma_k,0.05)$; estimates below 0.05 trigger a warning and clamp. Every $\beta_i$ and $b$ starts at this chain-scale vector; $W=I_K$ and $\sigma_d^2=a_d=1$.

If $n_j$ is the number of choices of inside alternative $j$ and $n_o$ the outside count,

$$\delta_j^{(0)}=\log\frac{n_j+1/2}{n_o+1/2}.$$

$\theta^{(0)}$ is the QR least-squares solution regressing $\delta^{(0)}$ on $Z$. These are starting values, not pseudo-count priors.

7. Retained quantities and diagnostics

For each chain the kernel retains thinned draws of $b$, vech$(W)$, $\delta$, $\theta$, $\sigma_d^2$, and the log likelihood. Chain 1 supplies the top-level summaries and object$draws; all chains are retained in object$chains. keep_beta_i="means" uses Welford updates for chain-1 respondent posterior means and SDs. "draws" additionally keeps the chain-1 $K\times N\times R_{keep}$ cube returned on the object, although all chains’ cubes coexist temporarily during sequential fitting. The memory guard estimates $1.9\times8KNR_{keep}$ bytes per chain and multiplies by the requested number of chains.

summary() reports rank-normalized folded split-$\widehat R$, bulk ESS, tail ESS, and MCSE for $b$, $\theta$, and $\sigma_d^2$, plus the worst $\widehat R$ among all $J$ columns of $\delta$. Fit-time warnings test every $\delta_j$ as well as the small blocks at $\widehat R>1.01$ or bulk ESS below 400. vech(W) is retained in every chain but is not included in this automatic table or warning; covariance heterogeneity is often among the slowest and weakest-identified quantities, so researchers should pass the per-chain w_vech matrices to rhat()/ess() explicitly and inspect sensitivity. The built-in thresholds are screening rules, not proofs of convergence. Multiple chains, traces, Monte Carlo error, prior sensitivity, and posterior predictive checks remain necessary.

8. Post-estimation mathematics and limits

All methods use evenly selected chain-1 draws. At population level, each posterior draw $r$ is paired with one simulation

$$\beta^{(r)}\sim N(b^{(r)},W^{(r)}),$$

so integration over taste heterogeneity is a one-draw-per-posterior-draw Monte Carlo average, not exact quadrature. At individual level, stored $\beta_i^{(r)}$ draws are used when available; otherwise posterior-mean $\beta_i$ is a plug-in. New alternative $a$ receives an independent draw

$$\delta_a^{(r)}\sim N(z_a'\theta^{(r)},\sigma_d^{2(r)}).$$

These post-estimation simulations use R’s RNG, so set.seed() controls their reproducibility.

For a log-normal coordinate, the posterior-draw population mean coefficient used by wtp() is

$$E(\gamma_k\mid b^{(r)},W^{(r)})=exp\{b_k^{(r)}+W_{kk}^{(r)}/2\};$$

normal coordinates use $b_k^{(r)}$. WTP draw ratios are $-\bar\gamma_a^{(r)}/\bar\gamma_p^{(r)}$ and are summarized by medians and equal-tailed intervals.

For HMNL the conditional expected maximum, apart from the Euler constant common to states, is

$$L_{it}^{(r)}=\log\left(1+\sum_{j\in C_{it}}e^{v_{ijt}^{(r)}}\right).$$

The implementation evaluates this without forming an unsafe $e^{v}$ at the original scale. With

$$m_{it}^{(r)}=\max\left\{0,\max_{j\in C_{it}}v_{ijt}^{(r)}\right\},$$

it computes

$$L_{it}^{(r)}=m_{it}^{(r)}+\log\left[ e^{-m_{it}^{(r)}}+\sum_{j\in C_{it}} e^{v_{ijt}^{(r)}-m_{it}^{(r)}}\right],$$

and obtains the inside and outside probabilities from the same scaled denominator. Thus R-side prediction and welfare remain finite even when a counterfactual utility is outside the direct numerical range of exp().

logsum() averages these simulated draws. consumer_surplus() divides by $-\bar\gamma_{price}^{(r)}$ and reports posterior medians and intervals; counterfactual CV uses the same simulated $\beta$ path for baseline and policy data. For optional non-negative task weights $q_{it}$, the aggregate CV draw is

$$CV^{(r)}=\frac{\sum_{i,t}q_{it} \{L_{it,1}^{(r)}-L_{it,0}^{(r)}\}} {-\bar\gamma_{price}^{(r)}}.$$

Equal task weights are the default. This calculation is appropriate only when the chosen price coefficient has a defensible, resolved sign and units; the weights define aggregation and do not repair a misspecified sampling or price process. The code requires the policy and baseline panels to contain exactly the same choice situations, matches them by (person_col, id_col) rather than by row position, and then applies unnamed weights in the baseline prepared data’s sorted task order. Alternative rows and complete tasks may therefore be reordered without changing the welfare contrast; adding or dropping tasks is rejected because it changes the aggregation population.

elasticities() and diversion_ratios() are aggregate finite-difference objects, not analytic individual derivatives. For each inside alternative $j$, the selected covariate is multiplied by $1+\epsilon$, probabilities are recomputed with common random numbers, and

$$E_{kj}=\frac{\{s_k^{pert(j)}-s_k\}/s_k}{\epsilon},\qquad DR(j\to k)=\frac{s_k^{pert(j)}-s_k}{-(s_j^{pert(j)}-s_j)}.$$

The outside option is included as a receiving alternative. ppc_shares() compares observed take rates with posterior predictive share intervals; matching shares alone is not a comprehensive model check.

9. Computational contract

The chain runs inside one persistent OpenMP region. Cache rebuilds and respondent proposals are work-shared; the $\delta$ sweep, hierarchy draws, recording, and interrupt polling are master-only. Worker code uses hand-written Cholesky and triangular solves and does not call the R API or threaded BLAS. Sufficient statistics are accumulated in fixed respondent or alternative order.

At iteration $r$, make_stream(seed,r,tag) assigns tags $i$ to respondent proposals, $2N+j$ to $\delta_j$, and $2N+J+\{0,1,2,3\}$ to $b,W,\theta,\sigma_d^2$. The unused $N+i$ block is reserved for HMNP coefficient draws. Chains run sequentially with seeds seed, seed+1, and so on. Draws are reproducible for a fixed seed and thread configuration; source comments conservatively allow cross-thread-count differences at floating-point reduction order, so users should not treat byte identity across toolchains or thread counts as an API promise.

References

• Berry, S., Levinsohn, J., & Pakes, A. (1995). Automobile prices in market equilibrium. Econometrica, 63(4), 841–890.
• Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3), 515–534.
• Makalic, E., & Schmidt, D. F. (2016). A simple sampler for the horseshoe estimator. IEEE Signal Processing Letters, 23(1), 179–182.
• Petrin, A., & Train, K. (2010). A control function approach to endogeneity in consumer choice models. Journal of Marketing Research, 47(1), 3–13.
• Roberts, G. O., Gelman, A., & Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. Annals of Applied Probability, 7(1), 110–120.
• Rossi, P. E., Allenby, G. M., & McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley.
• Train, K. E. (2009). Discrete Choice Methods with Simulation (2nd ed.). Cambridge University Press.
• Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved $\widehat R$ for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718.