Estimates a multinomial logit model via maximum likelihood.
Usage
run_mnlogit(
data = NULL,
id_col = NULL,
alt_col = NULL,
choice_col = NULL,
covariate_cols = NULL,
input_data = NULL,
optimizer = NULL,
control = list(),
weights = NULL,
weights_col = NULL,
outside_opt_label = NULL,
include_outside_option = FALSE,
use_asc = TRUE,
keep_data = TRUE,
scale_vars = c("bhhh", "none", "sd", "mad", "iqr"),
se_method = c("hessian", "bhhh", "sandwich", "cluster"),
cluster_col = NULL,
nloptr_opts = NULL
)Arguments
- data
Data frame containing choice data (convenience workflow). Mutually exclusive with
input_data.- id_col
Name of the column identifying choice situations (individuals).
- alt_col
Name of the column identifying alternatives.
- choice_col
Name of the column indicating chosen alternative (1 = chosen, 0 = not chosen).
- covariate_cols
Vector of names of columns to be used as covariates. With
use_asc = TRUE, none may be named like a constant,ASC_<label>.- input_data
List output from
prepare_mnl_data(advanced workflow). Mutually exclusive withdata.- optimizer
Optimizer to use:
"nloptr"(default),"optim", or a custom function with signaturef(theta_init, eval_f, lower, upper, control)whereeval_f(theta)returnslist(objective, gradient). Must return a list withpar/value(orsolution/objective). If the custom function acceptscontrolor..., thecontrolargument is forwarded; otherwise it is silently ignored. Withscale_varsother than"none", it receivestheta_init,lower,upperandeval_fin the scaled coordinates, and itsparis mapped back.- control
List of optimizer-specific control parameters passed to the chosen optimizer (e.g.,
list(maxeval = 2000)for nloptr). Withscale_varsother than"none", settings in parameter units apply to the scaled coordinates.- weights
Optional vector of weights for each choice situation. If
NULL, equal weights are used. All weights must be finite and strictly positive.- weights_col
Optional name of a column in
dataholding per-row weights (convenience workflow only). The column must be constant within eachid_col(one weight per choice situation) and is collapsed accordingly. Mutually exclusive withweights. All weights must be finite and strictly positive. Used for choice-based / WESML weighting; pair withse_method = "sandwich"for valid inference.- outside_opt_label
Label for the outside option (if any). If
NULL, no outside option is assumed.- include_outside_option
Logical indicating whether to include an outside option in the model.
- use_asc
Logical indicating whether to include alternative-specific constants (ASCs) in the model.
- keep_data
Logical. If
TRUE(default), stores prepared data in the returned object forpredict()and post-estimation functions.- scale_vars
How the optimizer's coordinates are scaled. The choice does not change the estimator; it changes how fast, and whether, the optimizer reaches the maximum.
"bhhh"(default) scales each parameter by the inverse square root of the diagonal of the BHHH (outer product of gradients) matrix at the start values (all zeros), rounded to a power of two, from one extra gradient pass, so that the BHHH diagonal, an estimate of the curvature, is near one at the start values in every coordinate, the constants' included (theirs follow the alternatives' shares). Where the diagonal is unusable (no finite, positive entry), the parameters are left unscaled, with a message."none"hands the optimizer the parameters as they are (the default in choicer 0.2.x). With"sd"(sample standard deviation),"mad"(stats::mad) or"iqr"(stats::IQR(x) / 1.349), the optimizer works on the coefficientsXwould have with every column divided by the chosen scale (each coefficient times its column's scale); the constants keep their scale, the data are not divided, and the fit holds no scaled copy of them. With any choice but"none"(so by default), the optimizer, a custom one included, works in the scaled coordinates, so settings in parameter units incontrol(nloptr'sxtol_abs, optim'sparscale) apply to the scaled coordinates;"none"keeps the parameters' units. Coefficients are reported in natural units and the standard errors computed in them (for any choice but"none", with the information matrix equilibrated before it is inverted), so, where the likelihood is well identified and every choice reaches its maximum, reported quantities do not depend on this choice beyond the optimizer's tolerance.- se_method
Method for computing standard errors:
"hessian"(default, analytical Hessian),"bhhh"(outer product of gradients),"sandwich"(robust Huber–White / WESML variance \(A^{-1} B A^{-1}\)), or"cluster"(cluster-robust sandwich; requirescluster_color a preparedinput_datawith aclusterfield). Use"sandwich"under choice-based / WESML weighting. Any of these can also be recomputed post hoc viavcov(fit, type = ).- cluster_col
Optional name of a column in
dataholding cluster labels for cluster-robust standard errors (e.g. a person id when the same decision maker contributes several choice situations). Must be constant within eachid_col. Supplyingcluster_colwithout an explicitse_methodselectsse_method = "cluster".- nloptr_opts
Deprecated. Use
optimizerandcontrolinstead.
Value
A choicer_mnl object (inherits from choicer_fit).
Standard S3 methods available: summary(), coef(), vcov(),
logLik(), AIC(), BIC(), nobs(), predict().
Details
Two workflows are supported:
- Convenience (default)
Supply
dataand column names. Data preparation (prepare_mnl_data) is handled automatically.- Advanced
Call
prepare_mnl_datayourself and pass the result viainput_data.
Column names
No column the model uses (the id, alternative, choice, covariate, weight,
cluster, nest or decision-maker column) may have a name starting with
.choicer_, which choicer reserves for its working columns, and the
alternative column may not be named alt_int, N_OBS, N_CHOICES,
TAKE_RATE or MKT_SHARE, the names of the other columns of
alt_mapping. Either is an error that names the column. Columns the model
does not use may have any name.
Examples
# \donttest{
library(data.table)
set.seed(42)
N <- 100; J <- 3; beta_true <- c(1.0, -0.5)
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#> id alt x1 x2
#> <int> <int> <num> <num>
#> 1: 1 1 1.37095845 -0.004620768
#> 2: 1 2 -0.56469817 0.760242168
#> 3: 1 3 0.36312841 0.038990913
#> 4: 2 1 0.63286260 0.735072142
#> 5: 2 2 0.40426832 -0.146472627
#> ---
#> 296: 99 2 -0.47733551 0.160327395
#> 297: 99 3 -0.16626149 -0.433641942
#> 298: 100 1 0.86256338 1.537412419
#> 299: 100 2 0.09734049 -2.170246577
#> 300: 100 3 -1.62561674 1.027004619
dt[, V := drop(as.matrix(.SD) %*% beta_true), .SDcols = c("x1","x2")]
#> id alt x1 x2 V
#> <int> <int> <num> <num> <num>
#> 1: 1 1 1.37095845 -0.004620768 1.37326883
#> 2: 1 2 -0.56469817 0.760242168 -0.94481926
#> 3: 1 3 0.36312841 0.038990913 0.34363295
#> 4: 2 1 0.63286260 0.735072142 0.26532653
#> 5: 2 2 0.40426832 -0.146472627 0.47750464
#> ---
#> 296: 99 2 -0.47733551 0.160327395 -0.55749920
#> 297: 99 3 -0.16626149 -0.433641942 0.05055948
#> 298: 100 1 0.86256338 1.537412419 0.09385717
#> 299: 100 2 0.09734049 -2.170246577 1.18246377
#> 300: 100 3 -1.62561674 1.027004619 -2.13911905
dt[, prob := exp(V) / sum(exp(V)), by = id]
#> id alt x1 x2 V prob
#> <int> <int> <num> <num> <num> <num>
#> 1: 1 1 1.37095845 -0.004620768 1.37326883 0.68700257
#> 2: 1 2 -0.56469817 0.760242168 -0.94481926 0.06764341
#> 3: 1 3 0.36312841 0.038990913 0.34363295 0.24535402
#> 4: 2 1 0.63286260 0.735072142 0.26532653 0.33940231
#> 5: 2 2 0.40426832 -0.146472627 0.47750464 0.41962617
#> ---
#> 296: 99 2 -0.47733551 0.160327395 -0.55749920 0.21814327
#> 297: 99 3 -0.16626149 -0.433641942 0.05055948 0.40069908
#> 298: 100 1 0.86256338 1.537412419 0.09385717 0.24525785
#> 299: 100 2 0.09734049 -2.170246577 1.18246377 0.72844833
#> 300: 100 3 -1.62561674 1.027004619 -2.13911905 0.02629382
dt[, choice := as.integer(alt == sample(alt, 1, prob = prob)), by = id]
#> id alt x1 x2 V prob choice
#> <int> <int> <num> <num> <num> <num> <int>
#> 1: 1 1 1.37095845 -0.004620768 1.37326883 0.68700257 1
#> 2: 1 2 -0.56469817 0.760242168 -0.94481926 0.06764341 0
#> 3: 1 3 0.36312841 0.038990913 0.34363295 0.24535402 0
#> 4: 2 1 0.63286260 0.735072142 0.26532653 0.33940231 1
#> 5: 2 2 0.40426832 -0.146472627 0.47750464 0.41962617 0
#> ---
#> 296: 99 2 -0.47733551 0.160327395 -0.55749920 0.21814327 0
#> 297: 99 3 -0.16626149 -0.433641942 0.05055948 0.40069908 0
#> 298: 100 1 0.86256338 1.537412419 0.09385717 0.24525785 0
#> 299: 100 2 0.09734049 -2.170246577 1.18246377 0.72844833 1
#> 300: 100 3 -1.62561674 1.027004619 -2.13911905 0.02629382 0
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
#> Optimization run time 0h:0m:0.01s
summary(fit)
#> Multinomial Logit (MNL) model
#>
#> Parameter Estimate Std.Error z-value Pr(>|z|)
#> x1 1.091325 0.192975 5.6553 1.56e-08 ***
#> x2 -0.511225 0.152893 -3.3437 8.27e-04 ***
#> ASC_2 -0.061376 0.288667 -0.2126 8.32e-01
#> ASC_3 -0.074843 0.289651 -0.2584 7.96e-01
#> ---
#> Signif. codes: '***' 0.001 '**' 0.01 '*' 0.05
#>
#> Std. Errors: Analytical Hessian
#> Log-likelihood: -81.2408
#> AIC: 170.482 | BIC: 180.902
#> McFadden R2: 0.261 (adj: 0.224) | Hit rate: 0.590
#> N: 100 | Parameters: 4
#> Optimization time: 0.01 s
#> Convergence: 1 ( NLOPT_SUCCESS: Generic success return value. )
coef(fit)
#> x1 x2 ASC_2 ASC_3
#> 1.09132487 -0.51122457 -0.06137617 -0.07484335
AIC(fit)
#> [1] 170.4816
predict(fit, type = "shares")
#> [,1]
#> [1,] 0.34
#> [2,] 0.32
#> [3,] 0.34
# }