Package {choicer}

choicer: Discrete Choice Models for Economic Applications

Description

Fast estimation of discrete-choice models for applied economics. Frequentist likelihoods, analytical gradients, and Hessians are implemented in C++ with 'OpenMP' parallelism, scaling efficiently to specifications with many alternative-specific constants. Compiled Gibbs samplers provide Bayesian multinomial probit and hierarchical models. Post-estimation routines cover predicted shares, own- and cross-price elasticities, diversion ratios, willingness to pay, and welfare counterfactuals. Supports multinomial logit ('MNL'), mixed logit ('MXL'), nested logit ('NL'), Bayesian multinomial probit ('MNP'), and hierarchical Bayesian multinomial logit and probit ('HMNL', 'HMNP').

Author(s)

Maintainer: Fernando Cordeiro fernandolpcordeiro@gmail.com [copyright holder]

Authors:

• Fernando Cordeiro fernandolpcordeiro@gmail.com [copyright holder]

See Also

Useful links:

• https://github.com/fpcordeiro/choicer

• https://fpcordeiro.github.io/choicer/

• Report bugs at https://github.com/fpcordeiro/choicer/issues

BLP contraction mapping

Description

Finds the ASC (delta) parameters such that predicted market shares match target shares, using the contraction mapping of Berry, Levinsohn, and Pakes (1995) doi:10.2307/2171802.

Usage

blp(object, target_shares, ...)

Arguments

Value

Converged delta (ASC) vector.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
blp(fit, target_shares = rep(1/J, J))

BLP contraction mapping for multinomial logit model

Description

BLP contraction mapping for multinomial logit model

Usage

## S3 method for class 'choicer_mnl'
blp(
  object,
  target_shares,
  delta_init = NULL,
  tol = 1e-08,
  max_iter = 1000,
  ...
)

Arguments

Value

Converged delta (ASC) vector.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
blp(fit, target_shares = rep(1/J, J))

BLP contraction mapping for mixed logit model

Description

BLP contraction mapping for mixed logit model

Usage

## S3 method for class 'choicer_mxl'
blp(
  object,
  target_shares,
  delta_init = NULL,
  tol = 1e-08,
  max_iter = 1000,
  ...
)

Arguments

Value

Converged delta (ASC) vector.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mxlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = "x1", random_var_cols = "w1", S = 50L
)
blp(fit, target_shares = rep(1/J, J))

BLP contraction mapping for nested logit model

Description

BLP contraction mapping for nested logit model

Usage

## S3 method for class 'choicer_nl'
blp(
  object,
  target_shares,
  delta_init = NULL,
  damping = 1,
  tol = 1e-08,
  max_iter = 1000,
  ...
)

Arguments

Value

Converged delta (ASC) vector.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, nest := rep(c(1L, 1L, 2L, 2L), N)]
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
blp(fit, target_shares = rep(1/J, J))

BLP95 contraction mapping to find delta given target shares

Description

BLP95 contraction mapping to find delta given target shares

Usage

blp_contraction(
  delta,
  target_shares,
  X,
  beta,
  alt_idx,
  M,
  weights,
  include_outside_option = FALSE,
  tol = 1e-08,
  max_iter = 1000L
)

Arguments

Value

vector with contraction's delta (ASCs) output

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
beta <- coef(fit)[fit$param_map$beta]
delta <- blp_contraction(rep(0, J), rep(1/J, J), fit$data$X,
  beta, fit$data$alt_idx, fit$data$M, fit$data$weights)
delta

Reconstruct variance matrix L from L_params

Description

Reconstruct variance matrix L from L_params

Usage

build_var_mat(L_params, K_w, rc_correlation)

Arguments

Value

matrix equal to LL', where L is the choleski decomposition of random coefficient matrix

Examples

L_params <- c(log(1.0), 0.3, log(0.5))
Sigma <- choicer:::build_var_mat(L_params, K_w = 2, rc_correlation = TRUE)
Sigma  # 2x2 covariance matrix

Extract coefficients from a choicer_fit object

Description

Extract coefficients from a choicer_fit object

Usage

## S3 method for class 'choicer_fit'
coef(object, ...)

Arguments

Value

Named numeric vector of estimated coefficients.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
coef(fit)

Extract posterior means from a hierarchical Bayes fit

Description

Extract posterior means from a hierarchical Bayes fit

Usage

## S3 method for class 'choicer_hb'
coef(object, component = c("beta", "theta", "delta", "xi"), ...)

Arguments

Value

Named numeric vector of posterior means.

Examples

sim <- simulate_hmnl_data(N = 50, T = 2, J = 3, seed = 42)
fit <- suppressWarnings(run_hmnlogit(sim$data, "task", "alt", "choice", c("x1", "x2"),
                    person_col = "pid",
                    mcmc = list(R = 300, burn = 100)))
coef(fit)
coef(fit, component = "delta")

Extract coefficients from a choicer_mnp object

Description

Returns the posterior means of the identified coefficients (\beta / \sqrt{\sigma_{11}}, computed per draw).

Usage

## S3 method for class 'choicer_mnp'
coef(object, ...)

Arguments

Value

Named numeric vector of posterior means.

Examples

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), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnprobit(dt, "id", "alt", "choice", c("x1", "x2"),
                    mcmc = list(R = 300, burn = 100))
coef(fit)

Expected consumer surplus

Description

Computes the expected consumer surplus per choice situation (Train 2009, Ch. 3):

E[CS_i] = \frac{logsum_i}{-\alpha},

where logsum_i is the expected maximum utility (see logsum) and \alpha is the (fixed) price coefficient, so that -\alpha is the marginal utility of income. The formula assumes no income effects: utility is linear in price, and the marginal utility of income is constant across the price changes considered.

Usage

## S3 method for class 'choicer_hmnl'
consumer_surplus(
  object,
  price_var,
  newdata = NULL,
  level = 0.95,
  weights = NULL,
  n_draws = 200L,
  ...
)

## S3 method for class 'choicer_hmnp'
consumer_surplus(
  object,
  price_var,
  newdata = NULL,
  level = 0.95,
  weights = NULL,
  ...
)

consumer_surplus(
  object,
  price_var,
  newdata = NULL,
  level = 0.95,
  weights = NULL,
  ...
)

## S3 method for class 'choicer_mnl'
consumer_surplus(
  object,
  price_var,
  newdata = NULL,
  level = 0.95,
  weights = NULL,
  ...
)

## S3 method for class 'choicer_mxl'
consumer_surplus(
  object,
  price_var,
  newdata = NULL,
  level = 0.95,
  weights = NULL,
  ...
)

## S3 method for class 'choicer_nl'
consumer_surplus(
  object,
  price_var,
  newdata = NULL,
  level = 0.95,
  weights = NULL,
  ...
)

Arguments

Details

Consumer surplus levels inherit the additive utility normalization (in particular the ASC normalization), so the level is only defined up to a constant; differences in CS between scenarios — e.g. consumer_surplus(fit, "price", newdata = scenario) minus the baseline — are the economically meaningful quantity.

For MNL fits, a delta-method standard error of the weighted mean CS is reported (weights are the stored fit weights, or the resolved newdata weights). For MXL and NL fits only point estimates are returned (se_mean_cs = NA): the delta method for the simulated MXL logsum and the nested logsum is deferred; simulation-based intervals (Krinsky-Robb: resample coefficients from their asymptotic distribution and recompute the mean CS) are a practical alternative.

The price variable must have a fixed coefficient. For mixed logit a random price coefficient is rejected (as in wtp): with a random denominator 1/(-\alpha) generally has no finite moments.

Value

A choicer_cs object: a list with cs (per-choice- situation surplus, length N), mean_cs (weighted mean), se_mean_cs (delta-method SE; NA for MXL/NL or when the variance-covariance matrix is unavailable), ci (confidence interval for the mean), price_var, level, and n.

Methods (by class)

• consumer_surplus(choicer_hmnl): Posterior consumer surplus for the hierarchical logit: per-task logsum divided by the (positive) marginal utility of income -\bar\gamma_{price}, per posterior draw. With newdata, the return also carries the compensating variation against the estimation data (attr(, "cv")), i.e. the posterior of (\mathrm{logsum}_{new} - \mathrm{logsum}_{base}) / (-\bar\gamma_{price}) summed over tasks. Requires a fixed-sign price coefficient; the posterior-median ratio discipline of wtp.choicer_hb() applies. When newdata is supplied, weights is an optional non-negative length-n_tasks vector used in the aggregate CV; equal task weights are the default. The counterfactual must contain the same choice situations as the estimation data, identified by (person_col, id_col) (or id_col when person_col = NULL); rows and tasks may be reordered. Unnamed weights follow the baseline tasks' sorted (person_col, id_col) order used by the prepared data.

• consumer_surplus(choicer_hmnp): Not available for the probit (see logsum.choicer_hmnp()); roadmapped via simulated Emax.

References

Train, K. (2009). Discrete Choice Methods with Simulation, 2nd ed., Ch. 3. Cambridge University Press.

See Also

logsum, wtp

Examples

library(data.table)
sim <- simulate_mnl_data(N = 1000, J = 3, beta = c(0.8, -0.6), seed = 123,
                         outside_option = FALSE, vary_choice_set = FALSE)
fit <- run_mnlogit(sim$data, "id", "alt", "choice", c("x1", "x2"))

# treat x2 as the price variable
cs0 <- consumer_surplus(fit, price_var = "x2")
cs0

# Change in consumer surplus from a price increase on alternative 2:
# levels depend on the ASC normalization, differences do not.
dt_cf <- copy(sim$data)[alt == 2, x2 := x2 + 0.5]
cs1 <- consumer_surplus(fit, price_var = "x2", newdata = dt_cf)
delta_cs <- cs1$mean_cs - cs0$mean_cs
delta_cs  # negative: the price increase lowers expected surplus

Compute aggregate diversion ratios

Description

Computes a J x J matrix of diversion ratios. Entry (i, j) is the fraction of demand lost by alternative j that is captured by alternative i when alternative j becomes less attractive.

Usage

## S3 method for class 'choicer_hb'
diversion_ratios(object, elast_var, eps = 0.01, n_draws = 100L, ...)

diversion_ratios(object, ...)

Arguments

Value

A J x J diversion ratio matrix with alternative labels.

Methods (by class)

• diversion_ratios(choicer_hb): Posterior-mean diversion ratios for hierarchical Bayes fits, from the same perturbation engine as elasticities.choicer_hb(): DR(j \to k) is the fraction of the share alternative j loses (when its elast_var worsens) that flows to k — including the outside option. Columns are the perturbed alternative j, rows the receiving alternative k; the diagonal is 0.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
diversion_ratios(fit)

Diversion ratios for multinomial logit model

Description

Diversion ratios for multinomial logit model

Usage

## S3 method for class 'choicer_mnl'
diversion_ratios(object, ...)

Arguments

Value

A J x J diversion ratio matrix with alternative labels.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
diversion_ratios(fit)

Diversion ratios for mixed logit model

Description

Computes the attribute-based diversion ratio matrix. Entry (k, j) is the fraction of demand lost by alternative j that is captured by alternative k when a marginal change in alternative j's wrt_var attribute reduces s_j.

Usage

## S3 method for class 'choicer_mxl'
diversion_ratios(object, wrt_var, is_random_coef = FALSE, ...)

Arguments

Details

Unlike MNL, the MXL diversion ratio depends on which variable is perturbed: the realised coefficient \beta_{ik}^s varies across individuals and draws and does not cancel in the ratio. For a variable with a fixed coefficient the result is independent of the variable (\beta cancels); for a random-coefficient variable it is not.

Value

A J x J diversion ratio matrix with alternative labels. Cross-products are averaged across simulation draws inside the integration to avoid Jensen-style bias.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mxlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = "x1", random_var_cols = "w1", S = 50L
)
diversion_ratios(fit, "x1")
diversion_ratios(fit, "w1", is_random_coef = TRUE)

Diversion ratios for nested logit model

Description

Diversion ratios for nested logit model

Usage

## S3 method for class 'choicer_nl'
diversion_ratios(object, ...)

Arguments

Value

A J x J diversion ratio matrix with alternative labels.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, nest := rep(c(1L, 1L, 2L, 2L), N)]
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
diversion_ratios(fit)

Compute aggregate elasticities

Description

Computes a J x J matrix of aggregate elasticities. Entry (i, j) is the percentage change in the probability of choosing alternative i when the attribute of alternative j changes by 1\

Usage

## S3 method for class 'choicer_hb'
elasticities(object, elast_var, eps = 0.01, n_draws = 100L, ...)

elasticities(object, ...)

Arguments

Value

A J x J elasticity matrix with alternative labels.

Methods (by class)

• elasticities(choicer_hb): Posterior-mean aggregate arc elasticities for hierarchical Bayes fits: each inside alternative's elast_var is perturbed by eps (default 1%) and shares are recomputed per posterior draw with a common random-coefficient path, giving E_{jk} = (\Delta s_j / s_j) / \epsilon. Rows are responding alternatives (including the outside option), columns the perturbed alternative.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
elasticities(fit, "x1")

Elasticities for multinomial logit model

Description

Elasticities for multinomial logit model

Usage

## S3 method for class 'choicer_mnl'
elasticities(object, elast_var, ...)

Arguments

Value

A J x J elasticity matrix with alternative labels.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
elasticities(fit, "x1")

Elasticities for mixed logit model

Description

Elasticities for mixed logit model

Usage

## S3 method for class 'choicer_mxl'
elasticities(object, elast_var, is_random_coef = FALSE, ...)

