---
title: "Hierarchical Bayes: panel tastes, product effects, and entry"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Hierarchical Bayes: panel tastes, product effects, and entry}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4)
options(digits = 4)
```

choicer's hierarchical Bayesian models put random effects at the two levels
where applied demand work actually needs them. For person $i$ facing
alternative $j$ in choice situation $t$,

$$U_{ijt} = x_{ijt}'\beta_i + \delta_j + \varepsilon_{ijt},
\qquad \beta_i \sim N(b, W),
\qquad \delta_j = z_j'\theta + \xi_j, \quad \xi_j \sim N(0, \sigma_d^2).$$

- **Respondent level.** Each person carries one taste vector $\beta_i$ across
  all of their choice situations. This is the genuine *panel* mixed logit: by
  contrast, `run_mxlogit()` maximizes a cross-sectional simulated likelihood in
  which every choice situation is integrated over tastes separately, with no
  draw held fixed across a person's likelihood contributions (see the
  [mixed logit vignette](mxl.html)). When people are observed repeatedly —
  survey panels, scanner data, repeated hospital admissions — the repetition,
  together with within-person attribute and menu variation, helps distinguish
  persistent taste heterogeneity $W$ from choice-level noise, and
  the hierarchical model is the one that uses it.

- **Alternative level.** Each alternative carries a mean-utility effect
  $\delta_j$, partially pooled toward a characteristics-based mean $z_j'\theta$
  with one scalar variance $\sigma_d^2$ regardless of $J$. The deviation
  $\xi_j$ is the micro-data counterpart of the unobserved product quality in
  Berry, Levinsohn and Pakes (1995). Partial pooling matters when per-alternative
  data are thin — with hundreds of hospitals or schools, fixed
  alternative-specific constants drown in incidental-parameter noise, while
  $\delta_j$ shrinks each alternative toward what its observed characteristics
  predict. And because $\theta$ and $\sigma_d^2$ describe alternatives *as a
  population*, the model has a posterior predictive for alternatives it has
  never seen — which is what makes entry counterfactuals possible below.

Both models require the implicit outside option (`include_outside_option =
TRUE`): the outside good, with systematic utility zero, anchors the *location*
of $\delta$, so $\delta_j$ is mean utility relative to not choosing any inside
alternative. Estimation is by Gibbs sampling in C++ — an adaptive
Metropolis-within-Gibbs for the logit (`run_hmnlogit()`), a fully conjugate
Albert-Chib sampler for the probit (`run_hmnprobit()`).

The hierarchy is shared; the utility-shock model is not:

| | HMNL (`run_hmnlogit`) | HMNP (`run_hmnprobit`) |
|---|---|---|
| Choice shock | iid Type-I extreme value | iid normal in utility levels, including a stochastic outside option |
| Persistent tastes | Normal or lognormal coordinates in `beta_i` | Normal coordinates only |
| Scale | Logit scale fixed by the EV1 distribution | Free expanded `sigma2`; every reported draw is scale-normalized |
| Sampler | Adaptive random-walk Metropolis-within-Gibbs | Fully conjugate Albert-Chib augmentation |
| Welfare | Logsum and consumer surplus | No expected-maximum/logsum welfare implementation |

In particular, HMNP is not the hierarchical analogue of the unrestricted
utility-difference covariance estimated by `run_mnprobit()`. HMNP gains a
conjugate, scalable panel sampler by imposing iid utility-level normal shocks;
its flexible substitution comes from persistent person tastes and pooled
alternative effects, not an unrestricted shock covariance.

```{r setup}
library(choicer)
set_num_threads(2)
```

## Simulate a panel

`simulate_hmnl_data()` draws a panel with known parameters: `N` people, `T`
choice situations each, `J` inside alternatives, structural covariates `x1`
and `x2`, and one alternative-level covariate `z1` feeding the $\delta$ mean
function.

```{r sim}
sim <- simulate_hmnl_data(N = 250, T = 6, J = 6, seed = 42)
head(sim$data, 8)
```

The two-level structure is visible in the identifiers: `pid` indexes people
(who share a taste vector) and `task` indexes choice situations. A task whose
inside rows are all `choice = 0` is one where the outside option was chosen.

## Fit

Point `person_col` at the person identifier — that is what groups choice
situations into a panel. With `person_col = NULL` every situation would be
its own respondent, the cross-sectional limit of the same model. We run two
chains; they feed the convergence diagnostics that `summary()` reports.

```{r fit}
set.seed(99)
fit <- run_hmnlogit(
  data               = sim$data,
  id_col             = "task",
  alt_col            = "alt",
  choice_col         = "choice",
  covariate_cols     = c("x1", "x2"),
  person_col         = "pid",
  alt_covariate_cols = "z1",
  chains             = 2,
  mcmc               = list(R = 24000, burn = 4000, thin = 8)
)
summary(fit)
```

Reading the output from the top: the population mean tastes `b`, the $\delta$
mean function `theta`, and the alternative-effect variance $\sigma_d^2$, each
summarized by their posterior. The *quality ladder* is the per-alternative
posterior of $\delta_j$ and $\xi_j$ — which alternatives deliver more mean
utility than their characteristics predict. The convergence table and
acceptance rates are discussed next.

One storage contract matters when `chains > 1`: the top-level `fit$draws`,
coefficient summaries, quality ladder, and post-estimation methods use chain 1.
`fit$chains` retains the hierarchical draws from every requested chain, and the
convergence table uses all of them. Thus additional chains diagnose whether the
reported chain-1 posterior is reproducible; choicer v0.2.0 does not pool chains
automatically for posterior summaries or policy calculations.

## Convergence is part of the result

MCMC output is only evidence about the posterior if the chains have mixed. The
consolidated table in `summary()` reports, per parameter, the rank-normalized
split R-hat, bulk and tail effective sample sizes, and the Monte Carlo
standard error of the posterior mean (Vehtari et al. 2021), plus one
worst-case row spanning all $J$ alternative effects. The $\delta$ block is
the one to watch: it is updated by a serial random-walk Metropolis sweep (its
full conditionals are coupled through the softmax denominators), so it mixes
more slowly than the conjugate blocks. If its R-hat or ESS looks poor, run
longer — the fit also warns at estimation time when any tracked parameter
fails the check.

Trace plots make the same point visually:

```{r trace-b}
traceplot(fit, block = "b")
```

```{r trace-delta}
traceplot(fit, block = "delta")
```

The same diagnostics are available programmatically, on any block, from the
retained per-chain draws in `fit$chains`:

```{r diagnostics}
b_chains <- lapply(fit$chains, function(ch) ch$b)
rhat(b_chains, rank = TRUE)
ess(b_chains)
mcse(b_chains)
```

Two chains keep this vignette build manageable and permit between-chain checks;
for a serious empirical run, four or more chains are often a better default.
choicer offsets the RNG seed across chains but currently gives them the same
data-driven initialization, so this is not an overdispersed-start diagnostic.
Report the number of people, tasks per person and alternatives;
all prior scales; HMNL acceptance rates; rank-normalized R-hat; bulk and tail
ESS; MCSE; trace plots; and posterior-predictive shares. Show prior sensitivity
for `W` and `sigma_d2`. For entry, defend why the entrant's residual quality is
exchangeable with incumbents rather than treating the posterior predictive as a
data-free ASC.

## Did we recover the truth?

Because the data are simulated, we can line the posterior up against the
generating parameters:

```{r recovery}
recovery_table(fit, sim)
```

The population means `b`, the taste variances `diag(W)`, and the realized
$\delta_j$ ladder recover tightly. Note the pattern in the `theta` and
`sigma_d` rows: with only $J = 6$ alternatives, the mean-function regression
of $\delta$ on $z$ has six observations, so those posteriors are wide and
lean on the prior — more alternatives is what sharpens them. The *level* of
$\delta$ has its own identification story: it is pinned by the outside-option
share, so with a small outside share the level posterior is diffuse while the
cross-alternative *contrasts* stay tight.

## Willingness to pay and predicted shares

Treating `x2` as the price-like attribute, `wtp()` forms the per-draw ratio of
population-mean utility coefficients and summarizes its posterior — a median
and quantile interval rather than a delta-method approximation, so the ratio's
skewness is carried through honestly:

```{r wtp}
wtp(fit, price_var = "x2")
```

`predict()` integrates shares over both the taste distribution and the
chain-1 posterior, returning posterior-predictive intervals (and the outside
share):

```{r shares}
set.seed(7)
predict(fit, n_draws = 200)
```

A posterior-predictive check compares observed choice shares with the model's
predictive distribution — a first-pass reality check on fit:

```{r ppc}
ppc_shares(fit, n_draws = 200)
```

Individual-level tastes are available too: `fit$beta_i` stores per-person
posterior summaries (or full draws with `keep_beta_i = "draws"`), and
`predict(fit, level = "individual")` conditions on them.

## A policy counterfactual

Cut the price `x2` of alternative 1 by 0.25 and re-predict — no refitting.
Reusing the RNG seed makes the baseline and counterfactual share calculations
use the same posterior and taste draws. This common-random-number pairing reduces
Monte Carlo noise in the contrast; it does not make a finite-draw integral exact:

```{r counterfactual}
cf <- sim$data
cf$x2[cf$alt == 1] <- cf$x2[cf$alt == 1] - 0.25

