Skip to content

Instantly share code, notes, and snippets.

@andrewheiss
Created April 16, 2024 02:11
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save andrewheiss/18a979eafcabe636ea0cac7afe2229a1 to your computer and use it in GitHub Desktop.
Save andrewheiss/18a979eafcabe636ea0cac7afe2229a1 to your computer and use it in GitHub Desktop.
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
)


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

preds_chocolate_mlogit %>% 
  group_by(dark) %>% 
  summarize(avg_pred = mean(estimate))
#> # A tibble: 2 × 2
#>   dark  avg_pred
#>   <fct>    <dbl>
#> 1 Milk    0.0500
#> 2 Dark    0.200


# Automatic predictions and averaging
by <- data.frame(dark = c("Milk", "Dark"), by = c("Milk", "Dark"))

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

This used to work though and would show these results:

## 
##  Estimate Std. Error    z Pr(>|z|)  2.5 % 97.5 %   By
##      0.20     0.0316 6.32   <0.001  0.138  0.262 Dark
##      0.05     0.0316 1.58    0.114 -0.012  0.112 Milk
## 
## Columns: estimate, std.error, statistic, p.value, conf.low, conf.high, by
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment