Type:	Package
Title:	Bayesian Probit Choice Modeling
Version:	1.2.0
Description:	Bayes estimation of probit choice models in cross-sectional and panel settings. The package can analyze binary, multivariate, ordered, and ranked choices, as well as heterogeneity of choice behavior among deciders. The main functionality includes model fitting via Gibbs sampling, tools for convergence diagnostic, choice data simulation, in-sample and out-of-sample choice prediction, and model selection using information criteria and Bayes factors. The latent class model extension facilitates preference-based decider classification, where the number of latent classes can be inferred via the Dirichlet process or a weight-based updating heuristic. This allows for flexible modeling of choice behavior without the need to impose structural constraints. For a reference on the method, see Oelschlaeger and Bauer (2021) https://trid.trb.org/view/1759753.
URL:	https://loelschlaeger.de/RprobitB/, https://github.com/loelschlaeger/RprobitB
BugReports:	https://github.com/loelschlaeger/RprobitB/issues
License:	GPL-3
Encoding:	UTF-8
Imports:	checkmate, cli, crayon, doSNOW, foreach, ggplot2, graphics, gridExtra, MASS, mixtools, oeli (≥ 0.7.5), parallel, plotROC, progress, Rcpp, Rdpack, rlang, stats, utils, viridis
LinkingTo:	oeli (≥ 0.7.5), Rcpp, RcppArmadillo, testthat
Suggests:	knitr, mlogit, testthat (≥ 3.0.0), xml2
Depends:	R (≥ 4.1.0)
RoxygenNote:	7.3.2
RdMacros:	Rdpack
VignetteBuilder:	knitr
LazyData:	true
LazyDataCompression:	xz
Config/testthat/edition:	3
NeedsCompilation:	yes
Packaged:	2025-08-25 17:51:29 UTC; Lennart Oelschläger
Author:	Lennart Oelschläger [aut, cre], Dietmar Bauer [ctb]
Maintainer:	Lennart Oelschläger <oelschlaeger.lennart@gmail.com>
Repository:	CRAN
Date/Publication:	2025-08-25 18:30:12 UTC

RprobitB: Bayesian Probit Choice Modeling

Description

Bayes estimation of probit choice models in cross-sectional and panel settings. The package can analyze binary, multivariate, ordered, and ranked choices, as well as heterogeneity of choice behavior among deciders. The main functionality includes model fitting via Gibbs sampling, tools for convergence diagnostic, choice data simulation, in-sample and out-of-sample choice prediction, and model selection using information criteria and Bayes factors. The latent class model extension facilitates preference-based decider classification, where the number of latent classes can be inferred via the Dirichlet process or a weight-based updating heuristic. This allows for flexible modeling of choice behavior without the need to impose structural constraints. For a reference on the method, see Oelschlaeger and Bauer (2021) https://trid.trb.org/view/1759753.

Author(s)

Maintainer: Lennart Oelschläger oelschlaeger.lennart@gmail.com (ORCID)

Other contributors:

• Dietmar Bauer dietmar.bauer@uni-bielefeld.de (ORCID) [contributor]

See Also

Useful links:

• https://loelschlaeger.de/RprobitB/

• https://github.com/loelschlaeger/RprobitB

• Report bugs at https://github.com/loelschlaeger/RprobitB/issues

Compute Gelman-Rubin statistic

Description

This function computes the Gelman-Rubin statistic R_hat.

Usage

R_hat(samples, parts = 2)

Arguments

Details

NA values in samples are ignored. The degenerate case is indicated by NA. The Gelman-Rubin statistic is bounded by 1 from below. Values close to 1 indicate reasonable convergence.

Value

The Gelman-Rubin statistic.

Examples

no_chains <- 2
length_chains <- 1e3
samples <- matrix(NA_real_, length_chains, no_chains)
samples[1, ] <- 1
Gamma <- matrix(c(0.8, 0.1, 0.2, 0.9), 2, 2)
for (c in 1:no_chains) {
  for (t in 2:length_chains) {
    samples[t, c] <- sample(1:2, 1, prob = Gamma[samples[t - 1, c], ])
  }
}
R_hat(samples)

Create object of class RprobitB_data

Description

This function constructs an object of class RprobitB_data.

Usage

RprobitB_data(
  data,
  choice_data,
  N,
  T,
  J,
  P_f,
  P_r,
  alternatives,
  ordered,
  ranked,
  base,
  form,
  re,
  ASC,
  effects,
  standardize,
  simulated,
  choice_available,
  true_parameter,
  res_var_names
)

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

## S3 method for class 'RprobitB_data'
summary(object, ...)

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

## S3 method for class 'RprobitB_data'
plot(x, by_choice = FALSE, alpha = 1, position = "dodge", ...)

Arguments

Value

An object of class RprobitB_data.

Examples

data <- simulate_choices(
  form = choice ~ cost | 0,
  N = 100,
  T = 10,
  J = 2,
  alternatives = c("bus", "car"),
  true_parameter = list("alpha" = -1)
)
plot(data, by_choice = TRUE)

Create object of class RprobitB_fit

Description

This function creates an object of class RprobitB_fit.

Usage

RprobitB_fit(
  data,
  scale,
  level,
  normalization,
  R,
  B,
  Q,
  latent_classes,
  prior,
  gibbs_samples,
  class_sequence,
  comp_time
)

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

## S3 method for class 'RprobitB_fit'
summary(object, FUN = c(mean = mean, sd = stats::sd, `R^` = R_hat), ...)

## S3 method for class 'summary.RprobitB_fit'
print(x, digits = 2, ...)

Arguments

Value

An object of class RprobitB_fit.

Create object of class RprobitB_gibbs_samples_statistics

Description

This function creates an object of class RprobitB_gibbs_samples_statistics.

Usage

RprobitB_gibbs_samples_statistics(gibbs_samples, FUN = list(mean = mean))

## S3 method for class 'RprobitB_gibbs_samples_statistics'
print(x, true = NULL, digits = 2, ...)

Arguments

Value

An object of class RprobitB_gibbs_samples_statistics, which is a list of statistics from gibbs_samples obtained by applying the elements of FUN.

Create object of class RprobitB_latent_classes

Description

This function creates an object of class RprobitB_latent_classes which defines the number of latent classes and their updating scheme.

Usage

RprobitB_latent_classes(latent_classes = NULL)

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

Arguments

Value

An object of class RprobitB_latent_classes.

Utility normalization

Description

This function creates an object of class RprobitB_normalization, which defines the utility scale and level.

Usage

RprobitB_normalization(
  level,
  scale = "Sigma_1,1 := 1",
  form,
  re = NULL,
  alternatives,
  base,
  ordered = FALSE
)

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

Arguments

Details

Utility models require normalization with respect to level and scale.

• For level normalization, {RprobitB} takes utility differences with respect to one alternative. For the ordered model where only one utility is modelled, {RprobitB} fixes the first utility threshold to 0.

• For scale normalization, {RprobitB} fixes one model parameter. Per default, the first error-term variance is fixed to 1. This is specified via scale = "Sigma_1,1 := 1". Alternatively, any error-term variance or any non-random coefficient can be fixed.

Value

An object of class RprobitB_normalization, which is a list of

• level, a list with the elements level (the number of the alternative specified by the input level) and name (the name of the alternative, i.e. the input level), or alternatively NA in the ordered probit case,

• and scale, a list with the elements parameter (either "s" for an element of Sigma or "a" for an element of alpha), the parameter index, and the fixed value. If parameter = "a", also the name of the fixed effect.

Define probit model parameter

Description

This function creates an object of class RprobitB_parameter, which contains the parameters of a probit model.

If sample = TRUE, missing parameters are sampled. All parameters are checked against the values of P_f, P_r, J, and N.

Note that parameters are automatically ordered with respect to a non-ascending s for class identifiability.

Usage

RprobitB_parameter(
  P_f,
  P_r,
  J,
  N,
  C = 1,
  ordered = FALSE,
  alpha = NULL,
  s = NULL,
  b = NULL,
  Omega = NULL,
  Sigma = NULL,
  Sigma_full = NULL,
  beta = NULL,
  z = NULL,
  d = NULL,
  sample = TRUE
)

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

Arguments

Value

An object of class RprobitB_parameter, which is a named list with the model parameters.

Examples

RprobitB_parameter(P_f = 1, P_r = 2, J = 3, N = 10, C = 2)

Compute WAIC value

Description

This function computes the WAIC value of an RprobitB_fit object.

Usage

WAIC(x)

## S3 method for class 'RprobitB_waic'
print(x, digits = 2, ...)

## S3 method for class 'RprobitB_waic'
plot(x, ...)

Arguments

Details

WAIC is short for Widely Applicable (or Watanabe-Akaike) Information Criterion. As for AIC and BIC, the smaller the WAIC value the better the model. Its definition is

WAIC = -2 \cdot lppd + 2 \cdot p_{WAIC},

where lppd stands for log pointwise predictive density and p_{WAIC} is a penalty term proportional to the variance in the posterior distribution that is sometimes called effective number of parameters. The lppd is approximated as follows. Let

p_{is} = \Pr(y_i\mid \theta_s)

be the probability of observation y_i given the sth set \theta_s of parameter samples from the posterior. Then

lppd = \sum_i \log S^{-1} \sum_s p_{si}.

The penalty term is computed as the sum over the variances in log-probability for each observation:

p_{WAIC} = \sum_i V_{\theta} \left[ \log p_{si} \right].

Value

A numeric, the WAIC value, with the following attributes:

• se_waic, the standard error of the WAIC value,

• lppd, the log pointwise predictive density,

• p_waic, the effective number of parameters,

• p_waic_vec, the vector of summands of p_waic,

• p_si, the output of compute_p_si.

Re-label alternative specific covariates

Description

In {RprobitB}, alternative specific covariates must be named in the format "<covariate>_<alternative>". This helper function generates the format for a given choice_data set.

Usage

as_cov_names(choice_data, cov, alternatives)

Arguments

Value

The choice_data input with updated column names.

Examples

data("Electricity", package = "mlogit")
cov <- c("pf", "cl", "loc", "wk", "tod", "seas")
alternatives <- 1:4
colnames(Electricity)
Electricity <- as_cov_names(Electricity, cov, alternatives)
colnames(Electricity)

Check model formula

Description

This function checks the input form.

Usage

check_form(form, re = NULL, ordered = FALSE)

Arguments

Value

A list that contains the following elements:

• The input form.

• The name choice of the dependent variable in form.

• The input re.

• A list vars of three character vectors of covariate names of the three covariate types.

• A boolean ASC, determining whether the model has ASCs.

See Also

overview_effects() for an overview of the model effects

Check prior parameters

Description

This function checks the compatibility of submitted parameters for the prior distributions and sets missing values to default values.

Usage

check_prior(
  P_f,
  P_r,
  J,
  ordered = FALSE,
  mu_alpha_0 = numeric(P_f),
  Sigma_alpha_0 = 10 * diag(P_f),
  delta = 1,
  mu_b_0 = numeric(P_r),
  Sigma_b_0 = 10 * diag(P_r),
  n_Omega_0 = P_r + 2,
  V_Omega_0 = diag(P_r),
  n_Sigma_0 = J + 1,
  V_Sigma_0 = diag(J - 1),
  mu_d_0 = numeric(J - 2),
  Sigma_d_0 = diag(J - 2)
)

Arguments

Details

A priori-distributions:

• \alpha \sim N(\mu_{\alpha_0}, \Sigma_{\alpha_0})

• s \sim Dir(\delta)

