The mixed (random-coefficients)
logit lets tastes vary across decision units. Instead of a single
coefficient on each random attribute, choicer estimates a
distribution of coefficients. Substitution then reflects both
observed covariates and the estimated distribution of tastes. Estimation
is by simulated maximum likelihood using Halton draws, with the
likelihood, gradient and Hessian evaluated in parallel C++. By default
run_mxlogit() estimates a cross-sectional
simulated likelihood: each choice situation’s probability is integrated
over the taste distribution separately. When the same decision maker is
observed in several choice situations, person_col switches
to the panel likelihood instead (Revelt and Train
1998): one taste draw is held fixed across that person’s situations and
their conditional probabilities are integrated jointly, so repeated
choices discipline the taste distribution directly — see Panel data:
one taste draw per decision maker below.
A useful way to see the mechanism is to write the mixed-logit probability as a logit kernel averaged over the taste distribution, \(P_{ij} = \int L_{ij}(\beta)\, f(\beta)\, d\beta\), where \(L_{ij}(\beta) = \exp(X_{ij}\beta) / \sum_k \exp(X_{ik}\beta)\). Conditional on a draw \(\beta\), the kernel is an ordinary logit. The empirical content of the mixed logit is in the averaging: people with different \(\beta\)’s place different values on the same attributes, so demand leaving an alternative need not go to the same destinations for all consumers. For counterfactual work, the distinction between observed heterogeneity and the estimated mixing distribution is central: flexibility that is not disciplined by the available variation is supplied by the maintained distribution \(f(\beta)\).
A robust recipe for mixed logit: for a large model, warm-start from a
plain MNL (as inst/simulations/mxl_simulation.R does); let
the optimizer work in coordinates scaled by the curvature at the start
values (the default, scale_vars = "bhhh"); and use enough
Halton draws. Here we estimate a full (correlated) covariance of the
random coefficients.
fit <- run_mxlogit(
data = sim$data,
id_col = "id",
alt_col = "alt",
choice_col = "choice",
covariate_cols = c("x1", "x2"), # fixed coefficients
random_var_cols = c("w1", "w2"), # random coefficients
rc_correlation = TRUE, # estimate their full covariance
S = 100L, # Halton draws per choice situation (no person_col here)
draws = "generate", # generate draws on the fly (low memory)
seed = 7L,
scale_vars = "bhhh", # the optimizer's coordinates (the default)
se_method = "bhhh"
)
#> Optimization run time 0h:0m:0.24s
summary(fit)
#> Mixed Logit (MXL) model
#>
#> Parameter Estimate Std.Error z-value Pr(>|z|)
#> x1 0.776408 0.070138 11.0697 0.00e+00 ***
#> x2 -0.675643 0.068318 -9.8897 0.00e+00 ***
#> Sigma_11 0.905183 0.451785 2.0036 4.51e-02 *
#> Sigma_21 0.894700 0.310679 2.8798 3.98e-03 **
#> Sigma_22 1.884033 0.575835 3.2718 1.07e-03 **
#> ASC_1 0.536159 0.078972 6.7893 1.13e-11 ***
#> ASC_2 -0.477738 0.105075 -4.5466 5.45e-06 ***
#> ASC_3 0.591089 0.077545 7.6225 2.49e-14 ***
#> ASC_4 -0.513880 0.105365 -4.8771 1.08e-06 ***
#> ---
#> Signif. codes: '***' 0.001 '**' 0.01 '*' 0.05
#>
#> Random coefficient covariance (Sigma):
#> w1 w2
#> w1 0.9052 0.8947
#> w2 0.8947 1.8840
#>
#> Std. Errors: BHHH (OPG)
#> Log-likelihood: -2435.84
#> AIC: 4889.67 | BIC: 4940.08
#> McFadden R2: 0.106 (adj: 0.102) | Hit rate: 0.452
#> N: 2000 | Parameters: 9
#> Optimization time: 0.24 s
#> Convergence: 1 ( NLOPT_SUCCESS: Generic success return value. )This vignette uses se_method = "bhhh" because the
outer-product calculation is fast and keeps package-build time short.
For final empirical work, compare it with the default analytical-Hessian
standard errors; when the data come from a choice-based or otherwise
weighted sample, use se_method = "sandwich" so the reported
covariance is the robust WESML sandwich rather than the inverse weighted
Hessian. If WESML observations are also dependent within people or
markets, use the weighted cluster sandwich
(se_method = "cluster", with cluster_col) and
have enough independent clusters for its asymptotics.
Tip. For real applications increase the number of draws (
S) until your estimates are stable, and keep the defaultscale_vars = "bhhh". Stability inSis a statistical requirement, not just numerical hygiene: simulated ML with a fixed number of draws is biased — the log of an unbiased probability simulator is not unbiased — and the classical asymptotics requireSto grow with the sample (see the math companion for the conditions). If the solver struggles, pass an explicittheta_init: place the MNL slope and ASC estimates in the corresponding MXL parameter blocks and initialize the additional mean and Cholesky coordinates deliberately. Bounds can keep Cholesky diagonals away from numerically pathological regions. Seeinst/simulations/mxl_simulation.Rfor a fully hardened example.
recovery_table(fit, sim$true_params)
#> <choicer_recovery> model=choicer_mxl level=0.95
#> parameter group true estimate se bias rel_bias_pct z_vs_true
#> <char> <char> <num> <num> <num> <num> <num> <num>
#> 1: x1 beta 0.8000 0.7764 0.0701 -0.0236 -2.949 -0.3364
#> 2: x2 beta -0.6000 -0.6756 0.0683 -0.0756 12.607 -1.1072
#> 3: L_11 sigma 0.0000 -0.0498 0.2496 -0.0498 NA -0.1996
#> 4: L_21 sigma 0.5000 0.9404 0.3057 0.4404 88.078 1.4406
#> 5: L_22 sigma 0.1116 -0.0002 0.3566 -0.1117 -100.136 -0.3133
#> 6: ASC_1 asc 0.5000 0.5362 0.0790 0.0362 7.232 0.4579
#> 7: ASC_2 asc -0.5000 -0.4777 0.1051 0.0223 -4.452 0.2119
#> 8: ASC_3 asc 0.5000 0.5911 0.0775 0.0911 18.218 1.1747
#> 9: ASC_4 asc -0.5000 -0.5139 0.1054 -0.0139 2.776 -0.1317
#> lower_ci upper_ci covers
#> <num> <num> <lgcl>
#> 1: 0.6389 0.9139 TRUE
#> 2: -0.8095 -0.5417 TRUE
#> 3: -0.5389 0.4393 TRUE
#> 4: 0.3412 1.5396 TRUE
#> 5: -0.6991 0.6988 TRUE
#> 6: 0.3814 0.6909 TRUE
#> 7: -0.6837 -0.2718 TRUE
#> 8: 0.4391 0.7431 TRUE
#> 9: -0.7204 -0.3074 TRUEThe beta rows are the fixed coefficients, the
sigma rows describe the covariance of the random
coefficients (its Cholesky elements), and the asc rows are
the alternative-specific constants.
With random coefficients, diversion is mediated by the distribution
of tastes. If the estimated mixing distribution captures economically
meaningful heterogeneity, people who value one alternative tend to value
nearby substitutes on that latent margin. Diversion can therefore depend
on which attribute is changing, so
diversion_ratios() takes a wrt_var:
elasticities(fit, elast_var = "x2")
#> 0 1 2 3 4
#> 0 0 -0.015656 -0.009793 -0.025914 -0.009414
#> 1 0 -0.023929 -0.005506 -0.017254 -0.005718
#> 2 0 -0.010723 -0.016649 -0.017014 -0.007051
#> 3 0 -0.008498 -0.006243 -0.006971 -0.006645
#> 4 0 -0.011956 -0.006111 -0.021578 -0.020221
diversion_ratios(fit, wrt_var = "x2")
#> 0 1 2 3 4
#> 0 0.0000 0.3435 0.3183 0.3533 0.3177
#> 1 0.3112 0.0000 0.2616 0.3195 0.2629
#> 2 0.1789 0.1623 0.0000 0.1658 0.1434
#> 3 0.3342 0.3336 0.2790 0.0000 0.2760
#> 4 0.1757 0.1605 0.1412 0.1614 0.0000
# For a random-coefficient attribute the perturbation coordinate matters.
elasticities(fit, elast_var = "w2", is_random_coef = TRUE)
#> 0 1 2 3 4
#> 0 0 -0.01809 -0.02890 -0.01330 -0.02617
#> 1 0 0.05861 -0.02036 -0.02040 -0.01905
#> 2 0 -0.02383 0.13609 -0.02394 -0.02000
#> 3 0 -0.02050 -0.01863 0.04831 -0.01820
#> 4 0 -0.02493 -0.02164 -0.02110 0.13556
diversion_ratios(fit, wrt_var = "w2", is_random_coef = TRUE)
#> 0 1 2 3 4
#> 0 0.0000 0.5513 -6.311 0.53192 -0.4087
#> 1 -2.0961 0.0000 8.744 0.62491 1.8795
#> 2 1.3819 0.5036 0.000 -0.08724 -0.7105
#> 3 1.4015 -0.4331 1.050 0.00000 0.2397
#> 4 0.3127 0.3782 -2.482 -0.06959 0.0000The rest of the toolkit — predict(), wtp(),
consumer_surplus(), blp() — uses the same
fitted object as in the getting-started
vignette. Prediction and share inversion integrate over the fitted
taste distribution; wtp() respects choicer’s supported
coefficient parameterizations. In v0.2.0,
consumer_surplus() requires a fixed price coefficient
because inverse moments of a random price coefficient need not
exist.
The mixed logit is a genuine generalization of the MNL — if tastes are in fact homogeneous, the estimator simply returns a near-zero variance and you are back to a logit. The issue is not that random coefficients are intrinsically fragile. The issue is identification: the additional substitution structure is carried by the mixing distribution \(f(\beta)\), and \(f(\beta)\) is often hardest to pin down where it matters most for welfare and diversion — in the tails.
Two consequences are worth keeping in front of you:
The tails drive the economics you report. A lognormal price coefficient puts a slice of the population at near-zero price sensitivity, which can give a willingness-to-pay distribution with no finite mean and explosive welfare numbers. An unbounded normal coefficient implies a fraction of consumers with the wrong sign (who prefer paying more). These artifacts come from the assumed shape of \(f\), not from the data, and the estimator will happily contort a tail to match an aggregate moment.
\(f(\beta)\) is hard to
identify and estimate in practice. A single cross-section of
choices — one decision per person, a fixed menu — carries little
information about the spread of tastes. In mixed-logit research,
reliable estimation typically uses repeated choices from the
same individual (panel data), substantial variation in
choice sets or attributes across markets (the BLP setting), or
the rich, designed attribute variation of a stated-preference
experiment. Without one of these, the random-coefficient variances are
weakly identified and the estimates can be fragile.
run_mxlogit() exploits the first source when
person_col is supplied: one taste draw is then held fixed
across a person’s choice situations, so it is their joint
probability that gets integrated over \(f(\beta)\) (Revelt and Train 1998) — see
Panel data: one taste draw per decision maker below. Without
person_col, the default cross-sectional likelihood
integrates each situation separately and does not use that within-person
link; cluster-robust standard errors (cluster_col=, or
vcov(fit, type = "cluster"); see Which standard errors, and when) then repair
the inference for within-person dependence without changing
what the point estimates target. The hierarchical Bayesian logit, run_hmnlogit(), is the Bayesian
alternative to the panel likelihood: the same taste-sharing idea, plus
partial pooling of alternative-level effects and posterior rather than
asymptotic inference.
Estimand. For decision maker \(n\) observed in choice situations \(t = 1, \dots, T_n\), the panel likelihood holds one taste draw \(\beta_n \sim f(\beta \mid \theta)\) fixed across all of \(n\)’s situations and integrates their joint probability (Revelt and Train 1998): \[L_n(\theta) = \int \prod_{t=1}^{T_n} P_{nt}(j_{nt} \mid \beta)\, f(\beta \mid \theta)\, d\beta.\] The target is the same taste distribution \(f(\beta \mid \theta)\); what differs is the maintained model of a person’s repeated choices. The default likelihood integrates every situation on its own — the full likelihood when tastes are redrawn each situation or each person is observed once.
Why repeated choices identify \(\Sigma\). A person who keeps picking the same kind of alternative, situation after situation, is revealing that their tastes sit away from the population average and stay there — that persistence is exactly the variation that identifies \(\Sigma\) most sharply (the marginals identify it too, but only through the shape of the mixture). The cross-sectional likelihood cannot see it: integrating each situation on its own treats a person’s five consistent choices the same as five different people who each chose once. It still targets the same population mixing distribution — situation by situation, each probability is correctly specified, so it stays consistent (a composite, or pseudo-, likelihood) — but it is statistically inefficient, and the loss falls squarely on \(\Sigma\).
Fitting the panel likelihood. Point
person_col at the respondent identifier;
id_col still identifies choice situations, and must stay
unique across the whole data set (situation ids do not restart within
person). S becomes draws per decision maker rather than per
choice situation, since one draw block is now shared across all of a
person’s situations. A product of \(T_n\) per-situation probabilities is more
sharply peaked in \(\beta\) than a
single choice probability and needs more draws to simulate accurately —
the recovery script in
inst/simulations/mxl_panel_simulation.R uses
S = 1000 at T = 10; this vignette keeps
S = 100 for speed.
T = 6 gives each of the 300 decision makers six choice
situations, sharing one realized taste vector
(panel_sim$true_params$gamma_i[, n]) across them.
panel_fit <- run_mxlogit(
data = panel_sim$data,
id_col = "id",
alt_col = "alt",
choice_col = "choice",
covariate_cols = c("x1", "x2"),
random_var_cols = c("w1", "w2"),
rc_correlation = TRUE,
person_col = "pid", # one taste draw per pid, shared across its T tasks
S = 100L, # Halton draws per decision maker
draws = "generate",
seed = 7L
)
#> Optimization run time 0h:0m:0.09s
summary(panel_fit)
#> Mixed Logit (MXL) model
#>
#> Parameter Estimate Std.Error z-value Pr(>|z|)
#> x1 0.867807 0.067285 12.8974 0.00e+00 ***
#> x2 -0.723011 0.069041 -10.4721 0.00e+00 ***
#> Sigma_11 1.241486 0.251382 4.9386 7.87e-07 ***
#> Sigma_21 0.532926 0.187874 2.8366 4.56e-03 **
#> Sigma_22 1.693020 0.301606 5.6134 1.98e-08 ***
#> ASC_1 0.537943 0.080112 6.7149 1.88e-11 ***
#> ASC_2 -0.414377 0.094427 -4.3883 1.14e-05 ***
#> ASC_3 0.485624 0.081056 5.9912 2.08e-09 ***
#> ASC_4 -0.590391 0.098788 -5.9763 2.28e-09 ***
#> ---
#> Signif. codes: '***' 0.001 '**' 0.01 '*' 0.05
#>
#> Random coefficient covariance (Sigma):
#> w1 w2
#> w1 1.2415 0.5329
#> w2 0.5329 1.6930
#>
#> Std. Errors: Analytical Hessian
#> Log-likelihood: -2126.78
#> AIC: 4271.57 | BIC: 4321.03
#> McFadden R2: 0.132 (adj: 0.128) | Hit rate: 0.436
#> N: 1800 | Respondents: 300 | Parameters: 9
#> Optimization time: 0.09 s
#> Convergence: 3 ( NLOPT_FTOL_REACHED: Optimization stopped because ftol_rel or ftol_abs (above) was reached. )recovery_table(panel_fit, panel_sim$true_params)
#> <choicer_recovery> model=choicer_mxl level=0.95
#> parameter group true estimate se bias rel_bias_pct z_vs_true
#> <char> <char> <num> <num> <num> <num> <num> <num>
#> 1: x1 beta 0.8000 0.8678 0.0673 0.0678 8.476 1.0078
#> 2: x2 beta -0.6000 -0.7230 0.0690 -0.1230 20.502 -1.7817
#> 3: L_11 sigma 0.0000 0.1082 0.1012 0.1082 NA 1.0683
#> 4: L_21 sigma 0.5000 0.4783 0.1600 -0.0217 -4.341 -0.1356
#> 5: L_22 sigma 0.1116 0.1907 0.1010 0.0791 70.897 0.7829
#> 6: ASC_1 asc 0.5000 0.5379 0.0801 0.0379 7.589 0.4736
#> 7: ASC_2 asc -0.5000 -0.4144 0.0944 0.0856 -17.125 0.9068
#> 8: ASC_3 asc 0.5000 0.4856 0.0811 -0.0144 -2.875 -0.1774
#> 9: ASC_4 asc -0.5000 -0.5904 0.0988 -0.0904 18.078 -0.9150
#> lower_ci upper_ci covers
#> <num> <num> <lgcl>
#> 1: 0.7359 0.9997 TRUE
#> 2: -0.8583 -0.5877 TRUE
#> 3: -0.0903 0.3066 TRUE
#> 4: 0.1646 0.7920 TRUE
#> 5: -0.0073 0.3887 TRUE
#> 6: 0.3809 0.6950 TRUE
#> 7: -0.5995 -0.2293 TRUE
#> 8: 0.3268 0.6445 TRUE
#> 9: -0.7840 -0.3968 TRUEWhat changes in inference. With
person_col, scores are already per decision maker, so
se_method = "sandwich" (equivalently
vcov(panel_fit, type = "robust")) is robust to
within-person dependence out of the box — there is nothing left within a
person to cluster on. cluster_col is for coarser groups
that nest decision makers (households, markets); it must be constant
within each person. Weights, if supplied, are decision-maker weights,
constant within person, since the objective is \(\sum_n w_n \log L_n\); WESML (choice-based)
weights are not supported with person_col, because a
choice-based weight varies with the alternative a person chose, which
can differ across their situations. nobs() still counts
choice situations; a panel fit additionally reports
n_persons, printed as “Respondents” above.
As a check on the efficiency argument above: refitting the same data
cross-sectionally, clustering on pid for valid inference,
repairs the standard errors’ validity but not their size.
cs_fit <- run_mxlogit(
data = panel_sim$data,
id_col = "id",
alt_col = "alt",
choice_col = "choice",
covariate_cols = c("x1", "x2"),
random_var_cols = c("w1", "w2"),
rc_correlation = TRUE,
cluster_col = "pid", # inference fix only: still one draw per situation
S = 100L,
draws = "generate",
seed = 7L
)
#> Optimization run time 0h:0m:0.17s
se_ratio <- sqrt(diag(vcov(cs_fit, type = "cluster"))) / sqrt(diag(vcov(panel_fit)))
round(se_ratio, 2)
#> x1 x2 L_11 L_21 L_22 ASC_1 ASC_2 ASC_3 ASC_4
#> 1.08 1.11 1.88 1.68 1.64 1.03 1.13 1.05 1.18The ratio sits close to one on the mean-utility parameters
(beta, ASCs — about 1.0-1.2 here) but visibly above one on
the Cholesky block — the cost of never sharing a taste draw within a
person, exactly where a composite likelihood’s efficiency loss
concentrates. (cs_fit’s numerator is the cluster-robust
sandwich; panel_fit’s denominator is the default
analytical-Hessian SE, so the ratio also reflects that choice of
estimator, not only the likelihood.) Because it targets a different
(composite) likelihood, the cross-sectional fit’s log-likelihood, AIC
and BIC are not comparable to the panel fit’s.
Individual-level tastes.
conditional_tastes() computes, for each decision maker, the
mean and SD of \(\beta_n\) conditional
on the choices they made (Revelt and Train 2000; Train 2009, ch. 11) —
the frequentist analogue of the hierarchical Bayesian
beta_i summaries of run_hmnlogit():
ct <- conditional_tastes(panel_fit)
ct
#> Conditional tastes E[beta | observed choices] for 300 decision makers (S = 100 draws each)
#> Unit-level means and SDs in $mean and $sd.
#> Population distribution vs conditional means (revealed + residual should be close to 1):
#> model_mean model_sd mean_cond sd_cond revealed residual
#> w1 0 1.114 0.064 0.767 0.473 0.497
#> w2 0 1.301 0.143 0.952 0.535 0.438Read this carefully: ct$mean[, n] is not
decision maker \(n\)’s own \(\beta_n\). It is the average taste among
people who would make \(n\)’s choices,
so it shrinks toward the population mean — more so for people with fewer
situations — and its spread across people understates the population
spread in \(\Sigma\). The printed
revealed and residual columns decompose
population variance by the law of total variance, \(\mathrm{Var}(\beta) = \mathrm{Var}(E[\beta \mid
y]) + E[\mathrm{Var}(\beta \mid y)]\); they should sum to about
one. This is a numerical check, not a specification test: the
random-coefficient covariance \(\Sigma\) is always freely estimated, so
this identity is the likelihood’s first-order condition for \(\Sigma\) and holds at a converged fit
whatever the true taste distribution — a departure points to too few
draws S or an unconverged fit. (Here
rc_mean = FALSE fixes the random-coefficient means at zero
rather than estimating them freely; a mean_cond far from
zero below would flag that restriction.) In this example the conditional
means correlate 0.67 (w1) and 0.74 (w2) with
the realized tastes in panel_sim$true_params$gamma_i. At
the true \(\theta\), shrinkage sets
this correlation exactly: \(\mathrm{cor}(E[\beta \mid y], \beta) = sd(E[\beta
\mid y]) / sd(\beta) =
\sqrt{\texttt{revealed}}\), so the same shrinkage that keeps
revealed below one keeps the correlation below one too. A
single choice reveals little about a person’s tastes, so fitting without
person_col leaves revealed much smaller — 0.17
vs. 0.47 for w1, 0.2 vs. 0.54 for w2 — than
with repeated choices; conditional means only earn their keep with
repeated choices.
Using conditional means as regressors in a second-stage analysis is
worse than merely ignoring their sampling and simulation error: the draw
weights behind each conditional mean are not exported, so regressing on
ct$mean directly measures tastes with error and biases the
second-stage coefficient toward zero (attenuation), on top of
understating its standard error. Prefer entering the observables in
utility instead — interactions with the random-coefficient columns, or
rc_mean shifts — or, if a second-stage regression is
unavoidable, propagate the uncertainty with a decision-maker-level
bootstrap rather than treating conditional means as if they were
data.
For a reportable run_mxlogit() application:
theta_init values. The default start moves with
the covariates’ units (a normal random coefficient’s standard deviation
starts at one over its covariate’s typical step: 1 for a dummy, the
standard deviation for a dense covariate), so the optimizer’s path
barely depends on them (exactly so, in exact arithmetic, under
scale_vars = "sd"); a shifted log-normal coefficient’s
random part, though, has median one in its covariate’s units: put such a
covariate in units where the coefficient is of order one.S until coefficients, covariance elements,
standard errors, WTP, elasticities and diversion ratios are stable—not
merely until the optimizer reports convergence.person_col
when tastes can be defended as constant within a decision maker’s
situations, and keep the clustered cross-sectional fit when they cannot
— it relies only on each situation’s marginal and pays in precision,
mostly on \(\Sigma\). Add
run_hmnlogit() when partially-pooled alternative effects or
posterior uncertainty are wanted too.Bounded/censored mixing distributions, WTP-space estimation, and
latent-class logit can be useful alternative modeling strategies, but
they are not implemented in choicer v0.2.0. They
require another implementation today and are roadmap candidates here;
they should not be presented as options to run_mxlogit().
The broader tradeoff is laid out in Choosing among choice
models.
McFadden, D. and Train, K. (2000). Mixed MNL models for discrete response. Journal of Applied Econometrics, 15(5), 447-470.
Revelt, D. and Train, K. (1998). Mixed logit with repeated choices: households’ choices of appliance efficiency level. Review of Economics and Statistics, 80(4), 647-657.
Revelt, D. and Train, K. (2000). Customer-specific taste parameters and mixed logit: households’ choice of electricity supplier. Working Paper E00-274, Department of Economics, University of California, Berkeley.
Train, K. E. (2009). Discrete Choice Methods with Simulation (2nd ed.). Cambridge University Press.
Train, K. and Weeks, M. (2005). Discrete choice models in preference space and willingness-to-pay space. In R. Scarpa and A. Alberini (Eds.), Applications of Simulation Methods in Environmental and Resource Economics. Springer.