```
// Index values, hyperprior values, observations, and covariates.
data {
int<lower = 1> R; // Number of respondents.
int<lower = 1> S; // Number of choice tasks.
int<lower = 2> A; // Number of choice alternatives.
int<lower = 1> I; // Number of observation-level covariates.
int<lower = 1> J; // Number of population-level covariates.
real Gamma_mean; // Mean of population-level means.
real<lower=0> Gamma_scale; // Scale of population-level means.
real<lower=0> Omega_shape; // Shape of population-level scale.
real tau_mean; // Mean of population-level scale.
real<lower=0> tau_scale; // Scale of population-level scale.
int<lower = 1, upper = A> Y[R, S]; // Matrix of observations.
matrix[A, I] X[R, S]; // Array of observation-level covariates.
matrix[R, J] Z; // Matrix of population-level covariates.
}
// Parameters and hyperparameters.
parameters {
matrix[J, I] Gamma; // Matrix of population-level hyperparameters.
corr_matrix[I] Omega; // Population model correlation matrix hyperparameters.
vector<lower = 0>[I] tau; // Population model vector of scale hyperparameters.
matrix[R, I] Beta; // Matrix of observation-level parameters.
}
// Hierarchical multinomial logit.
model {
// Hyperpriors.
to_vector(Gamma) ~ normal(Gamma_mean, Gamma_scale);
Omega ~ lkj_corr(Omega_shape);
tau ~ normal(tau_mean, tau_scale);
// Population model and likelihood.
for (r in 1:R) {
Beta[r,] ~ multi_normal(Z[r,] * Gamma, quad_form_diag(Omega, tau));for (s in 1:S) {
Y[r, s] ~ categorical_logit(X[r, s] * Beta[r,]');
}
} }
```

Choice models are common in marketing and other applications where researchers are interested in understanding both the drivers and trade-offs of choice. Since choice is typically manifest as a non-binary discrete outcome and we care about modeling consumer heterogeneity, a hierarchical Bayesian multinomial logit model is our default specification.

In marketing, choice models are often employed in conjunction with conjoint experiments, a survey-based approach to eliciting preferences where respondents choose from sets of hypothetical product alternatives. Because the conjoint experiment produces repeat observations at the respondent level, the groups in the hierarchical model are the respondents themselves.

In this post we will specify a choice model and seek to answer the following question: Does it matter if we use a centered or non-centered parameterization for a hierarchical multinomial logit model? We have previously explored using a non-centered parameterization for hierarchical regression and even argued that a non-centered parameterization should be our default approach in most applications. However, choice models used in practice employ a centered parameterization almost exclusively.

We will specify and compare performance of a centered parameterization and a non-centered parameterization of a hierarchical multinomial logit model estimated using Hamilton Monte Carlo (HMC) via Stan. Stan automates the required computation for many models and HMC is particularly verbose regarding model diagnostics. I am indebted to early tutorials on this topic by Elea Feit and Kevin Van Horn and Jim Savage.

## Centered parameterization

Let’s start with a centered parameterization of the hierarchical multinomial logit.

We’ll save this Stan script as `hmnl_centered.stan`

. In the `data`

block, we have `R`

respondents that each go through `S`

choice tasks where each choice task has `A`

choice alternatives to choose from. Thus there are `R * S`

total observations where each observation takes on a value from `1`

to `A`

, the chosen alternative. At the observation level, the choice alternatives are defined by `I`

covariates. These covariates are the attributes of the choice alternatives under consideration. At the population or group level, the respondents are defined by `J`

covariates. These covariates are used to explain preference heterogeneity across respondents in the population and improve our ability to predict preferences.

The `data`

block also includes the hyperprior values. These define our hyperpriors on the hyperparameters in the population-level model and can be more easily evaluated when specified as part of the `data`

block rather than hard-coded in the `model`

block. Finally, the observations are stored as a matrix `Y`

where each of the `R`

