如何基于dlply拆分的多lme模型，用ggeffects绘制边际效应图

针对你用dlply按Year拆分Soybean数据集并拟合含Variety随机截距的lme模型的场景，以下是用ggeffects或ggemmeans绘制边际效应/调整预测图的具体步骤：

1. 生成各年份模型的预测数据 #

首先遍历每个年份的模型，生成对应的边际效应或调整预测结果，并添加年份标识：

用ggeffects生成边际效应数据 #

# 遍历模型列表，生成Time的边际效应（默认固定效应层面）
pred_list <- llply(MxM1, function(mod) {
  # 若需包含随机截距的条件效应，添加参数 type = "random"
  ggeffect(mod, terms = "Time")
})

# 合并数据并添加Year列
pred_df <- ldply(names(pred_list), function(year) {
  cbind(pred_list[[year]], Year = as.integer(year))
})

用ggemmeans生成调整预测数据 #

如果需要基于最小二乘均值的调整预测，可使用ggemmeans：

pred_list_emm <- llply(MxM1, function(mod) {
  ggemmeans(mod, terms = "Time")
})

pred_df_emm <- ldply(names(pred_list_emm), function(year) {
  cbind(pred_list_emm[[year]], Year = as.integer(year))
})

2. 绘制可视化图 #

使用ggplot2对合并后的预测数据进行绘图，这里以分面展示各年份的趋势为例：

边际效应图（ggeffects结果） #

ggplot(pred_df, aes(x = x, y = predicted)) +
  geom_line(linewidth = 1) +
  # 添加置信区间
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  # 按年份分面对比
  facet_wrap(~Year) +
  labs(x = "Time", y = "Predicted Weight", title = "Marginal Effect of Time on Weight by Year") +
  theme_bw()

调整预测图（ggemmeans结果） #

ggplot(pred_df_emm, aes(x = x, y = predicted)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  facet_wrap(~Year) +
  labs(x = "Time", y = "Predicted Weight", title = "Adjusted Predictions of Weight by Time and Year") +
  theme_bw()

关键说明 #

• ggeffect(mod, terms = "Time")默认计算边际效应：仅基于固定效应的预测，不考虑Variety随机截距的个体变异，适合展示整体趋势。
• 若需展示包含随机截距的条件效应（即每个Variety的预测线），可修改为ggeffect(mod, terms = "Time", type = "random")，后续绘图时可添加aes(color = group)来区分品种。
• llply和ldply保持与你现有代码的plyr风格一致，也可替换为purrr包的map和map2_dfr函数实现相同功能。

内容的提问来源于stack exchange，提问作者Ahir Bhairav Orai