Package 'RprobitB'

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)

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

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)

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)

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.

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.

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.

Description

This helper 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.

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)

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

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)

Description

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

Usage

get_cov(x, id = NULL, idc = NULL, id_label = NULL, idc_label = NULL)

Arguments

Value

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

Examples

data <- simulate_choices(
  form = product ~ price,
  N = 10,
  T = 1:10,
  J = 3,
  ranked = TRUE
)
get_cov(data, id = 2)

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.

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)

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.

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

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.

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.

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

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.

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)

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.

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.

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 = 100, latent_classes = list(C = 2, "dp_update" = TRUE)
)
plot(model, type = "mixture")
plot(model, type = "acf", ignore = c("s", "Omega", "Sigma"))
plot(model, type = "trace", ignore = c("s", "Omega", "Sigma"))
plot(model, type = "class_seq")

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)

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.

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
)

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
)

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)

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)

Description

Sample allocation

Usage

sample_allocation(prob)

Arguments

Value

An integer 1:C.

Examples

sample_allocation(c(0.5, 0.3, 0.2))

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
)

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

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)

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.

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

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

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

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)

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)

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)

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

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

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

Description

Update class weight vector

Usage

update_s(delta, m)

Arguments

Value

An update for s.

Examples

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

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)

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

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)

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

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.