为何emmeans()需硬编码公式，无法从动态线性模型提取？

问题场景 #

在处理含异方差的模型事后检验时，若从动态模型中提取公式传入gls()，再用emmeans()执行事后检验会触发报错；但硬编码公式则能正常运行。以下分Shiny和非Shiny场景给出最小可复现示例：

Shiny场景MRE #

模拟大型应用中用户指定模型、ANOVA显著后执行事后检验的流程，需构建按效应组合设置不同方差的gls模型：

library(shiny)
library(nlme)
library(emmeans)

# 模拟数据
set.seed(123)
dat <- expand.grid(percent = factor(c("5%", "10%", "15%")),
                   cycle = factor(c("1", "2", "3")),
                   rep = 1:5)
dat$value <- rnorm(nrow(dat), 
                   mean = with(dat, as.numeric(percent) + as.numeric(cycle)),
                   sd = with(dat, as.numeric(percent)*0.5))

ui <- fluidPage(
  actionButton("run", "Run Post Hoc")
)

server <- function(input, output, session) {
  observeEvent(input$run, {
    # 模拟用户动态选择的线性模型
    lm_model <- lm(value ~ percent*cycle, data = dat)
    
    # 提取公式用于gls模型
    gls_formula <- formula(lm_model)
    
    # 构建带异方差的gls模型
    gls_model <- gls(gls_formula, 
                     data = dat,
                     weights = varIdent(form = ~1 | percent:cycle))
    
    # 动态公式传入emmeans（触发报错）
    results <- emmeans(gls_model, pairwise ~ percent*cycle)
    print(results)
    
    # 硬编码公式（正常运行）
    results2 <- emmeans(gls_model, pairwise ~ percent*cycle)
    print(results2)
  })
}

shinyApp(ui, server)

报错信息：Error in eval(call$model) : object 'gls_formula' not found

非Shiny场景MRE #

library(nlme)
library(emmeans)

set.seed(123)
dat <- expand.grid(percent = factor(c("5%", "10%", "15%")),
                   cycle = factor(c("1", "2", "3")),
                   rep = 1:5)
dat$value <- rnorm(nrow(dat), 
                   mean = with(dat, as.numeric(percent) + as.numeric(cycle)),
                   sd = with(dat, as.numeric(percent)*0.5))

# 初始模型
lm_model <- lm(value ~ percent*cycle, data = dat)
gls_formula <- formula(lm_model)

# 构建gls模型
gls_model <- gls(gls_formula, 
                 data = dat,
                 weights = varIdent(form = ~1 | percent:cycle))

# 执行事后检验（触发报错）
results <- emmeans(gls_model, pairwise ~ percent*cycle)

报错信息：Error in call$model[[2]] : object of type 'symbol' is not subsettable

问题根源 #

gls()会保存调用时的符号引用而非实际公式表达式。当你传入gls_formula这个符号对象时，gls的调用记录里只保存了gls_formula这个名称，而非value ~ percent*cycle的实际公式。emmeans()在解析模型时会尝试评估这个符号，导致找不到对象或无法子集化符号，从而报错。

解决方法 #

两种方式都能让gls()保存实际的公式表达式，避免符号引用问题：

方法1：在gls调用中强制解析公式 #

用eval()直接解析公式对象，让gls记录实际的公式表达式：

# 替换原gls模型构建代码
gls_model <- gls(eval(gls_formula), 
                 data = dat,
                 weights = varIdent(form = ~1 | percent:cycle))

方法2：提取公式的原始表达式 #

直接从原模型中提取公式的左右部分，重新构建公式表达式：

# 提取公式的响应变量和预测变量部分
gls_formula_expr <- formula(lm_model)[[2]] ~ formula(lm_model)[[3]]

# 构建gls模型
gls_model <- gls(gls_formula_expr, 
                 data = dat,
                 weights = varIdent(form = ~1 | percent:cycle))

两种方法都能让emmeans()正确解析模型，顺利执行事后检验。

内容的提问来源于stack exchange，提问作者Steven Ouellette