I have a multinomial model and am trying to estimate marginal effects for my model predictors, i.e., "how does the probability of choosing option C over A change when X increases 1 unit, or when Z is either group 1 or 2."

My question is first, if there's a best way to do this using the marginaleffects package on a multinom() object, and is it perhaps just with the avg_slopes() function?

I'm wondering this because the case example on the marginaleffects website doesn't use avg_slopes() which I thought would have been easiest. That makes me hesitant to trust my avg_slopes() blindly hence my question. It also uses mlogit() which I don't use so maybe that's why? I'm also a bit new to the marginaleffects package so maybe I'm probably missing something key.

Another layer of confusion though is that, I saw from this old stackoverflow post that the function marginaleffects() could be used--this is now deprecated and the warning message says I should use slopes() instead? There is still the warning that "The standard errors estimated by marginaleffects do not match those produced by Stata for nnet::multinom models. Please be very careful when interpreting the results." which again, makes me hesitant to rely on avg_slopes().

Any help and suggestions would be much appreciated!

Here's a short example from mtcars:

library(nnet)
df <- mtcars %>% 
  mutate(carb = factor(carb),
         gear = factor(gear),
         vs = factor(vs))
m1 <- multinom(gear ~ vs + carb, data = df)
m1_me <- avg_slopes(m1) # are these values reliable?