set.seed(7)
base_shares <- predict(fit, n_draws = 200)
set.seed(7)
cf_shares <- predict(fit, newdata = cf, n_draws = 200)

data.frame(
  alternative    = base_shares$alternative,
  baseline       = round(base_shares$share, 3),
  counterfactual = round(cf_shares$share, 3)
)
```

The welfare change is the posterior of the compensating variation — the
logsum difference divided by the marginal utility of income, draw by draw:

```{r welfare}
set.seed(7)
cs <- consumer_surplus(fit, price_var = "x2", newdata = cf, n_draws = 200)
attr(cs, "cv")
```

The three numbers are the lower quantile, posterior median, and upper quantile
of the **sum** of compensating variation over the prediction tasks. Supply
`weights = rep(1 / fit$nobs, fit$nobs)` for an equally weighted mean, or a
substantively justified task-weight vector for another aggregate. The function
does not normalize user-supplied weights.

Unlike a delta-method standard error, this interval carries posterior
uncertainty in the logsum and in the population-mean marginal utility of income
$-\bar\gamma_{price}$ used as the denominator. This is a population-mean money
metric, not the posterior distribution of person-specific compensating
variation when marginal utilities of income differ across people. The taste
distribution still enters the integrated logsum, and for a lognormal price
coordinate it also enters $\bar\gamma_{price} = \exp(b + W_{kk}/2)$. Report
that aggregation choice explicitly in applications with price heterogeneity;
an aggregate fixed-sign denominator does not by itself rule out individual
price coefficients near or across zero.
(`logsum()` and `consumer_surplus()` are available for the hierarchical
*logit* only; the probit expected-maximum counterpart is on the roadmap.)

## An entry counterfactual

The distinctive payoff of the BLP-style alternative level. A model with fixed
alternative-specific constants is silent about an alternative it has never
seen — it has no ASC for the entrant and no principled way to invent one. Here
the entrant is a draw from the estimated population of alternatives: given its
characteristics $z_{\text{new}}$, each posterior draw assigns it
$\delta_{\text{new}} \sim N(z_{\text{new}}'\theta_r,\; \sigma_{d,r}^2)$.

Add the entrant's rows to the data — same layout, a new `alt` label, its `z1`
value, `choice = 0` — and predict:

```{r entry}
entrant <- sim$data[sim$data$alt == 1, ]
entrant$alt <- 99L
entrant$z1 <- 0.4
entrant$choice <- 0L
entry_data <- rbind(sim$data, entrant)

