Which standard errors, and when • choicer

A demand estimate is only as reportable as its standard error. The default maximum-likelihood variance — the inverse Hessian — is correct when the model is correctly specified and choice situations are independent draws. Applied work routinely violates one or both: the same decision maker contributes several choice situations (panel and survey data), the sample was drawn by outcome (choice-based sampling), or the likelihood is at best a useful approximation. choicer treats the variance estimator as a choice you make — and can revise after estimation, without refitting.

Every MNL, MXL and NL fit exposes four estimators through one interface:

`type`	Estimator	When it is the right choice
`"hessian"`	Inverse analytical Hessian	Correct specification, independent choice situations. The classical default.
`"bhhh"`	Inverse outer product of scores (BHHH/OPG)	Same asymptotic justification as the Hessian; cheap; a common companion for simulated likelihoods.
`"robust"`	Huber-White sandwich $A^{-1} B A^{-1}$	Quasi-ML under misspecification; also the valid variance under WESML / choice-based weighting (the meat is the weight-squared score outer product).
`"cluster"`	Sandwich over within-cluster sums of scores	Dependence across choice situations from the same decision maker, market, or sampling unit.

Here the unit $i$ is the choice situation — or, for a mixed logit fitted with person_col (see below), the decision maker: several situations already collapse to one scored unit, so there is nothing left within a person to cluster on. Otherwise, if the same person supplies several situations, cluster on the person. No small-sample correction is applied—not a degrees-of-freedom multiplier, leverage correction, or few-cluster reference distribution—so the number of clusters should be reasonably large. With few clusters, the asymptotic cluster sandwich can be seriously anti-conservative; choicer does not currently implement CR2, a wild-cluster bootstrap, or randomization inference. Those designs require an external small-sample procedure or a research design with a defensible higher-level sampling argument.

library(choicer)
set_num_threads(2)

Dependence you should expect: panel data

The cleanest way to see what clustering buys is to build data with genuine within-person dependence. simulate_hmnl_data() draws a panel in which each person has their own taste vector $\beta_i \sim N(b, W)$ and makes several choices with it — a stylized version of the persistent within-person dependence in a survey or scanner panel.

sim <- simulate_hmnl_data(N = 500, T = 8, J = 4, seed = 3)

We deliberately fit a plain multinomial logit. Because a homogeneous logit is misspecified for this random-coefficients DGP, its point estimates target the pseudo-true homogeneous-logit projection, not generally the population mean $b$ . Its likelihood also treats all 500 * 8 = 4000 choice situations as independent — which they are not, because each person’s taste deviation persists across their eight tasks. Passing cluster_col supplies the person labels at fit time and selects cluster-robust standard errors for summary():

fit <- run_mnlogit(
  data                   = sim$data,
  id_col                 = "task",
  alt_col                = "alt",
  choice_col             = "choice",
  covariate_cols         = c("x1", "x2"),
  cluster_col            = "pid",
  include_outside_option = TRUE
)
#> Optimization run time 0h:0m:0.05s
summary(fit)
#> Multinomial Logit (MNL) model
#> 
#> Parameter    Estimate  Std.Error  z-value  Pr(>|z|)  
#> x1           0.740484   0.044384  16.6835  0.00e+00  ***
#> x2          -0.537569   0.044763 -12.0092  0.00e+00  ***
#> ASC_1        0.943559   0.054241  17.3956  0.00e+00  ***
#> ASC_2       -0.198141   0.068933  -2.8744  4.05e-03  **
#> ASC_3        0.766699   0.054849  13.9783  0.00e+00  ***
#> ASC_4        0.708742   0.055895  12.6798  0.00e+00  ***
#> ---
#> Signif. codes:  '***' 0.001 '**' 0.01 '*' 0.05
#> 
#> Std. Errors: Cluster-robust sandwich 
#> Log-likelihood: -5748.24 
#> AIC: 11508.5  | BIC: 11546.2 
#> McFadden R2: 0.107 (adj: 0.106) | Hit rate: 0.400 
#> N: 4000  | Parameters: 6 
#> Optimization time: 0.05 s
#> Convergence: 1 ( NLOPT_SUCCESS: Generic success return value. )

How much did it matter? vcov(fit, type = ) recomputes any of the four estimators from the data stored with the fit:

se_of <- function(V) sqrt(diag(V))
tab <- cbind(
  hessian = se_of(vcov(fit, type = "hessian")),
  bhhh    = se_of(vcov(fit, type = "bhhh")),
  robust  = se_of(vcov(fit, type = "robust")),
  cluster = se_of(vcov(fit, type = "cluster"))
)
round(tab, 4)
#>       hessian   bhhh robust cluster
#> x1     0.0356 0.0349 0.0364  0.0444
#> x2     0.0352 0.0342 0.0363  0.0448
#> ASC_1  0.0565 0.0566 0.0564  0.0542
#> ASC_2  0.0701 0.0701 0.0700  0.0689
#> ASC_3  0.0579 0.0579 0.0579  0.0548
#> ASC_4  0.0584 0.0583 0.0585  0.0559
round(tab[, "cluster"] / tab[, "hessian"], 2)
#>    x1    x2 ASC_1 ASC_2 ASC_3 ASC_4 
#>  1.25  1.27  0.96  0.98  0.95  0.96

The clustered standard errors on the slope coefficients are about 25% wider than the Hessian ones; those on the alternative-specific constants differ little. That pattern is the fingerprint of the data-generating process: the persistent component is each person’s taste deviation, so it is the slope scores that are correlated within person. Note also that "robust" alone recovers almost none of this — the Huber-White sandwich guards against misspecification of each situation’s likelihood, not against dependence across situations. If the same decision maker appears repeatedly, robust-without-clustering is the wrong comfort (the exception is a mixed logit fitted with person_col: there the decision maker is the likelihood unit, so "robust" is already clustered by construction — see below).

Post hoc, without refitting

Nothing above required deciding at fit time. Any fit stored with keep_data = TRUE (the default) can be re-inferenced later — useful when the clustering level is itself a robustness question:

fit0 <- run_mnlogit(
  data                   = sim$data,
  id_col                 = "task",
  alt_col                = "alt",
  choice_col             = "choice",
  covariate_cols         = c("x1", "x2"),
  include_outside_option = TRUE
)
#> Optimization run time 0h:0m:0.05s

# one cluster label per choice situation, named by choice-situation id
task_person <- unique(sim$data[, c("task", "pid")])
cl <- setNames(task_person$pid, task_person$task)
se_of(vcov(fit0, type = "cluster", cluster = cl))
#>      x1      x2   ASC_1   ASC_2   ASC_3   ASC_4 
#> 0.04438 0.04476 0.05424 0.06893 0.05485 0.05590

The cluster vector has one entry per choice situation — here, one pid per task — and is named by the choice-situation id. choicer uses those names to realign labels safely to its prepared order.

Clustering repairs inference, not efficiency

A caveat that matters for the mixed logit. Without person_col, run_mxlogit() maximizes a cross-sectional simulated likelihood: each choice situation’s probability is integrated over the mixing distribution separately, rather than holding one simulation draw fixed across a person’s tasks and integrating their joint probability. On panel data this is a composite (pseudo-) likelihood: every situation’s marginal probability is still correctly specified, so the point estimates stay consistent for the population mixing distribution, but the fit is statistically inefficient — it never uses the fact that the same tastes generated a person’s repeated choices, which is exactly the within-person variation that disciplines $\Sigma$ . type = "cluster" (or cluster_col at fit time) makes the standard errors valid under that within-person dependence, but leaves the point estimates, and their inefficiency, unchanged: clustering repairs the inference, not the likelihood. And because it is a different likelihood, its log-likelihood, AIC and BIC are not comparable to a panel fit’s.

When the panel structure is the object of interest — a taste distribution disciplined by repeated choices, not just valid standard errors around a composite fit — estimate the panel likelihood instead. Two options hold one taste draw fixed within each person (Revelt and Train 1998): run_mxlogit(..., person_col = ), covered in the mixed logit vignette, whose scores are then already per decision maker, so type = "robust" is clustered by construction; or the hierarchical Bayesian logit, run_hmnlogit(..., person_col = ), covered in the hierarchical Bayes vignette, which adds partial pooling of alternative-level effects and posterior rather than asymptotic uncertainty.

Choice-based sampling

When the sample was drawn by outcome (interviewing travellers at the terminal of the mode they chose, oversampling rare hospitals), the weights are not optional and the variance choice is forced: under WESML weighting the inverse weighted Hessian and the ordinary BHHH variance are both invalid, and type = "robust" — whose meat carries the squared weights — is the correct estimator, conditional on the supplied population shares Q. If weighted observations are also dependent within sampling units, use the weighted cluster sandwich instead; if Q is estimated, its uncertainty requires an additional resampling or sensitivity layer. The choice-based sampling vignette develops the full workflow, from population shares to weighted fit.

References

Berndt, E., Hall, B., Hall, R. and Hausman, J. (1974). Estimation and inference in nonlinear structural models. Annals of Economic and Social Measurement, 3(4), 653-665.

Cameron, A. C. and Miller, D. L. (2015). A practitioner’s guide to cluster-robust inference. Journal of Human Resources, 50(2), 317-372.

Manski, C. F. and Lerman, S. R. (1977). The estimation of choice probabilities from choice based samples. Econometrica, 45(8), 1977-1988.

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.

White, H. (1982). Maximum likelihood estimation of misspecified models. Econometrica, 50(1), 1-25.