含交互项by参数的mgcv GAM模型提取emmeans报错的解决方法

问题描述 #

我使用mgcv包拟合了带交互项分组平滑的GAM模型，代码如下：

gam(weight ~ s(Time, by = Diet:Treatment), data = dat)

但尝试用emmeans提取边际均值时触发错误：

Error in `[.data.frame`(tbl, , vars, drop = FALSE) : 
  undefined columns selected

完整可复现代码（基于datasets::ChickWeight数据集）：

library(mgcv)
library(emmeans)

dat <- as.data.frame(ChickWeight)
dat <- dat[dat$Chick %in% c(1:2, 21:22), ]
dat$Treatment <- factor(ifelse(dat$Chick %in% c(1, 21), "Control", "Treatment"))

fit <- gam(weight ~ s(Time, by = Diet:Treatment), data = dat)
emmeans(fit, specs = "Time", by = c("Diet", "Treatment"),
        at = list(Time = 0:20))
## 返回错误: undefined columns selected
emmeans::ref_grid(fit, at = list(Time = 0:20))
## 返回错误: undefined columns selected

错误原因 #

这个错误的核心是emmeans的recover_data.gam()函数无法正确解析GAM模型中by = Diet:Treatment这种直接使用因子交互项的语法。emmeans在恢复模型数据时，会尝试从原始数据中查找名为Diet:Treatment的列，但该交互项是模型拟合时动态生成的，原始数据中并没有这一列，因此触发了“未定义列选中”的错误。这属于emmeans对mgcv模型支持的局限性，而非严格意义上的bug——mgcv允许by参数直接使用交互项表达式，但emmeans的恢复逻辑更依赖原始数据中存在的显式变量。

解决办法 #

方法1：显式创建交互因子（临时方案优化） #

在数据中创建显式的交互因子，再将其作为by参数传入GAM模型，这样emmeans就能正常识别变量：

dat$Diet_Treatment <- interaction(dat$Diet, dat$Treatment)
fit2 <- gam(weight ~ s(Time, by = Diet_Treatment), data = dat)
emmeans(fit2, specs = "Time", by = "Diet_Treatment",
        at = list(Time = 0:20))

如果需要后续做特定对比，可通过emmeans的pairs()函数结合交互因子的水平实现，比如：

# 对比Control组的Diet=1 vs Diet=2，Treatment组的Diet=1 vs Diet=2
pairs(emmeans(fit2, specs = "Diet_Treatment"), 
      contrast = list(
        "Control_D1 vs Control_D2" = c(1, -1, 0, 0),
        "Treatment_D1 vs Treatment_D2" = c(0, 0, 1, -1)
      ))

方法2：手动构造参考网格（更灵活方案） #

如果不想修改原始模型，可以手动构造包含Diet、Treatment和Time的参考数据框，用predict()计算均值后包装成emmGrid对象：

# 构造参考数据
ref_dat <- expand.grid(
  Diet = levels(dat$Diet),
  Treatment = levels(dat$Treatment),
  Time = 0:20
)
# 预测均值
ref_dat$pred <- predict(fit, newdata = ref_dat)
# 转换为emmGrid对象
emm_obj <- emmeans::as.emmGrid(ref_dat, 
                               specs = ~ Time | Diet * Treatment,
                               resp = "pred")
# 查看结果
emm_obj

这种方法无需修改模型，就能直接按Diet和Treatment分组获取Time的边际均值，方便后续做任意对比。

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