set.seed(7)
predict(fit, newdata = entry_data, n_draws = 200)
```

The entrant takes share from every incumbent and from the outside good, and
its credible interval is typically wider than the incumbents' — appropriately so. The
model knows the entrant's observed characteristics but not its $\xi$, so the
prediction integrates over $\xi_{\text{new}} \sim N(0, \sigma_d^2)$: the
uncertainty about unobserved quality can be a dominant uncertainty about an
entrant, and the posterior predictive exposes it rather than hiding it. The
maintained assumption is exchangeability — the entrant's unobserved quality is
a draw from the same population as the incumbents'.

## Price endogeneity

If a price-like covariate is correlated with the unobserved quality $\xi_j$ —
the classic demand-estimation concern — the estimates are exogenous only
conditional on $Z$. The data preparations accept a control-function residual
(`cf_residual_col`, Petrin and Train 2010): regress price on instruments
outside the package, and pass the first-stage residual so it enters utility as
an ordinary covariate. The package does not run the first stage, and posterior
uncertainty does not propagate first-stage estimation error; joint Bayesian IV
is on the roadmap.

## The probit sibling

`run_hmnprobit()` estimates the same two-level structure with iid normal utility
shocks instead of extreme-value ones. The sampler is fully conjugate
(Albert-Chib data augmentation — no Metropolis step, no acceptance rates to
tune), and because it works in un-differenced utility space, unbalanced choice
sets pose no problem. The probit has a free scale, handled by parameter
expansion: a non-identified $\sigma^2$ chain wanders by design, and every
reported quantity is normalized per draw, so all summaries are on the
identified scale. Do not diagnose convergence from the wandering, unidentified
raw `sigma2` chain alone; diagnose the normalized structural and predictive
quantities that enter the economic conclusions.

```{r hmnp}
simp <- simulate_hmnp_data(N = 250, T = 5, J = 6, seed = 42)

