Hierarchical Bayes: panel tastes, product effects, and entry

In choicer: Discrete Choice Models for Economic Applications

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4)
options(digits = 4)

choicer's hierarchical Bayesian models put random effects at the two levels where applied demand work actually needs them. For person $i$ facing alternative $j$ in choice situation $t$,

$$U_{ijt} = x_{ijt}'\beta_i + \delta_j + \varepsilon_{ijt}, \qquad \beta_i \sim N(b, W), \qquad \delta_j = z_j'\theta + \xi_j, \quad \xi_j \sim N(0, \sigma_d^2).$$

• Respondent level. Each person carries one taste vector $\beta_i$ across all of their choice situations. This is the genuine panel mixed logit: by contrast, run_mxlogit() maximizes a cross-sectional simulated likelihood in which every choice situation is integrated over tastes separately, with no draw held fixed across a person's likelihood contributions (see the mixed logit vignette). When people are observed repeatedly — survey panels, scanner data, repeated hospital admissions — the repetition, together with within-person attribute and menu variation, helps distinguish persistent taste heterogeneity $W$ from choice-level noise, and the hierarchical model is the one that uses it.

• Alternative level. Each alternative carries a mean-utility effect $\delta_j$, partially pooled toward a characteristics-based mean $z_j'\theta$ with one scalar variance $\sigma_d^2$ regardless of $J$. The deviation $\xi_j$ is the micro-data counterpart of the unobserved product quality in Berry, Levinsohn and Pakes (1995). Partial pooling matters when per-alternative data are thin — with hundreds of hospitals or schools, fixed alternative-specific constants drown in incidental-parameter noise, while $\delta_j$ shrinks each alternative toward what its observed characteristics predict. And because $\theta$ and $\sigma_d^2$ describe alternatives as a population, the model has a posterior predictive for alternatives it has never seen — which is what makes entry counterfactuals possible below.

Both models require the implicit outside option (include_outside_option = TRUE): the outside good, with systematic utility zero, anchors the location of $\delta$, so $\delta_j$ is mean utility relative to not choosing any inside alternative. Estimation is by Gibbs sampling in C++ — an adaptive Metropolis-within-Gibbs for the logit (run_hmnlogit()), a fully conjugate Albert-Chib sampler for the probit (run_hmnprobit()).

The hierarchy is shared; the utility-shock model is not:

| | HMNL (run_hmnlogit) | HMNP (run_hmnprobit) | |---|---|---| | Choice shock | iid Type-I extreme value | iid normal in utility levels, including a stochastic outside option | | Persistent tastes | Normal or lognormal coordinates in beta_i | Normal coordinates only | | Scale | Logit scale fixed by the EV1 distribution | Free expanded sigma2; every reported draw is scale-normalized | | Sampler | Adaptive random-walk Metropolis-within-Gibbs | Fully conjugate Albert-Chib augmentation | | Welfare | Logsum and consumer surplus | No expected-maximum/logsum welfare implementation |

In particular, HMNP is not the hierarchical analogue of the unrestricted utility-difference covariance estimated by run_mnprobit(). HMNP gains a conjugate, scalable panel sampler by imposing iid utility-level normal shocks; its flexible substitution comes from persistent person tastes and pooled alternative effects, not an unrestricted shock covariance.

library(choicer)
set_num_threads(2)

Simulate a panel

simulate_hmnl_data() draws a panel with known parameters: N people, T choice situations each, J inside alternatives, structural covariates x1 and x2, and one alternative-level covariate z1 feeding the $\delta$ mean function.

sim <- simulate_hmnl_data(N = 250, T = 6, J = 6, seed = 42)
head(sim$data, 8)

The two-level structure is visible in the identifiers: pid indexes people (who share a taste vector) and task indexes choice situations. A task whose inside rows are all choice = 0 is one where the outside option was chosen.

Fit

Point person_col at the person identifier — that is what groups choice situations into a panel. With person_col = NULL every situation would be its own respondent, the cross-sectional limit of the same model. We run two chains; they feed the convergence diagnostics that summary() reports.

set.seed(99)
fit <- run_hmnlogit(
  data               = sim$data,
  id_col             = "task",
  alt_col            = "alt",
  choice_col         = "choice",
  covariate_cols     = c("x1", "x2"),
  person_col         = "pid",
  alt_covariate_cols = "z1",
  chains             = 2,
  mcmc               = list(R = 24000, burn = 4000, thin = 8)
)
summary(fit)

Reading the output from the top: the population mean tastes b, the $\delta$ mean function theta, and the alternative-effect variance $\sigma_d^2$, each summarized by their posterior. The quality ladder is the per-alternative posterior of $\delta_j$ and $\xi_j$ — which alternatives deliver more mean utility than their characteristics predict. The convergence table and acceptance rates are discussed next.