• b_c \sim N(\mu_{b_0}, \Sigma_{b_0}) for all c

• \Omega_c \sim IW(n_{\Omega_0}, V_{\Omega_0}) for all c

• \Sigma \sim IW(n_{\Sigma_0}, V_{\Sigma_0})

• d \sim N(\mu_{d_0}, \Sigma_{d_0})

Value

An object of class RprobitB_prior, which is a list containing all prior parameters.

Examples

check_prior(P_f = 1, P_r = 2, J = 3, ordered = TRUE)

Compute choice probabilities

Description

This function returns the choice probabilities of an RprobitB_fit object.

Usage

choice_probabilities(x, data = NULL, par_set = mean)

Arguments

Value

A data frame of choice probabilities with choice situations in rows and alternatives in columns. The first two columns are the decider identifier "id" and the choice situation identifier "idc".

Examples

data <- simulate_choices(form = choice ~ covariate, N = 10, T = 10, J = 2)
x <- fit_model(data)
choice_probabilities(x)

Preference-based classification of deciders

Description

This function classifies the deciders based on their allocation to the components of the mixing distribution.

Usage

classification(x, add_true = FALSE)

Arguments

Details

The relative frequencies of which each decider was allocated to the classes during the Gibbs sampling are displayed. Only the thinned samples after the burn-in period are considered.

Value

A data.frame.

The row names are the decider identifiers.

The first C columns contain the relative frequencies with which the ⁠deciders are allocated to classes. Next, the column⁠est' contains the estimated class of the decider based on the highest allocation frequency.

If add_true = TRUE, the next column true contains the true class memberships.

See Also

update_z() for the updating function of the class allocation vector.

Extract model effects

Description

This function extracts the estimated model effects.

Usage

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

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

## S3 method for class 'RprobitB_coef'
plot(x, sd = 1, het = FALSE, ...)

Arguments

Value

An object of class RprobitB_coef.

Compute probit choice probabilities

Description

This is a helper function for choice_probabilities and computes the probit choice probabilities for a single choice situation with J alternatives.

Usage

compute_choice_probabilities(X, alternatives, parameter, ordered = FALSE)

Arguments

Value

A probability vector of length length(alternatives).

Compute choice probabilities at posterior samples

Description

This function computes the probability for each observed choice at the (normalized, burned and thinned) samples from the posterior.

These probabilities are required to compute the WAIC and the marginal model likelihood mml.

Usage

compute_p_si(x, ncores = parallel::detectCores() - 1, recompute = FALSE)

Arguments

Value

The object x, including the object p_si, which is a matrix of probabilities, observations in rows and posterior samples in columns.

Extract estimated covariance matrix of mixing distribution

Description

This convenience function returns the estimated covariance matrix of the mixing distribution.

Usage

cov_mix(x, cor = FALSE)

Arguments

Value

The estimated covariance matrix of the mixing distribution. In case of multiple classes, a list of matrices for each class.

Create lagged choice covariates

Description

This function creates lagged choice covariates from the data.frame choice_data, which is assumed to be sorted by the choice occasions.

Usage

create_lagged_cov(choice_data, column = character(), k = 1, id = "id")

Arguments

Details

Say that choice_data contains the column column. Then, the function call

create_lagged_cov(choice_data, column, k, id)

returns the input choice_data which includes a new column named column.k. This column contains for each decider (based on id) and each choice occasion the covariate faced before k choice occasions. If this data point is not available, it is set to NA. In particular, the first k values of column.k will be NA (initial condition problem).

Value

The input choice_data with the additional columns named column.k for each element column and each number k containing the lagged covariates.

Examples

choice_data <- data.frame(id = rep(1:2, each = 3), cov = LETTERS[1:6])
create_lagged_cov(choice_data, column = "cov", k = 1:2)

Transform increments to thresholds

Description

Transform increments to thresholds

Usage

d_to_gamma(d)

Arguments

Value

The threshold vector of length J + 1.

Examples

d_to_gamma(c(0, 0, 0))

Sample from prior distributions

Description

This function returns a sample from each parameter's prior distribution.

Usage

draw_from_prior(prior, C = 1)

Arguments

Value

A list of draws for alpha, s, b, Omega, and Sigma (if specified for the model).

Filter Gibbs samples

Description

This is a helper function that filters Gibbs samples.

Usage

filter_gibbs_samples(
  x,
  P_f,
  P_r,
  J,
  C,
  cov_sym,
  ordered = FALSE,
  keep_par = c("s", "alpha", "b", "Omega", "Sigma", "d"),
  drop_par = NULL
)

Arguments

Value

An object of class RprobitB_gibbs_samples filtered by the labels of parameter_labels.

Fit probit model to choice data

Description

This function performs MCMC simulation for fitting different types of probit models (binary, multivariate, mixed, latent class, ordered, ranked) to discrete choice data.

Usage

fit_model(
  data,
  scale = "Sigma_1,1 := 1",
  R = 1000,
  B = R/2,
  Q = 1,
  print_progress = getOption("RprobitB_progress", default = TRUE),
  prior = NULL,
  latent_classes = NULL,
  fixed_parameter = list(),
  save_beta_draws = FALSE
)

Arguments

Value

An object of class RprobitB_fit.

See Also

• prepare_data() and simulate_choices() for building an RprobitB_data object

• update() for estimating nested models

• transform() for transforming a fitted model

Examples

set.seed(1)
form <- choice ~ var | 0
data <- simulate_choices(form = form, N = 100, T = 10, J = 3, re = "var")
model <- fit_model(data = data, R = 1000)
summary(model)

Extract covariates of choice occasion

Description

This convenience function returns the covariates and the choices of specific choice occasions.

Usage

get_cov(x, id, idc, idc_label)

Arguments

Value

A subset of the choice_data data frame specified in prepare_data().

Gibbs sampler for probit models

Description

Gibbs sampler for probit models

Usage

gibbs_sampler(
  sufficient_statistics,
  prior,
  latent_classes,
  fixed_parameter,
  R,
  B,
  print_progress,
  ordered,
  ranked,
  save_beta_draws = FALSE
)

Arguments

Details

This function is not supposed to be called directly, but rather via fit_model.

Value

A list of Gibbs samples for

• Sigma,

• alpha (only if P_f > 0),

• s, z, b, Omega (only if P_r > 0),

• d (only if ordered = TRUE),

and a vector class_sequence of length R, where the r-th entry is the number of classes after iteration r.

Compute ordered probit log-likelihood

Description

Compute ordered probit log-likelihood

Usage

ll_ordered(d, y, sys, Tvec)

Arguments

Value

The ordered probit log-likelihood value.

Examples

d <- c(0, 0, 0)
y <- matrix(c(1, 2, 1, NA), ncol = 2)
sys <- matrix(c(0, 0, 0, NA), ncol = 2)
Tvec <- c(2, 1)
ll_ordered(d = d, y = y, sys = sys, Tvec = Tvec)

Handle missing covariates

Description

This function checks for and replaces missing covariate entries in choice_data.

Usage

missing_covariates(
  choice_data,
  impute = "complete_cases",
  col_ignore = character()
)

Arguments

Value

The input choice_data, in which missing covariates are addressed.

Approximate marginal model likelihood

Description

This function approximates the model's marginal likelihood.

Usage

mml(x, S = 0, ncores = parallel::detectCores() - 1, recompute = FALSE)

## S3 method for class 'RprobitB_mml'
print(x, log = FALSE, ...)

## S3 method for class 'RprobitB_mml'
plot(x, log = FALSE, ...)

Arguments

Details

The model's marginal likelihood p(y\mid M) for a model M and data y is required for the computation of Bayes factors. In general, the term has no closed form and must be approximated numerically.

This function uses the posterior Gibbs samples to approximate the model's marginal likelihood via the posterior harmonic mean estimator. To check the convergence, call plot(x$mml), where x is the output of this function. If the estimation does not seem to have converged, you can improve the approximation by combining the value with the prior arithmetic mean estimator. The final approximation of the model's marginal likelihood than is a weighted sum of the posterior harmonic mean estimate and the prior arithmetic mean estimate, where the weights are determined by the sample sizes.

Value

The object x, including the object mml, which is the model's approximated marginal likelihood value.

Gibbs sample mode

Description

This function approximates the Gibbs sample mode.

Usage

mode_approx(samples)

Arguments

Value

The (approximated) mode.

Examples

samples <- oeli::rmixnorm(
  n = 1000, mean = matrix(c(-2, 2), ncol = 2),
  Sigma = matrix(c(1, 1), ncol = 2), proportions = c(0.7, 0.3)
)
hist(samples)
mean(samples) # expected: 0.7 * (-2) + 0.3 * 2 = -0.8
mode_approx(samples) # expected: -2

Compare fitted models

Description

This function returns a table with several criteria for model comparison.

Usage

model_selection(
  ...,
  criteria = c("npar", "LL", "AIC", "BIC"),
  add_form = FALSE
)

## S3 method for class 'RprobitB_model_selection'
print(x, digits = 2, ...)

Arguments

Details

See the vignette on model selection for more details.

Value

A data.frame, criteria in columns, models in rows.

Extract number of model parameters

Description

This function extracts the number of model parameters of an RprobitB_fit object.

Usage

npar(object, ...)

## S3 method for class 'RprobitB_fit'
npar(object, ...)

Arguments

Value

Either a numeric value (if just one object is provided) or a numeric vector.

Print effect overview

Description

This function gives an overview of the effect names, whether the covariate is alternative-specific, whether the coefficient is alternative-specific, and whether it is a random effect.

Usage

overview_effects(
  form,
  re = NULL,
  alternatives,
  base = tail(alternatives, 1),
  ordered = FALSE
)

Arguments

Value

A data.frame, each row is a effect, columns are the effect name "effect", and booleans whether the covariate is alternative-specific "as_value", whether the coefficient is alternative-specific "as_coef", and whether it is a random effect "random".

See Also

check_form() for checking the model formula specification.

Examples

overview_effects(
  form = choice ~ price + time + comfort + change | 1,
  re = c("price", "time"),
  alternatives = c("A", "B"),
  base = "A"
)

Create parameters labels

Description

This function creates model parameter labels.

Usage

parameter_labels(
  P_f,
  P_r,
  J,
  C,
  cov_sym,
  ordered = FALSE,
  keep_par = c("s", "alpha", "b", "Omega", "Sigma", "d"),
  drop_par = NULL
)

create_labels_s(P_r, C)

create_labels_alpha(P_f)

create_labels_b(P_r, C)

create_labels_Omega(P_r, C, cov_sym)

create_labels_Sigma(J, cov_sym, ordered = FALSE)

create_labels_d(J, ordered)

Arguments

Value

A list of labels for the selected model parameters.

Visualize fitted probit model

Description

This function is the plot method for an object of class RprobitB_fit.

Usage

## S3 method for class 'RprobitB_fit'
plot(x, type, ignore = NULL, ...)

Arguments

Value

No return value. Draws a plot to the current device.

Autocorrelation plot of Gibbs samples

Description

This function plots the autocorrelation of the Gibbs samples. The plots include the total Gibbs sample size TSS and the effective sample size ESS, see the details.

Usage

plot_acf(gibbs_samples, par_labels)

Arguments

Details

The effective sample size is the value

TSS / \sqrt{1 + 2\sum_{k\geq 1} \rho_k}

, where \rho_k is the auto correlation between the chain offset by k positions. The auto correlations are estimated via spec.ar.

Value

No return value. Draws a plot to the current device.

Plot class allocation (for P_r = 2 only)

Description

This function plots the allocation of decision-maker specific coefficient vectors beta given the allocation vector z, the class means b, and the class covariance matrices Omega.

Usage

plot_class_allocation(beta, z, b, Omega, ...)

Arguments

Details

Only applicable in the two-dimensional case, i.e. only if P_r = 2.

Value

No return value. Draws a plot to the current device.

Visualizing the number of classes during Gibbs sampling

Description

This function plots the number of latent Glasses during Gibbs sampling to visualize the class updating.

Usage

plot_class_seq(class_sequence, B)

Arguments

Value

No return value. Draws a plot to the current device.

Plot bivariate contour of mixing distributions

Description

This function plots an estimated bivariate contour mixing distributions.

Usage

plot_mixture_contour(means, covs, weights, names)

Arguments

Value

An object of class ggplot.

Examples

means <- list(c(0, 0), c(2, 2))
covs <- list(diag(2), 0.5 * diag(2))
weights <- c(0.7, 0.3)
names <- c("A", "B")
plot_mixture_contour(means, covs, weights, names)

Plot marginal mixing distributions

Description

This function plots an estimated marginal mixing distributions.

Usage

plot_mixture_marginal(mean, cov, weights, name)

Arguments

Value

An object of class ggplot.

Plot ROC curve

Description

This function draws receiver operating characteristic (ROC) curves.

Usage

plot_roc(..., reference = NULL)

Arguments

Value

No return value. Draws a plot to the current device.

Visualizing the trace of Gibbs samples.

Description

This function plots traces of the Gibbs samples.

Usage

plot_trace(gibbs_samples, par_labels)

Arguments

Value

No return value. Draws a plot to the current device.

Compute point estimates

Description

This function computes the point estimates of an RprobitB_fit. Per default, the mean of the Gibbs samples is used as a point estimate. However, any statistic that computes a single numeric value out of a vector of Gibbs samples can be specified for FUN.

Usage

point_estimates(x, FUN = mean)

Arguments

Value

An object of class RprobitB_parameter.

Examples

data <- simulate_choices(form = choice ~ covariate, N = 10, T = 10, J = 2)
model <- fit_model(data)
point_estimates(model)
point_estimates(model, FUN = median)

Parameter sets from posterior samples

Description

This function builds parameter sets from the normalized, burned and thinned posterior samples.

Usage

posterior_pars(x)

Arguments

Value

A list of RprobitB_parameter objects.

Compute prediction accuracy

Description

This function computes the prediction accuracy of an RprobitB_fit object. Prediction accuracy means the share of choices that are correctly predicted by the model, where prediction is based on the maximum choice probability.

Usage

pred_acc(x, ...)

Arguments

Value

A numeric.

Predict choices

Description

This function predicts the discrete choice behaviour.

Usage

## S3 method for class 'RprobitB_fit'
predict(object, data = NULL, overview = TRUE, digits = 2, ...)

Arguments

Details

Predictions are made based on the maximum predicted probability for each choice alternative.

See the vignette on choice prediction for a demonstration on how to visualize the model's sensitivity and specificity by means of a receiver operating characteristic (ROC) curve.

Value

Either a table if overview = TRUE or a data frame otherwise.

Examples

set.seed(1)
data <- simulate_choices(form = choice ~ cov, N = 10, T = 10, J = 2)
data <- train_test(data, test_proportion = 0.5)
model <- fit_model(data$train)

predict(model)
predict(model, overview = FALSE)
predict(model, data = data$test)
predict(
  model,
  data = data.frame("cov_A" = c(1, 1, NA, NA), "cov_B" = c(1, NA, 1, NA)),
  overview = FALSE
)

Check for flip in preferences after change in model scale.

Description

This function checks if a change in the model scale flipped the preferences.

Usage

preference_flip(model_old, model_new)

Arguments

Value

No return value, called for side-effects.

Prepare choice data for estimation

Description

This function prepares choice data for estimation.

Usage

prepare_data(
  form,
  choice_data,
  re = NULL,
  alternatives = NULL,
  ordered = FALSE,
  ranked = FALSE,
  base = NULL,
  id = "id",
  idc = NULL,
  standardize = NULL,
  impute = "complete_cases"
)

Arguments

Details

Requirements for the data.frame choice_data:

• It must contain a column named id which contains unique identifier for each decision maker.

• It can contain a column named idc which contains unique identifier for each choice situation of each decision maker. If this information is missing, these identifier are generated automatically by the appearance of the choices in the data set.

• It can contain a column named choice with the observed choices, where choice must match the name of the dependent variable in form. Such a column is required for model fitting but not for prediction.

• It must contain a numeric column named p_j for each alternative specific covariate p in form and each choice alternative j in alternatives.

• It must contain a numeric column named q for each covariate q in form that is constant across alternatives.

In the ordered case (ordered = TRUE), the column choice must contain the full ranking of the alternatives in each choice occasion as a character, where the alternatives are separated by commas, see the examples.

See the vignette on choice data for more details.

Value

An object of class RprobitB_data.

See Also

• check_form() for checking the model formula

• overview_effects() for an overview of the model effects

• create_lagged_cov() for creating lagged covariates

• as_cov_names() for re-labeling alternative-specific covariates

• simulate_choices() for simulating choice data

• train_test() for splitting choice data into a train and test subset

Examples

data <- prepare_data(
  form = choice ~ price + time + comfort + change | 0,
  choice_data = train_choice,
  re = c("price", "time"),
  id = "deciderID",
  idc = "occasionID",
  standardize = c("price", "time")
)

### ranked case
choice_data <- data.frame(
  "id" = 1:3, "choice" = c("A,B,C", "A,C,B", "B,C,A"), "cov" = 1
)
data <- prepare_data(
  form = choice ~ 0 | cov + 0,
  choice_data = choice_data,
  ranked = TRUE
)

