Runs mixed logit estimation

Estimates a mixed logit model via simulated maximum likelihood.

Usage

run_mxlogit(
  data = NULL,
  id_col = NULL,
  alt_col = NULL,
  choice_col = NULL,
  covariate_cols = NULL,
  random_var_cols = NULL,
  input_data = NULL,
  eta_draws = NULL,
  S = 100L,
  rc_dist = NULL,
  rc_mean = FALSE,
  rc_correlation = FALSE,
  use_asc = TRUE,
  theta_init = NULL,
  lower = NULL,
  upper = NULL,
  optimizer = NULL,
  control = list(),
  se_method = c("hessian", "bhhh", "sandwich", "cluster"),
  scale_vars = c("bhhh", "none", "sd", "mad", "iqr"),
  weights = NULL,
  outside_opt_label = NULL,
  include_outside_option = FALSE,
  draws = c("store", "generate"),
  seed = NULL,
  scramble = c("permuted", "none", "owen"),
  keep_data = TRUE,
  nloptr_opts = NULL,
  weights_col = NULL,
  cluster_col = NULL,
  person_col = NULL
)

Arguments

data

Data frame containing choice data (convenience workflow). Mutually exclusive with input_data.

id_col

Name of the column identifying choice situations.

alt_col

Name of the column identifying alternatives.

choice_col

Name of the column indicating chosen alternative (1/0).

covariate_cols

Vector of column names for fixed covariates. None may take a name the model gives its parameters, or prints for them in summary(): ASC_<label>, L_<i><j> and Sigma_<i><j>, and with rc_mean = TRUE Mu_<variable> and, for a log-normal coefficient, exp(Mu_<variable>).

random_var_cols

Vector of column names for random coefficients.

input_data

List output from prepare_mxl_data (advanced workflow). Mutually exclusive with data.

eta_draws

Array of shape K_w x S x U with standard normal draws, one K_w x S block per likelihood unit: U is the number of decision makers when input_data was prepared with person_col (length(input_data$Ti)), and the number of choice situations otherwise. Required for the advanced workflow; auto-generated from S in the convenience workflow. Post-hoc methods (vcov(fit, type = ), wesml_vcov, conditional_tastes, prediction) regenerate the Halton draws get_halton_normals builds from the recorded draw count, so build eta_draws with it for them to reproduce the estimation draws.

S

Integer number of Halton draws per decision maker, or per choice situation without person_col (convenience workflow only). Default 100.

rc_dist

Integer vector indicating distribution of random coefficients (0 = normal, 1 = log-normal). Default: all normal.

rc_mean

Logical indicating whether to estimate means for random coefficients.

rc_correlation

Logical indicating whether random coefficients are correlated (convenience workflow). Ignored when input_data is used (taken from the prepared data).

use_asc

Logical indicating whether to include alternative-specific constants.

theta_init

Initial parameter vector in natural-scale units. If NULL (default), the start is in the columns' units: zero for the \(\beta\), \(\mu\) and ASC blocks and the Cholesky factor's off-diagonal entries; on its diagonal, \(L_{pp} = 1 / h_p\) for a normal random coefficient, where \(h_p\) is the typical step of its column of \(W\): the distance between its values for a column with two (a dummy: 1), the mean distance from a value that fills more than half of the column over the other rows (a mostly-zero column), and the standard deviation otherwise; a column with no variation starts at \(L_{pp} = 1\). At the start, a step of \(h_p\) in the column then moves the utility by a random amount with standard deviation 1, whatever the column's units of measurement. A log-normal coefficient starts at \(L_{pp} = 0.5\), the standard deviation of the logarithm of its random part, which has no units (that random part has median 1 in the column's units, so a log-normal covariate belongs in units where its coefficient is of order one). Where the default start falls outside lower/upper, it is moved to the nearest bound, with a message. A supplied theta_init is used as given.

lower, upper

