使用nlme、ggeffects、sjplot绘制重复测量lme模型的群体预测图

一、ggemmeans实现（支持自定义边际数据集输出） #

报错原因说明 #

你之前运行报错的核心问题有两个：

• 错误使用type = "random"：边际效应仅需计算固定效应部分，该参数会纳入随机效应预测，无法匹配到固定效应中的交互项
• terms参数格式错误：双变量交互需传入两个变量的向量，而非冒号连接的字符串

代码实现 #

# 加载所需包
library(ggeffects)
library(ggplot2)

# 计算fu_time与年龄的交互边际效应
me_age <- ggemmeans(
  model = ris,
  terms = c("fu_time [all]", "age [quart2]"), # 第一个变量为X轴变量，第二个为分组变量，quart2表示取年龄25/50/75分位数
  type = "fixed", # 仅计算固定效应边际均值
  back.transform = FALSE, # 关闭自动反变换，手动计算避免识别误差
  wt.nuis = "proportional" # 分类协变量按样本实际占比加权调整，连续协变量默认取样本均值
)

# 手动对自然对数转换的结果做反变换，回到原始尺度
me_age$predicted <- exp(me_age$predicted)
me_age$conf.low <- exp(me_age$conf.low)
me_age$conf.high <- exp(me_age$conf.high)

输出的me_age即为符合要求的边际数据集，可直接用于自定义绘图，绘图示例如下：

ggplot(me_age, aes(x = x, y = predicted, color = group)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  labs(x = "随访时间（年）", y = "结局指标（原始尺度）", color = "年龄分位", fill = "年龄分位") +
  theme_bw()

若要计算分类变量（如高血压）和fu_time的交互，仅需修改terms参数为c("fu_time [all]", "高血压")即可。

二、sjPlot快速出图修正方法 #

针对plot_model输出为对数尺度、无法指定交互项的问题，可通过指定交互项+反变换参数解决：

library(sjPlot)

plot_model(
  ris,
  type = "int",
  terms = c("fu_time", "age"), # 按需指定要展示的交互项，顺序为X轴/分组
  transform = "exp", # 对自然对数转换的结果做指数反变换，回到原始尺度
  ci.lvl = 0.95
)

若去掉terms参数，会自动输出模型中所有交互项的原始尺度效应图。

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