因子水平组合缺失时，如何获取平均效应量估计值与置信区间

问题场景 #

使用metafor拟合包含调节变量交互项的混合效应模型时，因某一因子水平组合无对应数据，调用orchaRd的mod_results函数提取平均效应量及置信区间时出现如下报错：

Error in ref_grid(result, ...) : Something went wrong:
Non-conformable elements in reference grid.

手动补全系数和协方差矩阵后，结果与summary()输出完全一致，无法得到完整的各组合效应量估计。

用户原始代码：

library(orchaRd)
data(fish)
warm_dat <- fish
warm_dat$factor2<-as.factor(seq(c("A","B","C"),133) )

mod_fit <- metafor::rma.mv(yi = lnrr, V = lnrr_vi,
random = list(~1 | group_ID, ~1 | es_ID),
mods = ~ trait.type*factor2,
method = "REML", test = "t",
control=list(optimizer="optim", optmethod="Nelder-Mead"), data = warm_dat)

# 调用mod_results报错
results1 <- mod_results(mod_fit, group = "group_ID", mod = "trait.type", by = "Measurement_cat")
results2 <- mod_results(mod_fit, group = "group_ID", mod = "trait.type")

可行解决方案 #

方案1：添加虚拟数据补全因子组合 #

通过向数据集添加一行缺失组合的虚拟NA数据，让模型保留所有交互项参数，从而使orchaRd::mod_results正常工作：

# 构造缺失组合的虚拟行（根据实际缺失的水平调整）
missing_row <- data.frame(
  lnrr = NA,
  lnrr_vi = NA,
  group_ID = NA,
  es_ID = NA,
  trait.type = levels(warm_dat$trait.type)[which(!levels(warm_dat$trait.type) %in% unique(warm_dat$trait.type[warm_dat$factor2 == "C"]))],
  factor2 = "C",
  Measurement_cat = NA
)

# 合并数据
warm_dat_full <- rbind(warm_dat, missing_row)

# 重新拟合模型
mod_fit_full <- metafor::rma.mv(
  yi = lnrr, V = lnrr_vi,
  random = list(~1 | group_ID, ~1 | es_ID),
  mods = ~ trait.type*factor2,
  method = "REML", test = "t",
  control=list(optimizer="optim", optmethod="Nelder-Mead"), 
  data = warm_dat_full
)

# 正常调用mod_results
results <- orchaRd::mod_results(mod_fit_full, group = "group_ID", mod = "trait.type")
print(results)

方案2：使用metafor自带predict函数手动计算 #

直接通过predict.rma.mv结合expand.grid生成所有因子组合的效应量、置信区间和预测区间，无需依赖emmeans：

# 生成所有因子水平的完整组合
newdat <- expand.grid(
  trait.type = levels(warm_dat$trait.type),
  factor2 = levels(warm_dat$factor2),
  stringsAsFactors = TRUE
)

# 计算效应量、置信区间和预测区间
preds <- predict(mod_fit, newdata = newdat, 
                 interval = "confidence",  # 置信区间
                 predinterval = TRUE,      # 预测区间
                 addx = TRUE)

# 整理结果为数据框
result_df <- cbind(newdat, preds[, c("pred", "ci.lb", "ci.ub", "pi.lb", "pi.ub")])
print(result_df)

方案3：修正emmeans的qdrg调用 #

手动指定完整的因子水平，让emmeans正确构建参考网格：

# 构建包含所有水平的参考网格
grid <- emmeans::qdrg(
  formula = stats::formula(mod_fit),
  data = warm_dat,
  coef = mod_fit$b,
  vcov = stats::vcov(mod_fit),
  df = mod_fit$k - mod_fit$p,  # 使用模型正确的自由度
  levels = list(
    trait.type = levels(warm_dat$trait.type),
    factor2 = levels(warm_dat$factor2)
  ),
  allow.new.levels = TRUE
)

# 提取交互项的平均效应量
mm <- emmeans::emmeans(grid, specs = ~ trait.type * factor2)
print(mm)

# 计算预测区间（使用orchaRd函数）
mm_pi <- orchaRd::pred_interval_esmeans(mod_fit, mm, mod = "trait.type")
print(mm_pi)

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