如何绘制R中metafor包rma.mv模型的交互效应图？

我用metafor包的rma.mv()函数构建了包含Seedmass*Treat交互项的混合效应元分析模型，想绘制交互效应图。尝试用sjPlot包的plot_model()函数时，报错：Error: No interaction term found in model.，但这个函数对lm模型能正常工作。我有两个问题：

1. 若仅关注交互效应，用lm模型绘图是否可行？
2. 如何针对RMA/RMA.MV类模型绘制交互效应？

附数据集及代码示例：

数据集 #

LRRs.long<-structure(list(fullname = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 
3L, 3L, 3L), .Label = c("Argemone glauca", "Bacopa monnieri", 
"Chenopodium oahuense"), class = "factor"), Treat = structure(c(1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("10", "20", "35"), class = "factor"), 
    LRR = c(-0.954455598177961, -3.43398720448515, -3.43398720448515, 
    -0.206947246533427, -1.33340726127878, -3.59401011726068, 
    -0.0387710980425032, -0.7914266867076, -4.47267908945549), 
    LRR_var = c(0.336731047497513, 0.00631364800243193, 0.00631364800243193, 
    0.0106766272234543, 0.0245128677435665, 0.230755613454606, 
    0.00381238394178319, 0.0248419470726316, 0.00260356604424191
    ), Seedmass = c(0.0026402, 0.0026402, 0.0026402, 3.3e-05, 
    3.3e-05, 3.3e-05, 0.000234175, 0.000234175, 0.000234175)), row.names = c(1L, 
2L, 3L, 4L, 5L, 6L, 10L, 11L, 12L), class = "data.frame")

建模代码 #

library(metafor)
mass.initial.i<-rma.mv(LRR, LRR_var, mods=~Seedmass*Treat, random= list(~1|fullname), method="ML",digits=4,data=LRRs.long)
summary(mass.initial.i)
class(mass.initial.i)

mass.initial.Noi<-rma.mv(LRR, LRR_var, mods=~Seedmass + Treat, random= list(~1|fullname), method="ML",digits=4,data=LRRs.long)

anova(mass.initial.i, mass.initial.Noi)

library(sjPlot)
theme_set(theme_sjplot())
plot_model(mass.initial.i, type = "int")

问题解答 #

1. 用lm模型绘图是否可行？ #

绝对不建议这么做。rma.mv()是加权混合效应模型，会把LRR_var作为权重纳入计算，而lm是无权重的普通线性回归，两者的参数估计（包括交互项）会存在显著差异。如果用lm绘图，得到的是无权重下的交互关系，和你实际拟合的元分析模型结论完全不一致，参考价值极低。

2. 针对RMA/RMA.MV模型绘制交互效应的方法 #

方法1：手动计算预测值 + ggplot2绘图（最灵活）

通过手动生成预测数据集，从模型中提取预测值和置信区间，再用ggplot2绘图：

library(ggplot2)

# 1. 创建预测数据集：覆盖Seedmass的全范围，结合所有Treat水平
new_data <- expand.grid(
  Seedmass = seq(min(LRRs.long$Seedmass), max(LRRs.long$Seedmass), length.out = 100),
  Treat = unique(LRRs.long$Treat)
)

# 2. 从rma.mv模型中获取预测值和置信区间
preds <- predict(mass.initial.i, newmods = model.matrix(~Seedmass*Treat, data = new_data), addx = TRUE)
new_data$pred <- preds$pred
new_data$lower <- preds$ci.lb
new_data$upper <- preds$ci.ub

# 3. 绘制交互效应图
ggplot(new_data, aes(x = Seedmass, y = pred, color = Treat)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = Treat), alpha = 0.2, color = NA) +
  labs(x = "种子质量", y = "预测LRR值", color = "处理组", fill = "处理组") +
  theme_bw()

方法2：用emmeans包计算边际均值并绘图

emmeans支持rma.mv模型，可以直接计算不同处理组下种子质量的边际效应，再可视化：

library(emmeans)

# 计算不同Treat水平下，Seedmass的边际效应
emm <- emmeans(mass.initial.i, ~ Seedmass | Treat, at = list(Seedmass = seq(min(LRRs.long$Seedmass), max(LRRs.long$Seedmass), length.out = 100)))

# 转换为数据框并绘图
emm_df <- as.data.frame(emm)
ggplot(emm_df, aes(x = Seedmass, y = emmean, color = Treat)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = Treat), alpha = 0.2, color = NA) +
  labs(x = "种子质量", y = "预测LRR值", color = "处理组", fill = "处理组") +
  theme_bw()

为什么sjPlot的plot_model()会报错？ #

sjPlot的plot_model()目前对rma.mv模型的交互项识别支持不完善，无法正确解析模型中的交互项结构，因此会提示找不到交互项。上述两种手动方法可以完全替代它的功能，且更灵活可控。

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