Sample allocation

Description

Sample allocation

Usage

sample_allocation(prob)

Arguments

Value

An integer 1:C.

Examples

sample_allocation(c(0.5, 0.3, 0.2))

Simulate choice data

Description

This function simulates choice data from a probit model.

Usage

simulate_choices(
  form,
  N,
  T = 1,
  J,
  re = NULL,
  alternatives = NULL,
  ordered = FALSE,
  ranked = FALSE,
  base = NULL,
  covariates = NULL,
  true_parameter = list()
)

Arguments

Details

See the vignette on choice data for more details.

Value

An object of class RprobitB_data.

See Also

• check_form() for checking the model formula

• overview_effects() for an overview of the model effects

• create_lagged_cov() for creating lagged covariates

• as_cov_names() for re-labelling alternative-specific covariates

• prepare_data() for preparing empirical choice data

• train_test() for splitting choice data into a train and test subset

Examples

### simulate data from a binary probit model with two latent classes
data <- simulate_choices(
  form = choice ~ cost | income | time,
  N = 100,
  T = 10,
  J = 2,
  re = c("cost", "time"),
  alternatives = c("car", "bus"),
  true_parameter = list(
    "alpha" = c(-1, 1),
    "b" = matrix(c(-1, -1, -0.5, -1.5, 0, -1), ncol = 2),
    "C" = 2
  )
)

### simulate data from an ordered probit model
data <- simulate_choices(
  form = opinion ~ age + gender,
  N = 10,
  T = 1:10,
  J = 5,
  alternatives = c("very bad", "bad", "indifferent", "good", "very good"),
  ordered = TRUE,
  covariates = list(
    "gender" = rep(sample(c(0, 1), 10, replace = TRUE), times = 1:10)
  )
)

### simulate data from a ranked probit model
data <- simulate_choices(
  form = product ~ price,
  N = 10,
  T = 1:10,
  J = 3,
  alternatives = c("A", "B", "C"),
  ranked = TRUE
)

Compute sufficient statistics

Description

This function computes sufficient statistics from an RprobitB_data object for the Gibbs sampler to save computation time.

Usage

sufficient_statistics(data, normalization)

Arguments

Value

A list of sufficient statistics on the data for Gibbs sampling, containing

• the elements N, T, J, P_f and P_r from data,

• Tvec, the vector of choice occasions for each decider of length N,

• csTvec, a vector of length N with the cumulated sums of Tvec starting from 0,

• W, a list of design matrices differenced with respect to alternative number normalization$level$level for each decider in each choice occasion with covariates that are linked to a fixed coefficient (or NA if P_f = 0),

• X, a list of design matrices differenced with respect to alternative number normalization$level$level for each decider in each choice occasion with covariates that are linked to a random coefficient (or NA if P_r = 0),

• y, a matrix of dimension N x max(Tvec) with the observed choices of deciders in rows and choice occasions in columns, decoded to numeric values with respect to their appearance in data$alternatives, where rows are filled with NA in case of an unbalanced panel,

• WkW, a matrix of dimension P_f^2 x (J-1)^2, the sum over Kronecker products of each transposed element in W with itself,

• XkX, a list of length N, each element is constructed in the same way as WkW but with the elements in X and separately for each decider,

• rdiff (for the ranked case only), a list of matrices that reverse the base differencing and instead difference in such a way that the resulting utility vector is negative.

Stated Preferences for Train Traveling

Description

Data set of 2929 stated choices by 235 Dutch individuals deciding between two virtual train trip options "A" and "B" based on the price, the travel time, the number of rail-to-rail transfers (changes), and the level of comfort.

The data were obtained in 1987 by Hague Consulting Group for the National Dutch Railways. Prices were recorded in Dutch guilder and in this data set transformed to Euro at an exchange rate of 2.20371 guilders = 1 Euro.

Usage

train_choice

Format

A data.frame with 2929 rows and 11 columns:

deciderID

integer identifier for the decider

occasionID

integer identifier for the choice occasion

choice

character for the chosen alternative (either "A" or "B")

price_A

numeric price for alternative "A" in Euro

time_A

numeric travel time for alternative "A" in hours

change_A

integer number of changes for alternative "A"

comfort_A

integer comfort level (in decreasing comfort order) for alternative "A"

price_B

numeric price for alternative "B" in Euro

time_B

numeric travel time for alternative "B" in hours

change_B

integer number of changes for alternative "B"

comfort_B

integer comfort level (in decreasing comfort order) for alternative "B"

References

Ben-Akiva M, Bolduc D, Bradley M (1993). “Estimation of travel choice models with randomly distributed values of time.” Transportation Research Record, 1413.

Meijer E, Rouwendal J (2006). “Measuring welfare effects in models with random coefficients.” Journal of Applied Econometrics, 21(2).

Split choice data into train and test subset

Description

This function splits choice data into a train and a test part.

Usage

train_test(
  x,
  test_proportion = NULL,
  test_number = NULL,
  by = "N",
  random = FALSE
)

Arguments

Value

A list with two objects of class RprobitB_data, named "train" and "test".

Examples

### simulate choices for demonstration
x <- simulate_choices(form = choice ~ covariate, N = 10, T = 10, J = 2)

### 70% of deciders in the train subsample,
### 30% of deciders in the test subsample
train_test(x, test_proportion = 0.3, by = "N")

### 2 randomly chosen choice occasions per decider in the test subsample,
### the rest in the train subsample
train_test(x, test_number = 2, by = "T", random = TRUE)

Transform fitted probit model

Description

Given an object of class RprobitB_fit, this function can:

• change the length B of the burn-in period,

• change the the thinning factor Q of the Gibbs samples,

• change the utility scale.

Usage

## S3 method for class 'RprobitB_fit'
transform(
  `_data`,
  B = NULL,
  Q = NULL,
  scale = NULL,
  check_preference_flip = TRUE,
  ...
)

