系统发育广义最小二乘模型（PGLS）：如何可视化预测不确定性？

问题背景 #

使用phylolm::phylolm()构建模型平均的系统发育广义最小二乘（PGLS）模型，通过MuMIn完成模型平均流程（选择phylolm是看重其速度优势），需要可视化模型预测结果及不确定性，但遇到两个问题：

• ggeffects::ggpredict()处理模型平均对象时报错
• phylolm::phylolm()原生predict()方法不支持se.fit = TRUE参数

可复现代码 #

library(caper)
library(ggeffects)
library(phylolm)
library(MuMIn)

# 加载数据
data(shorebird)

# 拟合基础PGLS模型
mod <- phylolm::phylolm(Egg.Mass ~ M.Mass + F.Mass, data = shorebird.data, phy = shorebird.tree, model = "lambda")

# 模型平均流程
options(na.action = "na.fail")
mod.d <- MuMIn::dredge(mod, rank = "AICc")
mod.avg.fit <- MuMIn::model.avg(mod.d, revised.var = TRUE, fit = TRUE)

# ggpredict尝试失败
# plot( ggeffects::ggpredict(mod.avg.fit, terms = c("M.Mass")) )

# 原生predict无se.fit支持
# summary(predict(mod, se.fit = TRUE))

解决方案：手动计算预测值与不确定性并可视化 #

利用MuMIn为模型平均对象实现的predict()方法（支持se.fit = TRUE），手动计算置信区间后用ggplot2绘制：

步骤1：生成预测用新数据 #

固定其他协变量（这里用均值），生成目标变量的序列值：

library(ggplot2)

# 生成预测数据集：固定F.Mass为均值，M.Mass取全范围的100个点
new_data <- expand.grid(
  M.Mass = seq(min(shorebird.data$M.Mass), max(shorebird.data$M.Mass), length.out = 100),
  F.Mass = mean(shorebird.data$F.Mass, na.rm = TRUE)
)

步骤2：获取预测值与标准误 #

通过MuMIn的模型平均预测方法提取预测值和标准误：

# 从模型平均对象获取预测值和标准误
pred_results <- predict(mod.avg.fit, newdata = new_data, se.fit = TRUE)

# 将结果合并到新数据集中
new_data$predicted <- pred_results$fit
new_data$se <- pred_results$se.fit

步骤3：计算置信区间 #

以95%置信区间为例，计算上下限：

# 计算95%置信区间
new_data$lower_ci <- new_data$predicted - 1.96 * new_data$se
new_data$upper_ci <- new_data$predicted + 1.96 * new_data$se

步骤4：可视化结果 #

用ggplot2绘制趋势线、置信带和原始数据点：

ggplot(new_data, aes(x = M.Mass, y = predicted)) +
  # 绘制预测趋势线
  geom_line(color = "#2c3e50", linewidth = 1) +
  # 绘制置信带
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), alpha = 0.2, fill = "#3498db") +
  # 添加原始数据点
  geom_point(data = shorebird.data, aes(x = M.Mass, y = Egg.Mass), 
             alpha = 0.6, color = "#e74c3c") +
  # 设置标签和主题
  labs(x = "雄性体重", y = "卵重", title = "模型平均PGLS预测趋势（带95%置信区间）") +
  theme_minimal()

说明 #

• MuMIn的predict.model.avg()方法内置了对标准误的计算支持，无需依赖phylolm原生的predict方法
• 若需调整其他协变量的固定值（如中位数、分组水平），只需修改new_data的生成逻辑即可
• 该方法适用于所有MuMIn生成的模型平均对象，兼容性更强

内容的提问来源于stack exchange，提问作者M. Riera