Arguments

Value

A J x J elasticity matrix with alternative labels.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mxlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = "x1", random_var_cols = "w1", S = 50L
)
elasticities(fit, "x1")
elasticities(fit, "w1", is_random_coef = TRUE)

Elasticities for nested logit model

Description

Elasticities for nested logit model

Usage

## S3 method for class 'choicer_nl'
elasticities(object, elast_var, ...)

Arguments

Value

A J x J elasticity matrix with alternative labels.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, nest := rep(c(1L, 1L, 2L, 2L), N)]
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
elasticities(fit, "x1")

Rank-normalized effective sample size (bulk and tail)

Description

Computes the rank-normalized bulk and tail effective sample size (ESS) of Vehtari, Gelman, Simpson, Carpenter & Buerkner (2021, Bayesian Analysis) for each column of a matrix of posterior draws. ess_bulk measures the effective number of independent draws available for estimating the posterior mean/median; ess_tail is the minimum of the 5th- and 95th-percentile indicator ESS, relevant for tail-quantile / credible-interval precision. Autocovariance is computed via fft (Geyer's initial monotone sequence estimator of the integrated autocorrelation time).

Usage

ess(draws)

Arguments

Value

A numeric matrix, one row per parameter, two columns ("bulk", "tail"); NA for parameters with zero variance or a non-finite result.

Examples

set.seed(42)
draws <- matrix(rnorm(2000), ncol = 2,
                dimnames = list(NULL, c("a", "b")))
ess(draws)

Halton draws for mixed logit

Description

Create halton normal draws in appropriate format for mixed logit estimation

Usage

get_halton_normals(S, N, K_w)

Arguments

Value

K_w x S x N array with halton standard normal draws

Examples

draws <- get_halton_normals(S = 50, N = 10, K_w = 2)
dim(draws)  # 2 x 50 x 10

Goodness of fit for a fitted choice model

Description

Computes McFadden's pseudo R-squared (plain and adjusted) and the in-sample hit rate for a fitted model.

Usage

gof(object, null = c("equal_shares", "market_shares"), ...)

## S3 method for class 'choicer_fit'
gof(object, null = c("equal_shares", "market_shares"), ...)

Arguments

Details

Two null models are available for the pseudo R-squared R^2 = 1 - LL / LL_0 (adjusted: R^2_{adj} = 1 - (LL - K) / LL_0 with K the number of estimated parameters):

• "equal_shares" (default): every alternative in individual i's choice set is equally likely, so LL_0 = -\sum_i w_i \log(M_i + 1_{outside}). This is exact for unbalanced choice sets and arbitrary weights.

• "market_shares": the maximized log-likelihood of an ASC-only model, LL_0 = \sum_j N_j \log(s_j) with N_j the choice counts and s_j the observed market shares (including the outside option when present). This closed form is valid only for balanced choice sets and uniform weights; otherwise an error suggests refitting an ASC-only model.

The hit rate is the weighted share of individuals whose observed choice has the highest predicted probability. When the model includes an outside option, the outside good competes for the predicted maximum (its probability is 1 - \sum_j p_{ij}), and an individual predicted to choose the outside good is a hit when they actually did.

Both the null log-likelihood and the hit rate require the stored estimation data; models fitted with keep_data = FALSE return NA fields with a message.

Value

A choicer_gof object: a list with loglik, loglik_null, null, mcfadden_r2, mcfadden_r2_adj, hit_rate, nobs, and n_params.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
gof(fit)
gof(fit, null = "market_shares")

Gibbs sampler for the hierarchical Bayesian multinomial logit model

Description

Runs the adaptive random-walk Metropolis-within-Gibbs sampler for the hierarchical (random-coefficients, panel) multinomial logit with a BLP-style alternative-level random effect: inside utilities U_{ijt} = x_{ijt}'\gamma_i + \delta_j + EV1 against an implicit outside option with systematic utility 0, \beta_i \sim N(b, W) (\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).

Usage

hmnl_gibbs(
  X,
  Z,
  M,
  choice_pos,
  include_outside_option,
  alt_of_row,
  Ti,
  rc_dist,
  beta_pooled,
  delta_init,
  theta_init,
  b_bar,
  A,
  nu,
  V,
  theta_bar,
  A_theta,
  sd_prior,
  R,
  burn,
  thin,
  seed,
  keep_beta_i,
  s_init,
  accept_target,
  trace = 0L
)

Arguments

Details

The per-respondent \beta_i updates are parallelized with OpenMP; the \delta_j updates run as a strictly serial sweep (their conditionals are coupled through the shared softmax denominators). Each (iteration, unit) pair uses its own RNG stream, so draws are reproducible given the seed and a fixed thread count; across different thread counts they are invariant only up to floating-point reduction-order round-off (~1e-15), not bitwise (see set_num_threads()). This is the low-level engine behind run_hmnlogit, which handles initialization and post-processing.

Value

List with bdraw (R_keep x K), wdraw (R_keep x K(K+1)/2, lower triangle of W in row-major order), deltadraw (R_keep x J), thetadraw (R_keep x P), sigma_d2draw, loglik_trace, acceptance rates and final proposal scales (accept_rate_beta, accept_rate_delta, s_final, s_delta_final), posterior summaries beta_i_mean / beta_i_sd (K x N, NULL when keep_beta_i = 0), beta_i_draws (K x N x R_keep cube when keep_beta_i = 2), delta_mean / delta_sd / xi_mean / xi_sd (J x 1), and R_keep.

Examples

sim <- simulate_hmnl_data(N = 20, T = 2, J = 3, seed = 42)
d <- prepare_hmnl_data(sim$data, "task", "alt", "choice",
                       c("x1", "x2"), person_col = "pid")
out <- choicer:::hmnl_gibbs(d$X, d$Z, d$M, d$choice_pos, TRUE, d$alt_of_row, d$Ti,
  rc_dist = d$rc_dist, beta_pooled = rep(0, d$K_struct),
  delta_init = rep(0, d$J), theta_init = rep(0, d$P),
  b_bar = rep(0, d$K_struct), A = 0.01 * diag(d$K_struct),
  nu = d$K_struct + 3, V = (d$K_struct + 3) * diag(d$K_struct),
  theta_bar = rep(0, d$P), A_theta = 0.01 * diag(d$P),
  sd_prior = list(half_cauchy = TRUE, s_d = 1, c0 = 3, d0 = 3),
  R = 300, burn = 100, thin = 1, seed = 7, keep_beta_i = 1,
  s_init = 2.38 / sqrt(d$K_struct), accept_target = 0.234)
colMeans(out$bdraw)

Gibbs sampler for the hierarchical Bayesian multinomial probit model

Description

Runs the fully conjugate Albert-Chib Gibbs sampler for the hierarchical multinomial probit with iid N(0, \sigma^2) utility shocks in un-differenced utility space: inside utilities U_{ijt} = x_{ijt}'\beta_i + \delta_j + \epsilon against a stochastic implicit outside option U_{iot} = \epsilon, \beta_i \sim N(b, W), and \delta_j = z_j'\theta + \xi_j, \xi_j \sim N(0, \sigma_d^2). The chain runs on the non-identified parameterization (free \sigma^2); identified quantities are obtained by normalizing each draw by the matching power of \sigma (handled by run_hmnprobit).

Usage

hmnp_gibbs(
  X,
  Z,
  M,
  choice_pos,
  include_outside_option,
  alt_of_row,
  Ti,
  delta_init,
  theta_init,
  b_bar,
  A,
  nu,
  V,
  theta_bar,
  A_theta,
  sd_prior,
  a0,
  s0,
  R,
  burn,
  thin,
  seed,
  keep_beta_i,
  trace = 0L
)

Arguments

Details

The latent sweep and the \beta_i draws are parallelized with OpenMP across respondents; the \delta_j draws are parallelized across alternatives (conditionally independent given the augmented utilities — unlike the HMNL, whose delta sweep must be serial). Each (iteration, unit) pair uses its own RNG stream, so draws are reproducible given the seed and a fixed thread count; across different thread counts they are invariant only up to floating-point reduction-order round-off (~1e-15), not bitwise.

Value

List with RAW draw matrices bdraw, wdraw (lower triangle, row-major), deltadraw, thetadraw, sigma_d2draw, sigma2draw, identified-scale summaries beta_i_mean / beta_i_sd / beta_i_draws / delta_mean / delta_sd / xi_mean / xi_sd, and R_keep.

Examples

sim <- simulate_hmnp_data(N = 30, T = 2, J = 3, seed = 42)
d <- prepare_hmnp_data(sim$data, "task", "alt", "choice",
                       c("x1", "x2"), person_col = "pid")
out <- choicer:::hmnp_gibbs(d$X, d$Z, d$M, d$choice_pos, TRUE, d$alt_of_row, d$Ti,
  delta_init = rep(0, d$J), theta_init = rep(0, d$P),
  b_bar = rep(0, d$K_struct), A = 0.01 * diag(d$K_struct),
  nu = d$K_struct + 3, V = (d$K_struct + 3) * diag(d$K_struct),
  theta_bar = rep(0, d$P), A_theta = 0.01 * diag(d$P),
  sd_prior = list(half_cauchy = TRUE, s_d = 1, c0 = 3, d0 = 3),
  a0 = 3, s0 = 3, R = 300, burn = 100, thin = 1, seed = 7,
  keep_beta_i = 1)
colMeans(out$bdraw / sqrt(as.numeric(out$sigma2draw)))

Utility to compute analytical Jacobian of random coefficient matrix transformed by vech (dVech(Sigma) / dTheta)

Description

Utility to compute analytical Jacobian of random coefficient matrix transformed by vech (dVech(Sigma) / dTheta)

Usage

jacobian_vech_Sigma(L_params, K_w, rc_correlation = TRUE)

Arguments

Value

Jacobian (dVech(Sigma) / dTheta)

Examples

L_params <- c(log(0.8), 0.2, log(0.6))
J_mat <- choicer:::jacobian_vech_Sigma(L_params, K_w = 2, rc_correlation = TRUE)
dim(J_mat)  # 3 x 3 for K_w=2 correlated

Extract log-likelihood from a choicer_fit object

Description

Returns a logLik object, which enables AIC() and BIC() automatically.

Usage

## S3 method for class 'choicer_fit'
logLik(object, ...)

Arguments

Value

A logLik object with df and nobs attributes.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
logLik(fit)
AIC(fit)
BIC(fit)

Expected logsum (inclusive value) per choice situation

Description

Computes the expected maximum utility ("logsum" or inclusive value) for each choice situation, up to the additive constant of the extreme-value error:

• MNL: \log \sum_j \exp(V_{ij}).

• MXL: E_\beta[\log \sum_j \exp(V_{ij}(\beta))], simulated by averaging the log-sum-exp across the deterministic Halton draws (regenerated from object$draws_info). Taking the log-sum-exp of draw-averaged utilities would understate the expectation (Jensen's inequality), so a dedicated kernel (mxl_logsum) is used.

• NL: \log \sum_b \exp(\lambda_b I_{ib}) with the nest inclusive value I_{ib} = \log \sum_{j \in b} \exp(V_{ij}/\lambda_b) (singleton nests have \lambda_b = 1).

When the model includes an outside option, its normalized utility V = 0 contributes an \exp(0) term to the sum.

Usage

## S3 method for class 'choicer_hmnl'
logsum(object, newdata = NULL, n_draws = 200L, ...)

## S3 method for class 'choicer_hmnp'
logsum(object, newdata = NULL, ...)

logsum(object, newdata = NULL, ...)

## S3 method for class 'choicer_mnl'
logsum(object, newdata = NULL, ...)

## S3 method for class 'choicer_mxl'
logsum(object, newdata = NULL, ...)

## S3 method for class 'choicer_nl'
logsum(object, newdata = NULL, ...)

Arguments

Details

Logsum levels depend on the ASC normalization (and, more generally, on any additive utility normalization), so only logsum differences between scenarios (e.g. via newdata) are meaningful.

Value

Numeric vector with one logsum per choice situation. With a data.frame newdata, choice situations are ordered by id (as in predict()).

Methods (by class)

• logsum(choicer_hmnl): Posterior expected logsum for the hierarchical logit: per choice situation, \log(1 + \sum_j \exp V_j) against the outside-option anchor, averaged over posterior draws with one \beta \sim N(b_r, W_r) draw each. Returns the per-task posterior mean vector.

• logsum(choicer_hmnp): The probit expected maximum has no closed form; simulated-Emax surplus for the HMNP is on the roadmap.

See Also

consumer_surplus

Examples

library(data.table)
sim <- simulate_mnl_data(N = 500, J = 3, beta = c(0.8, -0.6), seed = 1,
                         outside_option = FALSE, vary_choice_set = FALSE)
fit <- run_mnlogit(sim$data, "id", "alt", "choice", c("x1", "x2"))
head(logsum(fit))

Asymptotic diagnostics for a Monte Carlo study

Description

Consumes a choicer_mc object and returns per-parameter asymptotic diagnostics: Monte Carlo bias (with MC standard error), empirical SD of the estimates, mean of the reported standard errors, SE-to-SD ratio (information-matrix-equality check), Wald coverage at nominal 90 / 95 / 99 percent with Wilson confidence bands, moments of the studentized statistic z = (theta_hat - theta_0) / se, and four normality tests on z (Shapiro-Wilk, Anderson-Darling via goftest::ad.test, a hand-coded Jarque-Bera statistic, and a one-sample Kolmogorov-Smirnov test against N(0, 1)).

Usage

mc_asymptotics(
  mc,
  level = 0.95,
  se_col = "se",
  conv_threshold = 0.99,
  se_ratio_threshold_floor = 0.1
)

Arguments

Details

Six logical pass / fail flags are attached to every parameter row: pass_bias requires ⁠|bias_mc_se| < 3⁠; pass_se_ratio requires ⁠|se_ratio - 1|⁠ to lie within max(se_ratio_threshold_floor, 3 * 1.4 / sqrt(R_used)) (a noise-aware band that widens at small R_used and tightens to the floor at large R_used); pass_cov95 requires the nominal 95 percent level to lie in the Wilson band for empirical coverage; pass_skew requires ⁠|skew_z| < 0.3⁠; pass_kurt requires excess kurtosis of z in ⁠[-0.5, 1.0]⁠; pass_convergence requires the per-parameter convergence rate (R_used / R_total) to meet conv_threshold.

Non-converged replications are excluded per parameter (reported in R_excluded). Winsorized (5 percent / 95 percent) versions of bias, sd_emp, and mean_se are reported in parallel columns (bias_w, sd_emp_w, mean_se_w) so silent outlier exclusion is transparent to the reader. Two robust SE-to-SD ratios accompany the Hessian-mean-based se_ratio: se_ratio_med (median SE divided by the empirical SD) and se_ratio_w (winsorized mean SE divided by the winsorized empirical SD); both stay near 1 when 1-2 replications produce near-singular Hessians that inflate mean_se. The companion se_med column reports the median per-replication SE used by se_ratio_med. Neither robust ratio drives a ⁠pass_*⁠ flag — they are purely informational.

Winsorized z-moment counterparts (mean_z_w, sd_z_w, skew_z_w, kurt_excess_z_w) are reported alongside the raw z-moments and feed an additional pass_z_w flag (Winsorized skew within the same band as pass_skew AND Winsorized excess kurtosis within the same band as pass_kurt). A companion pass_cov95_w flag is TRUE when either pass_cov95 is TRUE OR the per-rep Winsorized z-CI (the empirical 2.5 / 97.5 percentiles of the Winsorized z) covers truth-zero. These two flags are designed for boundary scenarios (e.g., near-zero variance components) where a small number of reps with vanishing SE inflate the raw z-moments without indicating an estimator defect.

Value

An object of class choicer_mc_asymptotics — a data.table with one row per unique parameter and columns documented above — with meta attached as an attribute (attr(x, "meta")).

Examples

sim_fun <- function(seed) simulate_mnl_data(N = 1000, J = 3, seed = seed)
fit_fun <- function(sim) run_mnlogit(
  data = sim$data, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = c("x1", "x2"), outside_opt_label = 0L,
  include_outside_option = FALSE, use_asc = TRUE,
  control = list(print_level = 0L)
)
mc <- monte_carlo(sim_fun, fit_fun, R = 50L, seed = 1L, progress = FALSE)
mc_asymptotics(mc)

Monte Carlo standard error of posterior summaries

Description

Monte Carlo standard error (MCSE) of the posterior mean or median, using the rank-normalized bulk/tail ESS from ess. kind = "mean": MCSE = SD / sqrt(ess_bulk). kind = "median": MCSE = sqrt(pi/2) * SD / sqrt(ess_tail) (the closed-form normal asymptotic-efficiency approximation; 2/pi is the asymptotic relative efficiency of the sample median vs. the mean under normality).

Usage

mcse(draws, kind = c("mean", "median"))

Arguments

Value

Named numeric vector, one value per parameter (NA when the underlying ESS or pooled SD is undefined).

Examples

set.seed(42)
draws <- matrix(rnorm(2000), ncol = 2,
                dimnames = list(NULL, c("a", "b")))
mcse(draws)
mcse(draws, kind = "median")

BHHH/OPG information matrix for multinomial logit model

Description

Computes the weighted outer product of per-individual scores \sum_i w_i\, s_i s_i^\top for the Multinomial Logit model. The per-individual score s_i is the (positive) gradient of individual i's log-likelihood contribution and is weight-free; the supplied weights enter only as the leading multiplier. Passing weights = w yields the ordinary weighted BHHH/OPG information; passing weights = w^2 yields the sandwich meat B = \sum_i w_i^2 s_i s_i^\top used for robust (WESML) inference.

Usage

mnl_bhhh_parallel(
  theta,
  X,
  alt_idx,
  choice_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

A symmetric positive-semidefinite information matrix \sum_i w_i\, s_i s_i^\top (same sign convention as the negated Hessian).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
B <- choicer:::mnl_bhhh_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
  fit$data$choice_idx, fit$data$M, fit$data$weights)
dim(B)

Compute MNL diversion ratios (parallelized over individuals)

Description

Computes the diversion ratio matrix DR(j->k), which measures the fraction of demand lost by alternative j that is captured by alternative k. For MNL: DR(j->k) = sum_n(w_n * P_nj * P_nk) / sum_n(w_n * P_nj * (1 - P_nj))

Usage

mnl_diversion_ratios_parallel(
  theta,
  X,
  alt_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

J x J matrix where entry (k, j) = DR(j->k). Diagonal is 0.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
dr <- choicer:::mnl_diversion_ratios_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
  fit$data$M, fit$data$weights)
dr

Compute aggregate elasticities for MNL model

Description

Computes the aggregate elasticity matrix (weighted average of individual elasticities) for the Multinomial Logit model.

Usage

mnl_elasticities_parallel(
  theta,
  X,
  alt_idx,
  choice_idx,
  M,
  weights,
  elast_var_idx,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

J x J matrix of aggregate elasticities

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
elas <- choicer:::mnl_elasticities_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
  fit$data$choice_idx, fit$data$M, fit$data$weights, elast_var_idx = 1L)
elas

Log-likelihood and gradient for multinomial logit model

Description

Computes the log-likelihood and its gradient for the Multinomial Logit model using OpenMP for parallelization. Allows for inclusion of alternative-specific constants, outside option, and observation weights.

Usage

mnl_loglik_gradient_parallel(
  theta,
  X,
  alt_idx,
  choice_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

List with loglikelihood and gradient evaluated at input arguments

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_mnl_data(dt, "id", "alt", "choice", c("x1", "x2"))
theta <- rep(0, ncol(d$X) + nrow(d$alt_mapping) - 1)
result <- choicer:::mnl_loglik_gradient_parallel(theta, d$X, d$alt_idx,
  d$choice_idx, d$M, d$weights)
result$objective  # negative log-likelihood

Hessian matrix for multinomial logit model

Description

Hessian matrix for multinomial logit model

Usage

mnl_loglik_hessian_parallel(
  theta,
  X,
  alt_idx,
  choice_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

Hessian matrix of the negative log-likelihood

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
H <- choicer:::mnl_loglik_hessian_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
  fit$data$choice_idx, fit$data$M, fit$data$weights)
dim(H)

Prediction of choice probabilities and utilities based on fitted model

Description

Prediction of choice probabilities and utilities based on fitted model

Usage

mnl_predict(
  theta,
  X,
  alt_idx,
  M,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

List with choice probability and utility for each choice situation evaluated at input arguments

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
pred <- choicer:::mnl_predict(coef(fit), fit$data$X, fit$data$alt_idx,
  fit$data$M, use_asc = TRUE)
head(pred$choice_prob)

Prediction of market shares based on fitted model

Description

Prediction of market shares based on fitted model

Usage

mnl_predict_shares(
  theta,
  X,
  alt_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

vector with predicted market shares for each alternative

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
shares <- choicer:::mnl_predict_shares(coef(fit), fit$data$X, fit$data$alt_idx,
  fit$data$M, fit$data$weights, use_asc = TRUE)
shares

Gibbs sampler for the Bayesian multinomial probit model

Description

Runs the McCulloch-Rossi (1994) Gibbs sampler with Albert-Chib data augmentation for the multinomial probit model in utility differences against a base alternative. The chain operates on the non-identified parameterization (unrestricted Sigma); identified quantities are obtained by normalizing each draw by sigma_11 (handled by run_mnprobit).

Usage

mnp_gibbs(X, y, p, beta_bar, A, nu, V, R, burn, thin, seed, trace = 0L)

Arguments

Details

The latent-utility sweep is parallelized with OpenMP across choice situations (they are conditionally independent given beta and Sigma). Each (iteration, observation) pair uses its own RNG stream, so draws are reproducible given the seed and a fixed thread count; across different thread counts they are invariant only up to floating-point reduction-order round-off (~1e-15), not bitwise (see set_num_threads()).

Value

List with betadraw (R_keep x K), sigmadraw (R_keep x p(p+1)/2, lower triangle of Sigma in row-major order), and R_keep.

Examples

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), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_mnp_data(dt, "id", "alt", "choice", c("x1", "x2"))
out <- choicer:::mnp_gibbs(d$X, d$y, d$p,
  beta_bar = rep(0, d$K), A = 0.01 * diag(d$K),
  nu = d$p + 3, V = (d$p + 3) * diag(d$p),
  R = 500, burn = 100, thin = 1, seed = 42)
colMeans(out$betadraw)

Intercity travel mode choice

Description

Stated choices of intercity travel mode for 210 travellers, each choosing among the same four modes: air, train, bus and car. This is the classic Greene & Hensher (1997) data set, reshaped into choicer's long layout (one row per traveller-by-alternative). It is a convenient, recognizable example for multinomial and nested logit models and the demand/welfare toolkit (elasticities, diversion ratios, willingness-to-pay, counterfactuals).

Usage

mode_choice

Format

A data frame with 840 rows (210 travellers x 4 modes) and 9 columns:

id

Integer traveller (choice situation) identifier, 1-210.

mode

Factor giving the travel mode: "air", "train", "bus" or "car". Use as the alternative column.

choice

Integer indicator, 1 for the chosen mode and 0 otherwise. Exactly one mode is chosen per traveller.

wait

Terminal waiting time in minutes (0 for car).

travel

In-vehicle travel time in minutes.

vcost

In-vehicle cost component, in currency units.

gcost

Generalized cost measure, in currency units.

income

Household income (traveller level, in thousands).

size

Size of the travelling party (traveller level).

Details

wait, travel, vcost and gcost vary across modes within a traveller, while income and size are traveller-level attributes that are constant across modes. A standard specification regresses the choice on wait, travel and vcost with alternative-specific constants; vcost then plays the role of price for willingness-to-pay and consumer-surplus calculations.

The sample is choice-based: the survey over-sampled the less popular modes (air, train, bus) and under-sampled car, so sample choice shares do not estimate population mode shares. With a full set of alternative-specific constants the slope coefficients remain consistently estimated under this design (Manski and Lerman, 1977), and willingness-to-pay ratios are unaffected; the constants, and any shares, elasticities or surplus levels computed from fitted probabilities, inherit the sampling design. To target population quantities, attach WESML weights with wesml_weights using external population shares; see vignette("wesml", package = "choicer").

Source

Greene, W. H. and Hensher, D. A. (1997). Reshaped from the TravelMode data distributed with the AER package (https://CRAN.R-project.org/package=AER). The same data appear in Greene's Econometric Analysis and in several other choice-modelling packages.

References

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

Examples

data(mode_choice)
head(mode_choice)
table(mode_choice$mode[mode_choice$choice == 1L])

Monte Carlo parameter recovery

Description

Replicates a (DGP -> fit) cycle R times with independent seeds and collects per-parameter estimates, standard errors, bias, and coverage. Returns a choicer_mc object; call summary() for aggregated statistics (mean estimate, bias, RMSE, coverage rate, convergence rate).

Usage

monte_carlo(
  sim_fun,
  fit_fun,
  R = 100,
  seed = 1L,
  parallel = FALSE,
  progress = TRUE,
  ...
)

Arguments

Details

Each iteration calls sim_fun(seed = seed + r - 1L), then fit_fun(sim). Write sim_fun as a closure that captures N, J, and other DGP settings and forwards seed. Write fit_fun as a closure that takes a choicer_sim and returns a fitted choicer_fit object, wrapping any data-preparation, draws, or optimizer-control setup.

Value

A choicer_mc object: a list with elements replications (a long data.table with one row per estimated parameter per replication) and meta (run metadata).

Examples

sim_fun <- function(seed) simulate_mnl_data(N = 1000, J = 4, seed = seed)
fit_fun <- function(sim) run_mnlogit(
  data = sim$data, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = c("x1", "x2"), outside_opt_label = 0L,
  include_outside_option = FALSE, use_asc = TRUE,
  control = list(print_level = 0L)
)
mc <- monte_carlo(sim_fun, fit_fun, R = 5, seed = 1L, progress = FALSE)
summary(mc)

BHHH (outer product of gradients) information matrix for Mixed Logit

Description

Computes the BHHH approximation to the observed information matrix for the Mixed Logit model: H_{BHHH} = \sum_i w_i \cdot s_i s_i^\top, where s_i is the per-individual score (gradient of \log \bar{P}_i). This outer product of gradients (OPG) estimator provides an alternative to the analytical Hessian for standard error computation that scales to large problems where the analytical Hessian is infeasible (e.g., many alternatives or simulation draws).

Usage

mxl_bhhh_parallel(
  theta,
  X,
  W,
  alt_idx,
  choice_idx,
  M,
  weights,
  eta_draws,
  rc_dist,
  rc_correlation = TRUE,
  rc_mean = FALSE,
  use_asc = TRUE,
  include_outside_option = FALSE,
  gen_seed = -1L,
  gen_scramble = 1L,
  gen_S = 0L
)

Arguments

Value

n_params x n_params PSD matrix representing the observed information matrix estimated by the outer product of gradients (same sign convention as the negated Hessian returned by mxl_hessian_parallel, so it can be inverted directly to obtain vcov).

Note

The BHHH/OPG estimator is only asymptotically equivalent to the Hessian-based information matrix at the true MLE. In finite samples it can underestimate standard errors, particularly when the model is mis-specified or away from the optimum.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
eta <- get_halton_normals(50, d$N, ncol(d$W))
theta <- rep(0, ncol(d$X) + ncol(d$W) + nrow(d$alt_mapping) - 1)
H <- choicer:::mxl_bhhh_parallel(theta, d$X, d$W, d$alt_idx, d$choice_idx,
  d$M, d$weights, eta, rc_dist = rep(0L, ncol(d$W)),
  rc_correlation = FALSE, rc_mean = FALSE)
dim(H)

BLP contraction mapping for mixed logit

Description

Finds the ASC (delta) parameters such that predicted market shares match target shares, using the contraction mapping of Berry, Levinsohn, and Pakes (1995).

Usage

mxl_blp_contraction(
  delta,
  target_shares,
  X,
  W,
  beta,
  mu,
  L_params,
  alt_idx,
  M,
  weights,
  eta_draws,
  rc_dist,
  rc_correlation = TRUE,
  rc_mean = FALSE,
  include_outside_option = FALSE,
  tol = 1e-08,
  max_iter = 1000L,
  gen_seed = -1L,
  gen_scramble = 1L,
  gen_S = 0L
)

Arguments

Value

vector with converged delta (ASC) values

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
eta <- get_halton_normals(50, d$N, ncol(d$W))
fit <- run_mxlogit(input_data = d, eta_draws = eta)
pm <- fit$param_map
delta <- mxl_blp_contraction(rep(0, J), rep(1/J, J), d$X, d$W,
  coef(fit)[pm$beta], rep(0, ncol(d$W)), coef(fit)[pm$sigma],
  d$alt_idx, d$M, d$weights, eta, rc_dist = rep(0L, ncol(d$W)),
  rc_correlation = FALSE, rc_mean = FALSE)
delta

Diversion ratios for Mixed Logit (simulated, derivative-based)

Description

Computes the matrix of attribute-based diversion ratios for a fitted Mixed Logit model. DR(k, j) is the fraction of demand lost by alternative j that is captured by alternative k when a marginal change in alternative j's elast_var attribute reduces s_j.

Usage

mxl_diversion_ratios_parallel(
  theta,
  X,
  W,
  alt_idx,
  M,
  weights,
  eta_draws,
  rc_dist,
  elast_var_idx,
  is_random_coef,
  rc_correlation = TRUE,
  rc_mean = FALSE,
  use_asc = TRUE,
  include_outside_option = FALSE,
  gen_seed = -1L,
  gen_scramble = 1L,
  gen_S = 0L
)

Arguments

Details

In MNL the per-draw realized coefficient is a constant, so it cancels in the ratio and the result is independent of the variable chosen. In MXL, the realized coefficient \beta_{ik}^s varies across individuals and draws, so the diversion ratio depends on which attribute is perturbed. For a variable with a fixed coefficient the dependence again vanishes (the constant cancels); for a random-coefficient variable it does not.

Value

J x J (or (J+1) x (J+1)) matrix of diversion ratios with zero diagonal.

Compute aggregate elasticities for mixed logit model

Description

Computes the aggregate elasticity matrix (weighted average of individual elasticities) for the Mixed Logit model. The elasticity E(i,j) represents the percentage change in the probability of choosing alternative i when the attribute of alternative j changes by 1%.

Usage

mxl_elasticities_parallel(
  theta,
  X,
  W,
  alt_idx,
  choice_idx,
  M,
  weights,
  eta_draws,
  rc_dist,
  elast_var_idx,
  is_random_coef,
  rc_correlation = TRUE,
  rc_mean = FALSE,
  use_asc = TRUE,
  include_outside_option = FALSE,
  gen_seed = -1L,
  gen_scramble = 1L,
  gen_S = 0L
)

Arguments

Value

J x J matrix of aggregate elasticities

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
eta <- get_halton_normals(50, d$N, ncol(d$W))
fit <- run_mxlogit(input_data = d, eta_draws = eta)
elas <- choicer:::mxl_elasticities_parallel(coef(fit), d$X, d$W, d$alt_idx,
  d$choice_idx, d$M, d$weights, eta, rc_dist = rep(0L, ncol(d$W)),
  elast_var_idx = 1L, is_random_coef = FALSE,
  rc_correlation = FALSE, rc_mean = FALSE)
elas

Analytical Hessian of the log-likelihood v2

Description

Computes the Hessian of the log-likelihood for the Mixed Logit model using OpenMP for parallelization. Mirrors the parameters of mxl_loglik_gradient_parallel.

Usage

mxl_hessian_parallel(
  theta,
  X,
  W,
  alt_idx,
  choice_idx,
  M,
  weights,
  eta_draws,
  rc_dist,
  rc_correlation = TRUE,
  rc_mean = FALSE,
  use_asc = TRUE,
  include_outside_option = FALSE,
  gen_seed = -1L,
  gen_scramble = 1L,
  gen_S = 0L
)

Arguments

Value

Hessian evaluated at input arguments

Note

For log-normal random coefficients (rc_dist=1) with rc_mean=TRUE, the distribution is a shifted log-normal: beta_k = exp(mu_k) + exp(L_k * eta), where exp(mu_k) shifts the location and exp(L_k * eta) ~ LogNormal(0, sigma_k^2). This differs from the textbook parameterization exp(mu_k + L_k * eta).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
eta <- get_halton_normals(50, d$N, ncol(d$W))
theta <- rep(0, ncol(d$X) + ncol(d$W) + nrow(d$alt_mapping) - 1)
H <- choicer:::mxl_hessian_parallel(theta, d$X, d$W, d$alt_idx, d$choice_idx,
  d$M, d$weights, eta, rc_dist = rep(0L, ncol(d$W)),
  rc_correlation = FALSE, rc_mean = FALSE)
dim(H)

Log-likelihood and gradient for Mixed Logit

Description

Computes the log-likelihood and its gradient for the Mixed Logit model using OpenMP for parallelization. Allows for inclusion of alternative-specific constants, outside option, observation weights, correlated random coefficients.

Usage

mxl_loglik_gradient_parallel(
  theta,
  X,
  W,
  alt_idx,
  choice_idx,
  M,
  weights,
  eta_draws,
  rc_dist,
  rc_correlation = TRUE,
  rc_mean = FALSE,
  use_asc = TRUE,
  include_outside_option = FALSE,
  gen_seed = -1L,
  gen_scramble = 1L,
  gen_S = 0L
)

Arguments

Value

List with loglikelihood and gradient evaluated at input arguments

Note

For log-normal random coefficients (rc_dist=1) with rc_mean=TRUE, the distribution is a shifted log-normal: beta_k = exp(mu_k) + exp(L_k * eta), where exp(mu_k) shifts the location and exp(L_k * eta) ~ LogNormal(0, sigma_k^2). This differs from the textbook parameterization exp(mu_k + L_k * eta).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
eta <- get_halton_normals(50, d$N, ncol(d$W))
K_x <- ncol(d$X); K_w <- ncol(d$W); J <- nrow(d$alt_mapping)
theta <- rep(0, K_x + K_w + J - 1)
result <- choicer:::mxl_loglik_gradient_parallel(theta, d$X, d$W, d$alt_idx,
  d$choice_idx, d$M, d$weights, eta, rc_dist = rep(0L, K_w),
  rc_correlation = FALSE, rc_mean = FALSE)
result$objective

Simulated expected logsum (inclusive value) for Mixed Logit

Description

Computes the simulated expected logsum (expected maximum utility, up to an additive constant) for each choice situation:

logsum_i = (1/S) \sum_s \log \sum_j \exp(V_{ij}^s),

where the inner sum runs over individual i's alternatives and includes the outside option's \exp(0) term when include_outside_option = TRUE. The log-sum-exp must be averaged across draws: applying log-sum-exp to the draw-averaged utilities returned by mxl_predict understates the expectation because log-sum-exp is convex (Jensen's inequality).

Usage

mxl_logsum(
  theta,
  X,
  W,
  alt_idx,
  M,
  eta_draws,
  rc_dist,
  rc_correlation = TRUE,
  rc_mean = FALSE,
  use_asc = TRUE,
  include_outside_option = FALSE,
  gen_seed = -1L,
  gen_scramble = 1L,
  gen_S = 0L
)

Arguments

Value

Vector of length N with the simulated expected logsum per choice situation.

Note

For log-normal random coefficients (rc_dist=1) with rc_mean=TRUE, the distribution is a shifted log-normal: beta_k = exp(mu_k) + exp(L_k * eta), where exp(mu_k) shifts the location and exp(L_k * eta) ~ LogNormal(0, sigma_k^2). This differs from the textbook parameterization exp(mu_k + L_k * eta).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
eta <- get_halton_normals(50, d$N, ncol(d$W))
fit <- run_mxlogit(input_data = d, eta_draws = eta)
ls <- choicer:::mxl_logsum(coef(fit), d$X, d$W, d$alt_idx, d$M, eta,
  rc_dist = rep(0L, ncol(d$W)), rc_correlation = FALSE, rc_mean = FALSE)
head(ls)

Per-observation simulated choice probabilities for Mixed Logit

Description

Returns the simulated choice probability for each (individual, alternative) row of X, averaged over the supplied Halton draws. Mirrors mnl_predict.

Usage

mxl_predict(
  theta,
  X,
  W,
  alt_idx,
  M,
  eta_draws,
  rc_dist,
  rc_correlation = TRUE,
  rc_mean = FALSE,
  use_asc = TRUE,
  include_outside_option = FALSE,
  gen_seed = -1L,
  gen_scramble = 1L,
  gen_S = 0L
)

Arguments

Value

List with choice_prob (length sum(M)), utility (length sum(M), simulated mean of the deterministic + W*gamma component), and, when include_outside_option = TRUE, choice_prob_outside (length N).

Predicted aggregate market shares for Mixed Logit

Description

Exported wrapper around the internal mxl_predict_shares_internal. Parses theta using the standard parameter ordering and returns the simulated weighted-average market shares.

Usage

mxl_predict_shares(
  theta,
  X,
  W,
  alt_idx,
  M,
  weights,
  eta_draws,
  rc_dist,
  rc_correlation = TRUE,
  rc_mean = FALSE,
  use_asc = TRUE,
  include_outside_option = FALSE,
  gen_seed = -1L,
  gen_scramble = 1L,
  gen_S = 0L
)

Arguments

Value

Vector of length J (or J+1 with outside option) of predicted shares.

Construct a choicer_sim object

Description

Wraps simulated data, true parameter values, and DGP settings into a classed list. Returned by simulate_mnl_data(), simulate_mxl_data(), and simulate_nl_data(), and consumed by recovery_table().

Usage

new_choicer_sim(data, true_params, settings, model)

Arguments

Value

A list of class choicer_sim.

BHHH/OPG information matrix for the Nested Logit model

Description

Computes the weighted outer product of per-individual scores \sum_i w_i\, s_i s_i^\top for the Nested Logit model. The per-individual score s_i (over the beta, lambda and delta/ASC blocks) is the (positive) gradient of individual i's log-likelihood contribution and is weight-free; the supplied weights enter only as the leading multiplier. Passing weights = w yields the ordinary weighted BHHH/OPG information; passing weights = w^2 yields the sandwich meat B = \sum_i w_i^2 s_i s_i^\top for robust (WESML) inference. Singleton-nest lambdas are fixed to 1 and contribute no score (mirroring the gradient kernel).

Usage

nl_bhhh_parallel(
  theta,
  X,
  alt_idx,
  choice_idx,
  nest_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

A symmetric positive-semidefinite information matrix \sum_i w_i\, s_i s_i^\top (same sign convention as the negated Hessian).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
K_x <- ncol(d$X); K_l <- length(unique(d$nest_idx))
theta <- c(rep(0, K_x), rep(0.5, K_l), rep(0, J - 1))
B <- choicer:::nl_bhhh_parallel(theta, d$X, d$alt_idx, d$choice_idx,
  d$nest_idx, d$M, d$weights)
dim(B)

BLP95 contraction mapping for the Nested Logit model

Description

Damped iterative fixed point recovering delta given target shares, using the NL probability structure. damping = 1 reproduces the plain BLP update.

Usage

nl_blp_contraction(
  delta,
  target_shares,
  X,
  beta,
  lambda,
  alt_idx,
  nest_idx,
  M,
  weights,
  include_outside_option = FALSE,
  damping = 1,
  tol = 1e-08,
  max_iter = 1000L
)

Arguments

Value

vector with contraction's delta (ASCs) output.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
beta <- coef(fit)[fit$param_map$beta]
lambda <- rep(1, length(unique(fit$data$nest_idx)))
lambda[as.integer(names(which(table(fit$data$nest_idx) > 1)))] <-
  coef(fit)[fit$param_map$lambda]
delta <- nl_blp_contraction(rep(0, J), rep(1/J, J), fit$data$X, beta, lambda,
  fit$data$alt_idx, fit$data$nest_idx, fit$data$M, fit$data$weights)
delta

Compute Nested Logit diversion ratios (parallelized over individuals)

Description

Computes the diversion ratio matrix DR(j->k) for the Nested Logit model. Entry (k, j) = fraction of demand lost by alternative j captured by k. Reduces to the MNL diversion ratios when all lambda = 1.

Usage

nl_diversion_ratios_parallel(
  theta,
  X,
  alt_idx,
  nest_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

J x J matrix where entry (k, j) = DR(j->k). Diagonal is 0.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
dr <- choicer:::nl_diversion_ratios_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
  fit$data$nest_idx, fit$data$M, fit$data$weights)
dr

Compute aggregate elasticities for the Nested Logit model

Description

Computes the aggregate (weighted-average) elasticity matrix for the Nested Logit model. Reduces to the MNL elasticities when all lambda = 1.

Usage

nl_elasticities_parallel(
  theta,
  X,
  alt_idx,
  choice_idx,
  nest_idx,
  M,
  weights,
  elast_var_idx,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

J x J matrix of aggregate elasticities (row = responding alt, col = perturbed alt).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
elas <- choicer:::nl_elasticities_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
  fit$data$choice_idx, fit$data$nest_idx, fit$data$M, fit$data$weights,
  elast_var_idx = 1L)
elas

Log-likelihood and gradient for Nested Logit model

Description

Computes the log-likelihood and its gradient for the Nested Logit model using OpenMP for parallelization. Especially handles singleton nests by fixing their lambda parameters to 1. Only non-singleton nests have a inclusive value coefficient estimated in theta.

Usage

nl_loglik_gradient_parallel(
  theta,
  X,
  alt_idx,
  choice_idx,
  nest_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

List with loglikelihood and gradient evaluated at input arguments

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
K_x <- ncol(d$X); K_l <- length(unique(d$nest_idx))
theta <- c(rep(0, K_x), rep(0.5, K_l), rep(0, J - 1))
result <- choicer:::nl_loglik_gradient_parallel(theta, d$X, d$alt_idx,
  d$choice_idx, d$nest_idx, d$M, d$weights)
result$objective

Analytical Hessian of the negated log-likelihood for the Nested Logit model

Description

Computes the exact (analytical) Hessian of the negated log-likelihood for the Nested Logit model using OpenMP parallelisation with thread-local accumulators. Covers all parameter blocks: beta-beta, beta-lambda, beta-delta, lambda-lambda, lambda-delta, and delta-delta. Singleton nests (lambda fixed to 1, not estimated) contribute no rows or columns to the lambda blocks.

Usage

nl_loglik_hessian_parallel(
  theta,
  X,
  alt_idx,
  choice_idx,
  nest_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

A symmetric (P x P) matrix: the Hessian of the negated log-likelihood evaluated at theta. Structurally identical to the output of nl_loglik_numeric_hessian; suitable for invert_hessian().

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
K_x <- ncol(d$X)
K_l <- sum(table(d$nest_idx) > 1)   # number of non-singleton nests (= 2)
theta <- c(rep(0, K_x), rep(0.8, K_l), rep(0, J - 1))
H <- choicer:::nl_loglik_hessian_parallel(theta, d$X, d$alt_idx, d$choice_idx,
  d$nest_idx, d$M, d$weights)
dim(H)

Numerical Hessian of the log-likelihood via finite differences

Description

Numerical Hessian of the log-likelihood via finite differences

Usage

nl_loglik_numeric_hessian(
  theta,
  X,
  alt_idx,
  choice_idx,
  nest_idx,
  M,
  weights,
  use_asc = TRUE,
  include_outside_option = FALSE,
  eps = 1e-06
)

Arguments

Value

Hessian evaluated at input arguments

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
K_x <- ncol(d$X); K_l <- length(unique(d$nest_idx))
theta <- c(rep(0, K_x), rep(0.5, K_l), rep(0, J - 1))
H <- choicer:::nl_loglik_numeric_hessian(theta, d$X, d$alt_idx, d$choice_idx,
  d$nest_idx, d$M, d$weights)
dim(H)

Prediction of choice probabilities and utilities for the Nested Logit model

Description

Prediction of choice probabilities and utilities for the Nested Logit model

Usage

nl_predict(
  theta,
  X,
  alt_idx,
  M,
  nest_idx,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

List with choice_prob (joint P_ij per stacked row) and utility (V_ij).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
pred <- choicer:::nl_predict(coef(fit), fit$data$X, fit$data$alt_idx, fit$data$M,
  fit$data$nest_idx, use_asc = TRUE)
head(pred$choice_prob)

Prediction of market shares for the Nested Logit model

Description

Prediction of market shares for the Nested Logit model

Usage

nl_predict_shares(
  theta,
  X,
  alt_idx,
  M,
  weights,
  nest_idx,
  use_asc = TRUE,
  include_outside_option = FALSE
)

Arguments

Value

vector with predicted market shares (outside-option share first when present).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
shares <- choicer:::nl_predict_shares(coef(fit), fit$data$X, fit$data$alt_idx,
  fit$data$M, fit$data$weights, fit$data$nest_idx, use_asc = TRUE)
shares

Extract number of observations from a choicer_fit object

Description

Extract number of observations from a choicer_fit object

Usage

## S3 method for class 'choicer_fit'
nobs(object, ...)

Arguments

Value

Integer number of choice situations.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
nobs(fit)

Number of choice situations behind a hierarchical Bayes fit

Description

Number of choice situations behind a hierarchical Bayes fit

Usage

## S3 method for class 'choicer_hb'
nobs(object, ...)

Arguments

Value

Integer count of choice situations (tasks).

Extract number of observations from a choicer_mnp object

Description

Extract number of observations from a choicer_mnp object

Usage

## S3 method for class 'choicer_mnp'
nobs(object, ...)

Arguments

Value

Integer number of choice situations.

Examples

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), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnprobit(dt, "id", "alt", "choice", c("x1", "x2"),
                    mcmc = list(R = 300, burn = 100))
nobs(fit)

Posterior-predictive share check for hierarchical Bayes fits

Description

Compares each alternative's observed take rate (share of choice situations in which it was chosen, including the outside option) with its posterior-predictive share from predict.choicer_hb(). Large systematic gaps indicate model misfit — e.g. a missing covariate or an outside-option share the delta level cannot rationalize.

Usage

ppc_shares(object, n_draws = 200L)

Arguments

Value

A data.table with columns alternative, observed, predicted, lower, upper (95% posterior-predictive interval), and covered (is the observed share inside the interval).

Examples

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",
                    mcmc = list(R = 500, burn = 200)))
ppc_shares(fit)

Posterior choice probabilities and shares for hierarchical Bayes fits

Description

Computes counterfactual choice probabilities, integrating over the posterior draws and (at the population level) over the random-coefficient distribution: for each kept draw (b_r, W_r, \delta_r, \ldots) one \beta \sim N(b_r, W_r) is drawn and the model probabilities are averaged. Alternatives in newdata that were not in the estimation sample receive a posterior-predictive \delta_{new} \sim N(z_{new}'\theta_r, \sigma_{d,r}^2) — the entry counterfactual unlocked by the random-effects \delta. Price or subsidy counterfactuals are just modified covariate columns in newdata.

Usage

## S3 method for class 'choicer_hb'
predict(
  object,
  newdata = NULL,
  level = c("population", "individual"),
  n_draws = 200L,
  aggregate = TRUE,
  ...
)

Arguments

Details

HMNL probabilities are closed-form logit; HMNP probabilities use the 1-D Gauss-Hermite representation of the iid-probit integral P(j) = \int \phi(u) \prod_{k \ne j} \Phi(V_j - V_k + u) du.

Value

With aggregate = TRUE, a data.table with columns alternative, share (posterior mean), sd, lower, upper (95% equal-tailed interval); the posterior share draws are attached as attr(, "draws"). With aggregate = FALSE, a numeric vector of posterior-mean choice probabilities, one per prediction row.

Examples

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)))
predict(fit)                       # posterior shares, estimation data
cf <- sim$data
cf$x1 <- cf$x1 + 0.5               # a counterfactual attribute change
predict(fit, newdata = cf)

Predict from a multinomial logit model

Description

Computes choice probabilities or aggregate market shares, either for the data used at fit time (default) or for counterfactual newdata.

Usage

## S3 method for class 'choicer_mnl'
predict(
  object,
  type = c("probabilities", "shares"),
  newdata = NULL,
  weights = NULL,
  ...
)

Arguments

Value

For "probabilities": a list with choice_prob and utility vectors. For "shares": a named numeric vector of market shares per alternative. With a data.frame newdata, rows are ordered by id, then by fit-time alternative code (alt_int in object$alt_mapping).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
predict(fit, type = "shares")
predict(fit, type = "probabilities")

# Counterfactual: increase x1 for alternative 2
dt_cf <- copy(dt)[alt == 2, x1 := x1 + 1]
predict(fit, type = "shares", newdata = dt_cf)

Predict from a mixed logit model

Description

Computes simulated choice probabilities or aggregate market shares using deterministic Halton draws, either for the data used at fit time (default) or for counterfactual newdata.

Usage

## S3 method for class 'choicer_mxl'
predict(
  object,
  type = c("probabilities", "shares"),
  newdata = NULL,
  weights = NULL,
  ...
)

Arguments

Value

For "probabilities": a list with choice_prob and utility vectors averaged across simulation draws. For "shares": a named numeric vector of simulated market shares per alternative. With a data.frame newdata, rows are ordered by id, then by fit-time alternative code (alt_int in object$alt_mapping).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mxlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = "x1", random_var_cols = "w1", S = 50L
)
predict(fit, type = "shares")
predict(fit, type = "probabilities")

Predict from a nested logit model

Description

Computes choice probabilities or aggregate market shares, either for the data used at fit time (default) or for counterfactual newdata.

Usage

## S3 method for class 'choicer_nl'
predict(
  object,
  type = c("probabilities", "shares"),
  newdata = NULL,
  weights = NULL,
  ...
)

Arguments

Value

For "probabilities": a list with choice_prob and utility vectors. For "shares": a named numeric vector of market shares per alternative. With a data.frame newdata, rows are ordered by id, then by fit-time alternative code (alt_int in object$alt_mapping).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, nest := rep(c(1L, 1L, 2L, 2L), N)]
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
predict(fit, type = "shares")
predict(fit, type = "probabilities")

Prepare inputs for hierarchical multinomial logit estimation

Description

Prepares and validates panel (or cross-sectional) choice data for the hierarchical Bayesian multinomial logit. The model has two random-effect levels: respondent-level structural tastes \beta_i \sim N(b, W) over the covariate_cols, and a global alternative-level effect \delta_j = z_j'\theta + \xi_j, \xi_j \sim N(0, \sigma_d^2), with mean-function design z_j built from alt_covariate_cols.

Usage

prepare_hmnl_data(
  data,
  id_col,
  alt_col,
  choice_col,
  covariate_cols,
  person_col = NULL,
  alt_covariate_cols = NULL,
  outside_opt_label = NULL,
  cf_residual_col = NULL,
  include_outside_option = TRUE,
  rc_dist = NULL
)

Arguments

Details

Structure. The design matrix X carries structural covariates only — no alternative-specific-constant dummies. The alternative effect \delta_j is indexed by alt_of_row (integer codes ⁠1..J⁠), so memory and compute scale with the number of rows, not with J extra design columns.

Outside option. With include_outside_option = TRUE (the default) the outside good is modelled implicitly, following the prepare_mnl_data() convention: physical outside rows (identified by outside_opt_label) are removed, the estimation kernels add the outside term (systematic utility 0), and a choice situation whose inside rows are all 0 in choice_col is coded as "outside chosen" (choice_pos = 0). The outside option anchors the location of \delta (mean utility relative to the outside good).

Cross-section vs panel. person_col groups choice situations into respondents sharing one \beta_i. With person_col = NULL (default) every choice situation is its own respondent (Ti all 1) — the cross-sectional random-coefficients mode.

Control function. cf_residual_col (a user-supplied first-stage residual, Petrin & Train 2010) is appended to X as an ordinary covariate; its provenance is recorded in data_spec. The first stage is NOT run here — supplying a valid residual is the user's responsibility.

Value

A list of class c("choicer_data_hmnl", "list") containing:

• X: Structural design matrix (total_rows x K_struct), no ASC columns; cf_residual_col last when supplied.

• alt_of_row: Integer alternative code per row (⁠1..J⁠).

• alt_idx: Alias of alt_of_row for the pooled-MLE init.

• Z: Alternative-level design (J x P), intercept first.

• M: Inside alternatives per choice situation.

• choice_pos: 1-based within-task position of the chosen row; 0 = outside option chosen.

• Ti: Choice situations per respondent.

• person_ids, N_persons, n_tasks, J, K_struct, P.

• include_outside_option: Logical flag.

• alt_mapping: Data.table mapping alternatives to summary statistics (outside option is alt_int = 0).

• param_map: Named list of index vectors (beta, theta), robust to collinearity drops.

• rc_dist: Integer vector aligned with the columns of X.

• dropped_cols, dropped_z_cols: Dropped column names, if any.

• data_spec: Column-name metadata (incl. person_col, outside_opt_label, cf_residual_col, alt_covariate_cols).

See Also

prepare_hmnp_data() for the hierarchical probit counterpart.

Examples

library(data.table)
set.seed(42)
N <- 20; T <- 3; J <- 4
dt <- data.table(
  pid  = rep(1:N, each = T * J),
  task = rep(seq_len(N * T), each = J),
  alt  = rep(1:J, N * T)
)
dt[, `:=`(x1 = rnorm(.N), x2 = runif(.N, -1, 1))]
dt[, quality := 0.1 * alt]                # alternative-level covariate
dt[, choice := 0L]
# leave some tasks all-zero: outside option chosen
dt[, choice := if (runif(1) < 0.8) sample(c(1L, rep(0L, J - 1))) else 0L,
   by = task]
input <- prepare_hmnl_data(dt, "task", "alt", "choice", c("x1", "x2"),
                           person_col = "pid",
                           alt_covariate_cols = "quality")
str(input$Z)
input$alt_mapping

Prepare inputs for hierarchical multinomial probit estimation

Description

Prepares and validates panel (or cross-sectional) choice data for the hierarchical Bayesian multinomial probit with iid N(0, \sigma^2) utility shocks. The model shares its two-level random-effect structure with prepare_hmnl_data(): respondent-level structural tastes \beta_i \sim N(b, W) over the covariate_cols (normal only — the probit keeps full conjugacy), and a global alternative-level effect \delta_j = z_j'\theta + \xi_j, \xi_j \sim N(0, \sigma_d^2).

Usage

prepare_hmnp_data(
  data,
  id_col,
  alt_col,
  choice_col,
  covariate_cols,
  person_col = NULL,
  alt_covariate_cols = NULL,
  outside_opt_label = NULL,
  cf_residual_col = NULL,
  include_outside_option = TRUE
)

Arguments

Details

The returned structure is identical to prepare_hmnl_data() (both preps share one internal engine), except there is no rc_dist field. Unlike prepare_mnp_data(), utilities are NOT differenced against a base alternative: the iid-shock model works in un-differenced utility space, so unbalanced choice sets are supported and the outside option is implicit (its latent utility is a stochastic N(0, \sigma^2) draw in the kernel, systematic utility 0).

Value

A list of class c("choicer_data_hmnp", "list") with the same components as prepare_hmnl_data() (minus rc_dist).

See Also

prepare_hmnl_data() for the component-by-component description.

Examples

library(data.table)
set.seed(42)
N <- 20; T <- 3; J <- 4
dt <- data.table(
  pid  = rep(1:N, each = T * J),
  task = rep(seq_len(N * T), each = J),
  alt  = rep(1:J, N * T)
)
dt[, `:=`(x1 = rnorm(.N), x2 = runif(.N, -1, 1))]
dt[, choice := 0L]
dt[, choice := if (runif(1) < 0.8) sample(c(1L, rep(0L, J - 1))) else 0L,
   by = task]
input <- prepare_hmnp_data(dt, "task", "alt", "choice", c("x1", "x2"),
                           person_col = "pid")
input$Ti[1:5]
input$alt_mapping

Prepare inputs for multinomial logit estimation

Description

Prepares and validates inputs for multinomial logit estimation routine.

Usage

prepare_mnl_data(
  data,
  id_col,
  alt_col,
  choice_col,
  covariate_cols,
  weights = NULL,
  outside_opt_label = NULL,
  include_outside_option = FALSE,
  weights_col = NULL,
  cluster_col = NULL
)

Arguments

Value

A list containing:

• X: Design matrix (sum(M) x K).

• alt_idx: Integer vector of alternative indices.

• choice_idx: Integer vector of chosen alternative indices.

• M: Integer vector with number of alternatives per choice situation.

• N: Number of choice situations.

• weights: Vector of weights.

• cluster: Vector of cluster labels (or NULL).

• situation_ids: Choice-situation ids in prepared (sorted) order.

• include_outside_option: Logical flag.

• alt_mapping: Data.table mapping alternatives to summary statistics.

• dropped_cols: Names of columns dropped due to collinearity, if any.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
input <- prepare_mnl_data(dt, "id", "alt", "choice", c("x1", "x2"))
str(input$X)
input$alt_mapping

Prepare inputs for Bayesian multinomial probit estimation

Description

Prepares and validates inputs for Bayesian multinomial probit estimation. Covariates are differenced against the base alternative, so the design matrix has one row per (choice situation, non-base alternative) pair. Balanced choice sets are required: every choice situation must contain the same J alternatives.

Usage

prepare_mnp_data(
  data,
  id_col,
  alt_col,
  choice_col,
  covariate_cols,
  base_alt = NULL,
  use_asc = TRUE
)

Arguments

Value

A list containing:

• X: Stacked differenced design matrix ((N * p) x K), covariate columns first, then ASC columns when use_asc = TRUE.

• y: Integer vector of choices (0 = base alternative, j in 1..p for the j-th non-base alternative), one per choice situation.

• p: Number of utility differences (J - 1).

• J: Number of alternatives.

• N: Number of choice situations.

• K: Number of columns of X.

• alt_mapping: Data.table mapping alternatives to summary statistics (the base alternative is alt_int = 1).

• base_alt: Resolved label of the base alternative.

• param_map: Named list of integer index vectors (beta, asc).

• use_asc: Logical flag.

• dropped_cols: Names of columns dropped due to collinearity, if any.

• data_spec: List with column name metadata.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
input <- prepare_mnp_data(dt, "id", "alt", "choice", c("x1", "x2"))
str(input$X)
input$alt_mapping

Prepare inputs for mixed logit estimation

Description

Prepares and validates inputs for mixed logit estimation routine.

Usage

prepare_mxl_data(
  data,
  id_col,
  alt_col,
  choice_col,
  covariate_cols,
  random_var_cols,
  weights = NULL,
  outside_opt_label = NULL,
  include_outside_option = FALSE,
  rc_correlation = FALSE,
  weights_col = NULL,
  cluster_col = NULL
)

Arguments

Value

A choicer_data_mxl object (list) containing:

• X: Fixed-coefficient design matrix (sum(M) x K_x).

• W: Random-coefficient design matrix (sum(M) x K_w).

• alt_idx: Integer vector of alternative indices.

• choice_idx: Integer vector of chosen alternative indices.

• M: Integer vector with number of alternatives per choice situation.

• N: Number of choice situations.

• weights: Vector of weights.

• cluster: Vector of cluster labels (or NULL).

• situation_ids: Choice-situation ids in prepared (sorted) order.

• include_outside_option: Logical flag.

• rc_correlation: Logical flag.

• alt_mapping: data.table mapping alternatives to summary statistics.

• dropped_cols: Names of columns dropped due to collinearity, if any.

• data_spec: List with column-name metadata.

Examples

library(data.table)
set.seed(42)
N <- 50; 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))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
input <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", c("w1", "w2"))
str(input$X)
str(input$W)

Prepare inputs for nested logit estimation

Description

Validates inputs, builds design matrices, and constructs nest structure for nested logit estimation. Calls prepare_mnl_data internally for base data preparation, then adds nest-specific fields.

Usage

prepare_nl_data(
  data,
  id_col,
  alt_col,
  choice_col,
  covariate_cols,
  nest_col,
  weights = NULL,
  outside_opt_label = NULL,
  include_outside_option = FALSE,
  weights_col = NULL,
  cluster_col = NULL
)

Arguments

Value

A choicer_data_nl object (list) containing:

• All fields from prepare_mnl_data (X, alt_idx, choice_idx, M, N, weights, cluster, situation_ids, include_outside_option, alt_mapping, dropped_cols).

• nest_idx: Integer vector of length J mapping each alternative (in alt_mapping row order) to its nest.

• data_spec: List with column name metadata including nest_col.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
input <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
input$nest_idx
input$alt_mapping

Print a consumer surplus summary

Description

Print a consumer surplus summary

Usage

## S3 method for class 'choicer_cs'
print(x, digits = 4, ...)

Arguments

Value

The object invisibly.

Print a choicer_fit object

Description

Prints a brief summary of the fitted model.

Usage

## S3 method for class 'choicer_fit'
print(x, ...)

Arguments

Value

The object invisibly.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
print(fit)

Print goodness-of-fit measures

Description

Print goodness-of-fit measures

Usage

## S3 method for class 'choicer_gof'
print(x, ...)

Arguments

Value

The object invisibly.

Print a hierarchical Bayes fit

Description

Print a hierarchical Bayes fit

Usage

## S3 method for class 'choicer_hb'
print(x, ...)

Arguments

Value

The object invisibly.

Examples

sim <- simulate_hmnl_data(N = 50, T = 2, J = 3, seed = 42)
fit <- suppressWarnings(run_hmnlogit(sim$data, "task", "alt", "choice", c("x1", "x2"),
                    person_col = "pid",
                    mcmc = list(R = 300, burn = 100)))
print(fit)

Print a choicer_mnp object

Description

Prints a brief summary of the fitted Bayesian multinomial probit model.

Usage

## S3 method for class 'choicer_mnp'
print(x, ...)

Arguments

Value

The object invisibly.

Examples

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), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnprobit(dt, "id", "alt", "choice", c("x1", "x2"),
                    mcmc = list(R = 300, burn = 100))
print(fit)

Print a WTP table

Description

Print a WTP table

Usage

## S3 method for class 'choicer_wtp'
print(x, digits = 4, ...)

Arguments

Value

The object invisibly.

Print the summary of a hierarchical Bayes fit

Description

Print the summary of a hierarchical Bayes fit

Usage

## S3 method for class 'summary.choicer_hb'
print(x, ...)

Arguments

Value

The object invisibly.

Print summary for multinomial logit model

Description

Print summary for multinomial logit model

Usage

## S3 method for class 'summary.choicer_mnl'
print(x, ...)

Arguments

Value

The object invisibly.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
print(summary(fit))

Print summary for Bayesian multinomial probit model

Description

Print summary for Bayesian multinomial probit model

Usage

## S3 method for class 'summary.choicer_mnp'
print(x, ...)

Arguments

Value

The object invisibly.

Examples

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), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnprobit(dt, "id", "alt", "choice", c("x1", "x2"),
                    mcmc = list(R = 300, burn = 100))
print(summary(fit))

Print summary for mixed logit model

Description

Print summary for mixed logit model

Usage

## S3 method for class 'summary.choicer_mxl'
print(x, ...)

Arguments

Value

The object invisibly.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mxlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = "x1", random_var_cols = "w1", S = 50L
)
print(summary(fit))

