Obtaining predicted values of the outcome variable using margins command

Question

I have fit a multinomial logistic model as shown below. Using the margins command, I would like to obtain the predicted values for my outcome variable while setting the predictors to some specific values. In other words, I would like to replicate the Stata code below. Chatgpt gave me code that I pasted below as well but it returns an error of "Error in find_terms_in_model.default(model, variables = variables): Some values in 'variables' are not in the model terms."

my multinomial logist code

library(nnet)
model <- multinom(obesity_E ~ age_100 + I(age_100^2) + obesity + age_100:obesity +
                  covid:race + age_100:race + education + education:race +
                  rabplace_5, data = female_98_2020, maxit=1000)

Code suggested by ChatGPT

margins_model <- margins(model, variables = "obesity:race:rabplace_5",
                         at = list(age_100 = seq(0, 25, 1)), atMethod = "mean",
                         method = "probs", force = TRUE, noSe = TRUE,
                         save = "tran_point_F", replace = TRUE)

My Stata code

margins  , at (age_100=(0 (1) 25) obesity=(1 (1) 4) race=(0 (1) 3) rabplace_5=(1 (1) 5)) atmeans force nose saving(tran_point_F, replace ) 

Answer 1

The margins package only computes slopes, not predictions. Also, that package is not being actively maintained and developed anymore (perhaps only critical bug fixes).

You can use the marginaleffects package instead, which is a newer package I developed as a more flexible successor to margins (note: conflict of interest).

See this vignette for a syntax comparison with Stata: https://vincentarelbundock.github.io/marginaleffects/articles/alternative_software.html

The code you need will probably look similar to this:

library(nnet)
library(marginaleffects)

m <- multinom(Species ~ Petal.Length + Petal.Width, data = iris, trace = FALSE)

avg_predictions(m, newdata = datagrid(Petal.Length = c(1, 4, 6)))
# 
#       Group Estimate Std. Error    z Pr(>|z|)    S  2.5 % 97.5 % Petal.Width Petal.Length
#  setosa        0.333     0.0479 6.95  < 0.001 38.0 0.2389  0.427         1.2            1
#  versicolor    0.382     0.1195 3.20  0.00139  9.5 0.1477  0.616         1.2            4
#  virginica     0.285     0.1192 2.39  0.01672  5.9 0.0516  0.519         1.2            6
# 
# Columns: rowid, group, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, Species, Petal.Width, Petal.Length