Package 'choicer'

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

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

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

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,
  keep_draws_bytes = 2^32,
  ...
)

Arguments

Details

The contraction integrates each choice situation over the same draws as predict() (for a cross-sectional fit, the estimation draws; a panel fit's estimation draws are per decision maker), formed on the fly, in parallel, and kept across iterations within keep_draws_bytes. With stored draws (draws = "store") they are the points of the cube get_halton_normals would build, formed without it; with more than 128 random coefficients, the cube itself, read at every iteration (and nothing more is kept).

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

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

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)

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

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)

Description

For each decision maker, the mean and standard deviation of the random coefficients conditional on the choices they were observed to make, $E[\beta_n \mid y_n, x_n; \hat\theta]$ (Revelt and Train 2000; Train 2009, ch. 11): the frequentist counterpart of the individual-level beta_i summaries of run_hmnlogit.

Usage

conditional_tastes(object, ...)

## Default S3 method:
conditional_tastes(object, ...)

## S3 method for class 'choicer_mxl'
conditional_tastes(object, ...)

Arguments

Details

Estimand. The fit estimates the population distribution of tastes, $f(\beta \mid \theta)$. Among decision makers who face $x_n$ and make the choices $y_n$, tastes follow the conditional distribution

$h(\beta \mid y_n, x_n; \theta) = \frac{P(y_n \mid x_n, \beta) \, f(\beta \mid \theta)} {\int P(y_n \mid x_n, \beta) \, f(\beta \mid \theta) \, d\beta},$

where $P(y_n \mid x_n, \beta) = \prod_t P_{nt}(j_{nt} \mid \beta)$ is the product of the decision maker's logit probabilities. Its mean and standard deviation are simulated at $\hat\theta$ with the estimation draws, as weighted averages over draws $\beta_{ns}$ with weights $\omega_{ns} \propto P(y_n \mid x_n, \beta_{ns})$. Coefficients are on the scale they enter utility: $\mu_k + (L\eta)_k$ for a normal coefficient and $\exp(\mu_k) + \exp((L\eta)_k)$ for the shifted log-normal, with the $\mu$ term present only when rc_mean = TRUE. If a variable's mean is carried instead by the same variable in covariate_cols (with rc_mean = FALSE), only its random part is reported. This is the utility scale; the beta_i summaries of run_hmnlogit() are on the scale of the underlying normal, i.e. the log of the coefficient for log-normal coordinates.

What it is not. The conditional mean is not the decision maker's own $\beta_n$. It is the average taste among people who would make $n$'s choices, so it shrinks toward the population mean, the more so the fewer choice situations $n$ contributes, and its dispersion across decision makers understates the population dispersion. By the law of iterated expectations the conditional means average, in large samples, to the population mean (mean_cond vs model_mean in population), and by the law of total variance $Var(\beta) = Var(E[\beta \mid y]) + E[Var(\beta \mid y)]$: the share of taste variance revealed by choices (revealed) and the share left unresolved (residual) should sum to about one. These are numerical checks, not specification tests: for normal random coefficients with free means, both identities are the likelihood's first-order conditions for $\mu$ and $\Sigma$, so they hold at a converged maximum whatever the true taste distribution. A departure points to too few draws S or an unconverged $\hat\theta$; when a random coefficient's mean is fixed at zero (rc_mean = FALSE and no fixed coefficient on the same variable), a mean_cond far from zero can also flag that restriction. For log-normal coefficients the variance identity holds on the underlying normal scale, not the utility scale reported here, so expect departures of a few percent in either direction at moderate S. The sd is the dispersion of the conditional distribution of tastes, not a standard error: it ignores the sampling error in $\hat\theta$.

Cross-section. Without person_col each choice situation is its own decision maker. A single choice reveals little about tastes, so the conditional means stay close to the population mean (revealed small); fit with person_col to learn about tastes from repeated choices.

Choice probabilities are handled in log space, so every unit has a finite conditional distribution however improbable its choices. Only a unit whose utilities overflow at the estimates (non-finite log-likelihood) gets NA, and it is left out of the population summaries.

Value

An object of class choicer_tastes: a list with

mean, sd

$K_w \times U$ matrices of conditional means and standard deviations; rows are the random coefficients (random_var_cols), columns the $U$ likelihood units, named by decision-maker id (person_col) or by choice-situation id.

unit

"decision maker" or "choice situation".

n_units

$U$.

S

Simulation draws per unit.

weights

Unit weights used in population.

population

