Extract posteriors and compute odds ratio in multinomial regression

Hi! I’am fitting a multinomial logistic regression with a within-subject factor (mimicry) as predictor. My response variable is a 3 level factor (perc) and I would like to see if there is an effect on each perc probability as a function of mimicry.
This is my dataset:

library(tidyverse)
library(tidybayes)
library(brms)

perc <- c("happy", "mixed", "neutral")
mimicry <- c("blocked", "free")
subject <- 1:10
trial <- 1:10

dat <- expand_grid(subject, trial, perc, mimicry)

dat$resp <- sample(perc, nrow(dat), replace = T)

Then I fit the model:

fit_ex <- brm(resp ~ mimicry + (1|subject),
              data = dat,
              family = categorical(link = "logit"))

summary(fit_exp)

Now my model parameters are not directly what I’m looking for because I would like to have the odds ratio between let’s say happy (from perc) in the mimicry|blocked condition and the mimicry|free condition. So using the tidybayes package I have extracted the posteriors and computed the odds ratio as \frac{\frac{p(happy|free)}{1-p(happy|free)}}{\frac{p(happy|blocked)}{1-p(happy|blocked)}}

new_dat <- expand_grid(perc, mimicry)

post <- add_fitted_draws(new_dat, fit_ex, re_formula = NA)

post %>% 
    group_by(mimicry, .category, .draw) %>% 
    summarise(.value = mean(.value)) %>% 
    ungroup() %>% 
    pivot_wider(names_from = mimicry, values_from = .value) %>% 
    mutate(odds_ratio = (blocked/(1-blocked))/(free/(1-free))) %>% 
    group_by(.category) %>% 
    median_qi(odds_ratio)

Is that a correct approach using the posterior distributions?
Thanks

Hi there,

I initially thought something was off with your code but I was thinking of multinomial rather than categorical models (they can be confusing). But for categorical I think you are correct. The key is to compute your derived variable of interest for each draw and then summarise at the end, which is what you did.

I think the only thing I’d change is that your new_dat doesn’t need perc. You can just do

new_dat <- tibble(mimicry=mimicry)
post <- add_fitted_draws(new_dat, fit_ex, re_formula = NA)

This way there’s also no need to call summarise . Depending on your ultimate goal you may also choose to include the uncertainty from the subject effects with something like

new_dat <- expand_grid(mimicry=mimicry, subject=1000) # any new subject level
post <- add_fitted_draws(new_dat, fit_ex, re_formul = NULL, allow_new_levels=TRUE)

Hugo

EDIT: PS, tagging the post with ‘brms’ may be a good idea next time :)

2 Likes

Thanks Hugo! Yes I have realized that multinomial and categorical are a distinct type of stuff. Thanks also for the code improvement suggestion! I would like to ask you a clarification about this part

new_dat <- expand_grid(mimicry=mimicry, subject=1000) # any new subject level

What’s the role of the subject inside the newdata object?

If you wanted to ‘predict’, in a sense, or see how things may look for a new subject you can add a new ‘level’ to the grouping variable and then use one of the sample_new_levels options to account for this uncertainty . In your simulated data it hardly makes a difference but in real data it will likely widen your 95% HDI.

A nice discussion here: https://github.com/paul-buerkner/brms/issues/82#issuecomment-231440994

And a recent update: