multinomialemmeansr-marginaleffects

Estimate pairwise comparisons in marginaleffects for a multinomial model with an intercept only


There's something that's trivial to do in emmeans but it eludes me in marginaleffects.

Imagine you want to estimate a multinomial logit model with nnet::multinom().

library(tidyverse)
library(nnet)
library(marginaleffects)
library(emmeans)
library(magrittr)

d1 <- tibble(
  dv = letters[1:3] %>% 
    extract(
      1:3 %>% 
        rep(
          c(600, 300, 100)
          )
      ) %>% 
    factor(
      letters[1:3]
    )
  )

mn1 <- multinom(
  dv ~ 1, data = d1
  )

marginaleffects will return the probabilities for each of the categories

mn1  %>% 
  avg_predictions()
 Group Estimate Std. Error    z Pr(>|z|)     S  2.5 % 97.5 %
 a      0.6    0.01549 38.7   <0.001   Inf 0.5696  0.630
 b      0.3    0.01449 20.7   <0.001 313.9 0.2716  0.328
 c      0.1    0.00949 10.5   <0.001  83.9 0.0814  0.119

but not the comparisons

mn1  %>% 
  avg_comparisons()
Error: There is no valid predictor variable. Please change the `variables` 
argument or supply an alternative data frame to the `newdata` argument.

But it's trivial in emmeans

emm_p <- emmeans(mn1, ~ dv, type = "response")
pairs(emm_p)
 contrast estimate     SE df t.ratio p.value
 a - b         0.3 0.0285  2  10.541  0.0162
 a - c         0.5 0.0212  2  23.570  0.0027
 b - c         0.2 0.0190  2  10.541  0.0162

How do I get marginaleffects to report the comparisons for an intercept only multinomial model?


Solution

  • The avg_comparisons() is only appropriate when you are trying to estimate how the predicted outcome changes when one (or more) predictors in your model change. In your model, there are no predictors at all, so you should not use that function.

    Instead, you can use avg_predictions() and its hypothesis argument. This is documented in the manual and in the hypothesis testing chapter of the free online book.

    avg_predictions(mn1, hypothesis = ~revpairwise)