Print summary for nested logit model

Description

Print summary for nested logit model

Usage

## S3 method for class 'summary.choicer_nl'
print(x, ...)

Arguments

Value

The object invisibly.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = c("x1", "x2"), nest_col = "nest"
)
print(summary(fit))

Parameter recovery table

Description

Compares fitted coefficients to a set of true parameter values on the same scale as the estimator's internal parameterization. Returns one row per estimated parameter with true value, estimate, standard error, bias, relative bias (%), z-score against the truth, Wald CI, and a coverage indicator.

Usage

recovery_table(object, truth = NULL, level = 0.95, ...)

## S3 method for class 'choicer_fit'
recovery_table(object, truth = NULL, level = 0.95, ...)

## S3 method for class 'choicer_mnp'
recovery_table(object, truth = NULL, level = 0.95, ...)

## S3 method for class 'choicer_mc'
recovery_table(object, truth = NULL, level = 0.95, ...)

## S3 method for class 'choicer_hb'
recovery_table(object, truth = NULL, level = 0.95, ...)

Arguments

Details

For MXL fits the sigma block compares the raw Cholesky parameters (L_params), not the reconstructed covariance matrix. For log-normal random-coefficient means the raw mu estimate is compared directly; callers who want recovery on the DGP scale (exp(mu)) should transform both sides before calling.

