Model specification and variants
In RprobitB: Bayesian Probit Choice Modeling

knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4,
  fig.align = "center"
)

This vignette describes how to specify a probit choice model with fit(). It covers the normalization of the utility scale and the prior distribution, which every fit involves, and then the model variants: the three types of covariates and the alternative-specific constants, choice sets that differ between occasions, and ordered and ranked responses. Preference heterogeneity between deciders is the subject of the vignette Modeling preference heterogeneity, and @Oelschlaeger2026c gives the methodological background. From the covariate types onward, each section first fits the variant to simulated data and then to a data set of the AER package [@Kleiber2008], the MASS package [@VenablesRipley2002], or the mlogit package [@Croissant2020]. Without data, fit() simulates the requested model before estimating it, and summary() prints the true values in a dgp column beside the posterior summaries, converted to the normalization of the fit.

library(RprobitB)
set.seed(1)

Normalization

The probit model introduced in the vignette Get started with RprobitB assigns every alternative $j$ at occasion $t$ of decider $n$ the latent utility $U_{ntj} = X_{ntj}^\top \beta_n + \epsilon_{ntj}$ with jointly normal errors $\epsilon_{nt} \sim \mathrm{N}(0, \Sigma)$, and the decider chooses the alternative with the largest utility. The choice probability of an alternative is the probability that its utility exceeds the utilities of all others. These probabilities are invariant to adding a constant to all utilities and to multiplying all utilities by a positive number, so neither the level nor the scale of the utilities is identified. RprobitB therefore works with utility differences relative to a base alternative, by default the first in alphabetical order, and fixes one parameter through scale:

scale = NULL fixes the error variance of the first utility difference to one. This is the default. The covariance matrix of the error differences is reported as Sigma[<alternative>,<alternative>], with rows and columns named by the alternatives other than the base. For three alternatives A, B, and C with base A, its entries are Sigma[B,B], Sigma[C,B], and Sigma[C,C].
scale = c(z = -1) fixes the coefficient of z to -1 instead, so that all other coefficients are measured in units of it. With a price coefficient fixed in this way, the other coefficients are willingness-to-pay values.

The sampler draws all parameters without restriction and rescales every retained draw afterwards. Variables that the normalization fixes remain in the draws but are omitted from summary(), coef(), and vcov().

The following demonstration shows the effect of the normalization on the reported values. We simulate data with coefficients 1 for x and -0.5 for z and fit them under the default scale. Simulated alternatives are labeled with capital letters unless alternatives names them, here A and B, and A is the base. update() then refits the same simulated data with the coefficient of z fixed to -1, so only the normalization differs between the two fits.

scale_default <- fit(
  choice ~ x + z | 0,
  dgp_parameters = list(beta = c(x = 1, z = -0.5)),
  n_deciders = 300,
  chains = 1
)
summary(scale_default)
scale_z <- update(scale_default, scale = c(z = -1))
summary(scale_z)

Under the default scale, the dgp column shows the coefficients as specified. Sigma[B,B], the error variance of the utility difference between alternative B and the base A, is fixed to one and therefore not listed. Under the coefficient normalization, beta[z] is fixed to -1 instead of its true -0.5, so all utilities are multiplied by 2: the true beta[x] becomes 2 and the true error variance 4. summary() converts the true values to the normalization of the fit, so the dgp column remains comparable.

The convergence diagnostics of the two tables differ because update() runs the sampler again and because the rescaled variables are different functions of the draws: beta[x] is now a ratio of two coefficients, which mixes differently than a single coefficient.

Prior distribution

RprobitB estimates every model with a Gibbs sampler, which draws each block of parameters in turn from its conditional posterior distribution given the other parameters and the latent utilities. The prior is conjugate for every block, so each conditional posterior belongs to the same family as the prior. Fixed coefficients and class means have normal priors, covariance matrices have inverse Wishart priors, class weights have a Dirichlet prior, and the log-increments between ordered thresholds have a normal prior. The section "Prior distribution" of ?fit lists all components with their defaults, which are weakly informative for coefficients and covariances. Single components are overridden through a named list, and the complete prior of a fit is stored in its prior component.