Arguments

Details

See the vignette "Model fitting" for more details: vignette("v03_model_fitting", package = "RprobitB").

Value

The transformed RprobitB_fit object.

Transformation of Gibbs samples

Description

This function normalizes, burns and thins the Gibbs samples.

Usage

transform_gibbs_samples(gibbs_samples, R, B, Q, normalization)

Arguments

Value

A list, the first element gibbs_sampes_raw is the input gibbs_samples, the second element is the normalized, burned, and thinned version of gibbs_samples called gibbs_samples_nbt. The list gets the class RprobitB_gibbs_samples.

Transformation of parameter values

Description

This function transforms parameter values based on normalization.

Usage

transform_parameter(parameter, normalization, ordered = FALSE)

Arguments

Value

An object of class RprobitB_parameter.

Update and re-fit probit model

Description

This function estimates a nested probit model based on a given RprobitB_fit object.

Usage

## S3 method for class 'RprobitB_fit'
update(
  object,
  form,
  re,
  alternatives,
  id,
  idc,
  standardize,
  impute,
  scale,
  R,
  B,
  Q,
  print_progress,
  prior,
  latent_classes,
  ...
)

Arguments

Details

All parameters (except for object) are optional and if not specified retrieved from the specification for object.

Value

An object of class RprobitB_fit.

Update class covariances

Description

Update class covariances

Usage

update_Omega(beta, b, z, m, n_Omega_0, V_Omega_0)

Arguments

Value

A matrix of updated covariance matrices for each class in columns.

Examples

N <- 100
b <- cbind(c(0, 0), c(1, 1))
Omega <- matrix(c(1, 0.3, 0.3, 0.5, 1, -0.3, -0.3, 0.8), ncol = 2)
z <- c(rep(1, N / 2), rep(2, N / 2))
m <- as.numeric(table(z))
beta <- sapply(
  z, function(z) oeli::rmvnorm(n = 1, b[, z], matrix(Omega[, z], 2, 2))
)
update_Omega(
  beta = beta, b = b, z = z, m = m,
  n_Omega_0 = 4, V_Omega_0 = diag(2)
)

Update covariance of a single class

Description

Update covariance of a single class

Usage

update_Omega_c(S_c, m_c, n_Omega_0, V_Omega_0)

Arguments

Value

An update for Omega_c.

Examples

update_Omega_c(S_c = diag(2), m_c = 10, n_Omega_0 = 4, V_Omega_0 = diag(2))

Update error covariance matrix

Description

Update error covariance matrix

Usage

update_Sigma(n_Sigma_0, V_Sigma_0, N, S)

Arguments

Value

An update for Sigma.

Examples

(Sigma_true <- matrix(c(1, 0.5, 0.2, 0.5, 1, 0.2, 0.2, 0.2, 2), ncol = 3))
beta <- matrix(c(-1, 1), ncol = 1)
N <- 100
X <- replicate(N, matrix(rnorm(6), ncol = 2), simplify = FALSE)
eps <- replicate(
  N, oeli::rmvnorm(n = 1, mean = c(0, 0, 0), Sigma = Sigma_true),
  simplify = FALSE
)
U <- mapply(function(X, eps) X %*% beta + eps, X, eps, SIMPLIFY = FALSE)
n_Sigma_0 <- 4
V_Sigma_0 <- diag(3)
outer_prod <- function(X, U) (U - X %*% beta) %*% t(U - X %*% beta)
S <- Reduce(`+`, mapply(
  function(X, U) (U - X %*% beta) %*% t(U - X %*% beta), X, U,
  SIMPLIFY = FALSE
))
Sigma_draws <- replicate(100, update_Sigma(n_Sigma_0, V_Sigma_0, N, S))
apply(Sigma_draws, 1:2, mean)

Update utility vector

Description

Update utility vector

Usage

update_U(U, y, sys, Sigma_inv)

Arguments

Value

An update for (a single) U.

Examples

U <- sys <- c(0, 0, 0)
Sigma_inv <- diag(3)
lapply(1:4, function(y) update_U(U, y, sys, Sigma_inv))

Update ranked utility vector

Description

Update ranked utility vector

Usage

update_U_ranked(U, sys, Sigma_inv)

Arguments

Value

An update for (a single) ranked U.

Examples

U <- sys <- c(0, 0)
Sigma_inv <- diag(2)
update_U_ranked(U, sys, Sigma_inv)

Update class means

Description

Update class means

Usage

update_b(beta, Omega, z, m, Sigma_b_0_inv, mu_b_0)

Arguments

Value

A matrix of updated means for each class in columns.

Examples

N <- 100
b <- cbind(c(0, 0), c(1, 1))
Omega <- matrix(c(1, 0.3, 0.3, 0.5, 1, -0.3, -0.3, 0.8), ncol = 2)
z <- c(rep(1, N / 2), rep(2, N / 2))
m <- as.numeric(table(z))
beta <- sapply(
  z, function(z) oeli::rmvnorm(n = 1, b[, z], matrix(Omega[, z], 2, 2))
)
update_b(
  beta = beta, Omega = Omega, z = z, m = m,
  Sigma_b_0_inv = diag(2), mu_b_0 = c(0, 0)
)

Update mean of a single class

Description

Update mean of a single class

Usage

update_b_c(bar_b_c, Omega_c, m_c, Sigma_b_0_inv, mu_b_0)

Arguments

Value

An update for b_c.

Examples

update_b_c(
  bar_b_c = c(0, 0), Omega_c = diag(2), m_c = 10,
  Sigma_b_0_inv = diag(2), mu_b_0 = c(0, 0)
)

Dirichlet process class updates

Description

Dirichlet process class updates

Usage