When the estimator has normalized the first inside alternative's ASC to zero (which happens for MNL/MXL with include_outside_option = FALSE and no outside option baked into the fit), the first entry of truth$delta is dropped before the comparison so lengths match.

Value

See class-specific methods.

Methods (by class)

• recovery_table(choicer_fit): Returns a choicer_recovery object (a data.table) with columns parameter, group, true, estimate, se, bias, rel_bias_pct, z_vs_true, lower_ci, upper_ci, covers.

• recovery_table(choicer_mnp): Method for Bayesian MNP fits (choicer_mnp). The estimate column holds posterior means of the identified draws and se holds their posterior standard deviations, so lower_ci / upper_ci are normal-approximation credible intervals. In addition to the beta and asc blocks, a sigma block compares the identified covariance of the utility differences (lower triangle in the estimator's row-major Sigma_ij order); its first row is the sigma_11 = 1 normalization and is exact by construction. truth must be on the identified scale, as returned by simulate_mnp_data().

• recovery_table(choicer_mc): For a choicer_mc object, delegates to summary(object, level) and returns a choicer_mc_summary. Inspect object$replications directly for per-rep detail.

• recovery_table(choicer_hb): Method for hierarchical Bayes fits (choicer_hmnl / choicer_hmnp). estimate holds posterior means and se posterior standard deviations, so lower_ci / upper_ci are normal-approximation credible intervals. Blocks: beta (population means b vs truth$beta), w (diag(W) vs diag(truth$W)), theta (delta mean function), sigma_d (the SD, compared through the sqrt of the sigma_d^2 draws), and delta (per-alternative effects vs the realized truth$delta). For an HMNP fit, truth must be on the identified scale, as returned by simulate_hmnp_data().

Examples

sim <- simulate_mnl_data(N = 2000, J = 4, seed = 123)
fit <- run_mnlogit(
  data = sim$data, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = c("x1", "x2"),
  outside_opt_label = 0L, include_outside_option = FALSE, use_asc = TRUE
)
recovery_table(fit, sim)

Split-\widehat{R} convergence diagnostic

Description

Computes the split-\widehat{R} (potential scale reduction factor) of Gelman et al. for each column of a matrix of posterior draws. Every chain is split in half, so the diagnostic detects non-stationarity within a single chain as well as disagreement across chains; values near 1 indicate convergence, and values above roughly 1.05 warrant a longer run.

Usage

rhat(draws, rank = FALSE)

Arguments

Value

Named numeric vector with one \widehat{R} per parameter (NA for parameters with zero variance).

Examples

set.seed(42)
draws <- matrix(rnorm(2000), ncol = 2,
                dimnames = list(NULL, c("a", "b")))
rhat(draws)          # ~1: white noise is stationary
drifting <- cbind(a = cumsum(rnorm(1000)))
rhat(drifting)       # >> 1: a random walk is not
rhat(draws, rank = TRUE)   # rank-normalized variant

Fit a hierarchical Bayesian multinomial logit (HMNL)

Description

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

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".

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.

See Also

prepare_hmnl_data(), simulate_hmnl_data(), recovery_table(), rhat(), ess(), mcse(), traceplot()

Examples

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)))
summary(fit)
coef(fit, component = "delta")

