多水平元分析的Egger检验与漏斗图绘制技术问询

我正在用metafor包的rma.mv函数结合clubSandwich工具，针对多水平元分析做Egger检验来检测发表偏倚。网上关于多水平元数据适配Egger检验及生成漏斗图的资料很少，不确定自己写的代码是否合理，希望有相关经验的人帮忙检查代码正确性。

以下是我编写的代码：

egger_multi <- rma.mv(Effectsize, Effectsize_Va^2, random = ~ 1 |
    Study_Id / Effectsize_Id, mods = ~ Effectsize_Va, data = Data_DSMIV)

coef_test(egger_multi, vcov = "CR2")

#create a funnel plot with the results of the Egger's test.

library(metafor)

dat <- Data_DSMIV
res <- rma.mv(Effectsize, Effectsize_Va, random = ~ 1 | Study_Id/Effectsize_Id, data=dat)
res
funnel(res)

cols <- palette.colors(length(unique(dat$Study_Id)), palette="Alphabet")
cols <- cols[as.numeric(factor(dat$Study_Id))]
funnel(res, col=cols)

代码问题分析与修正建议 #

Egger检验核心问题 #

Egger检验的核心逻辑是用**效应量的标准误（或其转换形式）**作为预测变量，同时权重应为效应量方差的倒数（即1/SE²）。你的代码里存在两处关键错误：

1. 权重用了Effectsize_Va^2（方差的平方），这完全不符合权重的设定逻辑，权重应该是效应量的方差（如果Effectsize_Va是方差）或者方差的倒数（如果Effectsize_Va是标准误）。
2. 预测变量用了Effectsize_Va，但如果Effectsize_Va是方差，正确的预测变量应该是标准误（sqrt(Effectsize_Va)）；如果是标准误，直接用它作为预测变量即可。

修正后的Egger检验代码

假设Effectsize_Va是效应量的方差，正确代码如下：

# 先计算标准误
Data_DSMIV$SE <- sqrt(Data_DSMIV$Effectsize_Va)

# 构建Egger检验的多水平模型：权重为1/SE²（即方差的倒数），预测变量为SE
egger_multi <- rma.mv(Effectsize, 1/SE^2, 
                      random = ~ 1 | Study_Id/Effectsize_Id, 
                      mods = ~ SE, 
                      data = Data_DSMIV)

# 用CR2方差估计做系数检验，适配多水平结构
coef_test(egger_multi, vcov = "CR2")

如果Effectsize_Va是标准误，则代码调整为：

# 权重为1/SE²，预测变量为SE
egger_multi <- rma.mv(Effectsize, 1/Effectsize_Va^2, 
                      random = ~ 1 | Study_Id/Effectsize_Id, 
                      mods = ~ Effectsize_Va, 
                      data = Data_DSMIV)

coef_test(egger_multi, vcov = "CR2")

漏斗图部分优化 #

你的漏斗图代码逻辑没问题，用Study_Id分组着色能清晰区分不同研究的效应量。可以额外添加Egger检验的回归直线，让偏倚趋势更直观：

library(metafor)

dat <- Data_DSMIV
res <- rma.mv(Effectsize, Effectsize_Va, 
              random = ~ 1 | Study_Id/Effectsize_Id, 
              data=dat)

cols <- palette.colors(length(unique(dat$Study_Id)), palette="Alphabet")
cols <- cols[as.numeric(factor(dat$Study_Id))]

# 绘制带分组颜色的漏斗图，并添加Egger回归直线
funnel(res, col=cols)
abline(egger_multi$beta[1], egger_multi$beta[2], lwd=2, col="red")

其他注意事项 #

• 确认Effectsize_Va的定义：是方差还是标准误，这直接影响模型的权重和预测变量设定，必须对应正确。
• 你使用的~1 | Study_Id/Effectsize_Id随机效应结构是合理的，适配研究-效应量的两层嵌套结构。
• 多水平Egger检验的结果解读和普通Egger检验一致：关注mods项的系数显著性，若显著则提示存在发表偏倚。

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