Nothing
Modeling preference heterogeneity
In RprobitB: Bayesian Probit Choice Modeling
knitr::opts_chunk$set(
collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4,
fig.align = "center"
)
The probit model assigns alternative $j$ at occasion $t$ of
decider $n$ the latent utility
$U_{ntj} = X_{ntj}^\top \beta_n + \epsilon_{ntj}$ with the coefficient
vector $\beta_n$. A model with fixed coefficients sets $\beta_n = \beta$
for all deciders and thereby assumes that all deciders weigh price, time,
and comfort in the same way. This is rarely true: some travelers react mainly to
the fare, others to travel time, and some households would pay a premium for
a local electricity supplier while others would not. RprobitB lets
coefficients differ between deciders in two ways, which can be combined:
random coefficients follow a continuous distribution over the population,
and latent classes divide the deciders into groups. Both are requested
through arguments of fit(). @Oelschlaeger2026c treats the methodological
background. Each variant is first estimated on simulated data, where the
estimates can be compared with the parameters that generated them, and then
applied to data of the mlogit package [@Croissant2020].
library(RprobitB)
set.seed(1)
Random coefficients
Each term named in random_effects receives a separate coefficient for every
decider, drawn from a population distribution whose parameters the sampler
estimates together with the other parameters. Such hierarchical models
capture continuous preference heterogeneity and allow inference about the
coefficients of individual deciders [@Allenby1998]. The population
distribution is normal on a latent scale: the random coefficients of decider
$n$ are $\beta_n \sim \mathrm{N}(\mu, \Omega)$ with the mean vector $\mu$ and
the covariance matrix $\Omega$, reported as mu[<effect>] and
Omega[<effect>,<effect>]. The mixing distribution of an effect decides how
its latent normal variable enters the utility. An unnamed vector requests
correlated normal effects. A named vector selects the mixing distribution of
each term as "<covariate>" = "<distribution>", where "ASC" names the
alternative-specific constants:
| Value | Distribution | Coefficient |
|:--|:--|:--|
| "cn" | correlated normal | any sign |
| "n" | uncorrelated normal | any sign |
| "cln" | correlated log-normal | positive |
| "ln" | uncorrelated log-normal | positive |
| "cln-" | correlated log-normal | negative |
| "ln-" | uncorrelated log-normal | negative |
The six values combine two choices. The first is the shape of the
distribution. A normal coefficient can take any value, which suits an
attribute that some deciders like and others dislike. A log-normal
coefficient is the exponential of a normal variable and therefore always
positive, or, with the trailing minus, always negative. It suits an attribute
whose sign is fixed by theory: a higher price does not raise the utility of
any decider. Log-normal effects are estimated on the latent normal scale, so
mu, Omega, and the individual draws refer to the normal variable whose
exponential enters the utility. The second choice is whether an effect is
correlated with the other correlated effects. Deciders who value travel time
may also value comfort, and the c prefix adds the covariance between such
effects to the model. An uncorrelated effect has its own variance but no
covariance with any other effect, which appears as zeros in Omega.
The following demonstration combines both choices. The price coefficient is
negative log-normal and uncorrelated. Travel time and comfort
receive correlated normal effects with a positive covariance.
mixing <- fit(
choice ~ price + time + comfort | 0,
random_effects = c(price = "ln-", time = "cn", comfort = "cn"),
n_deciders = 100,
n_occasions = 10,
dgp_parameters = list(
beta = c(price = -1, time = -0.8, comfort = 0.5),
Omega = rbind(c(0.25, 0, 0), c(0, 0.4, 0.2), c(0, 0.2, 0.3))
),
iterations = 4000,
warmup = 2000,
thin = 2,
chains = 1,
save_individual_draws = TRUE
)
summary(mixing)
The dgp column lists the parameters that generated the data. For the
price coefficient, mu[price] is the mean of the latent normal variable, so
the coefficient that enters the utility is minus its exponential and
negative for every decider. time and comfort share the estimated
covariance Omega[time,comfort], while price has no covariance entry with
either of them.
interpret() converts the coefficients into trade-offs. For a random effect
it uses the coefficient of the median decider, which for the log-normal
price is minus the exponential of mu[price]:
interpret(mixing, reference = "price")
The dgp_parameters set the latent mean of the price effect to -1, so
the true median price coefficient is -exp(-1), about -0.37. Dividing the
true time and comfort coefficients by it gives the true trade-offs, with
which the posterior means above can be compared.
The panel structure makes individual coefficients estimable.
coef(level = "individual") returns the posterior mean coefficient of every
decider, on the latent normal scale for log-normal effects. The histogram
shows the price coefficient of each of the 100 deciders, all negative as the
log-normal specification enforces, with the left tail containing the
deciders who react most strongly to a higher price.
individual <- coef(mixing, level = "individual")
head(individual)
price <- -exp(individual[, "price"])
hist(
price,
breaks = 30, col = "grey85", border = "white",
main = "", xlab = "price coefficient of a decider"
)
Willingness to pay for electricity supplier attributes
The Electricity data of the mlogit package come from a stated choice
experiment in which 361 US households chose 8 to 12 times among four
hypothetical suppliers. The suppliers differ in price (pf), contract
length (cl), whether the supplier is local (loc) or well-known (wk),
and whether it offers time-of-day (tod) or seasonal (seas) rates
[@Huber2001]. The attribute columns end in the supplier number without a
delimiter, pf1 to pf4, so an underscore is inserted first, and the
choice occasions of a household are numbered in the order of the rows.
Contract length and locality receive correlated normal random coefficients,
the other attributes one coefficient for all households. Fixing the price
coefficient to -1 identifies the scale and expresses every other
coefficient in cents per kWh, that is, directly as a willingness to pay.
The fit uses the first 100 households to keep the computation short.
data("Electricity", package = "mlogit")
names(Electricity) <- sub("([a-z]+)([1-4])$", "\\1_\\2", names(Electricity))
Electricity$occasion <- ave(Electricity$id, Electricity$id, FUN = seq_along)
households <- Electricity[Electricity$id %in% unique(Electricity$id)[1:100], ]
electricity <- fit(
choice ~ pf + cl + loc + wk + tod + seas | 0,
data = households,
random_effects = c("cl", "loc"),
column_decider = "id",
column_occasion = "occasion",
scale = c(pf = -1),
iterations = 4000,
warmup = 2000,
thin = 2,
chains = 1,
save_individual_draws = TRUE,
progress = FALSE
)
summary(electricity)
Because the price coefficient is fixed, interpret() reports the
coefficients directly as willingness to pay in cents per kWh, with credible
intervals:
interpret(electricity, effects = c("cl", "loc", "wk"))
On average, households would accept a price about
r round(coef(electricity)[["mu[loc]"]], 1) cents per kWh higher for a
local supplier, and would need a price about
r round(abs(coef(electricity)[["mu[cl]"]]), 2) cents lower for each
additional year of contract length.
Latent classes
A single normal distribution is a strong assumption about the distribution
of preferences in a population. There may be commuters and leisure
travelers, or households that focus on price and households that focus on
service. latent_class_effects names the effects that differ between
classes latent classes. Every decider belongs to exactly one class. The
class weights $w_1, \dots, w_K$ are the probabilities of membership and sum
to one, and the class allocation of every decider is a latent variable that
the sampler draws together with the parameters. Naming a random effect in
latent_class_effects replaces its normal distribution by a finite mixture
of normals, which yields the latent-class mixed multinomial probit model of
@Oelschlaeger2021. Random effects that are not named keep one distribution
for all deciders, and a named coefficient that is not random takes one value
per class and does not vary within it, as in the classical latent class
model [@Kamakura1989; @Greene2003]. class_update selects how the number
of classes is treated:
class_update = "fixed" when exactly $K$ substantively meaningful classes
are assumed;class_update = "sparse" when classes is a generous upper bound and
redundant classes should empty;class_update = "dirichlet_process" for a Dirichlet-process mixture whose
number of occupied classes is itself random;class_update = "weight_based" for the split, remove, and merge heuristic.
All four updates are demonstrated on one simulated data set with two
classes. The fixed-class fit simulates the data; the other three are refits
through update(), which reuses the simulated data.
A fixed number of classes
For class_update = "fixed", classes = K fixes the number of classes.
The class weights have the prior
$(w_1,\ldots,w_K)\sim\operatorname{Dirichlet}(\delta,\ldots,\delta)$, where
the default class_concentration = 1 is uniform on the weight simplex. The
class labels are arbitrary: swapping them leaves the likelihood unchanged,
so the sampler may swap them during a run. RprobitB therefore relabels
the retained draws after sampling, so that every draw uses the same labels
[@Dahl2006; @Papastamoulis2010], and numbers the classes by decreasing
weight.
The following demonstration uses eight occasions per decider and two
well-separated classes: a majority with coefficients centered at -1 and a
minority centered at 2.
mixture <- fit(
choice ~ x | 0,
random_effects = "x",
latent_class_effects = "x",
classes = 2,
n_deciders = 100,
n_occasions = 8,
save_individual_draws = TRUE,
dgp_parameters = list(
beta = list(c(x = -1), c(x = 2)),
Omega = list(matrix(0.2), matrix(0.2)),
weights = c(0.6, 0.4)
),
iterations = 1500,
chains = 2,
progress = FALSE
)
summary(mixture, variables = c(
"weight[1]", "weight[2]", "mu[x,1]", "mu[x,2]",
"Omega[x,x,1]", "Omega[x,x,2]"
))
latent_class_diagnostics() returns the posterior distribution of the
number of occupied classes, the membership probabilities of the deciders
after relabeling, and the co-clustering matrix. Its entries are the
posterior probabilities that two deciders belong to the same class:
class_diagnostics <- latent_class_diagnostics(mixture)
class_diagnostics$occupancy
class_diagnostics$membership[1:6, ]
class_diagnostics$co_clustering[1:6, 1:6]
Class-specific coefficients
Do the train travelers of the vignette Get started with RprobitB all
trade time against money at the same rate, or are there classes with
different values of time? The fit below uses the first 60 travelers, fixes
the price coefficient to -1, and gives the time coefficient two
class-specific values through latent_class_effects without a random
effect, so each class has its own value of travel time and the remaining
coefficients are common to both classes.
data("Train", package = "mlogit")
Train$price_A <- Train$price_A / 100 / 2.20371
Train$price_B <- Train$price_B / 100 / 2.20371
Train$time_A <- Train$time_A / 60
Train$time_B <- Train$time_B / 60
train_small <- Train[Train$id %in% unique(Train$id)[1:60], ]
train_classes <- fit(
choice ~ price + time + change + factor(comfort) | 0,
data = train_small,
latent_class_effects = "time",
classes = 2,
column_decider = "id",
column_occasion = "choiceid",
scale = c(price = -1),
iterations = 1500,
chains = 2,
progress = FALSE
)
summary(train_classes, variables = c(
"weight[1]", "weight[2]", "beta[time,1]", "beta[time,2]"
))
With the price coefficient fixed at -1, the class-specific time
coefficients are values of travel time in euro per hour, which interpret()
reports class by class:
time_by_class <- interpret(train_classes, effects = "time")
time_by_class
The larger class values an hour at about
r round(abs(time_by_class$mean[1])) euro, the smaller class at about
r round(abs(time_by_class$mean[2])) euro. The smaller class chooses the
faster trip almost regardless of its price.
Weight-based class updates
@Oelschlaeger2021 presents a weight-based update scheme for latent class
analysis. Every buffer warmup iterations, it removes the
smallest class if its weight is below epsmin, splits the largest class if
its weight is above epsmax, or merges the closest pair of classes if the
distance of their means is below deltamin, at most one operation in this
order. weight_based_control overrides the defaults of these constants,
and max_classes bounds the splitting. These dimension changes correspond
to no prior on the number of classes, and they stop after warmup, so the
reported n_classes is the outcome of a search, not a posterior
distribution.
This refit and the two in the following subsections change the class update,
so summary() has no dgp column for them. The helper recover_classes()
provides the comparison with the true values instead. The classes are
numbered by decreasing weight, so the first class should be the majority
with weight 0.6 and mean -1, and the second class the minority with
weight 0.4 and mean 2. The helper puts the posterior means of these four
variables and the most probable number of occupied classes beside the true
values. Applied to the fit with two fixed classes, it gives the reference
for the refits:
recover_classes <- function(x) {
variables <- c("weight[1]", "weight[2]", "mu[x,1]", "mu[x,2]")
occupancy <- latent_class_diagnostics(x)$occupancy
data.frame(
variable = c("n_classes", variables),
dgp = c(2, 0.6, 0.4, -1, 2),
estimate = round(c(
occupancy$n_classes[which.max(occupancy$probability)],
coef(x)[variables]
), 2),
row.names = NULL
)
}
recover_classes(mixture)
The weight-based refit is compared in the same way:
weight_based <- update(mixture, class_update = "weight_based")
recover_classes(weight_based)
The run ends with the two classes that generated the data, and their weights
and means agree closely with those of the fixed fit.
Sparse finite mixtures
A sparse finite mixture fixes a generous upper bound $K$, permits empty
classes, and places the symmetric Dirichlet prior with a small concentration
$e_0$ on the weights, which favors emptying redundant classes
[@Rousseau2011; @FruehwirthSchnatter2019]. The default
class_concentration = c(shape = 1, rate = 200) is the gamma hyperprior
$e_0\sim\operatorname{Gamma}(1,200)$ with mean 0.005. Under a fixed
$e_0$, the prior expected number of occupied classes among $N$ deciders is
[
K\left[1-
\frac{\Gamma(Ke_0)\,\Gamma((K-1)e_0+N)}
{\Gamma((K-1)e_0)\,\Gamma(Ke_0+N)}\right],
]
which translates $e_0$ into a statement about the number of classes. At the
prior mean $e_0 = 0.005$, with $K = 6$ and the 100 deciders of the simulated
data, the prior expects close to a single occupied class. The refit sets the
upper bound to six classes for the two-class data.
sparse <- update(mixture, classes = 6, class_update = "sparse")
summary(sparse)
latent_class_diagnostics(sparse)$occupancy
recover_classes(sparse)
Although the prior favors a single class, the posterior concentrates on two
occupied classes, whose weights and means are close to the true values.
When the number of classes varies, a class
exists only in part of the draws, and the occupied column of summary()
reports this share.
Dirichlet-process mixtures
For class_update = "dirichlet_process", the precision $\alpha$ controls
the prior tendency to open new classes, with the default
$\alpha\sim\operatorname{Gamma}(2,4)$ of mean 0.5. Conditional on a fixed
$\alpha$, the prior expects $\sum_{i=1}^{N}\alpha/(\alpha+i-1)$ occupied
classes among $N$ deciders, which for $\alpha = 0.5$ and 100 deciders is
about three. RprobitB
updates $\alpha$ with the augmentation of @Escobar1995 and the allocations
with the algorithm of @Neal2000. max_classes caps the number of classes
the sampler may open, and fit() warns when the draws reach the cap.
dynamic <- update(
mixture, class_update = "dirichlet_process", max_classes = 15,
iterations = 1000
)
summary(dynamic)
latent_class_diagnostics(dynamic)$occupancy
recover_classes(dynamic)
The occupancy distribution is wider than under the sparse finite prior and
assigns most of its probability to three or more classes. This reflects the
prior, which expects about three occupied classes for 100 deciders. The two
largest classes nevertheless match the two true groups, because the
superfluous classes contain only a few deciders.
Ordered and ranked responses
Both mechanisms also work for ordered and ranked responses, which the
vignette Model specification and variants introduces. The following
demonstration uses a simulated panel of rankings:
100 deciders order three alternatives five times each, with a coefficient
that varies normally around -1.
ranked_random <- fit(
rank ~ x | 0,
choice_type = "ranked",
random_effects = "x",
n_deciders = 100,
n_occasions = 5,
n_alternatives = 3,
dgp_parameters = list(beta = c(x = -1), Omega = matrix(0.3)),
iterations = 2000,
chains = 1,
progress = FALSE
)
summary(ranked_random)
The summary shows the population mean and variance beside their true values.
Ordered responses accept the same arguments, and latent classes are
requested in the same way as above, with latent_class_effects and
classes.
Further reading
The vignette Posterior prediction uses the individual coefficients to
compute choice probabilities for each decider in the fit. The vignette
Bayesian model evaluation shows how to decide whether random
coefficients improve a model.
References
Try the RprobitB package in your browser
Any scripts or data that you put into this service are public.
RprobitB documentation built on Sept. 20, 2026, 9:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4, fig.align = "center" )
The probit model assigns alternative $j$ at occasion $t$ of
decider $n$ the latent utility
$U_{ntj} = X_{ntj}^\top \beta_n + \epsilon_{ntj}$ with the coefficient
vector $\beta_n$. A model with fixed coefficients sets $\beta_n = \beta$
for all deciders and thereby assumes that all deciders weigh price, time,
and comfort in the same way. This is rarely true: some travelers react mainly to
the fare, others to travel time, and some households would pay a premium for
a local electricity supplier while others would not. RprobitB lets
coefficients differ between deciders in two ways, which can be combined:
random coefficients follow a continuous distribution over the population,
and latent classes divide the deciders into groups. Both are requested
through arguments of fit(). @Oelschlaeger2026c treats the methodological
background. Each variant is first estimated on simulated data, where the
estimates can be compared with the parameters that generated them, and then
applied to data of the mlogit package [@Croissant2020].
library(RprobitB) set.seed(1)
Random coefficients
Each term named in random_effects receives a separate coefficient for every
decider, drawn from a population distribution whose parameters the sampler
estimates together with the other parameters. Such hierarchical models
capture continuous preference heterogeneity and allow inference about the
coefficients of individual deciders [@Allenby1998]. The population
distribution is normal on a latent scale: the random coefficients of decider
$n$ are $\beta_n \sim \mathrm{N}(\mu, \Omega)$ with the mean vector $\mu$ and
the covariance matrix $\Omega$, reported as mu[<effect>] and
Omega[<effect>,<effect>]. The mixing distribution of an effect decides how
its latent normal variable enters the utility. An unnamed vector requests
correlated normal effects. A named vector selects the mixing distribution of
each term as "<covariate>" = "<distribution>", where "ASC" names the
alternative-specific constants:
| Value | Distribution | Coefficient |
|:--|:--|:--|
| "cn" | correlated normal | any sign |
| "n" | uncorrelated normal | any sign |
| "cln" | correlated log-normal | positive |
| "ln" | uncorrelated log-normal | positive |
| "cln-" | correlated log-normal | negative |
| "ln-" | uncorrelated log-normal | negative |
The six values combine two choices. The first is the shape of the
distribution. A normal coefficient can take any value, which suits an
attribute that some deciders like and others dislike. A log-normal
coefficient is the exponential of a normal variable and therefore always
positive, or, with the trailing minus, always negative. It suits an attribute
whose sign is fixed by theory: a higher price does not raise the utility of
any decider. Log-normal effects are estimated on the latent normal scale, so
mu, Omega, and the individual draws refer to the normal variable whose
exponential enters the utility. The second choice is whether an effect is
correlated with the other correlated effects. Deciders who value travel time
may also value comfort, and the c prefix adds the covariance between such
effects to the model. An uncorrelated effect has its own variance but no
covariance with any other effect, which appears as zeros in Omega.
The following demonstration combines both choices. The price coefficient is negative log-normal and uncorrelated. Travel time and comfort receive correlated normal effects with a positive covariance.
mixing <- fit( choice ~ price + time + comfort | 0, random_effects = c(price = "ln-", time = "cn", comfort = "cn"), n_deciders = 100, n_occasions = 10, dgp_parameters = list( beta = c(price = -1, time = -0.8, comfort = 0.5), Omega = rbind(c(0.25, 0, 0), c(0, 0.4, 0.2), c(0, 0.2, 0.3)) ), iterations = 4000, warmup = 2000, thin = 2, chains = 1, save_individual_draws = TRUE ) summary(mixing)
The dgp column lists the parameters that generated the data. For the
price coefficient, mu[price] is the mean of the latent normal variable, so
the coefficient that enters the utility is minus its exponential and
negative for every decider. time and comfort share the estimated
covariance Omega[time,comfort], while price has no covariance entry with
either of them.
interpret() converts the coefficients into trade-offs. For a random effect
it uses the coefficient of the median decider, which for the log-normal
price is minus the exponential of mu[price]:
interpret(mixing, reference = "price")
The dgp_parameters set the latent mean of the price effect to -1, so
the true median price coefficient is -exp(-1), about -0.37. Dividing the
true time and comfort coefficients by it gives the true trade-offs, with
which the posterior means above can be compared.
The panel structure makes individual coefficients estimable.
coef(level = "individual") returns the posterior mean coefficient of every
decider, on the latent normal scale for log-normal effects. The histogram
shows the price coefficient of each of the 100 deciders, all negative as the
log-normal specification enforces, with the left tail containing the
deciders who react most strongly to a higher price.
individual <- coef(mixing, level = "individual") head(individual) price <- -exp(individual[, "price"]) hist( price, breaks = 30, col = "grey85", border = "white", main = "", xlab = "price coefficient of a decider" )
Willingness to pay for electricity supplier attributes
The Electricity data of the mlogit package come from a stated choice
experiment in which 361 US households chose 8 to 12 times among four
hypothetical suppliers. The suppliers differ in price (pf), contract
length (cl), whether the supplier is local (loc) or well-known (wk),
and whether it offers time-of-day (tod) or seasonal (seas) rates
[@Huber2001]. The attribute columns end in the supplier number without a
delimiter, pf1 to pf4, so an underscore is inserted first, and the
choice occasions of a household are numbered in the order of the rows.
Contract length and locality receive correlated normal random coefficients,
the other attributes one coefficient for all households. Fixing the price
coefficient to -1 identifies the scale and expresses every other
coefficient in cents per kWh, that is, directly as a willingness to pay.
The fit uses the first 100 households to keep the computation short.
data("Electricity", package = "mlogit") names(Electricity) <- sub("([a-z]+)([1-4])$", "\\1_\\2", names(Electricity)) Electricity$occasion <- ave(Electricity$id, Electricity$id, FUN = seq_along) households <- Electricity[Electricity$id %in% unique(Electricity$id)[1:100], ] electricity <- fit( choice ~ pf + cl + loc + wk + tod + seas | 0, data = households, random_effects = c("cl", "loc"), column_decider = "id", column_occasion = "occasion", scale = c(pf = -1), iterations = 4000, warmup = 2000, thin = 2, chains = 1, save_individual_draws = TRUE, progress = FALSE ) summary(electricity)
Because the price coefficient is fixed, interpret() reports the
coefficients directly as willingness to pay in cents per kWh, with credible
intervals:
interpret(electricity, effects = c("cl", "loc", "wk"))
On average, households would accept a price about
r round(coef(electricity)[["mu[loc]"]], 1) cents per kWh higher for a
local supplier, and would need a price about
r round(abs(coef(electricity)[["mu[cl]"]]), 2) cents lower for each
additional year of contract length.
Latent classes
A single normal distribution is a strong assumption about the distribution
of preferences in a population. There may be commuters and leisure
travelers, or households that focus on price and households that focus on
service. latent_class_effects names the effects that differ between
classes latent classes. Every decider belongs to exactly one class. The
class weights $w_1, \dots, w_K$ are the probabilities of membership and sum
to one, and the class allocation of every decider is a latent variable that
the sampler draws together with the parameters. Naming a random effect in
latent_class_effects replaces its normal distribution by a finite mixture
of normals, which yields the latent-class mixed multinomial probit model of
@Oelschlaeger2021. Random effects that are not named keep one distribution
for all deciders, and a named coefficient that is not random takes one value
per class and does not vary within it, as in the classical latent class
model [@Kamakura1989; @Greene2003]. class_update selects how the number
of classes is treated:
class_update = "fixed"when exactly $K$ substantively meaningful classes are assumed;class_update = "sparse"whenclassesis a generous upper bound and redundant classes should empty;class_update = "dirichlet_process"for a Dirichlet-process mixture whose number of occupied classes is itself random;class_update = "weight_based"for the split, remove, and merge heuristic.
All four updates are demonstrated on one simulated data set with two
classes. The fixed-class fit simulates the data; the other three are refits
through update(), which reuses the simulated data.
A fixed number of classes
For class_update = "fixed", classes = K fixes the number of classes.
The class weights have the prior
$(w_1,\ldots,w_K)\sim\operatorname{Dirichlet}(\delta,\ldots,\delta)$, where
the default class_concentration = 1 is uniform on the weight simplex. The
class labels are arbitrary: swapping them leaves the likelihood unchanged,
so the sampler may swap them during a run. RprobitB therefore relabels
the retained draws after sampling, so that every draw uses the same labels
[@Dahl2006; @Papastamoulis2010], and numbers the classes by decreasing
weight.
The following demonstration uses eight occasions per decider and two
well-separated classes: a majority with coefficients centered at -1 and a
minority centered at 2.
mixture <- fit( choice ~ x | 0, random_effects = "x", latent_class_effects = "x", classes = 2, n_deciders = 100, n_occasions = 8, save_individual_draws = TRUE, dgp_parameters = list( beta = list(c(x = -1), c(x = 2)), Omega = list(matrix(0.2), matrix(0.2)), weights = c(0.6, 0.4) ), iterations = 1500, chains = 2, progress = FALSE ) summary(mixture, variables = c( "weight[1]", "weight[2]", "mu[x,1]", "mu[x,2]", "Omega[x,x,1]", "Omega[x,x,2]" ))
latent_class_diagnostics() returns the posterior distribution of the
number of occupied classes, the membership probabilities of the deciders
after relabeling, and the co-clustering matrix. Its entries are the
posterior probabilities that two deciders belong to the same class:
class_diagnostics <- latent_class_diagnostics(mixture) class_diagnostics$occupancy class_diagnostics$membership[1:6, ] class_diagnostics$co_clustering[1:6, 1:6]
Class-specific coefficients
Do the train travelers of the vignette Get started with RprobitB all
trade time against money at the same rate, or are there classes with
different values of time? The fit below uses the first 60 travelers, fixes
the price coefficient to -1, and gives the time coefficient two
class-specific values through latent_class_effects without a random
effect, so each class has its own value of travel time and the remaining
coefficients are common to both classes.
data("Train", package = "mlogit") Train$price_A <- Train$price_A / 100 / 2.20371 Train$price_B <- Train$price_B / 100 / 2.20371 Train$time_A <- Train$time_A / 60 Train$time_B <- Train$time_B / 60 train_small <- Train[Train$id %in% unique(Train$id)[1:60], ] train_classes <- fit( choice ~ price + time + change + factor(comfort) | 0, data = train_small, latent_class_effects = "time", classes = 2, column_decider = "id", column_occasion = "choiceid", scale = c(price = -1), iterations = 1500, chains = 2, progress = FALSE ) summary(train_classes, variables = c( "weight[1]", "weight[2]", "beta[time,1]", "beta[time,2]" ))
With the price coefficient fixed at -1, the class-specific time
coefficients are values of travel time in euro per hour, which interpret()
reports class by class:
time_by_class <- interpret(train_classes, effects = "time") time_by_class
The larger class values an hour at about
r round(abs(time_by_class$mean[1])) euro, the smaller class at about
r round(abs(time_by_class$mean[2])) euro. The smaller class chooses the
faster trip almost regardless of its price.
Weight-based class updates
@Oelschlaeger2021 presents a weight-based update scheme for latent class
analysis. Every buffer warmup iterations, it removes the
smallest class if its weight is below epsmin, splits the largest class if
its weight is above epsmax, or merges the closest pair of classes if the
distance of their means is below deltamin, at most one operation in this
order. weight_based_control overrides the defaults of these constants,
and max_classes bounds the splitting. These dimension changes correspond
to no prior on the number of classes, and they stop after warmup, so the
reported n_classes is the outcome of a search, not a posterior
distribution.
This refit and the two in the following subsections change the class update,
so summary() has no dgp column for them. The helper recover_classes()
provides the comparison with the true values instead. The classes are
numbered by decreasing weight, so the first class should be the majority
with weight 0.6 and mean -1, and the second class the minority with
weight 0.4 and mean 2. The helper puts the posterior means of these four
variables and the most probable number of occupied classes beside the true
values. Applied to the fit with two fixed classes, it gives the reference
for the refits:
recover_classes <- function(x) { variables <- c("weight[1]", "weight[2]", "mu[x,1]", "mu[x,2]") occupancy <- latent_class_diagnostics(x)$occupancy data.frame( variable = c("n_classes", variables), dgp = c(2, 0.6, 0.4, -1, 2), estimate = round(c( occupancy$n_classes[which.max(occupancy$probability)], coef(x)[variables] ), 2), row.names = NULL ) } recover_classes(mixture)
The weight-based refit is compared in the same way:
weight_based <- update(mixture, class_update = "weight_based") recover_classes(weight_based)
The run ends with the two classes that generated the data, and their weights and means agree closely with those of the fixed fit.
Sparse finite mixtures
A sparse finite mixture fixes a generous upper bound $K$, permits empty
classes, and places the symmetric Dirichlet prior with a small concentration
$e_0$ on the weights, which favors emptying redundant classes
[@Rousseau2011; @FruehwirthSchnatter2019]. The default
class_concentration = c(shape = 1, rate = 200) is the gamma hyperprior
$e_0\sim\operatorname{Gamma}(1,200)$ with mean 0.005. Under a fixed
$e_0$, the prior expected number of occupied classes among $N$ deciders is
[ K\left[1- \frac{\Gamma(Ke_0)\,\Gamma((K-1)e_0+N)} {\Gamma((K-1)e_0)\,\Gamma(Ke_0+N)}\right], ]
which translates $e_0$ into a statement about the number of classes. At the prior mean $e_0 = 0.005$, with $K = 6$ and the 100 deciders of the simulated data, the prior expects close to a single occupied class. The refit sets the upper bound to six classes for the two-class data.
sparse <- update(mixture, classes = 6, class_update = "sparse") summary(sparse) latent_class_diagnostics(sparse)$occupancy recover_classes(sparse)
Although the prior favors a single class, the posterior concentrates on two
occupied classes, whose weights and means are close to the true values.
When the number of classes varies, a class
exists only in part of the draws, and the occupied column of summary()
reports this share.
Dirichlet-process mixtures
For class_update = "dirichlet_process", the precision $\alpha$ controls
the prior tendency to open new classes, with the default
$\alpha\sim\operatorname{Gamma}(2,4)$ of mean 0.5. Conditional on a fixed
$\alpha$, the prior expects $\sum_{i=1}^{N}\alpha/(\alpha+i-1)$ occupied
classes among $N$ deciders, which for $\alpha = 0.5$ and 100 deciders is
about three. RprobitB
updates $\alpha$ with the augmentation of @Escobar1995 and the allocations
with the algorithm of @Neal2000. max_classes caps the number of classes
the sampler may open, and fit() warns when the draws reach the cap.
dynamic <- update( mixture, class_update = "dirichlet_process", max_classes = 15, iterations = 1000 ) summary(dynamic) latent_class_diagnostics(dynamic)$occupancy recover_classes(dynamic)
The occupancy distribution is wider than under the sparse finite prior and assigns most of its probability to three or more classes. This reflects the prior, which expects about three occupied classes for 100 deciders. The two largest classes nevertheless match the two true groups, because the superfluous classes contain only a few deciders.
Ordered and ranked responses
Both mechanisms also work for ordered and ranked responses, which the
vignette Model specification and variants introduces. The following
demonstration uses a simulated panel of rankings:
100 deciders order three alternatives five times each, with a coefficient
that varies normally around -1.
ranked_random <- fit( rank ~ x | 0, choice_type = "ranked", random_effects = "x", n_deciders = 100, n_occasions = 5, n_alternatives = 3, dgp_parameters = list(beta = c(x = -1), Omega = matrix(0.3)), iterations = 2000, chains = 1, progress = FALSE ) summary(ranked_random)
The summary shows the population mean and variance beside their true values.
Ordered responses accept the same arguments, and latent classes are
requested in the same way as above, with latent_class_effects and
classes.
Further reading
The vignette Posterior prediction uses the individual coefficients to compute choice probabilities for each decider in the fit. The vignette Bayesian model evaluation shows how to decide whether random coefficients improve a model.
References
Try the RprobitB package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.