基于mids对象的lme混合效应模型生成交互图遇错求解

我通过mice包生成的mids对象（多重插补数据集）运行了包含三向连续变量交互项的线性混合效应模型，尝试用interactions::interact_plot绘制交互图时持续报错，推测是模型基于mids对象而非普通数据框导致的。

相关代码 #

# 基于多重插补数据拟合模型
MIDmod1 <- with(data = df.mids, exp = lmer(GC ~ Age + Sex + Edu + Stress*Time*HLI + (1|ID)))
summary(pool(MIDmod1))

# 尝试绘制三向交互图
interact_plot(
  model=MIDmod1, 
  pred = Time,
  modx=Stress,
  mod2=HLI,
  data = df.mids,
  interval=TRUE,
  y.label='Global cognition composite score',
  modx.labels=c('Low Baseline Stress (-1SD)','Moderate Baseline Stress (Mean)', 'High Baseline Stress (+1SD)'),
  mod2.labels=c('Low HLI (-1SD)', 'Moderate HLI (Mean)', 'High HLI (+1SD)'),
  legend.main='') + ylim(-2,2)

报错信息 #

• 传入data参数时：

Error in rep(1, times = nrow(data)) : invalid 'times' argument

• 不传入data参数时：

Error in formula.default(object, env = baseenv()) : invalid formula

注：当模型基于普通数据框时，可正常生成目标交互图。

interactions::interact_plot不支持直接处理mice包中with()返回的mira类模型对象，以下两种方法可解决该问题：

方法1：用单个插补数据集快速绘图 #

从mids对象中提取某一个插补完成的数据集，重新拟合模型后再调用interact_plot。这种方法适合快速验证交互趋势，若要整合多重插补结果，需手动合并不同插补数据集的绘图：

# 提取第一个插补后的数据集
df_imputed <- complete(df.mids, 1)

# 用单个插补数据集重新拟合模型
mod_single <- lmer(GC ~ Age + Sex + Edu + Stress*Time*HLI + (1|ID), data = df_imputed)

# 绘制交互图
interact_plot(
  model=mod_single, 
  pred = Time,
  modx=Stress,
  mod2=HLI,
  interval=TRUE,
  y.label='Global cognition composite score',
  modx.labels=c('Low Baseline Stress (-1SD)','Moderate Baseline Stress (Mean)', 'High Baseline Stress (+1SD)'),
  mod2.labels=c('Low HLI (-1SD)', 'Moderate HLI (Mean)', 'High HLI (+1SD)'),
  legend.main='') + ylim(-2,2)

方法2：基于多重插补合并结果生成严谨的边际效应图 #

这种方法先对每个插补模型的预测结果进行合并，再用ggplot2绘图，能整合多重插补的统计结果：

library(dplyr)
library(ggplot2)

# 定义变量的均值和标准差，用于分组
stress_mean <- mean(complete(df.mids, 1)$Stress, na.rm = TRUE)
stress_sd <- sd(complete(df.mids, 1)$Stress, na.rm = TRUE)
hli_mean <- mean(complete(df.mids, 1)$HLI, na.rm = TRUE)
hli_sd <- sd(complete(df.mids, 1)$HLI, na.rm = TRUE)

# 创建预测网格：覆盖Time的全范围，Stress和HLI取均值±1SD，控制协变量为均值/基准水平
pred_grid <- expand.grid(
  Time = seq(min(complete(df.mids,1)$Time), max(complete(df.mids,1)$Time), length.out = 100),
  Stress = c(stress_mean - stress_sd, stress_mean, stress_mean + stress_sd),
  HLI = c(hli_mean - hli_sd, hli_mean, hli_mean + hli_sd),
  Age = mean(complete(df.mids,1)$Age, na.rm = TRUE),
  Sex = levels(complete(df.mids,1)$Sex)[1],
  Edu = mean(complete(df.mids,1)$Edu, na.rm = TRUE),
  ID = complete(df.mids,1)$ID[1] # 固定随机效应，只预测固定效应部分
)

# 对每个插补数据集进行预测
pred_list <- lapply(1:df.mids$m, function(i) {
  mod_i <- lmer(GC ~ Age + Sex + Edu + Stress*Time*HLI + (1|ID), data = complete(df.mids, i))
  predict(mod_i, newdata = pred_grid, re.form = NA)
})

# 合并预测结果，计算均值和95%置信区间
pred_grid <- pred_grid %>%
  mutate(
    GC_mean = rowMeans(do.call(cbind, pred_list)),
    GC_lower = apply(do.call(cbind, pred_list), 1, quantile, 0.025),
    GC_upper = apply(do.call(cbind, pred_list), 1, quantile, 0.975),
    # 添加分组标签
    Stress_group = case_when(
      Stress == stress_mean - stress_sd ~ 'Low Baseline Stress (-1SD)',
      Stress == stress_mean ~ 'Moderate Baseline Stress (Mean)',
      Stress == stress_mean + stress_sd ~ 'High Baseline Stress (+1SD)'
    ),
    HLI_group = case_when(
      HLI == hli_mean - hli_sd ~ 'Low HLI (-1SD)',
      HLI == hli_mean ~ 'Moderate HLI (Mean)',
      HLI == hli_mean + hli_sd ~ 'High HLI (+1SD)'
    )
  )

# 绘制交互图
ggplot(pred_grid, aes(x = Time, y = GC_mean, color = Stress_group)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = GC_lower, ymax = GC_upper, fill = Stress_group), alpha = 0.2, color = NA) +
  facet_wrap(~HLI_group) +
  labs(x = 'Time', y = 'Global cognition composite score', color = '', fill = '') +
  ylim(-2, 2) +
  theme_minimal()

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