Fit a hierarchical Bayesian multinomial probit (HMNP)

Description

Runs the fully conjugate Albert-Chib Gibbs sampler for the hierarchical multinomial probit with iid normal utility shocks in un-differenced utility space:

U_{ijt} = x_{ijt}'\beta_i + \delta_j + \epsilon_{ijt}, \qquad U_{iot} = \epsilon_{iot}, \qquad \epsilon \sim N(0, \sigma^2),

choice by argmax within the task including the stochastic implicit outside option, \beta_i \sim N(b, W) (normal coordinates only — log-normal would break conjugacy), and \delta_j = z_j'\theta + \xi_j, \xi_j \sim N(0, \sigma_d^2).

Usage

run_hmnprobit(
  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,
  prior = list(),
  mcmc = list(),
  chains = 1,
  keep_beta_i = c("means", "draws", "none"),
  keep_data = TRUE
)

Arguments

Details

Identification. The probit likelihood is invariant to a common rescaling of utilities and \sigma, so the chain runs on the non-identified parameterization (free \sigma^2, better mixing via parameter expansion) and every kept draw is normalized by the matching power of the CURRENT \sigma: reported b/\sigma, W/\sigma^2, \delta/\sigma, \theta/\sigma, \sigma_d^2/\sigma^2. Raw chains are kept in ⁠draws$*_raw⁠. The outside option anchors the location of \delta exactly as in run_hmnlogit().

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".

