---
title: "Modeling preference heterogeneity"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Modeling preference heterogeneity}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: ../inst/REFERENCES.bib
link-citations: true
---

```{r, include=FALSE}
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].

```{r setup}
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.

```{r mixing}
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]`:

```{r mixing-interpret}
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.

```{r individual}
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.

```{r electricity}
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:

```{r electricity-interpret}
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`.

```{r classes}
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:

```{r class-diagnostics}
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][v01] 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.

```{r train-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:

```{r train-classes-interpret}
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:

```{r recover-classes}
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:

```{r weight-based}
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.

```{r sparse}
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.

```{r dirichlet}
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][v02] 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`.

```{r ranked-random}
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][v04] uses the individual coefficients to
compute choice probabilities for each decider in the fit. The vignette
[Bayesian model evaluation][v05] shows how to decide whether random
coefficients improve a model.

[v01]: https://loelschlaeger.de/RprobitB/articles/v01_get_started.html
[v02]: https://loelschlaeger.de/RprobitB/articles/v02_model_variants.html
[v04]: https://loelschlaeger.de/RprobitB/articles/v04_prediction.html
[v05]: https://loelschlaeger.de/RprobitB/articles/v05_model_evaluation.html

## References