One storage contract matters when chains > 1: the top-level fit$draws, coefficient summaries, quality ladder, and post-estimation methods use chain 1. fit$chains retains the hierarchical draws from every requested chain, and the convergence table uses all of them. Thus additional chains diagnose whether the reported chain-1 posterior is reproducible; choicer v0.2.0 does not pool chains automatically for posterior summaries or policy calculations.

Convergence is part of the result

MCMC output is only evidence about the posterior if the chains have mixed. The consolidated table in summary() reports, per parameter, the rank-normalized split R-hat, bulk and tail effective sample sizes, and the Monte Carlo standard error of the posterior mean (Vehtari et al. 2021), plus one worst-case row spanning all $J$ alternative effects. The $\delta$ block is the one to watch: it is updated by a serial random-walk Metropolis sweep (its full conditionals are coupled through the softmax denominators), so it mixes more slowly than the conjugate blocks. If its R-hat or ESS looks poor, run longer — the fit also warns at estimation time when any tracked parameter fails the check.

Trace plots make the same point visually:

traceplot(fit, block = "b")

traceplot(fit, block = "delta")

The same diagnostics are available programmatically, on any block, from the retained per-chain draws in fit$chains:

b_chains <- lapply(fit$chains, function(ch) ch$b)
rhat(b_chains, rank = TRUE)
ess(b_chains)
mcse(b_chains)

Two chains keep this vignette build manageable and permit between-chain checks; for a serious empirical run, four or more chains are often a better default. choicer offsets the RNG seed across chains but currently gives them the same data-driven initialization, so this is not an overdispersed-start diagnostic. Report the number of people, tasks per person and alternatives; all prior scales; HMNL acceptance rates; rank-normalized R-hat; bulk and tail ESS; MCSE; trace plots; and posterior-predictive shares. Show prior sensitivity for W and sigma_d2. For entry, defend why the entrant's residual quality is exchangeable with incumbents rather than treating the posterior predictive as a data-free ASC.

Did we recover the truth?

Because the data are simulated, we can line the posterior up against the generating parameters:

recovery_table(fit, sim)

The population means b, the taste variances diag(W), and the realized $\delta_j$ ladder recover tightly. Note the pattern in the theta and sigma_d rows: with only $J = 6$ alternatives, the mean-function regression of $\delta$ on $z$ has six observations, so those posteriors are wide and lean on the prior — more alternatives is what sharpens them. The level of $\delta$ has its own identification story: it is pinned by the outside-option share, so with a small outside share the level posterior is diffuse while the cross-alternative contrasts stay tight.

Willingness to pay and predicted shares

Treating x2 as the price-like attribute, wtp() forms the per-draw ratio of population-mean utility coefficients and summarizes its posterior — a median and quantile interval rather than a delta-method approximation, so the ratio's skewness is carried through honestly:

wtp(fit, price_var = "x2")

predict() integrates shares over both the taste distribution and the chain-1 posterior, returning posterior-predictive intervals (and the outside share):

set.seed(7)
predict(fit, n_draws = 200)

A posterior-predictive check compares observed choice shares with the model's predictive distribution — a first-pass reality check on fit:

ppc_shares(fit, n_draws = 200)

Individual-level tastes are available too: fit$beta_i stores per-person posterior summaries (or full draws with keep_beta_i = "draws"), and predict(fit, level = "individual") conditions on them.

A policy counterfactual

Cut the price x2 of alternative 1 by 0.25 and re-predict — no refitting. Reusing the RNG seed makes the baseline and counterfactual share calculations use the same posterior and taste draws. This common-random-number pairing reduces Monte Carlo noise in the contrast; it does not make a finite-draw integral exact:

cf <- sim$data
cf$x2[cf$alt == 1] <- cf$x2[cf$alt == 1] - 0.25

set.seed(7)
base_shares <- predict(fit, n_draws = 200)
set.seed(7)
cf_shares <- predict(fit, newdata = cf, n_draws = 200)

data.frame(
  alternative    = base_shares$alternative,
  baseline       = round(base_shares$share, 3),
  counterfactual = round(cf_shares$share, 3)
)

The welfare change is the posterior of the compensating variation — the logsum difference divided by the marginal utility of income, draw by draw:

set.seed(7)
cs <- consumer_surplus(fit, price_var = "x2", newdata = cf, n_draws = 200)
attr(cs, "cv")

The three numbers are the lower quantile, posterior median, and upper quantile of the sum of compensating variation over the prediction tasks. Supply weights = rep(1 / fit$nobs, fit$nobs) for an equally weighted mean, or a substantively justified task-weight vector for another aggregate. The function does not normalize user-supplied weights.