update_classes_dp(
  beta,
  z,
  b,
  Omega,
  delta,
  mu_b_0,
  Sigma_b_0,
  n_Omega_0,
  V_Omega_0,
  identify_classes = FALSE,
  Cmax = 10L
)

Arguments

Value

A list of updated values for z, b, Omega, and C.

Examples

set.seed(1)
z <- c(rep(1, 10),rep(2, 10))
b <- matrix(c(5, 5, 5, -5), ncol = 2)
Omega <- matrix(c(1, 0.3, 0.3, 0.5, 1, -0.3, -0.3, 0.8), ncol = 2)
beta <- sapply(
  z, function(z) oeli::rmvnorm(n = 1, b[, z], matrix(Omega[, z], 2, 2))
)
beta[, 1] <- c(-10, 10)
update_classes_dp(
  beta = beta, z = z, b = b, Omega = Omega,
  delta = 1, mu_b_0 = numeric(2), Sigma_b_0 = diag(2),
  n_Omega_0 = 4, V_Omega_0 = diag(2)
)

Weight-based class updates

Description

Weight-based class updates

Usage

update_classes_wb(
  s,
  b,
  Omega,
  epsmin = 0.01,
  epsmax = 0.7,
  deltamin = 0.1,
  deltashift = 0.5,
  identify_classes = FALSE,
  Cmax = 10L
)

Arguments

Details

The following updating rules apply:

• Class c is removed if s_c < \epsilon_{min}.

• Class c is split into two classes, if s_c > \epsilon_{max}.

• Two classes c_1 and c_2 are merged to one class, if ||b_{c_1} - b_{c_2}|| < \delta_{min}.

Value

A list of updated values for s, b, and Omega and the indicator update_type which signals the type of class update:

• 0: no update

• 1: removed class

• 2: split class

• 3: merged classes

Examples

s <- c(0.7, 0.3)
b <- matrix(c(1, 1, 1, -1), ncol = 2)
Omega <- matrix(c(0.5, 0.3, 0.3, 0.5, 1, -0.1, -0.1, 0.8), ncol = 2)

### no update
update_classes_wb(s = s, b = b, Omega = Omega)

### remove class 2
update_classes_wb(s = s, b = b, Omega = Omega, epsmin = 0.31)

### split class 1
update_classes_wb(s = s, b = b, Omega = Omega, epsmax = 0.69)

### merge classes 1 and 2
update_classes_wb(s = s, b = b, Omega = Omega, deltamin = 3)

Update coefficient vector

Description

Update coefficient vector

Usage

update_coefficient(mu_beta_0, Sigma_beta_0_inv, XSigX, XSigU)

Arguments

Value

An update for the coefficient vector.

Examples

beta_true <- matrix(c(-1, 1), ncol = 1)
Sigma <- matrix(c(1, 0.5, 0.2, 0.5, 1, 0.2, 0.2, 0.2, 2), ncol = 3)
N <- 100
X <- replicate(N, matrix(rnorm(6), ncol = 2), simplify = FALSE)
eps <- replicate(
  N, oeli::rmvnorm(n = 1, mean = c(0, 0, 0), Sigma = Sigma),
  simplify = FALSE
)
U <- mapply(
  function(X, eps) X %*% beta_true + eps, X, eps, SIMPLIFY = FALSE
)
mu_beta_0 <- c(0, 0)
Sigma_beta_0_inv <- diag(2)
XSigX <- Reduce(
  `+`, lapply(X, function(X) t(X) %*% solve(Sigma) %*% X)
)
XSigU <- Reduce(
  `+`, mapply(function(X, U) t(X) %*% solve(Sigma) %*% U, X, U,
  SIMPLIFY = FALSE)
)
R <- 10
beta_draws <- replicate(
  R, update_coefficient(mu_beta_0, Sigma_beta_0_inv, XSigX, XSigU),
  simplify = TRUE
)
rowMeans(beta_draws)

Update utility threshold increments

Description

Update utility threshold increments

Usage

update_d(d, y, sys, ll, mu_d_0, Sigma_d_0, Tvec, step_scale = 0.1)

Arguments

Value

An update for d.

Examples

set.seed(1)
N <- 1000
d_true <- rnorm(2)
gamma <- c(-Inf, 0, cumsum(exp(d_true)), Inf)
X <- matrix(rnorm(N * 2L), ncol = 2L)
beta <- c(0.8, -0.5)
mu <- matrix(as.vector(X %*% beta), ncol = 1L)
U <- rnorm(N, mean = mu[, 1], sd = 1)
yvec <- as.integer(cut(U, breaks = gamma, labels = FALSE))
y <- matrix(yvec, ncol = 1L)
Tvec <- rep(1, N)
mu_d_0 <- c(0, 0)
Sigma_d_0 <- diag(2) * 5
d <- rnorm(2)
ll <- -Inf
R <- 1000
for (iter in seq_len(R)) {
  ans <- update_d(
    d = d, y = y, sys = mu, ll = ll, mu_d_0 = mu_d_0, Sigma_d_0 = Sigma_d_0,
    Tvec = Tvec
  )
  d  <- ans$d
  ll <- ans$ll
}
cbind("true" = d_true, "sample" = d) |> round(2)

Update class sizes

Description

Update class sizes

Usage

update_m(C, z, non_zero = FALSE)

Arguments

Value

An update for m.

Examples

update_m(C = 4, z = c(1, 1, 1, 2, 2, 3))

Update class weight vector

Description

Update class weight vector

Usage

update_s(delta, m)

Arguments

Value

An update for s.

Examples

update_s(delta = 1, m = 4:1)

Update class allocation vector

Description

Update class allocation vector

Usage

update_z(s, beta, b, Omega)

Arguments

Value

An update for z.

Examples

update_z(
  s = c(0.6, 0.4), beta = matrix(c(-2, 0, 2), ncol = 3),
  b = cbind(0, 1), Omega = cbind(1, 1)
)