Value

A choicer_hmnp object (classed c("choicer_hmnp", "choicer_hb")); the same layout as run_hmnlogit()'s return, with all reported summaries on the identified scale, raw chains in ⁠draws$*_raw⁠, the non-identified draws$sigma2 trace, and all chains' retained draws (identified scale) in chains.

See Also

prepare_hmnp_data(), simulate_hmnp_data(), run_hmnlogit(), ess(), mcse(), traceplot()

Examples

sim <- simulate_hmnp_data(N = 100, T = 3, J = 4, seed = 42)
fit <- suppressWarnings(run_hmnprobit(sim$data, "task", "alt", "choice", c("x1", "x2"),
                     person_col = "pid", alt_covariate_cols = "z1",
                     mcmc = list(R = 500, burn = 200)))
summary(fit)

Runs multinomial logit estimation

Description

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("none", "sd", "mad", "iqr"),
  se_method = c("hessian", "bhhh", "sandwich", "cluster"),
  cluster_col = NULL,
  nloptr_opts = NULL
)

Arguments

Details

Two workflows are supported:

Convenience (default)

Supply data and column names. Data preparation (prepare_mnl_data) is handled automatically.

Advanced

Call prepare_mnl_data yourself and pass the result via input_data.

Value

A choicer_mnl object (inherits from choicer_fit). Standard S3 methods available: summary(), coef(), vcov(), logLik(), AIC(), BIC(), nobs(), predict().

Examples

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))]
dt[, V := drop(as.matrix(.SD) %*% beta_true), .SDcols = c("x1","x2")]
dt[, prob := exp(V) / sum(exp(V)), by = id]
dt[, choice := as.integer(alt == sample(alt, 1, prob = prob)), by = id]

fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
summary(fit)
coef(fit)
AIC(fit)
predict(fit, type = "shares")

Runs Bayesian multinomial probit estimation

Description

Estimates a multinomial probit model by Gibbs sampling with data augmentation (Albert & Chib 1993; McCulloch & Rossi 1994). The model is specified in utility differences against a base alternative: for choice situation i with J alternatives, w_i = X_i \beta + \epsilon_i with \epsilon_i \sim N_{J-1}(0, \Sigma).

Usage

run_mnprobit(
  data = NULL,
  id_col = NULL,
  alt_col = NULL,
  choice_col = NULL,
  covariate_cols = NULL,
  input_data = NULL,
  base_alt = NULL,
  use_asc = TRUE,
  prior = list(),
  mcmc = list(),
  keep_data = TRUE
)

Arguments

Details

Two workflows are supported:

Convenience (default)

Supply data and column names. Data preparation (prepare_mnp_data) is handled automatically.

Advanced

Call prepare_mnp_data yourself and pass the result via input_data.

Identification. The multinomial probit likelihood is invariant to a common rescaling (\beta, \Sigma) \to (c\beta, c^2\Sigma). The sampler runs on the non-identified parameterization (unrestricted \Sigma with an inverse-Wishart prior) and identified quantities are computed by normalizing each kept draw by \sigma_{11}: \beta / \sqrt{\sigma_{11}} and \Sigma / \sigma_{11}. This is the McCulloch & Rossi (1994) default, which keeps all Gibbs conditionals conjugate and mixes better than the fully identified sampler of McCulloch, Polson & Rossi (2000). Reported coefficients, standard deviations, and credible intervals are posterior summaries of the identified draws.

Reproducibility. The sampler uses its own thread-safe RNG with one stream per (iteration, observation), so results are reproducible independent of the number of OpenMP threads (see set_num_threads()). When mcmc$seed is not supplied, a master seed is drawn from R's RNG, so set.seed() controls the run.

Scope. Balanced choice sets are required: every choice situation must contain the same J alternatives. To model an outside option, include it as explicit rows with zero covariates and set base_alt to its label.

Value

A choicer_mnp object. S3 methods available: summary(), coef() (posterior means of identified coefficients), vcov() (posterior covariance of identified coefficient draws), nobs(). Posterior draws are stored in $draws (beta / sigma on the identified scale, beta_raw / sigma_raw unnormalized).

References

Albert, J. H., & Chib, S. (1993). Bayesian Analysis of Binary and Polychotomous Response Data. Journal of the American Statistical Association, 88(422), 669-679.

McCulloch, R., & Rossi, P. E. (1994). An exact likelihood analysis of the multinomial probit model. Journal of Econometrics, 64(1-2), 207-240.

Examples

library(data.table)
set.seed(42)
N <- 200; 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))]
dt[, U := drop(as.matrix(.SD) %*% beta_true) + rnorm(.N), .SDcols = c("x1", "x2")]
dt[, choice := as.integer(U == max(U)), by = id]

fit <- run_mnprobit(dt, "id", "alt", "choice", c("x1", "x2"),
                    mcmc = list(R = 500, burn = 100))
summary(fit)
coef(fit)

Runs mixed logit estimation

Description

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("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
)

Arguments

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.

Value

A choicer_mxl object (inherits from choicer_fit). Standard S3 methods available: summary(), coef(), vcov(), logLik(), AIC(), BIC(), nobs().

Examples

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))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]

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
)
summary(fit)

Runs nested logit estimation

Description

Estimates a nested logit model via maximum likelihood.

Usage

run_nestlogit(
  data = NULL,
  id_col = NULL,
  alt_col = NULL,
  choice_col = NULL,
  covariate_cols = NULL,
  nest_col = NULL,
  input_data = NULL,
  use_asc = TRUE,
  theta_init = NULL,
  param_names = NULL,
  optimizer = NULL,
  control = list(),
  weights = NULL,
  weights_col = NULL,
  outside_opt_label = NULL,
  include_outside_option = FALSE,
  keep_data = TRUE,
  se_method = c("hessian", "numeric", "bhhh", "sandwich", "cluster"),
  cluster_col = NULL,
  nloptr_opts = NULL
)

Arguments

Details

Two workflows are supported:

Convenience

Supply data and column names (including nest_col). Data preparation (prepare_nl_data) is handled automatically.

Advanced

Call prepare_nl_data (or build the input list manually) and pass it via input_data.

Value

A choicer_nl object (inherits from choicer_fit). Standard S3 methods available: summary(), coef(), vcov(), logLik(), AIC(), BIC(), nobs().

Examples

library(data.table)
set.seed(42)
N <- 100; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]

fit <- run_nestlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = c("x1", "x2"), nest_col = "nest"
)
summary(fit)

Draw a choice-based sample stratified by the chosen alternative

Description

Subsamples whole choice situations from a population data set according to fixed per-stratum quotas, where strata are defined by the chosen alternative. The input data set is treated as the population, so the population shares Q(j) are known exactly; the returned sample carries a ready-to-use WESML weight column (see wesml_weights).

Usage

sample_by_choice(
  data,
  id_col,
  alt_col,
  choice_col,
  n_per_alt = NULL,
  frac_per_alt = NULL,
  seed = NULL,
  weight_name = ".wesml_weight",
  outside_opt_label = NULL,
  include_outside_option = FALSE
)

Arguments

Details

Sampling is by choice situation (id), never by row: all alternative-rows of a sampled situation are kept together. Sampling is without replacement.

Value

A data.table subsample with the weight column appended and "Q", "H", and "choice_sampling" attributes (the last records the scheme, shares, quotas, and meat = "robust").

References

Manski, C. F. and Lerman, S. R. (1977). Econometrica 45(8), 1977-1988.

See Also

wesml_weights, run_mxlogit

Examples

library(data.table)
set.seed(1)
N <- 600L; J <- 3L
pop <- data.table(id = rep(seq_len(N), each = J), alt = rep(1:J, N))
pop[, x1 := rnorm(.N)]
pop[, w1 := rnorm(.N)]
pop[, choice := as.integer(seq_len(.N) == sample.int(.N, 1L)), by = id]

s <- sample_by_choice(pop, "id", "alt", "choice", n_per_alt = 50L, seed = 1L)
attr(s, "choice_sampling")$H   # realized sample shares
head(s[[".wesml_weight"]])

Set the number of OpenMP threads used by choicer

Description

Set the number of OpenMP threads used by choicer

Usage

set_num_threads(n_threads)

Arguments

Value

Invisibly returns NULL.

Simulate hierarchical multinomial logit data

Description

Generates synthetic panel choice data from the hierarchical (random coefficients + alternative-level random effects) logit DGP: respondents ⁠i = 1..N⁠ face T choice situations each, with utilities

U_{ijt} = x_{ijt}'\gamma_i + \delta_j + \epsilon_{ijt}, \qquad U_{iot} = \epsilon_{iot},

i.i.d. Gumbel shocks (including a shock on the outside option, whose systematic utility is 0), \beta_i \sim N(\beta, W) with \gamma_{ik} = \beta_{ik} or \exp(\beta_{ik}) per rc_dist, and \delta_j = z_j'\theta + \xi_j with \xi_j \sim N(0, \sigma_d^2). Covariates are Uniform(-1, 1); the alternative-level covariates ⁠z*⁠ are constant within each alternative.

Usage

simulate_hmnl_data(
  N = 500,
  T = 10,
  J = 4,
  beta = c(0.8, -0.6),
  W = NULL,
  theta = c(0.5, -0.4),
  sigma_d = 0.5,
  Z = NULL,
  rc_dist = NULL,
  include_outside = TRUE,
  seed = 123,
  vary_choice_set = FALSE
)

Arguments

Details

Log-normal coordinates are reported on the chain (log) scale in true_params$beta — the scale on which the estimator's hierarchy operates — while entering utility as exp(beta_ik).

Value

A choicer_sim object. true_params contains beta, W, theta, sigma_d, the realized delta and xi vectors, the full mean-function design Z (intercept first), and rc_dist.

Examples

sim <- simulate_hmnl_data(N = 100, T = 4, J = 4, seed = 123)
print(sim)
sim$true_params$delta

Simulate hierarchical multinomial probit data

Description

Generates synthetic panel choice data from the hierarchical probit DGP with iid normal utility shocks:

U_{ijt} = x_{ijt}'\beta_i + \delta_j + \epsilon_{ijt}, \qquad U_{iot} = \epsilon_{iot}, \qquad \epsilon \sim N(0, \sigma^2),

choice by argmax within the task. The outside option is stochastic — it carries its own N(0, \sigma^2) shock on top of systematic utility 0, exactly as in the estimator. \beta_i \sim N(\beta, W) (normal only) and \delta_j = z_j'\theta + \xi_j, \xi_j \sim N(0, \sigma_d^2), as in simulate_hmnl_data().

Usage

simulate_hmnp_data(
  N = 500,
  T = 10,
  J = 4,
  beta = c(0.8, -0.6),
  W = NULL,
  theta = c(0.5, -0.4),
  sigma_d = 0.5,
  Z = NULL,
  include_outside = TRUE,
  seed = 123,
  vary_choice_set = FALSE,
  sigma = 1
)

Arguments

Details

The iid-probit likelihood identifies parameters only up to the common scale \sigma, so true_params is reported on the identified scale: beta = \beta/\sigma, W = W/\sigma^2, theta = \theta/\sigma, sigma_d = \sigma_d/\sigma, delta = \delta/\sigma, xi = \xi/\sigma. With the default sigma = 1 the DGP scale and the identified scale coincide.

Value