set.seed(99)
fitp <- run_hmnprobit(
  data               = simp$data,
  id_col             = "task",
  alt_col            = "alt",
  choice_col         = "choice",
  covariate_cols     = c("x1", "x2"),
  person_col         = "pid",
  alt_covariate_cols = "z1",
  chains             = 2,
  mcmc               = list(R = 30000, burn = 5000, thin = 10)
)
summary(fitp)
```

The identified probit coefficients live on a different scale from the logit's.
A common variance-matching rule notes that the EV1 shock has standard deviation
$\pi/\sqrt{6} \approx 1.28$ against the probit's 1, and therefore multiplies
probit coefficients by $\pi/\sqrt{6}$ for a rough comparison. This is not an
exact transformation: normal and Type-I extreme-value shocks have different
shapes, and coefficients also need matched utility specifications and data:

```{r scale}
rbind(
  probit         = coef(fitp),
  "logit scale"  = coef(fitp) * pi / sqrt(6)
)
```

The `choicer_hb` post-estimation suite — `predict()` (probabilities by a
deterministic fixed-node one-dimensional Gauss-Hermite approximation), `wtp()`,
`elasticities()`, `diversion_ratios()`, `ppc_shares()`, `recovery_table()` — is
available through the same interfaces on the probit fit. The exceptions are
`logsum()` and `consumer_surplus()`,
which are logit-only as noted above.

## Further reading

The full derivations — priors, the Gibbs sweeps, identification, and the
implementation contract — are in the math companions:
[hierarchical MNL](https://fpcordeiro.github.io/choicer/articles/hierarchical_mnl_math.html) and
[hierarchical MNP](https://fpcordeiro.github.io/choicer/articles/hierarchical_mnp_math.html). For where these
models sit among choicer's estimators, see
[Choosing among choice models](choicer.html#choosing-among-choice-models).

## References

Berry, S., Levinsohn, J. and Pakes, A. (1995). Automobile prices in market
equilibrium. *Econometrica*, 63(4), 841-890.

Petrin, A. and Train, K. (2010). A control function approach to endogeneity in
consumer choice models. *Journal of Marketing Research*, 47(1), 3-13.

Rossi, P. E., Allenby, G. M. and McCulloch, R. (2005). *Bayesian Statistics
and Marketing*. Wiley.

Train, K. E. (2009). *Discrete Choice Methods with Simulation* (2nd ed.).
Cambridge University Press.

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. and Bürkner, P.-C. (2021).
Rank-normalization, folding, and localization: An improved R-hat for assessing
convergence of MCMC. *Bayesian Analysis*, 16(2), 667-718.