Optional parameter bounds for the optimizer, in natural-scale units (mapped to the optimizer's coordinates when scale_vars != "none"); lower must not exceed upper, nor be Inf (nor upper -Inf). Each accepts three forms:

NULL

(default) Unbounded (-Inf/Inf).

Unnamed numeric vector of length n_params

Full-length vector ordered exactly like theta_init (the nloptr-native form).

Named numeric vector

Names must be a subset of the parameter names (\(\beta\) block: column names of X; \(\mu\) block: Mu_<col> (if rc_mean = TRUE); Cholesky block: L_<i><j> for \(i \ge j\); ASC block: ASC_<level>). Unlisted parameters default to \(\pm\infty\). This is the recommended form for typical use, e.g. lower = c(L_11 = -5, L_22 = -5) to clip Cholesky diagonals.

optimizer

Optimizer to use: "nloptr" (default), "optim", or a custom function. See run_mnlogit for details.

control

List of optimizer-specific control parameters. With scale_vars other than "none", settings in parameter units apply to the scaled coordinates.

se_method

Method for computing standard errors. One of "hessian" (default) for the analytical Hessian of the simulated log-likelihood, "bhhh" for the BHHH/outer-product-of-gradients (OPG) estimator, "sandwich" for the robust (Huber-White) variance \(V = A^{-1} B A^{-1}\) (bread \(A\) = weighted negated Hessian, meat \(B\) = weight-squared OPG), or "cluster" for the cluster-robust sandwich (requires cluster_col or a prepared input_data with a cluster field). Use "sandwich" for valid inference under choice-based / WESML weighting, where the inverse-Hessian and ordinary BHHH are invalid; it reduces to the usual robust variance under uniform weights. BHHH needs no second derivatives and so costs less than the Hessian. For both, the work per decision maker involves only the parameters that enter that decision maker's likelihood (the fixed coefficients, the random coefficients' means and Cholesky factor, and the constants of the alternatives they faced), so it does not grow with the number of alternatives in the data; the parameters-by-parameters matrix and its inversion do. BHHH may underestimate standard errors in finite samples or away from the optimum. Any of these variances can also be recomputed post hoc via vcov(fit, type = ). Without person_col, clustering repairs the inference, not the likelihood: the cross-sectional likelihood still treats each choice situation as an independent draw from the mixing distribution, and on panel data is a less efficient composite likelihood for the same taste distribution. With person_col the likelihood unit is the decision maker, so scores are per decision maker and "sandwich" is already robust to dependence across a decision maker's choice situations; "cluster" is then needed only for coarser groups (e.g. households or markets) that nest decision makers.

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, 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: it reaches the constants and the Cholesky factor, whose curvatures can differ by orders of magnitude that no column scale reaches (the constants of large and small alternatives, say). A start near the estimates (theta_init) gives it a better diagonal. Where the diagonal is unusable (no finite, positive entry, as when the likelihood overflows at the start values), 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, i.e. 1.4826 \(\times\) median absolute deviation; SD-equivalent under normality) or "iqr" (stats::IQR(x) / 1.349; also SD-equivalent under normality), the optimizer works on the parameters the design would have with every column of X and W divided by the chosen scale; the constants keep their scale, the data are not divided, and the fit holds no scaled copy of them. Robust scales ("mad"/"iqr") better capture the bulk for heavy-tailed columns where SD is dominated by outliers, but stats::mad can return zero when more than half of a column's entries are identical (e.g., a sparse 0/1 dummy) and will then trigger the same near-constant-column error as "sd". The parameters of log-normal random coefficients (rc_dist == 1) keep their scale under these three, since the shifted log-normal parameterization does not admit a closed-form rescaling under multiplicative column scaling. With any choice but "none" (so by default), the optimizer, a custom one included, works in the scaled coordinates (theta_init, lower and upper are given in natural units and mapped), so settings in parameter units in control (nloptr's xtol_abs; optim's parscale, and pgtol, which then applies to the scaled gradient) 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.

weights

Optional weight vector (convenience workflow), one weight per choice situation in ascending-id order (see prepare_mxl_data). With person_col, situations are reordered by decision maker, so a positional vector is accepted only when that leaves the id order unchanged; otherwise supply weights_col. If NULL, equal weights are used. All weights must be finite and strictly positive. With person_col they are decision-maker weights, constant within each decision maker: the objective is \(\sum_n w_n \log L_n\).

outside_opt_label

Label for the outside option (convenience workflow).

include_outside_option

Logical whether to include an outside option (convenience workflow).

draws