respondents’ `S`

choices are stored as a row while the observation-level covariates `X`

is an `R x S`

array where each element is the `A x I`

matrix of covariates specific to that respondent and choice task.

The `parameters`

block includes the population-level hyperparameters `Gamma`

or the coefficients associated with the population-level covariates `Z`

, the population-level hyperparameters `Omega`

and `tau`

which are the accompanying correlation matrix and scale of the population-level model, and the observation-level parameters `Beta`

. Again, because a conjoint experiment produces repeat observations at the respondent level, we get a set of these `Beta`

parameters for each respondent. These parameters are often referred to as part-worth utilities or preference parameters in the conjoint literature as they describe the preferences each respondent has for each of the attributes that define the choice alternatives. Thus the hyperparameters are estimates of population preference heterogeneity.

The `model`

block includes the specification of hyperpriors, using the hyperprior values specified in the `data`

block, and the hierarchy of models: the population model on the `Beta`

parameters and the observation model or likelihood on the observed choices `Y`

. We use the helper function `to_vector()`

to easily use a `normal()`

hyperprior on the `Gamma`

matrix. By decomposing the covariance matrix of the multivariate normal population model into a correlation matrix and a scale vector means we can use an LKJ hyperprior on the correlation matrix `Omega`

and a normal hyperprior on the scale vector `tau`

. We then reform the covariance matrix using `quad_form_diag()`

as part of the multivariate population model. Note that `categorical_logit()`

is the multinomial logit likelihood and that it must be used within nested for loops since it isn’t vectorized. Also note that since `Beta`

is a matrix where each row is a vector of respondent-specific parameters, the vector output of `multi_normal()`

is transposed as part of the likelihood.

Since Bayesian models are generative, we can translate this Stan script into `data`

and `generated quantities`

blocks and use Stan to generate data for us.

```
// Index values, hyperprior values, and covariates.
data {
int<lower = 1> R; // Number of respondents.
int<lower = 1> S; // Number of choice tasks.
int<lower = 2> A; // Number of choice alternatives.
int<lower = 1> I; // Number of observation-level covariates.
int<lower = 1> J; // Number of population-level covariates.
real Gamma_mean; // Mean of population-level means.
real<lower=0> Gamma_scale; // Scale of population-level means.
real<lower=0> Omega_shape; // Shape of population-level scale.
real tau_df; // Degrees of freedom of population-level scale.
matrix[A, I] X[R, S]; // Array of observation-level covariates.
matrix[R, J] Z; // Matrix of population-level covariates.
}
// Generate data according to the hierarchical multinomial logit.
generated quantities {
int<lower = 1, upper = A> Y[R, S]; // Matrix of observations.
matrix[J, I] Gamma; // Matrix of population-level hyperparameters.
corr_matrix[I] Omega; // Population model correlation matrix hyperparameters.
vector[I] tau; // Population model vector of scale hyperparameters.
matrix[R, I] Beta; // Matrix of observation-level parameters.
// Draw parameter values and generate data.
for (j in 1:J) {
for (i in 1:I) {
Gamma[j, i] = normal_rng(Gamma_mean, Gamma_scale);
}
}for (i in 1:I) {
tau[i] = chi_square_rng(tau_df);
}
Omega = lkj_corr_rng(I, Omega_shape);for (r in 1:R) {
Beta[r,] = multi_normal_rng(Z[r,] * Gamma, quad_form_diag(Omega, tau))';for (s in 1:S) {
Y[r, s] = categorical_logit_rng(X[r, s] * Beta[r,]');
}
} }
```

We’ll save this Stan script as `generate_data.stan`

. While providing the covariate matrices `X`

and `Z`

in `data`

, largely because this is more straightforward to generate in R, we are generating hyperparameter and parameters values in `generated quantities`

. This includes using the LKJ hyperprior to generate the values of the correlation matrix `Omega`

. However, instead of using a normal hyperprior to generate the scale vector hyperparameters `tau`

