在R中如何用点须图展示线性混合模型的全部预测值？

问题背景 #

我在R中构建了以sex（性别）、Light（昼夜）为预测变量，Displacement（深度位移）为响应变量的线性混合模型，代码如下：

displacement_lmm_hour <- lmer(Displacement~sex*Light + (1|Hour), data = avg_depth_df_hour)

想要用点须图展示预测变量所有组合的预测值，遇到以下问题：

• 使用dotwhisker包的dwplot函数：

  dwplot(displacement_lmm_hour, effects = "fixed")
  

  生成的图仅显示部分预测值（无雄性、白天数值），模型summary输出也仅包含这些结果。
• 使用plot_model函数：

  plot_model(displacement_lmm_hour, type = "pred", terms = c("sex","Light"), axis.title = c("Sex", "Displacement"))
  

  虽能拆分预测值，但认为误差条不准确。

解决方案 #

一、排查模型系数显示不全的原因 #

summary和dwplot只显示部分结果，通常是因为变量因子水平设置问题，或数据中存在某类变量组合无样本的情况：

1. 先检查所有变量组合的样本量：

   table(avg_depth_df_hour$sex, avg_depth_df_hour$Light)
   

2. 若数据齐全，重新设置因子水平确保所有类别被纳入：

   # 替换为你的变量实际类别
   avg_depth_df_hour$sex <- factor(avg_depth_df_hour$sex, levels = c("male", "female"))
   avg_depth_df_hour$Light <- factor(avg_depth_df_hour$Light, levels = c("day", "night"))
   # 重新拟合模型
   displacement_lmm_hour <- lmer(Displacement~sex*Light + (1|Hour), data = avg_depth_df_hour)
   

二、用dwplot绘制所有组合的预测值 #

dwplot默认展示模型固定效应系数，而非预测值。需手动生成预测数据集后绘图：

library(dotwhisker)
library(lme4)
library(ggplot2)

# 生成所有预测变量组合的数据集
newdata <- expand.grid(
  sex = unique(avg_depth_df_hour$sex),
  Light = unique(avg_depth_df_hour$Light),
  Hour = unique(avg_depth_df_hour$Hour) # 若要忽略Hour的随机效应，可设为mean(avg_depth_df_hour$Hour)
)

# 计算仅考虑固定效应的预测值与95%置信区间
preds <- predict(displacement_lmm_hour, newdata = newdata, interval = "confidence", re.form = NA)
newdata$fit <- preds[,1]
newdata$lower <- preds[,2]
newdata$upper <- preds[,3]

# 合并变量组合为分组标签
newdata$group <- paste(newdata$sex, newdata$Light, sep = "_")

# 绘制点须图
dwplot(newdata, x = "fit", y = "group", xmin = "lower", xmax = "upper") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Predicted Displacement", y = "Sex-Light Combination") +
  theme_bw()

三、替代方案：用ggplot2直接绘制 #

若需要更灵活的自定义，直接用ggplot2绘制：

ggplot(newdata, aes(x = group, y = fit)) +
  geom_point(size = 3, color = "#2c3e50") +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, color = "#2c3e50") +
  labs(x = "Sex-Light Combination", y = "Predicted Displacement") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

四、修正plot_model的误差条 #

plot_model(type = "pred")默认展示预测区间（包含随机效应变异），若需要置信区间，添加interval = "confidence"参数：

plot_model(displacement_lmm_hour, type = "pred", terms = c("sex","Light"), 
           axis.title = c("Sex", "Displacement"), interval = "confidence")

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