Fit a hierarchical Bayesian multinomial logit (HMNL)

Runs the adaptive RW-Metropolis-within-Gibbs sampler for the hierarchical (random-coefficients, panel or cross-sectional) multinomial logit with a BLP-style alternative-level random effect: $$U_{ijt} = x_{ijt}'\gamma_i + \delta_j + \epsilon_{ijt}, \qquad U_{iot} = \epsilon_{iot},$$ with i.i.d. Gumbel shocks (including on the implicit outside option, whose systematic utility is 0), $\beta_i \sim N(b, W)$ over the structural covariates ($\gamma_{ik} = \beta_{ik}$ or $\exp(\beta_{ik})$ per rc_dist), and $\delta_j = z_j'\theta + \xi_j$, $\xi_j \sim N(0, \sigma_d^2)$. Partial pooling shrinks each $\delta_j$ toward its characteristics-based mean $z_j'\theta$; the outside option anchors the level of $\delta$ (mean utility relative to the outside good), so no base alternative or sum-to-zero constraint is needed.

Usage

run_hmnlogit(
  data = NULL,
  id_col = NULL,
  alt_col = NULL,
  choice_col = NULL,
  covariate_cols = NULL,
  person_col = NULL,
  alt_covariate_cols = NULL,
  outside_opt_label = NULL,
  cf_residual_col = NULL,
  input_data = NULL,
  include_outside_option = TRUE,
  rc_dist = NULL,
  prior = list(),
  mcmc = list(),
  chains = 1,
  keep_beta_i = c("means", "draws", "none"),
  keep_data = TRUE
)

Arguments

data: Data frame (convenience pathway). Supply either data (with the column names) or input_data, not both.
id_col: Name of the column identifying choice situations (tasks). Task ids only need to be unique within a respondent.
alt_col: Name of the column identifying alternatives.
choice_col: Name of the column indicating the chosen alternative (1 = chosen, 0 = not chosen).
covariate_cols: Vector of names of structural covariate columns (the random-coefficient dimensions).
person_col: Name of the respondent column grouping choice situations. NULL (default) makes each choice situation its own respondent.
alt_covariate_cols: Names of alternative-level covariate columns (constant within each alternative) forming the $\delta$ mean function. NULL (default) gives an intercept-only design (P = 1).
outside_opt_label: Label of physical outside-option rows, removed when include_outside_option = TRUE (the outside good is implicit).
cf_residual_col: Name of a first-stage residual column (control function for an endogenous covariate), appended to X. Default NULL.
input_data: A choicer_data_hmnl object from prepare_hmnl_data() (advanced pathway).
include_outside_option: Logical; if TRUE (default) an implicit outside option with systematic utility 0 is part of every choice set.
rc_dist: Integer vector, one entry per column of covariate_cols: 0 for a normal random coefficient, 1 for log-normal (the coefficient enters utility as exp(beta_ik); hierarchy normal on the log scale). Default NULL is all-normal. Automatically aligned through dropped columns; a cf_residual_col coordinate is always normal.
prior: Named list overriding prior defaults: b_bar (0), A (0.01 I), nu (K + 3), V (nu I), theta_bar (0), A_theta (0.01 I), sd_prior (list(half_cauchy = TRUE, s_d = 1, c0 = 3, d0 = 3)).
mcmc: Named list overriding MCMC defaults: R (10000), burn (R %/% 5), thin (1), seed (drawn via sample.int() so set.seed() governs), trace (0), s_init (2.38 / sqrt(K)), accept_target (0.234).
chains: Number of independent chains (seeds offset by 1, run sequentially). Chain 1 provides the reported draws; all chains feed the rank-normalized split-R-hat table and the retained chains field (all per-chain b/w_vech/delta/theta/sigma_d2/loglik draws, consumed by ess(), mcse(), and traceplot()).
keep_beta_i: "means" (default) stores posterior means/SDs of the individual-level $\beta_i$; "draws" additionally stores the full (K, N, R_keep) draw cube (memory-guarded, budgeted per chain); "none" stores neither.
keep_data: Logical; keep the prepared data on the fit (default TRUE, needed by post-estimation methods).

Value