The following demonstration simulates 100 deciders with a true coefficient of -1 under the default prior and then refits the same data under two priors with mean 1: a moderate one with variance 0.5 and a tight one with variance 0.01.

default_prior <- fit(
  choice ~ x | 0,
  dgp_parameters = list(beta = c(x = -1)),
  chains = 1
)
default_prior$prior
summary(default_prior)
moderate_prior <- update(
  default_prior, prior = list(fixed_mean = 1, fixed_covariance = matrix(0.5))
)
tight_prior <- update(
  default_prior, prior = list(fixed_mean = 1, fixed_covariance = matrix(0.01))
)
data.frame(
  variable = "beta[x]", dgp = -1, default = coef(default_prior),
  moderate = coef(moderate_prior), tight = coef(tight_prior),
  row.names = NULL
)

The table compares the posterior means under the three priors with the true value. The moderate prior shifts the posterior mean only slightly towards the prior mean. The tight prior has a standard deviation of 0.1 around 1, outweighs the data, and pulls the posterior mean far away from the true value. A prior on a coefficient is unproblematic as long as its variance reflects the actual uncertainty about the coefficient.

Covariate types and alternative-specific constants

The formula choice ~ A | B | C distinguishes three kinds of covariates: attributes of the alternatives with one shared coefficient (A), characteristics of the decider with alternative-specific coefficients (B), and attributes of the alternatives with alternative-specific coefficients (C). Alternative-specific constants are included by default and removed with 0 in the second part. The following demonstration simulates all three types at once, and summary() shows the posterior summaries beside the true values.

covariate_types <- fit(
  choice ~ x | z | w,
  n_deciders = 500,
  dgp_parameters = list(beta = c(
    x = 0.5, z_B = -0.5, ASC_B = 0.25, w_A = -0.5, w_B = 0.5
  )),
  chains = 1
)
summary(covariate_types)

The summary lists the shared coefficient of x, the coefficient of z and the constant of alternative B, and the two alternative-specific coefficients of w, each beside its true value. base selects the alternative against which the alternative-specific coefficients are measured, by default the first in alphabetical order, here A. Switching the base to B changes only the parameterization, not the model:

base_b <- update(covariate_types, base = "B")
coef(base_b)[c("beta[x]", "beta[z_A]", "beta[ASC_A]")]

The shared coefficient of x is unchanged up to Monte Carlo error, while the coefficient of z and the constant now describe alternative A relative to B and therefore change sign.

In 1987, 210 travelers between Sydney and Melbourne reported which of four modes they had taken: air, train, bus, or car. The TravelMode data of the AER package [@Kleiber2008] are in long format with one row per mode and a "yes"/"no" indicator of the chosen mode, which fit() expects as a logical or 0/1 variable. Terminal waiting time (wait), in-vehicle cost (vcost), and travel time (travel) vary across modes and enter as type A. Household income and the size of the traveling party describe the traveler and enter as type B, with the alternative-specific constants included by default. Cost and income are converted from Australian dollars to euro.

travel_formula <- choice ~ wait + vcost + travel | income + size

The first alternative in alphabetical order, air, is the base, so the constants and the coefficients of income and party size describe the other modes relative to flying.

data("TravelMode", package = "AER")
TravelMode$choice <- TravelMode$choice == "yes"
TravelMode$vcost <- TravelMode$vcost / 1.6196
TravelMode$income <- TravelMode$income / 1.6196
travel <- fit(
  formula = travel_formula,
  data = TravelMode,
  format = "long",
  column_decider = "individual",
  column_alternative = "mode",
  iterations = 6000,
  warmup = 3000,
  thin = 30,
  chains = 2,
  progress = FALSE
)
summary(travel)

The three attribute coefficients are negative: waiting, cost, and travel time reduce the utility of a mode. The income and party size coefficients are alternative-specific and have no direct interpretation on the utility scale, so interpret(type = "mea") computes marginal effects, the derivatives of the choice probabilities with respect to a covariate, for a traveler with average covariates:

mode_effects <- interpret(travel, type = "mea")
mode_effects[mode_effects$covariate == "income", ]