A data.frame with one row per random coefficient: model_mean and model_sd, the moments of $f(\beta \mid \hat\theta)$; mean_cond and sd_cond, the weighted mean and standard deviation of the conditional means across units; revealed $= Var(E[\beta \mid y]) / Var(\beta)$; and residual $= E[Var(\beta \mid y)] / Var(\beta)$.

References

Revelt, D. and Train, K. (2000). Customer-specific taste parameters and mixed logit: households' choice of electricity supplier. Working Paper E00-274, Department of Economics, University of California, Berkeley.

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

See Also

run_mxlogit (person_col), simulate_mxl_data (true_params$gamma_i)

Examples

sim <- simulate_mxl_data(N = 150, T = 5, J = 3, seed = 1)
fit <- run_mxlogit(sim$data, "id", "alt", "choice", c("x1", "x2"),
                   c("w1", "w2"), outside_opt_label = 0L,
                   rc_correlation = TRUE, person_col = "pid", S = 50L)
ct <- conditional_tastes(fit)
ct
# conditional means track the realized tastes, with shrinkage
cor(ct$mean["w1", ], sim$true_params$gamma_i["w1", colnames(ct$mean)])

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, ordered as in logsum), 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

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)

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)

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)

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)

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

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

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)

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

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)

Description

Create halton normal draws in appropriate format for mixed logit estimation

Usage

get_halton_normals(S, N, K_w)

Arguments

Details

Draw unit $i$ receives the $S$ consecutive points $(i - 1)S + 1, \ldots, iS$ of the $K_w$-dimensional Halton sequence of randtoolbox::halton(), mapped to standard normals. The cube takes $8 K_w S N$ bytes and is filled a block of units at a time, so it never holds a second copy of the sequence. Store mode supports at most $2^{31} - 1$ points ($S \times N$), the largest starting index halton(start = ) accepts; draws = "generate" in run_mxlogit has no such limit and never materializes the cube. predict(), logsum(), consumer_surplus(), elasticities(), diversion_ratios(), blp() and gof() on a store-mode fit with up to 128 random coefficients do not build their cube (one block per choice situation) with this function: they form the same draws on the fly, from the same uniforms mapped by R's qnorm() (bit for bit where randtoolbox and choicer are compiled with the same floating-point contraction, as on CRAN's platforms).

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

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

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)

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 (by decision maker first for a panel mixed logit whose person_col is present in newdata), 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))

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)

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

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])

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)

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

Details

With gen_seed >= 0 the draws are formed again at every evaluation of the shares. blp() on a fitted model runs the same contraction but keeps them across iterations within its keep_draws_bytes (see blp.choicer_mxl).

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

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.

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

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)

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

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)

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)

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)

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, one-column matrices with one row per row of the data; with a data.frame newdata, rows are ordered by id, then by fit-time alternative code (alt_int in object$alt_mapping). For "shares": the market shares as a one-column matrix with one row per alternative of object$alt_mapping, in its order (the outside option's first when it is modelled), 0 for an alternative that newdata does not offer.

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)

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, one-column matrices averaged across simulation draws with one row per row of the data, and with the outside option choice_prob_outside, one row per choice situation; with a data.frame newdata, rows are ordered by id, then by fit-time alternative code (alt_int in object$alt_mapping), and for a panel fit whose person_col is present in newdata, by decision maker first, which reproduces the in-sample row order. For "shares": the simulated market shares as a one-column matrix with one row per alternative of object$alt_mapping, in its order (the outside option's first when it is modelled), 0 for an alternative that newdata does not offer. Predictions are unconditional: each choice situation integrates over the population taste distribution (see conditional_tastes() for tastes conditional on observed choices).

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

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, one-column matrices with one row per row of the data; with a data.frame newdata, rows are ordered by id, then by fit-time alternative code (alt_int in object$alt_mapping). For "shares": the market shares as a one-column matrix with one row per alternative of object$alt_mapping, in its order (the outside option's first when it is modelled), 0 for an alternative that newdata does not offer.

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

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, double), 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).

Column names

No column the model uses (the id, alternative, choice, covariate, weight, cluster, nest or decision-maker column) may have a name starting with .choicer_, which choicer reserves for its working columns, and the alternative column may not be named alt_int, N_OBS, N_CHOICES, TAKE_RATE or MKT_SHARE, the names of the other columns of alt_mapping. Either is an error that names the column. Columns the model does not use may have any name.

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

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

Column names

No column the model uses (the id, alternative, choice, covariate, weight, cluster, nest or decision-maker column) may have a name starting with .choicer_, which choicer reserves for its working columns, and the alternative column may not be named alt_int, N_OBS, N_CHOICES, TAKE_RATE or MKT_SHARE, the names of the other columns of alt_mapping. Either is an error that names the column. Columns the model does not use may have any name.

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