I am using the marginaleffects package in R to calculate marginal effects from a multinomial logistic regression model (MNL) that includes an interaction term. I want to calculate the Marginal Effects at the Mean (MEM) for each term in the MNL; the marginaleffects documentation refers me to the slopes function with the argument newdata = means. My MNL contains an interaction term, and the slopes function does not seem to be recognized correctly. Here’s a summary of my problem and what I’ve tried:

My model includes both main effects and an interaction term:

modelCO <- multinom(choice ~ risk + crop_loss + tech_loss + Policy + crop_loss:tech_loss, data = risk_simple_CO)

The interaction term on crop_loss and tech_loss is correctly specified in the model

modelCO[["coefnames"]]
[1] "(Intercept)"         "risk"                "crop_loss"          
[4] "tech_loss"           "Policy"              "crop_loss:tech_loss"

Issue with slopes Function

When I try to calculate the MEM using the slopes() function with newdata = “means”, I receive a warning that the interaction term is not found:

memCO <- slopes(modelCO, newdata = "means", variables = c("risk", "crop_loss", "tech_loss", "Policy", "crop_loss:tech_loss"))

Warning message:
These variables were not found: crop_loss:tech_loss. Try specifying the
  `newdata` argument explicitly and make sure the missing variable is
  included.

I tried manually created the newdata data frame to create the interaction term

`meansCO <- risk_simple_CO %>%
  summarize(
    risk = mean(risk, na.rm = TRUE),
    crop_loss = mean(crop_loss, na.rm = TRUE),
    tech_loss = mean(tech_loss, na.rm = TRUE),
    Policy = mean(Policy, na.rm = TRUE)
  )
meansCO <- meansCO %>%
  mutate(`crop_loss:tech_loss` = crop_loss * tech_loss)
newdata_manual <- as.data.frame(meansCO)
memCO <- slopes(modelCO, newdata = newdata_manual, variables = c("risk", "crop_loss", "tech_loss", "Policy", "crop_loss:tech_loss"))`

However, I still encountered issues with the interaction term not being correctly evaluated. Manually defining newdata does change the MEM estimate from slopes, though.

Why does the slopes function not handle the interaction term correctly when using newdata = “means”?

How can I correctly specify the newdata argument to ensure the interaction term is included and properly evaluated?