使用as.formula()传递公式至emmeans::emmeans()时nnet::multinom()报错

I’ve run into this exact issue before—let’s break down what’s happening and how to fix it with a clean wrapper function.

The Root Cause #

When you pass a formula built with as.formula() to nnet::multinom(), the resulting model object stores the formula as a character string in its call attribute, instead of a proper formula object. emmeans::emmeans() relies on parsing the formula directly from the model’s call to map out factor relationships, and it chokes on the string version. Hardcoding the formula works because the model retains the formula object intact.

Solution: A Wrapper Function for multinom() #

We can build a simple wrapper that fits the multinomial model, then replaces the string formula in the model’s call with the original formula object. This lets emmeans() work as expected.

First, let’s reproduce the error to confirm:

# Load required packages
library(nnet)
library(emmeans)

# Load your dataset (adjust this to match your data source)
dat <- read.csv("your_data.csv")

# Create a dynamic formula with as.formula()
resp_var <- "outcome"
pred_vars <- c("predictor1", "predictor2", "predictor1:predictor2")
form_str <- paste(resp_var, "~", paste(pred_vars, collapse = " + "))
my_formula <- as.formula(form_str)

# Fit model with dynamic formula
mod <- multinom(my_formula, data = dat, trace = FALSE)

# Try emmeans() - this will throw an error!
emmeans(mod, pairwise ~ predictor1 | predictor2)

Now here’s the wrapper function to fix this:

multinom_wrap <- function(formula, data, ...) {
  # Fit the multinomial model as usual
  mod <- multinom(formula = formula, data = data, ...)
  
  # Replace the string formula in the model's call with the original formula object
  mod$call$formula <- formula
  
  # Return the corrected model
  return(mod)
}

Test the Wrapper #

Use the wrapper instead of the raw multinom() function, and emmeans() will work correctly:

# Fit model with the wrapper
mod_fixed <- multinom_wrap(my_formula, data = dat, trace = FALSE)

# Now emmeans() runs without errors!
post_hoc_results <- emmeans(mod_fixed, pairwise ~ predictor1 | predictor2)
print(post_hoc_results)

Why This Works #

By explicitly setting mod$call$formula to the original formula object (not the string version), we’re giving emmeans() the structured formula it needs to correctly identify predictors, interactions, and their relationships. This bypasses the quirk where multinom() converts dynamic formulas to strings in the model call.

You can extend this wrapper to handle additional multinom() arguments (like weights, subset, etc.) by passing them through with the ... parameter, just like we did above.

内容的提问来源于stack exchange，提问作者Mark White