blog/2023/08/12/conjoint-multilevel-multinomial-guide/

[giscus[bot]]

blog/2023/08/12/conjoint-multilevel-multinomial-guide/

Learn how to use R, {brms}, {marginaleffects}, and {tidybayes} to analyze discrete choice conjoint data with fully specified hierarchical multilevel multinomial models

https://www.andrewheiss.com/blog/2023/08/12/conjoint-multilevel-multinomial-guide/

[joscani]

Wow. awesome post.
I'm thinking about it is possible use this kind of technique to model how effective is to offer differents retention offers to our customers in a telco sector.
Example. Offers only for customers who like to change from our company to other.
Choice 1. 20 % discount during 24 months
Choice 2. More tv channels by only 2$
Choice 3. Discount in phone terminal.

Etc..

[pablobernabeu]

Impressive article!

To complement, should any readers want to compare the lme4 and the brms results in the same plot, check out https://pablobernabeu.github.io/2022/why-can-t-we-be-friends-plotting-frequentist-lmertest-and-bayesian-brms-mixed-effects-models

[1531295]

Hi Andrew,

when i am trying to apply the your codes to analyze my own discrete choice experiment, i am struggling with the differences between the clogit() and the mlogit() function in R. After reading the helpfiles, I am still uncertain whether they are the same thing. Especially when it comes to a mixed logit model (choice characteristic and individual characteristics), or the full hierarchical multinomial logit (if i understood it correctly), which function is more appropriate?

My understanding so far is that both functions should deliver similar results and calculation of AMCEs should be unproblematic for both. But i am unsure. Do you have any suggestions or references that i could read? Thank you so much for your help!

[tupolevy]

Thank you for an excellent post , trying to replicate the example in R , however gets stuck and got an error when trying to predict with mlogit and feeding a new dataframe to it , any ideas n how to solve?

 library(mlogit)
> 
> # Creating 'by' data frame with levels of the 'dark' variable
> #by <- datagrid(dark = c("Milk", "Dark"))
> 
> # Using 'by' in predictions
> preds_chocolate_mlogit <- predictions(
+   model_chocolate_mlogit,
+   by = by
+ )
Error: Unable to compute predicted values with this model. You can try to supply a different dataset to the `newdata` argument. This error was
  also raised:
  
  <text>:2:0: unexpected end of input
1: . ~   . |  . |  . |  
   ^
  
  Bug Tracker: https://github.com/vincentarelbundock/marginaleffects/issues

[andrewheiss]

From this:

<text>:2:0: unexpected end of input
1: . ~   . |  . |  . |  

it looks like you might have some other code that's messing it up, likely in your model creation code. This should work:

library(tidyverse)
library(mlogit)
library(dfidx)
library(marginaleffects)

chocolate <- read_csv("https://www.andrewheiss.com/blog/2023/08/12/conjoint-multilevel-multinomial-guide/data/choco_candy.csv") %>% 
  mutate(
    dark = case_match(dark, 0 ~ "Milk", 1 ~ "Dark"),
    dark = factor(dark, levels = c("Milk", "Dark")),
    soft = case_match(soft, 0 ~ "Chewy", 1 ~ "Soft"),
    soft = factor(soft, levels = c("Chewy", "Soft")),
    nuts = case_match(nuts, 0 ~ "No nuts", 1 ~ "Nuts"),
    nuts = factor(nuts, levels = c("No nuts", "Nuts"))
  )

chocolate_idx <- dfidx(
  chocolate,
  idx = list("subj", "alt"),
  choice = "choice",
  shape = "long"
)

model_chocolate_mlogit <- mlogit(
  choice ~ dark + soft + nuts | 0 | 0, 
  data = chocolate_idx
)

# Using unique()
preds_chocolate_mlogit <- predictions(
  model_chocolate_mlogit, 
  newdata = datagrid(dark = unique, soft = unique, nuts = unique)
)

# Using specific levels
preds_chocolate_mlogit <- predictions(
  model_chocolate_mlogit, 
  newdata = datagrid(dark = c("Milk", "Dark"), soft = unique, nuts = unique)
)
preds_chocolate_mlogit
#> 
#>  Group dark  soft    nuts Estimate Std. Error     z Pr(>|z|)    S    2.5 %  97.5 %
#>      A Milk Chewy No nuts    0.054    0.04335 1.246    0.213  2.2 -0.03096  0.1390
#>      B Milk Chewy Nuts       0.126    0.08490 1.484    0.138  2.9 -0.04039  0.2924
#>      C Milk Soft  No nuts    0.006    0.00743 0.808    0.419  1.3 -0.00856  0.0206
#>      D Milk Soft  Nuts       0.014    0.01622 0.863    0.388  1.4 -0.01780  0.0458
#>      E Dark Chewy No nuts    0.216    0.11212 1.926    0.054  4.2 -0.00375  0.4358
#>      F Dark Chewy Nuts       0.504    0.14163 3.559   <0.001 11.4  0.22641  0.7816
#>      G Dark Soft  No nuts    0.024    0.02583 0.929    0.353  1.5 -0.02663  0.0746
#>      H Dark Soft  Nuts       0.056    0.05509 1.016    0.309  1.7 -0.05198  0.1640
#> Columns: rowid, group, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, choice, dark, soft, nuts 
#> Type:  response

[1531295]

Hi, thank you for the awesome post!
Unfortunately, i had the same problem as the user above when trying to replicate the results in R. Your response has worked - though the error occurs when i try to collapase the eight possible choices into those with milk chocolate and those with dark chocolate, and subsequently calculate the contrast between these two averages (categorical AMCEs). That is, the causal effect of moving from milk to dark, holding all other features constant. I reckon something is wrong with the command:

by <- data.frame(dark = c("Milk", "Dark"), by = c("Milk", "Dark"))

predictions(
model_chocolate_mlogit,
by = by
)

Do you have any hints to solve this issue?

[andrewheiss]

Ah, yeah, I think {marginaleffects} had some new updates regarding how it works with multinomial models. Let me check with Vincent about the new syntax

[andrewheiss]

Turns out it was a bug with {mlogit} and {marginaleffects} (vincentarelbundock/marginaleffects#1086), but it's fixed now! Install the latest development version with

install.packages("marginaleffects", repos = c("https://vincentarelbundock.r-universe.dev"))

…and everything should work correctly:

by <- data.frame(dark = c("Milk", "Dark"), by = c("Milk", "Dark"))

predictions(
  model_chocolate_mlogit,
  by = by
)

#>    By Estimate Std. Error    z Pr(>|z|)    S  2.5 % 97.5 %
#>  Milk     0.05     0.0316 1.58    0.114  3.1 -0.012  0.112
#>  Dark     0.20     0.0316 6.32   <0.001 31.9  0.138  0.262
#> 
#> Columns: estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, by 
#> Type:  response

[1531295]

Hi Andrew, thank you so much for the fast reply! I updated my R, installed the package, and now the codes work! :)

[WillemSleegers]

Thanks for the post Andrew! Going through it now myself to try and understand it all. While doing that I saw that rintro-chapter13conjoint.csv links to the chocolate data and not to the minivan data.

[andrewheiss]

Ah good catch! It's fixed now

[ChrisTwie]

Thanks so much for this wonderful guide, Andrew! Can you advise on how to do it for a two-alternative forced choice task? I.e., subjects made 40 choices each among two unlabeled alternatives that varied on 3 attributes. The alternatives A and B have no inherent meaning, only varied in attribute levels. Can this be be done with mlogit or do you know of a package (besides survival) that can be used for the frequentist approach?
For brms, bernoulli family should do the trick. For a model like the following (I'll withhold the priors etc for now), when I include the fixed intercept, I get a large intercept estimate and only a tiny coefficient for gendermale (no genderfemale coefficient at all), but when I remove the fixed intercept, I get large estimates for gendermale AND genderfemale. Does it even make sense to use the intercept, given that alternatives are unlabeled?
model_9b <- brm(
bf(choice~
# choice-level predictors that are nested within respondents
1 + return + risk + social *SVOangle + age + gender + finKnow + degree +
# with random respondent-specific slops for the nested choice-level predictors
(1 + return + risk + social | choiceid | id)),

[ab-hanifa]

Hey Andrew,
I am trying to model a multilevel multinomial model with a categorical response yij where yij has three category A, B and C.
My data has 20 clusters and each cluster contains 30 individuals and each individual can take any of response A, B and C.

A sample data

A tibble: 600 × 5
cluster cluster_size x1c x2b yij

1 1 1 0.587 1 A
2 1 2 1.04 1 B
3 1 3 0.601 1 A
4 1 4 -2.78 1 C
5 1 5 -4.13 1 C
6 1 6 -0.247 1 C
7 1 7 3.67 1 B
8 1 8 1.93 1 A
9 1 9 0.818 0 C
10 1 10 0.715 1 B

[tslumley]

SAS apparently fits these models with proportional hazard survival-style models, which feels weird, but there’s probably a mathematical or statistical reason for it.

There is! The partial likelihood is actually based on a discrete-choice process.

David Cox noticed that under a proportional hazards model, the number and timing of deaths gave you information about the baseline hazard but gave you very little information about the covariate effects. The partial likelihood takes advantage of this. It says "I'll tell you one person died at this time, but you can use the relative hazard $e^{z_i(t)\beta}$ to model how likely it was that Atropos picked person $i$ in her discrete-choice task"

Folks use Cox-model software to do the calculations because the biostatisticians actually got off their asses and wrote the routines in SAS and S-PLUS and R and Stata.

We can recreate the same proportional hazards model with coxph() from the {survival} package.

You can, but you can also use survival::clogit, which does some of the wrapping work for you, and makes it a little less obvious that you're just repurposing Cox-model software.

[philipshirk]

For what it's worth, fitting the model_minivans_mega_mlm_brms model fails for me. (The only change that I made to your code is that I'm using the default backend instead of backend = "cmdstanr", because cmdstanr is not available for R 4.4.0.) I get an error "anon model does not contain any samples". Adding init = 0 to the brm function allowed the model to fit.

[andrewheiss]

cmdstanr should work with 4.4.x (I'm on 4.4.2 and have been running it fine). The Stan team recently moved all its R development builds to r-universe.dev, so you have to install it from there, but that should be in the cmdstanr readme now: https://mc-stan.org/cmdstanr/

install.packages("cmdstanr", repos = c('https://stan-dev.r-universe.dev', getOption("repos")))