, we use a Chi-square to ensure positive values. Once again, since `Beta`

is a matrix where each row is a vector of respondent parameters, the vector output of `multi_normal_rng()`

is transposed. Note that while sampling statements like `~ normal()`

from the `model`

block of the estimation code are vectorized, none of the corresponding RNG statements like `= normal_rng()`

are, hence all of the additional for loops in the corresponding `generated quantities`

block in `generate_data.stan`

.

In an R script, let’s load the necessary packages, allow Stan to use as many cores as we have available and save compiled code, and specify the structure of the data and the hyperprior values.

```
# Load packages.
library(tidyverse)
library(rstan)
library(bayesplot)
library(tidybayes)
# Set Stan options.
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
# Specify data and hyperprior values.
<- list(
sim_values R = 500, # Number of respondents.
S = 10, # Number of choice tasks.
A = 4, # Number of choice alternatives.
L = c(3, 4, 5), # Number of levels in each discrete attribute.
I = 12, # Number of estimable attribute levels, including the brand intercept.
J = 3, # Number of population-level covariates.
Gamma_mean = 0, # Mean of population-level means.
Gamma_scale = 5, # Scale of population-level means.
Omega_shape = 2, # Shape of population-level scale.
tau_df = 2 # Degrees of freedom of population-level scale.
)
```

Here `sim_values`

is close to the `data`

input for using `generate_data.stan`

. These are typical dimensions for a conjoint experiment with `R = 500`

respondents, `S = 10`

choice tasks, `A = 4`

choice alternatives, and `I = 12`

estimable attribute levels. One important difference from the `data`

block is that we’ve included an additional variable `L`

that specifies the number of levels in each *discrete* attribute. In conjoint experiments, most of the attributes are discrete. We want to keep track of which of the simulated attributes are discrete, by design appearing before any continuous attributes, so we can easily impose dummy coding by dropping the first level of each attribute. The first attribute is typically the brand of the product and rather than including an intercept, we will allow for all brand levels to be present, thus creating a “brand intercept” where everything for a given product alternative not explained by the associated attributes is included in the brand.

Thus `I`

, the number of *estimable* attribute levels, will be the number of brand levels plus the number of levels for the remaining attribute levels minus the hold-out level for each attribute plus the number of any continuous attributes. Since `I = 12`

and there are three brands, `L[1] = 3`

, we have specified two continuous attributes, as follows.

```
# Array of observation-level covariates.
<- array(
X NA,
dim = c(sim_values$R, sim_values$S, sim_values$A, sim_values$I)
)for (r in 1:sim_values$R) {
for (s in 1:sim_values$S) {
# Discrete predictors.
<- NULL
X_s for (l in 1:length(sim_values$L)) {
<- NULL
X_l for (a in 1:sim_values$A) {
<- matrix(0, nrow = 1, ncol = sim_values$L[l])
X_a 1, sample(seq(1, sim_values$L[l]), 1)] <- 1
X_a[if (l == 1) X_l <- rbind(X_l, X_a)
if (l != 1) X_l <- rbind(X_l, X_a[, -1])
}<- cbind(X_s, X_l)
X_s
}# Continuous predictors.
<- sim_values$I - (sum(sim_values$L) - length(sim_values$L) + 1)
L_n if(L_n != 0) {
<- cbind(X_s, matrix(rnorm(sim_values$A * L_n), ncol = L_n))
X_s
}<- X_s
X[r, s, , ]
}
}$X <- X sim_values
```

We can check the design matrix for the first respondent to ensure we have 12 estimable attribute levels, including the brand intercept and two continuous variables.

`$X[1,1,,] sim_values`

```
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 1 0 0 0 0 1 1 0 0 0 1.4235293 -1.1998511
[2,] 1 0 0 0 0 1 0 0 0 0 -0.7737283 0.3324931
[3,] 0 0 1 1 0 0 0 0 0 0 -0.2326687 1.4514409
[4,] 1 0 0 0 1 0 0 0 0 1 -1.0825479 1.1937547
```

Now we can generate the covariate matrix `Z`

and then generate the remainder of our data according to the hierarchical multinomial logit by calling `generate_data.stan`

.

```
# Matrix of population-level covariates.
<- cbind(
Z rep(1, sim_values$R),
matrix(
runif(sim_values$R * (sim_values$J - 1), min = 2, max = 5),
nrow = sim_values$R
)
)$Z <- Z
sim_values
# Generate data.
<- stan(
sim_data file = here::here("Code", "generate_data.stan"),
data = sim_values,
iter = 1,
chains = 1,
seed = 42,
algorithm = "Fixed_param"
)
# Extract simulated data and parameters.
<- extract(sim_data)$Y[1,,]
sim_Y <- extract(sim_data)$Gamma[1,,]
sim_Gamma <- extract(sim_data)$Omega[1,,]
sim_Omega <- extract(sim_data)$tau[1,]
sim_tau <- extract(sim_data)$Beta[1,,] sim_Beta
```

Now let’s estimate our model by calling `hmnl_centered.stan`

.

```
<- list(
data R = sim_values$R, # Number of respondents.
S = sim_values$S, # Number of choice tasks.
A = sim_values$A, # Number of choice alternatives.
I = sim_values$I, # Number of observation-level covariates.
J = sim_values$J, # Number of population-level covariates.
Gamma_mean = 0, # Mean of population-level means.
Gamma_scale = 5, # Scale of population-level means.
Omega_shape = 2, # Shape of population-level scale.
tau_mean = 0, # Mean of population-level scale.
tau_scale = 5, # Scale of population-level scale.
Y = sim_Y, # Matrix of observations.
X = sim_values$X, # Array of observation-level covariates.
Z = sim_values$Z # Matrix of population-level covariates.
)
<- stan(
fit_centered file = here::here("Code", "hmnl_centered.stan"),
data = data,
iter = 4000,
thin = 2,
seed = 42
)
```

```
Chain 4: Elapsed Time: 21419.5 seconds (Warm-up)
Chain 4: 23963 seconds (Sampling)
Chain 4: 45382.5 seconds (Total)
Chain 4:
Warning messages:
1: There were 311 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
2: There were 3689 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
3: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
4: Examine the pairs() plot to diagnose sampling problems
5: The largest R-hat is NA, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
6: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess
7: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess
```

Estimating the model gives us more than 300 divergent transitions. Recall that a divergent transition or divergence is a unique HMC diagnostic that identifies problems navigating the posterior distribution. These difficulties with posterior geometry are true regardless of the sampler, but HMC makes the issue more transparent. We continue to see divergences even if we set `adapt_delta = 0.99`

as we’ve done with hierarchical models previously. There are also warnings suggesting we run the chains for longer. Even up to `iter = 10000`

, five times the default 2000 draws, we see the same problems.

Regardless, the fact that we run into divergences when using simulated data with typical conjoint experiment dimensions says that in order to produce a posterior geometry that can be navigated, we need to reparameterize our model.

## Non-centered parameterization

A non-centered parameterization re-expresses the population model and likelihood for a hierarchical model by including an intermediate, deterministic transformation on a set of the hyperparameters. The benefit is that difficult dependencies between the layers in the hierarchy are broken, producing a simpler posterior geometry. See a previous post for a walkthrough of the re-parameterization using a simplified population model as part of a hierarchical linear model.

Let’s see a non-centered parameterization of the hierarchical multinomial logit.