A higher household income lowers the probability of the train and raises the probabilities of the plane and the car.

Individual choice sets

Not every alternative is available at every occasion: a traveler without a car cannot drive, and a route without a rail link has no train option. Unordered choices can therefore be made from occasion-specific subsets of the alternatives. In long format, an occasion lists only the rows of its available alternatives, and no further argument is needed. The sampler imputes the latent utilities of unavailable alternatives without restriction, so they do not affect the choice, and predictions assign them probability zero.

Between Montreal and Toronto, travelers can fly, drive, or take the train, but not all modes are available on every trip. The ModeCanada data of the mlogit package cover 4324 trips in this corridor [@Bhat1995].

data("ModeCanada", package = "mlogit")
head(ModeCanada)

Cost (cost), in-vehicle time (ivt), out-of-vehicle time (ovt), and service frequency (freq) vary across modes; household income and the number of urban trip endpoints are trip-specific. The data also contain the bus, which was chosen on only r sum(ModeCanada$alt == "bus" & ModeCanada$choice == 1) trips, too few to identify its coefficients and error covariances, so the bus rows and the trips on which it was chosen are removed. Cost and income are converted from Canadian dollars to euro, and trips with a single remaining alternative are dropped. The first 1000 of the remaining trips keep the computation short.

ModeCanada$cost <- ModeCanada$cost / 1.6151
ModeCanada$income <- ModeCanada$income / 1.6151
bus_trips <- ModeCanada$case[ModeCanada$alt == "bus" & ModeCanada$choice == 1]
canada_data <- ModeCanada[
  ModeCanada$alt != "bus" & !(ModeCanada$case %in% bus_trips),
]
set_size <- table(canada_data$case)
canada_data <- canada_data[set_size[as.character(canada_data$case)] > 1, ]
canada_data <- canada_data[
  canada_data$case %in% unique(canada_data$case)[1:1000],
]
table(table(canada_data$case))

The choice sets are read from the rows of each trip. Every trip is one decider, and the model has the same structure as the travel mode model above: three attributes of type A, two trip characteristics of type B, and the constants, all relative to air.

canada <- fit(
  choice ~ cost + ivt + ovt + freq | income + urban,
  data = canada_data,
  format = "long",
  column_decider = "case",
  column_alternative = "alt",
  iterations = 6000,
  warmup = 3000,
  thin = 15,
  chains = 2,
  progress = FALSE
)
summary(canada)

The income coefficients of car and train are negative: a higher income increases the probability of flying.

Ordered responses

Some responses are ordered levels. An ordered model has one latent utility per occasion and compares it with increasing thresholds gamma; the level is the interval into which the utility falls. alternatives gives the response levels in increasing order. Latent-variable data augmentation provides a direct Bayesian treatment of ordered probit models [@Albert1993]. The following demonstration estimates a simulated three-category model and reports the coefficient and the free threshold beside their true values.

ordered_sim <- fit(
  choice ~ x | 0,
  alternatives = c("low", "middle", "high"),
  choice_type = "ordered",
  n_deciders = 500,
  dgp_parameters = list(beta = c(x = 1), gamma = c(0, 1)),
  chains = 1
)
summary(ordered_sim)

The summary lists the coefficient and the free threshold gamma[2] beside their true values. The first threshold is fixed to zero and the error variance to one, which identifies the level and the scale of the utility.

The survey data of the MASS package [@VenablesRipley2002] come from 237 statistics students at the University of Adelaide who reported how often they smoke, together with their age and how much they exercise. The smoking level Smoke is stored as a factor whose levels are in alphabetical order; alternatives puts them in their natural order from never to heavy. The other survey questions are not used; fit() ignores their missing values.

data("survey", package = "MASS")
levels(survey$Smoke)
smoking_levels <- c("Never", "Occas", "Regul", "Heavy")
smoking <- fit(
  Smoke ~ Age + Exer | 0,
  data = survey,
  alternatives = smoking_levels,
  choice_type = "ordered",
  column_decider = NULL,
  chains = 1
)
summary(smoking)