A choicer_sim object. true_params contains beta, W, theta, sigma_d, the realized delta and xi, and the full mean-function design Z — all on the identified scale (see Details).

Examples

sim <- simulate_hmnp_data(N = 100, T = 4, J = 4, seed = 123)
print(sim)
sim$true_params$delta

Simulate multinomial logit data

Description

Generates synthetic choice data with i.i.d. Gumbel errors, optionally with varying choice-set sizes and an outside option (alt = 0). Choices are determined by argmax of utility; covariates are drawn as Uniform(-1, 1).

Usage

simulate_mnl_data(
  N = 5000,
  J = 5,
  beta = c(0.8, -0.6),
  delta = NULL,
  seed = 123,
  outside_option = TRUE,
  vary_choice_set = TRUE
)

Arguments

Value

A choicer_sim object.

Examples

sim <- simulate_mnl_data(N = 1000, J = 5, seed = 123)
print(sim)

Simulate multinomial probit data

Description

Generates synthetic choice data from the MNP data-generating process estimated by run_mnprobit(): latent utility differences against the base alternative (alternative 1),

w_i = X_i \beta + \delta + \varepsilon_i, \qquad \varepsilon_i \sim N_{J-1}(0, \Sigma),

with alternative j > 1 chosen iff w_{ij} > \max(0, \max_{k \neq j} w_{ik}) and the base chosen iff all w_{ij} < 0. Covariates are Uniform(-1, 1). Choice sets are balanced (every individual faces all J alternatives), as the MNP estimator requires; there is no outside-option flag — model an outside good as a zero-covariate base alternative instead.

Usage

simulate_mnp_data(
  N = 5000,
  J = 3,
  beta = c(0.8, -0.6),
  delta = NULL,
  Sigma = matrix(c(1, 0.5, 0.5, 1.5), nrow = 2),
  seed = 123
)

Arguments

Details

The MNP likelihood only identifies parameters up to scale, so true_params is reported on the identified scale (normalized by \sigma_{11}): beta = \beta / \sqrt{\sigma_{11}}, delta = \delta / \sqrt{\sigma_{11}}, and Sigma = \Sigma / \sigma_{11} — the scale on which run_mnprobit() reports its posterior. With the default Sigma (\sigma_{11} = 1) the DGP scale and the identified scale coincide.

Value

A choicer_sim object. true_params contains beta, delta, and Sigma on the identified scale (see Details).

Examples

sim <- simulate_mnp_data(N = 1000, J = 3, seed = 123)
print(sim)

Simulate mixed logit data

Description

Generates synthetic choice data with random coefficients drawn from a multivariate normal (optionally log-normal per dimension) and an additional mean shifter mu. Random coefficients are parameterized via the lower Cholesky factor of Sigma. Covariates are Uniform(-1, 1) by default; columns named in price_cols are drawn as -Uniform(0.1, 3) to mimic strictly-negative price variables.

Usage

simulate_mxl_data(
  N = 5000,
  J = 4,
  beta = c(0.8, -0.6),
  delta = NULL,
  mu = NULL,
  Sigma = matrix(c(1, 0.5, 0.5, 1.5), nrow = 2),
  rc_dist = NULL,
  rc_correlation = NULL,
  price_cols = NULL,
  seed = 123,
  outside_option = TRUE,
  vary_choice_set = TRUE
)

Arguments

Details

Random coefficients are constructed to match the estimator's parameterization in src/mxlogit.cpp. For every dimension the raw draw is L %*% eta where eta ~ N(0, I). A normal random coefficient (rc_dist = 0) is then ⁠gamma_k = mu_k + (L %*% eta)_k⁠. A log-normal random coefficient (rc_dist = 1) follows the shifted log-normal ⁠beta_k = exp(mu_k) + exp((L %*% eta)_k)⁠ – not the textbook exp(mu_k + sigma_k * eta) – so mu_k in true_params$mu is on the same scale the estimator recovers and recovery_table() can compare like-for-like.

Value

A choicer_sim object. true_params includes beta, delta, Sigma, L_params (packed Cholesky parameters), mu, rc_dist, rc_correlation.

Examples

sim <- simulate_mxl_data(N = 1000, J = 4, seed = 123)
print(sim)

Simulate nested logit data

Description

Generates synthetic choice data with nested logit probabilities computed analytically (log-sum-exp over inclusive values), then samples choices from the implied multinomial. The outside option (j = 0) sits in a singleton nest with lambda = 1.

Usage

simulate_nl_data(
  N = 10000,
  beta = c(1.5, -0.8),
  delta = c(`1` = 0.5, `2` = 0.3, `3` = -0.2, `4` = -0.5, `5` = 0.4),
  nests = list(c(1, 2), c(3, 4, 5)),
  lambdas = c(0.8, 0.2),
  seed = 123
)

Arguments

Value

A choicer_sim object. true_params includes beta, delta, lambdas; settings includes the nest_structure. The returned data retains a nest column (integer, with 0L for the outside option) for convenient use with run_nestlogit().

Note

Unlike simulate_mnl_data() and simulate_mxl_data(), this function does not expose outside_option or vary_choice_set flags. The outside option (j = 0) is always present as a singleton nest with lambda = 1, and every individual faces the full set of inside alternatives. Add these flags if downstream use cases need them.

Examples

sim <- simulate_nl_data(N = 2000, seed = 123)
print(sim)

Summarize a hierarchical Bayes fit

Description

Posterior summaries (mean, SD, equal-tailed credible interval) for the population coefficients b, the mean-function coefficients \theta, the alternative-effect variance \sigma_d^2 (and, for the HMNP, the raw shock variance trace), plus the \delta_j / \xi_j quality ladder, acceptance diagnostics, and a consolidated convergence-diagnostic table (rank-normalized R-hat, ESS bulk/tail, MCSE) built from all retained chains (see rhat(), ess(), mcse()).

Usage

## S3 method for class 'choicer_hb'
summary(object, prob = 0.95, ...)

Arguments

Value

A summary.choicer_hb object.

Examples

sim <- simulate_hmnl_data(N = 50, T = 2, J = 3, seed = 42)
fit <- suppressWarnings(run_hmnlogit(sim$data, "task", "alt", "choice", c("x1", "x2"),
                    person_col = "pid",
                    mcmc = list(R = 300, burn = 100)))
summary(fit)

Summary for multinomial logit model

Description

Computes and returns a coefficient summary table with standard errors, z-values, p-values, and significance codes. Triggers lazy Hessian computation if standard errors have not been computed yet.

Usage

## S3 method for class 'choicer_mnl'
summary(object, gof = TRUE, ...)

Arguments

Value

A summary.choicer_mnl object (list with coefficients table and metadata, including a gof element with goodness-of-fit measures from gof; its fields are NA when the model was fitted with keep_data = FALSE).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
summary(fit)

Summary for Bayesian multinomial probit model

Description

Posterior summaries (mean, SD, equal-tailed credible interval) of the identified coefficient and covariance draws.

Usage

## S3 method for class 'choicer_mnp'
summary(object, prob = 0.95, ...)

Arguments

Value

A summary.choicer_mnp object (list with coefficient and Sigma posterior tables plus metadata).

Examples

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), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnprobit(dt, "id", "alt", "choice", c("x1", "x2"),
                    mcmc = list(R = 300, burn = 100))
summary(fit)

Summary for mixed logit model

Description

Computes coefficient summary with delta-method transformation for variance parameters (Cholesky to covariance scale) and log-normal mean parameters. Triggers lazy Hessian computation if standard errors have not been computed yet.

Usage

## S3 method for class 'choicer_mxl'
summary(object, gof = TRUE, ...)

Arguments

Value

A summary.choicer_mxl object (includes a gof element with goodness-of-fit measures from gof; its fields are NA when the model was fitted with keep_data = FALSE).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mxlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = "x1", random_var_cols = "w1", S = 50L
)
summary(fit)

Summary for nested logit model

Description

Triggers lazy Hessian computation if standard errors have not been computed yet.

Usage

## S3 method for class 'choicer_nl'
summary(object, gof = TRUE, ...)

Arguments

Value

A summary.choicer_nl object (includes a gof element with goodness-of-fit measures from gof; its fields are NA when the model was fitted with keep_data = FALSE).

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 4
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, nest := ifelse(alt <= 2, "A", "B")]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_nestlogit(
  data = dt, id_col = "id", alt_col = "alt", choice_col = "choice",
  covariate_cols = c("x1", "x2"), nest_col = "nest"
)
summary(fit)

Query choicer OpenMP thread settings

Description

Query choicer OpenMP thread settings

Usage

thread_info()

Value

A list with OpenMP availability, active/max thread settings, CPU thread capacity reported by OpenMP, thread limits, and relevant environment variables.

Traceplot for a hierarchical Bayes fit

Description

Generic dispatching on the fit's class. See traceplot.choicer_hb() for the choicer_hmnl / choicer_hmnp method.

Usage

traceplot(object, ...)

Arguments

Value

The object, invisibly.

Traceplot method for hierarchical Bayes fits

Description

Overlaid per-chain traceplots (base graphics, no new dependency) for a choicer_hmnl / choicer_hmnp fit's population coefficients (b), delta mean-function coefficients (theta), alternative-effect variance (sigma_d2), and, opt-in, a representative subset of the (potentially ~200-column) alternative effects (delta).

Usage

## S3 method for class 'choicer_hb'
traceplot(object, block = c("b", "theta", "sigma_d2"), which = NULL, ...)

Arguments

Value

object, invisibly.

Examples

sim <- simulate_hmnl_data(N = 60, T = 2, J = 4, seed = 1)
fit <- suppressWarnings(run_hmnlogit(sim$data, "task", "alt", "choice", c("x1", "x2"),
                    person_col = "pid",
                    mcmc = list(R = 300, burn = 100), chains = 2))
traceplot(fit)

Extract variance-covariance matrix from a choicer_fit object

Description

With no arguments, returns the variance-covariance matrix implied by the fit's own se_method (triggering lazy computation if needed). Passing type recomputes a different variance estimator post hoc from the stored data — no refit needed (requires keep_data = TRUE):

"hessian"

Inverse of the analytical negated Hessian.

"bhhh"

Inverse of the BHHH/OPG information \sum_i w_i s_i s_i'.

"robust"

Huber-White sandwich A^{-1} (\sum_i w_i^2 s_i s_i') A^{-1} — also the valid WESML variance under choice-based weighting.

"cluster"

Cluster-robust sandwich A^{-1} (\sum_g g_g g_g') A^{-1} with g_g = \sum_{i \in g} w_i s_i the within-cluster sum of weighted scores. Requires cluster (or a fit made with cluster_col). No small-sample correction is applied.

Here i indexes choice situations. For repeated choices by the same decision maker (panel data), cluster on the decision maker.

Usage

## S3 method for class 'choicer_fit'
vcov(object, type = NULL, cluster = NULL, ...)

Arguments

Details

Note (mixed logit): clustering repairs the inference, not the estimand. run_mxlogit() treats each choice situation as an independent draw from the mixing distribution (a cross-sectional MSL likelihood, not the panel product form), so on panel data the point estimates target that cross-sectional model; type = "cluster" makes their standard errors robust to within-person dependence but does not turn the fit into a panel mixed logit. For panel random coefficients use run_hmnlogit (person_col).

Value

Named variance-covariance matrix, or NULL if unavailable.

Examples

library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
dt[, person := rep(1:10, each = 5)[id]]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
vcov(fit)                          # as fitted (hessian)
vcov(fit, type = "robust")         # Huber-White, post hoc
# named by situation id -> safe regardless of order
cl <- dt[, person[1L], by = id]
vcov(fit, type = "cluster", cluster = setNames(cl$V1, cl$id))

Posterior covariance of the population coefficients

Description

Posterior covariance of the population coefficients

Usage

## S3 method for class 'choicer_hb'
vcov(object, ...)

Arguments

Value

K x K posterior covariance matrix of the b draws.

Extract variance-covariance matrix from a choicer_mnp object

Description

Returns the posterior covariance matrix of the identified coefficient draws (computed eagerly at fit time; no Hessian is involved).

Usage

## S3 method for class 'choicer_mnp'
vcov(object, ...)

Arguments

Value

Named posterior covariance matrix.

Examples

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), x2 = rnorm(.N))]
dt[, choice := 0L]
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
fit <- run_mnprobit(dt, "id", "alt", "choice", c("x1", "x2"),
                    mcmc = list(R = 300, burn = 100))
vcov(fit)

Robust (sandwich) variance for a weighted / choice-based logit fit

Description

Recomputes the robust Huber-White sandwich variance V = A^{-1} B A^{-1} for a fitted multinomial (MNL), mixed (MXL) or nested (NL) logit, where the bread A = \sum_i w_i (-H_i) is the weighted negated Hessian and the meat B = \sum_i w_i^2 s_i s_i' is the weight-squared outer product of the per-individual scores. This is the appropriate variance under choice-based (endogenous stratified) / WESML weighting, where the inverse-Hessian and the ordinary BHHH variance are invalid. It can be called on any fitted model (e.g. one estimated with se_method = "hessian") to obtain robust standard errors post hoc, without refitting.

Usage

wesml_vcov(object, ...)

## S3 method for class 'choicer_mxl'
wesml_vcov(object, type = c("vcov", "se"), ...)

## S3 method for class 'choicer_mnl'
wesml_vcov(object, type = c("vcov", "se"), ...)

## S3 method for class 'choicer_nl'
wesml_vcov(object, type = c("vcov", "se"), ...)

Arguments