含交互项的线性混合效应模型ggplot绘图异常求助

当拟合包含AGE_Switch、SEX与治疗时间双项交互的线性混合模型后，直接用geom_line或geom_smooth绘图容易出现混乱或结果不匹配，核心原因是ggplot的内置函数无法正确复用你已拟合的复杂模型，尤其是带自定义分段函数的双斜率结构。正确的做法是基于模型生成专门的预测数据集，再绘制预测值的趋势线，步骤如下：

1. 生成预测数据集 #

需要覆盖所有自变量的关键组合，确保能展示不同性别、年龄组在治疗调整前后的血压趋势：

# 确保自定义分段函数在环境中
Delai <- 0
DelaiNeg <- function(x, Delai) ifelse(x < Delai, x, 0)
DelaiPos <- function(x, Delai) ifelse(x < Delai, 0, x)

# 提取关键变量的代表性取值
# 年龄取均值及±1标准差，覆盖多数样本范围
age_ref <- with(IeDEA_BP_Switch, c(mean(AGE_Switch)-sd(AGE_Switch), mean(AGE_Switch), mean(AGE_Switch)+sd(AGE_Switch)))
# 性别取所有分类水平
sex_ref <- unique(IeDEA_BP_Switch$SEX)
# 治疗时间序列，覆盖观察到的全范围
delai_seq <- seq(min(IeDEA_BP_Switch$DELAI_SWITCH, na.rm=T), max(IeDEA_BP_Switch$DELAI_SWITCH, na.rm=T), length.out=100)

# 生成所有组合的预测数据集
pred_df <- expand.grid(
  DELAI_SWITCH = delai_seq,
  AGE_Switch = age_ref,
  SEX = sex_ref,
  PATIENT = sample(unique(IeDEA_BP_Switch$PATIENT), 1) # 仅用于占位，忽略随机效应时不影响
)

2. 从模型生成预测值 #

使用predict()函数基于已拟合的模型生成预测值，通过re.form参数控制是否包含随机效应：

# 生成固定效应的预测值（展示整体趋势，忽略个体随机效应）
pred_df$bp_pred <- predict(model_Diasto, newdata=pred_df, re.form=~0)

# 如果需要展示个体随机效应的差异，可改为：
# pred_df <- expand.grid(
#   DELAI_SWITCH = delai_seq,
#   AGE_Switch = age_ref,
#   SEX = sex_ref,
#   PATIENT = unique(IeDEA_BP_Switch$PATIENT)
# )
# pred_df$bp_pred <- predict(model_Diasto, newdata=pred_df, re.form=~(1|PATIENT))

3. 绘制趋势图 #

用ggplot绘制预测值的趋势线，可叠加原始数据散点作为参考：

library(ggplot2)

ggplot() +
  # 原始数据散点（透明度调低避免遮挡趋势线）
  geom_point(data=IeDEA_BP_Switch, aes(x=DELAI_SWITCH, y=DIASTO, color=SEX), alpha=0.2) +
  # 模型预测的趋势线，按性别分组、年龄区分线型
  geom_line(data=pred_df, aes(x=DELAI_SWITCH, y=bp_pred, color=SEX, linetype=factor(AGE_Switch)), linewidth=1) +
  # 标注治疗调整时间点（x=0）
  geom_vline(xintercept=0, linetype="dashed", color="gray50") +
  # 调整标签与主题
  labs(
    x = "Time since treatment adjustment",
    y = "Diastolic Blood Pressure (mmHg)",
    color = "Sex",
    linetype = "Age at adjustment"
  ) +
  theme_minimal()

关键注意事项 #

• 避免用geom_smooth(method="lmer")：该函数会重新拟合模型，无法识别你自定义的DelaiNeg/DelaiPos分段函数，导致拟合结果与原模型不一致。
• 变量一致性：确保预测数据集的变量名与模型中完全匹配，且Delai变量在当前环境中存在，否则predict()会报错。
• 随机效应控制：若要展示个体差异，需在预测数据集中包含所有患者ID，但这样图会较杂乱，优先展示固定效应的整体趋势更直观。

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