Wrapping nnet::multinom() + ggeffects::ggemmeans() in a custom function fails though regular code works: object of type 'symbol' is not subsettable

I want to fit a multinomial model with nnet::multinom() and get predictions with ggeffects::ggemmeans(). Whereas such procedure works in regular code, I fail to wrap this in a function.

Example

Data

library(dplyr)

my_mtcars <- 
  mtcars %>%
  mutate(across(c(vs, carb), as.factor)) %>%
  as_tibble()

Fitting and predicting works in the following way

library(nnet) # 7.3-15
library(emmeans) # 1.5.4
library(ggeffects) # 1.0.2

m <- multinom(carb ~ vs, data = my_mtcars)
ggemmeans(model = m, terms = "vs")

## # Predicted probabilities of carb
## # x = vs

## # Response Level = 1

## x | Predicted |       95% CI
## ----------------------------
## 0 |      0.00 | [0.00, 0.00]
## 1 |      0.50 | [0.43, 0.57]

## # Response Level = 2

## x | Predicted |       95% CI
## ----------------------------
## 0 |      0.28 | [0.24, 0.32]
## 1 |      0.36 | [0.30, 0.42]

## # Response Level = 3

## x | Predicted |       95% CI
## ----------------------------
## 0 |      0.17 | [0.14, 0.19]
## 1 |      0.00 | [0.00, 0.00]

## # Response Level = 4

## x | Predicted |       95% CI
## ----------------------------
## 0 |      0.44 | [0.39, 0.50]
## 1 |      0.14 | [0.12, 0.17]

## # Response Level = 6

## x | Predicted |       95% CI
## ----------------------------
## 0 |      0.06 | [0.05, 0.06]
## 1 |      0.00 | [0.00, 0.00]

## # Response Level = 8

## x | Predicted |       95% CI
## ----------------------------
## 0 |      0.06 | [0.05, 0.06]
## 1 |      0.00 | [0.00, 0.00]

But when I try to wrap this procedure in a custom function it fails

my_multinom <- function(dat, dv, expl) {
  
  frmla <- as.formula(paste0(dv, "~", expl))
  
  model_fit <- nnet::multinom(frmla, data = dat)
  ggemmeans(model = model_fit, terms = expl)
}

my_multinom(dat = my_mtcars, dv = "carb", expl = "vs")

Error in object$call$formula[[2]] :
object of type 'symbol' is not subsettable

Notably, it seems that the problem lies in the interaction between multinom() and ggemmeans(). If we omit ggemmeans() from my_multinom() then it seems to work OK:

my_multinom_no_ggemmeans <- function(dat, dv, expl) {
  
  frmla <- as.formula(paste0(dv, "~", expl))
  model_fit <- nnet::multinom(frmla, data = dat)
  model_fit
}

my_multinom_no_ggemmeans(dat = my_mtcars, dv = "carb", expl = "vs")

## # weights:  18 (10 variable)
## initial  value 57.336303 
## iter  10 value 38.192450
## iter  20 value 37.940409
## final  value 37.940164 
## converged
## Call:
## nnet::multinom(formula = frmla, data = dat)

## Coefficients:
##   (Intercept)       vs1
## 2    13.44961 -13.78607
## 3    12.93879 -33.99280
## 4    13.91961 -15.17237
## 6    11.84015 -23.96194
## 8    11.84015 -23.96194

## Residual Deviance: 75.88033 
## AIC: 95.88033 

Any idea why my_multinom() wrapper fails?

---

UPDATE

---

I may have found a solution but I don't understand why it works. Based on this github issue (a different package), I've adapted the following solution:

my_multinom_with_do.call <- function(dat, dv, expl) {

  frmla <- as.formula(paste0(dv, "~", expl))

  model_fit <- do.call(multinom, args = list(formula = frmla, data = dat))
  ggemmeans(model = model_fit, terms = expl)
}

And it works:

my_multinom_with_do.call(dat = my_mtcars, dv = "carb", expl = "vs")

But why this works whereas my original my_multinom() didn't?

It doesn't work because of lazy evaluation. The call member of model_fit has formula = frmla, unevaluated. The emmeans support for that model is expecting a formula there. It will work if you add a line to the original function:

my_multinom <- function(dat, dv, expl) {
    
    frmla <- as.formula(paste0(dv, "~", expl))
    
    model_fit <- nnet::multinom(frmla, data = dat)
    model_fit$call$formula <- frmla
    ggemmeans(model = model_fit, terms = expl)
}

The reason the do.call method does work is that frmla is evaluated when you create the list that is passed to do.call.