Unlike a delta-method standard error, this interval carries posterior uncertainty in the logsum and in the population-mean marginal utility of income $-\bar\gamma_{price}$ used as the denominator. This is a population-mean money metric, not the posterior distribution of person-specific compensating variation when marginal utilities of income differ across people. The taste distribution still enters the integrated logsum, and for a lognormal price coordinate it also enters $\bar\gamma_{price} = \exp(b + W_{kk}/2)$. Report that aggregation choice explicitly in applications with price heterogeneity; an aggregate fixed-sign denominator does not by itself rule out individual price coefficients near or across zero. (logsum() and consumer_surplus() are available for the hierarchical logit only; the probit expected-maximum counterpart is on the roadmap.)

An entry counterfactual

The distinctive payoff of the BLP-style alternative level. A model with fixed alternative-specific constants is silent about an alternative it has never seen — it has no ASC for the entrant and no principled way to invent one. Here the entrant is a draw from the estimated population of alternatives: given its characteristics $z_{\text{new}}$, each posterior draw assigns it $\delta_{\text{new}} \sim N(z_{\text{new}}'\theta_r,\; \sigma_{d,r}^2)$.

Add the entrant's rows to the data — same layout, a new alt label, its z1 value, choice = 0 — and predict:

entrant <- sim$data[sim$data$alt == 1, ]
entrant$alt <- 99L
entrant$z1 <- 0.4
entrant$choice <- 0L
entry_data <- rbind(sim$data, entrant)

set.seed(7)
predict(fit, newdata = entry_data, n_draws = 200)

The entrant takes share from every incumbent and from the outside good, and its credible interval is typically wider than the incumbents' — appropriately so. The model knows the entrant's observed characteristics but not its $\xi$, so the prediction integrates over $\xi_{\text{new}} \sim N(0, \sigma_d^2)$: the uncertainty about unobserved quality can be a dominant uncertainty about an entrant, and the posterior predictive exposes it rather than hiding it. The maintained assumption is exchangeability — the entrant's unobserved quality is a draw from the same population as the incumbents'.

Price endogeneity

If a price-like covariate is correlated with the unobserved quality $\xi_j$ — the classic demand-estimation concern — the estimates are exogenous only conditional on $Z$. The data preparations accept a control-function residual (cf_residual_col, Petrin and Train 2010): regress price on instruments outside the package, and pass the first-stage residual so it enters utility as an ordinary covariate. The package does not run the first stage, and posterior uncertainty does not propagate first-stage estimation error; joint Bayesian IV is on the roadmap.

The probit sibling

run_hmnprobit() estimates the same two-level structure with iid normal utility shocks instead of extreme-value ones. The sampler is fully conjugate (Albert-Chib data augmentation — no Metropolis step, no acceptance rates to tune), and because it works in un-differenced utility space, unbalanced choice sets pose no problem. The probit has a free scale, handled by parameter expansion: a non-identified $\sigma^2$ chain wanders by design, and every reported quantity is normalized per draw, so all summaries are on the identified scale. Do not diagnose convergence from the wandering, unidentified raw sigma2 chain alone; diagnose the normalized structural and predictive quantities that enter the economic conclusions.

simp <- simulate_hmnp_data(N = 250, T = 5, J = 6, seed = 42)

set.seed(99)
fitp <- run_hmnprobit(
  data               = simp$data,
  id_col             = "task",
  alt_col            = "alt",
  choice_col         = "choice",
  covariate_cols     = c("x1", "x2"),
  person_col         = "pid",
  alt_covariate_cols = "z1",
  chains             = 2,
  mcmc               = list(R = 30000, burn = 5000, thin = 10)
)
summary(fitp)

The identified probit coefficients live on a different scale from the logit's. A common variance-matching rule notes that the EV1 shock has standard deviation $\pi/\sqrt{6} \approx 1.28$ against the probit's 1, and therefore multiplies probit coefficients by $\pi/\sqrt{6}$ for a rough comparison. This is not an exact transformation: normal and Type-I extreme-value shocks have different shapes, and coefficients also need matched utility specifications and data:

rbind(
  probit         = coef(fitp),
  "logit scale"  = coef(fitp) * pi / sqrt(6)
)

The choicer_hb post-estimation suite — predict() (probabilities by a deterministic fixed-node one-dimensional Gauss-Hermite approximation), wtp(), elasticities(), diversion_ratios(), ppc_shares(), recovery_table() — is available through the same interfaces on the probit fit. The exceptions are logsum() and consumer_surplus(), which are logit-only as noted above.

Further reading

The full derivations — priors, the Gibbs sweeps, identification, and the implementation contract — are in the math companions: hierarchical MNL and hierarchical MNP. For where these models sit among choicer's estimators, see Choosing among choice models.

References

Berry, S., Levinsohn, J. and Pakes, A. (1995). Automobile prices in market equilibrium. Econometrica, 63(4), 841-890.

Petrin, A. and Train, K. (2010). A control function approach to endogeneity in consumer choice models. Journal of Marketing Research, 47(1), 3-13.

Rossi, P. E., Allenby, G. M. and 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. and Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667-718.

choicer documentation built on Sept. 5, 2026, 1:07 a.m.