基于MuMIn::model.avg()的GAMMs模型open_water平均偏效应评估问询

我用MuMIn::model.avg()生成了一组权重占所有候选模型95%的最优模型，基于该模型平均对象生成预测结果。3个最优模型如下：

catch ~ s(effort) + s(year) + s(month) + s(open_water) + s(pop_dens)
catch ~ s(effort) + s(year) + s(month) + s(forest)
catch ~ s(effort) + s(year) + s(month) + s(open_water)

其中open_water/forest是比例变量。我想绘制open_water的平均偏效应图，但这个变量仅出现在2个最优模型中，且各模型权重不同，有没有更直接的评估方法？

我当前的做法是：从每个包含目标变量的模型生成估计值，乘以对应模型权重后相加，得到跨模型的平均估计值，示例代码如下：

# 获取包含目标变量的模型
water_mods <- top_mod_list[grep(pattern = "water", x = names(top_mod_list))]
# 标准化权重
water_weights <-  model.sel(water_mods) %>% as.data.frame() %>%
  mutate(mod = rownames(.))

# 固定其他预测变量为均值，仅改变目标变量
pdat <- means %>% select(-water)
pdat$water <- seq(min(landing_dat$water),
                  max(landing_dat$water),
                  length.out = 200)

# 根据目标变量的变化预测响应值
water_predict <- 
  map(water_mods,
      ~pdat %>% 
        mutate(pred = predict(.x,
                              newdata = pdat, 
                              # type = "response",
                              newdata.guaranteed = TRUE),
               weighted_pred = pred * water_weights$weight))) 

不需要手动逐个模型计算加权预测值，利用MuMIn和mgcv的内置功能就能更高效地完成，同时保证结果严谨性：

1. 基于模型平均对象直接生成加权预测 #

如果你已经通过model.avg()得到了模型平均对象（假设命名为avg_mod），直接调用predict()函数就能自动按模型权重计算加权预测值，无需手动筛选模型。

示例代码：

# 构建预测数据集：其他变量固定为均值，仅生成open_water的序列值
pdat <- means %>% 
  select(-open_water) %>% 
  mutate(open_water = seq(min(landing_dat$open_water), 
                          max(landing_dat$open_water), 
                          length.out = 200))

# 从模型平均对象直接获取加权预测值
pdat$weighted_pred <- predict(avg_mod, newdata = pdat, 
                              type = "response", # 按需选择link或response类型
                              newdata.guaranteed = TRUE)

2. 直接计算平均偏效应 #

如果需要的是open_water的偏效应（控制其他变量后的边际影响），可以结合predict.gam()的type = "terms"参数提取目标变量的项，再按模型权重加权平均，同时自动处理不含open_water的模型（这类模型中该变量的偏效应视为0）。

示例代码：

# 提取所有最优模型的权重
all_weights <- model.sel(top_mod_list)$weight
names(all_weights) <- names(top_mod_list)

# 构建预测数据集
pdat <- means %>% 
  select(-open_water) %>% 
  mutate(open_water = seq(min(landing_dat$open_water), 
                          max(landing_dat$open_water), 
                          length.out = 200))

# 遍历所有模型，提取open_water的项并加权平均
water_effects <- map_dfr(names(top_mod_list), function(mod_name) {
  mod <- top_mod_list[[mod_name]]
  # 判断模型是否包含open_water
  if ("s(open_water)" %in% colnames(predict(mod, type = "terms"))) {
    term_pred <- predict(mod, newdata = pdat, type = "terms")[, "s(open_water)"]
    tibble(open_water = pdat$open_water,
           effect = term_pred,
           weight = all_weights[mod_name])
  } else {
    # 不含open_water的模型，偏效应为0
    tibble(open_water = pdat$open_water,
           effect = 0,
           weight = all_weights[mod_name])
  }
}) %>% 
  group_by(open_water) %>% 
  summarize(avg_effect = sum(effect * weight) / sum(weight)) # 加权平均计算最终偏效应

# 绘制平均偏效应图
ggplot(water_effects, aes(x = open_water, y = avg_effect)) +
  geom_line(linewidth = 1) +
  labs(x = "Open Water Proportion", y = "Average Partial Effect on Catch")

方法优势 #

• 自动覆盖所有95%权重内的模型，无需手动筛选，避免遗漏信息；
• 利用内置函数处理权重和预测，减少手动计算的出错概率；
• 直接生成偏效应结果，无需额外转换预测值。

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