With four levels, the thresholds gamma[2] and gamma[3] are estimated. Each student has one latent utility, normally distributed around its systematic part, and the thresholds partition it into the four levels. The area under the density between two thresholds is the probability of that level. The figure shows the density of a student whose systematic utility is zero.

thresholds <- coef(smoking)[c("gamma[2]", "gamma[3]")]
cuts <- c(-Inf, 0, thresholds, Inf)
shades <- grey(seq(0.45, 0.9, length.out = length(smoking_levels)))
utility <- seq(-3.5, 3.5, length.out = 400)
plot(
  utility, dnorm(utility),
  type = "n", axes = FALSE, ylab = "",
  xlab = "latent utility of a student"
)
for (k in seq_along(smoking_levels)) {
  inside <- utility >= cuts[k] & utility <= cuts[k + 1]
  polygon(
    c(max(cuts[k], -3.5), utility[inside], min(cuts[k + 1], 3.5)),
    c(0, dnorm(utility[inside]), 0),
    col = shades[k], border = NA
  )
}
lines(utility, dnorm(utility), lwd = 2)
axis(1, at = c(-3, 0, 3))
legend(
  "topright", legend = smoking_levels, fill = shades, border = NA, bty = "n"
)

A covariate shifts the density along the utility axis, so one coefficient per covariate describes its effect on all four levels. The age coefficient is negative: older students report smoking less. The factor Exer enters through its dummy variables relative to the students who exercise frequently. interpret(type = "ame") translates the age coefficient into probabilities: it differentiates the probability of each level with respect to age and averages the derivatives over the students. The vignette Posterior prediction explains marginal effects in more detail.

age_effects <- interpret(smoking, type = "ame")
age_effects

One more year of age raises the probability of never smoking by about r round(100 * age_effects$mean[age_effects$alternative == "Never"], 1) percentage points and lowers the probabilities of the other three levels.

Ranked responses

Ranked data record the complete order of the alternatives. In wide format, the response columns combine the response name with each alternative, for example rank_Xbox. The model is the same probit model as for unordered choices, but the likelihood uses the full ordering of the utilities. The following demonstration simulates rankings of three alternatives and compares the coefficient and the free covariance parameters with the true values.

ranked_sim <- fit(
  rank ~ x | 0,
  choice_type = "ranked",
  n_deciders = 300,
  dgp_parameters = list(
    beta = c(x = 1),
    Sigma = rbind(c(0, 0, 0), c(0, 1, 0.2), c(0, 0.2, 1))
  ),
  chains = 1
)
summary(ranked_sim, variables = c("beta[x]", "Sigma[C,B]", "Sigma[C,C]"))

The Game data of the mlogit package contain complete rankings of six gaming platforms by 91 Dutch respondents, together with whether they own each platform (own), their age, and their weekly gaming hours; the source study develops a rank-ordered choice model for these data [@Fok2012]. The ranks are stored in the columns ch.Xbox, ch.PlayStation, and so on, so the response in the formula is ch and delimiter = "." separates it from the alternative.

data("Game", package = "mlogit")
gaming <- fit(
  ch ~ own | age + hours,
  data = Game,
  alternatives = c(
    "Xbox", "PlayStation", "PSPortable", "GameCube", "GameBoy", "PC"
  ),
  choice_type = "ranked",
  delimiter = ".",
  column_decider = NULL,
  iterations = 1000,
  warmup = 500,
  thin = 20,
  chains = 2,
  progress = FALSE
)
coef(gaming)[1:6]

The coefficient of own is positive: owning a platform raises its rank. The alternative-specific constants and the coefficients of age and hours are relative to the base alternative GameBoy. Which platform gains from additional gaming hours? For a ranked model, the marginal effects of interpret() refer to the probability of being ranked first, here for a respondent of average age who plays the average number of hours:

platform_effects <- interpret(gaming, type = "mea")
platform_effects[platform_effects$covariate == "hours", ]

pc_hours <- platform_effects$mean[
  platform_effects$covariate == "hours" & platform_effects$alternative == "PC"
]

Heavy gamers prefer the PC: every additional weekly hour raises the probability of ranking it first by about r round(100 * pc_hours, 1) percentage points.