Draw storage mode. One of "store" (default) or "generate". "store" pre-materializes the full \(K_w \times S \times U\) Halton cube (\(8 K_w S U\) bytes), one block per likelihood unit (U decision makers with person_col, choice situations otherwise; existing behavior, exact reproducibility). It supports at most \(2^{31} - 1\) points (\(S \times U\)). predict(), logsum(), consumer_surplus(), elasticities(), diversion_ratios(), blp() and gof() (hence summary()'s fit statistics) integrate over one block per choice situation instead (see get_halton_normals). With up to 128 random coefficients they form each block on the fly, in parallel, from the points get_halton_normals(S, N, K_w) would hold for the \(N\) choice situations, mapped by R's qnorm(), so they never hold that cube and have no limit on the number of points; with more, they build the cube. Their draws, and so their results, are the cube's, bit for bit where randtoolbox and choicer are compiled with the same floating-point contraction, as on CRAN's platforms (otherwise some of the underlying uniforms differ in their last bit, which moves those draws by at most a few parts in \(10^9\), in the far tails); sums over choice situations are bit for bit at one thread, and at more threads vary in their last bits with the threads' order, as they always did. blp() may keep the draws it forms across its iterations, within its keep_draws_bytes (4 GiB by default), with either storage mode. In the convenience workflow, run_mxlogit() warns before it builds an estimation cube above 1 GiB: the cube is held while the model is fitted, and vcov(type = ), wesml_vcov and conditional_tastes rebuild it whole. "generate" computes each unit's draws on-the-fly in C++ from a stored seed, eliminating the O(U) cube; recommended for memory-constrained or large-N settings. With scramble = "permuted", each base-\(b\) digit position in each dimension receives a seeded permutation shared across sequence indices. This is not Owen's nested-uniform scramble and does not carry standard randomized-QMC unbiasedness or error-estimation guarantees. Only supported in the convenience workflow.

seed

Integer master seed for the on-the-fly generator. Used only when draws = "generate". If NULL (default), a seed is drawn from R's RNG at call time (so set.seed() governs reproducibility). Ignored when draws = "store".

scramble

Scrambling mode for on-the-fly Halton draws. One of "permuted" (default) for seeded position-wise digit permutations or "none" for plain Halton (identity permutations). The historical value "owen" is accepted with a deprecation warning as an alias for "permuted"; the implementation is not Owen's nested-uniform scramble. "none" uses the points of the randtoolbox sequence, mapped to normals by choicer's own inverse normal CDF, which is within about \(10^{-14}\) of qnorm(): its draws agree with store mode's to that precision, not bit for bit. Simulation-draw sensitivity should be assessed by increasing S and, for "permuted", varying seed. Used only when draws = "generate".

keep_data

Logical. If TRUE (default), stores prepared data in the returned object for post-estimation functions.

nloptr_opts

Deprecated. Use optimizer and control instead.

weights_col

Optional name of a column in data holding a per-row weight (constant within each choice situation, finite and strictly positive). Mutually exclusive with weights; the recommended way to pass WESML weights from sample_by_choice / wesml_weights, since alignment is by id rather than by position. Convenience workflow only. If data carries choice-based-sampling provenance (a "choice_sampling" attribute, as attached by sample_by_choice / wesml_weights) and neither weights nor weights_col is supplied, the recorded weight column is auto-detected and applied (with a message); if that column is absent the call errors rather than silently fitting an unweighted model under a WESML label. With person_col the weight must also be constant within each decision maker, and WESML is not supported: choice-based weights vary with the chosen alternative across a decision maker's situations, so they cannot weight the panel likelihood.

cluster_col

Optional name of a column in data holding cluster labels for cluster-robust standard errors (e.g. a person id when the same decision maker contributes several choice situations to a cross-sectional fit). Must be constant within each id_col, and within each decision maker when person_col is used (clusters must nest decision makers). Supplying cluster_col without an explicit se_method selects se_method = "cluster".

person_col

Optional name of the column in data identifying decision makers (respondents). When supplied, all choice situations of a decision maker share one draw of the random coefficients (panel likelihood); id_col must still identify choice situations uniquely across the data set. NULL (default) makes each choice situation its own decision maker (cross-sectional likelihood). Convenience workflow only; in the advanced workflow pass person_col to prepare_mxl_data.

Value

A choicer_mxl object (inherits from choicer_fit). Standard S3 methods available: summary(), coef(), vcov(), logLik(), AIC(), BIC(), nobs(). A panel fit (person_col) also carries n_persons, the number of decision makers; nobs() remains the number of choice situations, which is also the \(n\) in BIC's \(\log n\) penalty.

Details

Two workflows are supported:

Convenience

Supply data and column names. Data preparation (prepare_mxl_data) and Halton draw generation (get_halton_normals) are handled automatically.

Advanced

Call prepare_mxl_data and get_halton_normals yourself, then pass the results via input_data and eta_draws.

Cross-section vs panel. With person_col, each decision maker \(n\) draws one taste vector \(\beta_n\) from the mixing distribution \(f(\beta \mid \theta)\) and keeps it in all of their choice situations \(t = 1, \dots, T_n\). The likelihood integrates the product of their conditional logit probabilities (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,$$ simulated with S draws per decision maker, so repeated choices speak to each decision maker's tastes, and conditional_tastes recovers what they reveal. Without person_col (the default), each choice situation is an independent draw from the mixing distribution: \(L_i(\theta) = \int P_i(j_i \mid \beta) f(\beta \mid \theta) \, d\beta\). That is the right model when tastes are redrawn in every situation (or each decision maker chooses once). On panel data with stable tastes it is a composite (pseudo-) likelihood: it still targets the population mixing distribution, but it ignores the within-person correlation of choices, so it is less efficient, its standard errors must be clustered by decision maker (cluster_col), and its log-likelihood, AIC and BIC are not comparable with those of the panel fit.

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.

References

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.

See also

conditional_tastes, run_hmnlogit (the hierarchical Bayes counterpart)

Examples

# \donttest{
library(data.table)
set.seed(42)
N <- 100; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N), w2 = rnorm(.N))]
#>         id   alt          x1           w1         w2
#>      <int> <int>       <num>        <num>      <num>
#>   1:     1     1  1.37095845 -0.004620768 -0.2484829
#>   2:     1     2 -0.56469817  0.760242168  0.4223204
#>   3:     1     3  0.36312841  0.038990913  0.9876533
#>   4:     2     1  0.63286260  0.735072142  0.8355682
#>   5:     2     2  0.40426832 -0.146472627 -0.6605219
#>  ---                                                
#> 296:    99     2 -0.47733551  0.160327395  0.1704735
#> 297:    99     3 -0.16626149 -0.433641942  1.2006682
#> 298:   100     1  0.86256338  1.537412419 -0.1634059
#> 299:   100     2  0.09734049 -2.170246577  1.2824759
#> 300:   100     3 -1.62561674  1.027004619  2.7271964
dt[, choice := 0L]
#>         id   alt          x1           w1         w2 choice
#>      <int> <int>       <num>        <num>      <num>  <int>
#>   1:     1     1  1.37095845 -0.004620768 -0.2484829      0
#>   2:     1     2 -0.56469817  0.760242168  0.4223204      0
#>   3:     1     3  0.36312841  0.038990913  0.9876533      0
#>   4:     2     1  0.63286260  0.735072142  0.8355682      0
#>   5:     2     2  0.40426832 -0.146472627 -0.6605219      0
#>  ---                                                       
#> 296:    99     2 -0.47733551  0.160327395  0.1704735      0
#> 297:    99     3 -0.16626149 -0.433641942  1.2006682      0
#> 298:   100     1  0.86256338  1.537412419 -0.1634059      0
#> 299:   100     2  0.09734049 -2.170246577  1.2824759      0
#> 300:   100     3 -1.62561674  1.027004619  2.7271964      0
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#>         id   alt          x1           w1         w2 choice
#>      <int> <int>       <num>        <num>      <num>  <int>
#>   1:     1     1  1.37095845 -0.004620768 -0.2484829      0
#>   2:     1     2 -0.56469817  0.760242168  0.4223204      1
#>   3:     1     3  0.36312841  0.038990913  0.9876533      0
#>   4:     2     1  0.63286260  0.735072142  0.8355682      0
#>   5:     2     2  0.40426832 -0.146472627 -0.6605219      0
#>  ---                                                       
#> 296:    99     2 -0.47733551  0.160327395  0.1704735      1
#> 297:    99     3 -0.16626149 -0.433641942  1.2006682      0
#> 298:   100     1  0.86256338  1.537412419 -0.1634059      0
#> 299:   100     2  0.09734049 -2.170246577  1.2824759      0
#> 300:   100     3 -1.62561674  1.027004619  2.7271964      1

