Hypothesis testing with multinomial model

Interfaces brms
1

Hi,
I have looked across so many questions about multinomial models, but in the end the models in question just do not resemble the one I have at hand.
I understand that each coefficient in the model output is a difference in log-odds from my reference level. And I think my reference level is the response ‘alive_good’ at sex=1???
This is the model I am fitting, and it is a train data that I made from an existing dataset.

#data is from package 'survival'
lung_new <- lung %>%
  mutate(
    alive_good = ifelse(status == 1 & ph.karno >= 80, 1, 0),  # Alive with good performance status
    dead = ifelse(status == 2, 1, 0),                         # Dead
    alive_sick = ifelse(status == 1 & ph.karno < 80, 1, 0) ,
    time_c=as.factor(time),
    sex=as.factor(sex)# Alive with poor performance status
  ) %>%
  select(-status)%>%group_by(time) %>%
  mutate(unitid = row_number()) %>%
  ungroup()%>%arrange(unitid,time)%>%mutate(y=cbind(alive_good,dead,alive_sick))

brfit<-brm(y|trials(1)~1+sex*time+(1|inst),data=lung_new, family=multinomial())

summary(brfit)
 Family: multinomial 
  Links: mudead = logit; mualivesick = logit 
Formula: y | trials(1) ~ 1 + sex * time + (1 | inst) 
   Data: lung_new (Number of observations: 227) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~inst (Number of levels: 18) 
                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(mudead_Intercept)          0.29      0.23     0.01     0.86 1.00     1500
sd(mualivesick_Intercept)     1.39      0.91     0.09     3.56 1.01     1182
                          Tail_ESS
sd(mudead_Intercept)          1529
sd(mualivesick_Intercept)     1472

Regression Coefficients:
                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
mudead_Intercept          2.15      0.41     1.38     2.95 1.00     4962     3686
mualivesick_Intercept    -5.59      2.24   -10.67    -2.05 1.00     1920     1477
mudead_sex2              -1.35      0.61    -2.53    -0.17 1.00     5628     2945
mudead_time              -0.00      0.00    -0.00     0.00 1.00     4441     3516
mudead_sex2:time          0.00      0.00    -0.00     0.00 1.00     5069     3488
mualivesick_sex2          2.73      2.27    -1.21     7.74 1.00     2391     1879
mualivesick_time          0.00      0.00    -0.00     0.01 1.00     2321     1867
mualivesick_sex2:time    -0.00      0.00    -0.01     0.00 1.00     2681     2201

My first object is to get estimates for my response variable at the levels of my explanatory variables. Just list in Table from < Multinomial analysis of behavior: statistical methods | Behavioral Ecology and Sociobiology>

I think I am in right direction by computing predictions of the estimates. I also add posterior_draws() for plotting.

predictions(brfit1, newdata=datagrid(sex=c(1,2),time=250),
            by="sex", re_formula=NA)

My second objective here is to test hypothesis such as:
Q1- Is Alive different from Dead when at sex==1? (I am ignoring the interaction here, just to know how to do it)
Q2- Does sex==1 differ from sex==2 for the dead, alive_good, alive_sick responses? (I am ignoring the interaction here)
Q3- For the interaction, it would be good to be able to estimate the slopes for each level. And also, to plot the ‘distribution of the difference’ as in < 7 Interactions | Statistical rethinking with brms, ggplot2, and the tidyverse>. Or maybe even using hypothesis() just to test if the slopes are different from zero.

How does the hypothesis test work for this type of multinomial model? I do not know how to include the different response variables in the hypothesis testing.

Thank you for your help,

2

Yes, I think your predictions() approach is sound for the table.

For your Q1, I’d switch to a tidybayes-based approach, which gives you access to add_epred_draws(), compare_levels(), and median_qi(). Here’s what that might look like:

# Load
library(tidybayes)

# Compute
data.frame(sex = 1,
           time = 250) %>% 
  add_epred_draws(brfit, re_formula = NA) %>% 
  compare_levels(variable = .epred, by = .category) %>% 
  median_qi()

# A tibble: 3 × 9
    sex  time .category               .epred .lower  .upper .width .point .interval
  <dbl> <dbl> <chr>                    <dbl>  <dbl>   <dbl>  <dbl> <chr>  <chr>    
1     1   250 alive_sick - alive_good -0.154 -0.232 -0.0930   0.95 median qi       
2     1   250 alive_sick - dead       -0.834 -0.897 -0.754    0.95 median qi       
3     1   250 dead - alive_good        0.680  0.525  0.799    0.95 median qi

Here’s how you might use the tidybayes framework for your Q2.

data.frame(sex = 1:2,
           time = 250) %>% 
  add_epred_draws(brfit, re_formula = NA) %>% 
  compare_levels(variable = .epred, by = sex) %>% 
  median_qi()

# A tibble: 3 × 9
  sex    time .category   .epred    .lower  .upper .width .point .interval
  <chr> <dbl> <fct>        <dbl>     <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 2 - 1   250 alive_good  0.184   0.0588    0.319    0.95 median qi       
2 2 - 1   250 dead       -0.210  -0.347    -0.0815   0.95 median qi       
3 2 - 1   250 alive_sick  0.0219  0.000471  0.0884   0.95 median qi

As to your Q3, I’m not clear on what slopes you’re talking about. Maybe flesh that out a bit.

2 Likes
3

After going through this code and generally looking at more examples of analyses with brms.
My idea was to look at the interaction main effect, initially. As I would normally approach whether I would consider the interaction relevant or not, and then go to an ANCOVA approach.
However, in brms/ Bayesian style I think the problem is directly approached by testing against the model with no interaction. And if there is an interaction, the approach would be to look at the contrasts directly if the interaction was categorical*categorical. Am I right?
I am still trying to wrap my head around the whole Bayesian inference approach as well as this type of model specifically. Sorry if it is a pretty basic question.

I am still curious in how to test the slopes of a categoricalcontinuous, as for example 'alive_goodtime at sex2 =0’, is the slope different from zero?

Cheers,

4

You’ve referenced hypothesis tests a few times in this thread. You could continue working within that paradigm if you want, but IMO the Bayesian framework also allows you to drop it. An alternative is to just fit a model and then compute the various parameters and contrasts of interest. If your theory or your experimental procedure indicate there should be an interaction term, just include the interaction term and find meaningful ways to express that interaction as contrasts, slopes, and so on.

But sure, you could also fit a simpler model that does not include the interaction term and compare the two models by their WAIC or LOO values.

1 Like