```
// Index values, hyperprior values, observations, and covariates.
data {
int<lower = 1> R; // Number of respondents.
int<lower = 1> S; // Number of choice tasks.
int<lower = 2> A; // Number of choice alternatives.
int<lower = 1> I; // Number of observation-level covariates.
int<lower = 1> J; // Number of population-level covariates.
real Gamma_mean; // Mean of population-level means.
real<lower=0> Gamma_scale; // Scale of population-level means.
real<lower=0> Omega_shape; // Shape of population-level scale.
real tau_mean; // Mean of population-level scale.
real<lower=0> tau_scale; // Scale of population-level scale.
int<lower = 1, upper = A> Y[R, S]; // Matrix of observations.
matrix[A, I] X[R, S]; // Array of observation-level covariates.
matrix[R, J] Z; // Matrix of population-level covariates.
}
// Parameters and hyperparameters.
parameters {
matrix[J, I] Gamma; // Matrix of population-level hyperparameters.
corr_matrix[I] Omega; // Population model correlation matrix hyperparameters.
vector<lower = 0>[I] tau; // Population model vector of scale hyperparameters.
matrix[R, I] Delta; // Matrix of non-centered observation-level parameters.
}
// Deterministic transformation.
transformed parameters {
// Matrix of centered observation-level parameters.
matrix[R, I] Beta;
// Non-centered parameterization.
for (r in 1:R) {
Beta[r,] = Z[r,] * Gamma + Delta[r,] * quad_form_diag(Omega, tau);
}
}
// Hierarchical multinomial logit model.
model {
// Hyperpriors.
to_vector(Gamma) ~ normal(Gamma_mean, Gamma_scale);
Omega ~ lkj_corr(Omega_shape);
tau ~ normal(tau_mean, tau_scale);
// Non-centered population model and likelihood.
for (r in 1:R) {
0, 1);
Delta[r,] ~ normal(for (s in 1:S) {
Y[r, s] ~ categorical_logit(X[r, s] * Beta[r,]');
}
} }
```

We’ll save this Stan script as `hmnl_noncentered.stan`

. We now have the `Delta`

hyperparameters in place of `Beta`

in the `parameters`

block, a `transformed parameters`

block that introduces the `Beta`

hyperparameters and imposes the deterministic non-centered parameterization, and now a standard normal population model in the `model`

block.

The generated data doesn’t change, so we can go ahead and estimate our model by calling `hmnl_noncentered.stan`

.

```
<- list(
data R = sim_values$R, # Number of respondents.
S = sim_values$S, # Number of choice tasks.
A = sim_values$A, # Number of choice alternatives.
I = sim_values$I, # Number of observation-level covariates.
J = sim_values$J, # Number of population-level covariates.
Gamma_mean = 0, # Mean of population-level means.
Gamma_scale = 5, # Scale of population-level means.
Omega_shape = 2, # Shape of population-level scale.
tau_mean = 0, # Mean of population-level scale.
tau_scale = 5, # Scale of population-level scale.
Y = sim_Y, # Matrix of observations.
X = sim_values$X, # Array of observation-level covariates.
Z = sim_values$Z # Matrix of population-level covariates.
)
# Calibrate the model.
<- stan(
fit_noncentered file = here::here("Code", "hmnl_noncentered.stan"),
data = data,
iter = 4000,
thin = 2,
seed = 42
)
```

The model runs without divergences! There also aren’t any diagnostic warnings, beyond an initial suggestion for longer chains, hence we added `iter = 4000`

. For convenience, we also thin the posterior draws with `thin = 2`

to reduce the size of the model output and to keep the dimensions consistent with the default 2000 draws.

We can also visually check model performance by first considering the trace plots. Given the number of parameters (`500 x 12`

), we will only consider the hyperparameters.

```
# Check population model trace plots.
<- str_c("Gamma[", 1:data$J, ",", 1, "]")
gamma_string <- str_c("Omega[", 1:data$I, ",", 1, "]")
omega_string <- str_c("tau[", 1:data$I, "]")
tau_string for (i in 2:data$I) {
<- str_c("Gamma[", 1:data$J, ",", i, "]")
gamma_temp <- c(gamma_string, gamma_temp)
gamma_string <- str_c("Omega[", 1:data$I, ",", i, "]")
omega_temp <- c(omega_string, omega_temp)
omega_string
}
# Gamma.
%>%
fit_noncentered mcmc_trace(
pars = gamma_string,
n_warmup = 500,
facet_args = list(
ncol = 4,
labeller = label_parsed
) )
```

```
# Omega.
%>%
fit_noncentered mcmc_trace(
pars = omega_string,
n_warmup = 500,
facet_args = list(
ncol = 4,
labeller = label_parsed
) )
```

Note that the diagonal in a correlation matrix is fixed to 1.

```
# tau.
%>%
fit_noncentered mcmc_trace(
pars = tau_string,
n_warmup = 500,
facet_args = list(
ncol = 4,
labeller = label_parsed
) )
```

Now let’s demonstrate parameter recovery.

```
# Recover Gamma values.
<- tibble(
gamma_values j = sort(rep(1:(data$J), data$I)),
i = rep(1:(data$I), data$J),
.variable = str_c("Gamma", "_", j, "_", i),
values = as.vector(t(matrix(sim_Gamma, ncol = data$I)))
%>%
) select(.variable, values)
%>%
fit_noncentered gather_draws(Gamma[j, i]) %>%
unite(.variable, .variable, j, i) %>%
ggplot(aes(x = .value, y = .variable)) +
geom_halfeyeh(.width = .95) +
geom_vline(aes(xintercept = values), gamma_values, color = "red") +
facet_wrap(
~ .variable,
ncol = 4,
scales = "free"
)
```

```
# Recover Omega values.
<- tibble(
omega_values j = sort(rep(1:(data$I), data$I)),
i = rep(1:(data$I), data$I),
.variable = str_c("Omega", "_", j, "_", i),
values = as.vector(t(matrix(sim_Omega, ncol = data$I)))
%>%
) select(.variable, values)
%>%
fit_noncentered gather_draws(Omega[j, i]) %>%
unite(.variable, .variable, j, i) %>%
ggplot(aes(x = .value, y = .variable)) +
geom_halfeyeh(.width = .95) +
geom_vline(aes(xintercept = values), omega_values, color = "red") +
facet_wrap(
~ .variable,
ncol = 4,
scales = "free"
)
```

Once again, the diagonal in a correlation matrix is fixed at 1.

```
# Recover tau values.
<- tibble(
tau_values i = 1:(data$I),
.variable = str_c("tau", "_", i),
values = as.vector(sim_tau)
%>%
) select(.variable, values)
%>%
fit_noncentered gather_draws(tau[i]) %>%
unite(.variable, .variable, i) %>%
ggplot(aes(x = .value, y = .variable)) +
geom_halfeyeh(.width = .95) +
geom_vline(aes(xintercept = values), tau_values, color = "red") +
facet_wrap(
~ .variable,
ncol = 4,
scales = "free"
)
```

While we have some difficulty recovering the variance parameters in our covariance matrix decomposition, which is typical, we can see that beyond HMC model diagnostics we have both good convergence and parameter recovery for our *many* hyperparameters.

## Final thoughts

Does it matter if we use a centered or non-centered parameterization for a hierarchical multinomial logit model? Yes. Based on HMC model diagnostics and using simulated data typical of a conjoint experiment, the centered parameterization is not sufficient. Additionally, though our primary concern isn’t computational efficiency, with both models running for 4000 iterations, the centered parameterization took 12 and a half hours to run while the non-centered parameterization took 5 hours to run.

This difference in choice model parameterizations is a concern because the centered parameterization is used almost exclusively in practice. However, since HMC isn’t as widespread, the diagnostics accompanying other estimation paradigms might not indicate any problems. This post doesn’t go so far as to establish the edge cases of when either parameterization might work best or consider *how* wrong the centered parameterization is. Regardless, given the savings in computation time and the evidence presented using simulated data, the non-centered parameterization should be our default approach in choice modeling.