A choicer_hmnl object (classed c("choicer_hmnl", "choicer_hb")) with posterior summaries (coefficients, se, vcov for $b$; theta_summary; sigma_d2_summary; W_mean; delta and xi quality-ladder tables; beta_i), the raw thinned draws (chain 1), acceptance diagnostics in accept, the rank-normalized split-R-hat table in rhat, all chains' retained draws in chains, and sampler metadata.

Details

Initialization. $\beta_i$ start at the pooled MNL maximum likelihood estimate over the structural covariates (log-normal coordinates transformed to the chain scale with a warn-and-clamp at 0.05); $\delta$ starts at shrunk log choice-share contrasts against the outside option; $\theta$ at the OLS regression of the initial $\delta$ on Z.

Priors. $b \sim N(b\_bar, A^{-1})$, $W \sim IW(\nu, V)$, $\theta \sim N(\theta\_bar, A_\theta^{-1})$, and $\sigma_d \sim$ half-Cauchy$(0, s_d)$ via the Makalic-Schmidt scale mixture (set sd_prior$half_cauchy = FALSE for a plain $IG(c_0, d_0)$ on $\sigma_d^2$).

Endogeneity. If a price-like covariate is endogenous (correlated with $\xi_j$), supply a first-stage residual via cf_residual_col (Petrin & Train 2010); posterior uncertainty does NOT propagate first-stage estimation error.

For beta_i draws at very large scale beyond the memory guard's threshold, a future disk-streaming path (writing each kept slice to disk instead of retaining it in memory) is on the roadmap but not built in this phase; users needing per-respondent draws at that scale should reduce R, reduce chains, or use keep_beta_i = "means".

Examples

# \donttest{
sim <- simulate_hmnl_data(N = 100, T = 3, J = 4, seed = 42)
fit <- suppressWarnings(run_hmnlogit(sim$data, "task", "alt", "choice", c("x1", "x2"),
                    person_col = "pid", alt_covariate_cols = "z1",
                    mcmc = list(R = 500, burn = 200)))
#> MCMC run time 0h:0m:0.04s
summary(fit)
#> Hierarchical Bayesian Multinomial Logit (HMNL) model
#> 
#> Population coefficients b (posterior):
#> Parameter        Mean         SD       2.5%     Median      97.5%
#> x1           0.753790   0.161606   0.451668   0.734653   1.088809
#> x2          -0.597166   0.213346  -0.980699  -0.598832  -0.192214
#> 
#> Delta mean function theta (posterior):
#> Parameter          Mean         SD       2.5%     Median      97.5%
#> (Intercept)    0.414833   0.184127   0.070274   0.441049   0.711342
#> z1            -0.408093   0.257467  -0.969343  -0.429151   0.130759
#> 
#> Alternative-effect variance (posterior):
#> Parameter        Mean         SD       2.5%     Median      97.5%
#> sigma_d^2    0.071588   0.297311   0.000085   0.009473   0.599862
#> 
#> Quality ladder (delta = mean utility vs the outside option; xi = delta - z'theta):
#>  alternative delta_mean delta_sd xi_mean  xi_sd
#>            1     0.0568   0.1506 -0.0195 0.1403
#>            2     0.0780   0.1652  0.0199 0.1702
#>            3     0.5801   0.1510 -0.0092 0.2189
#>            4     0.1382   0.1493 -0.0069 0.1319
#> 
#> Convergence diagnostics (1 chain, 300 draws each)
#> Block                R-hat  ESS_bulk  ESS_tail  MCSE(mean)
#> b[x1]                1.022        15        36      0.0417
#> b[x2]                1.192         9        21      0.0695
#> theta[(Intercept)]   1.009        20       141      0.0410
#> theta[z1]            1.190         4        77      0.1297
#> sigma_d^2            1.206         5        17      0.1347
#> delta (J=4)         1.291*        3*       10*         —
#> *worst: delta[2]
#> Acceptance: beta 0.24, delta 0.51
#> 
#> Respondents: 100  Choice situations: 300  Alternatives: 4 
#> Draws kept: 300  Chains: 1 
#> MCMC run time 0h:0m:0.04s 
coef(fit, component = "delta")
#>          1          2          3          4 
#> 0.05678103 0.07797479 0.58013663 0.13823079 
# }