fit <- run_mxlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = "x1", random_var_cols = c("w1", "w2"), S = 50L
)
#> Optimization run time 0h:0m:0.01s
summary(fit)
#> Mixed Logit (MXL) model
#> 
#> Parameter    Estimate  Std.Error  z-value  Pr(>|z|)  
#> x1           0.005704   0.121525   0.0469  9.63e-01  
#> Sigma_11     0.000000   0.000000   0.0000  1.00e+00  
#> Sigma_22     0.000000   0.000000   0.0000  1.00e+00  
#> ASC_2        0.031083   0.248189   0.1252  9.00e-01  
#> ASC_3        0.089236   0.244715   0.3647  7.15e-01  
#> ---
#> Signif. codes:  '***' 0.001 '**' 0.01 '*' 0.05
#> 
#> Random coefficient covariance (Sigma):
#>              w1          w2
#> w1 3.768717e-20 0.00000e+00
#> w2 0.000000e+00 1.04647e-17
#> 
#> Std. Errors: Analytical Hessian 
#> Log-likelihood: -109.79 
#> AIC: 229.581  | BIC: 242.607 
#> McFadden R2: 0.001 (adj: -0.045) | Hit rate: 0.350 
#> N: 100  | Parameters: 5 
#> Optimization time: 0.01 s
#> Convergence: 3 ( NLOPT_FTOL_REACHED: Optimization stopped because ftol_rel or ftol_abs (above) was reached. )
# }