Metafor 2.5-101开发版多级元分析自由度调整问题咨询

问题描述 #

编辑：已修改代码，添加了test = "t"参数
我希望进一步了解Metafor 2.5-101开发版本如何调整多级模型的自由度，以防控I类错误。
我对该功能的认知来源于Nakagawa的预印本《生态与进化元分析发表偏倚检验方法》及其配套的Supplemental_Impleentation_Example.Rmd文件，参考了其中第133-142行的内容：

在进行相关校正操作前，用户需要知晓最新版metafor（2.5-101版本）现已支持降低I类错误风险。计算自由度时可不再使用效应量数量，转而采用总文献数。我们的模拟结果显示“文献数-1”的自由度设置效果较好。安装metafor开发版后可按如下方式实现：

mod_multilevel_pdf = rma.mv(yi = yi, V = vi, mods = ~1, 
                                    random=list(~1|study_id,~1|obs), 
                                    data=data, test="t", dfs = "contain")
summary(mod_multilevel_pdf)

可以看到模型的df从149变为29，p值也对应完成调整。
按照我的理解，此处模型显示的df为29，对应30篇原始文献减1，而非30篇文献每篇5个效应量共150个减1。
我将该方法适配到自有代码中，我的数据集包含18篇文献、共24个效应量，按预期使用上述代码后自由度应该调整为17（18篇文献减1），但我最初误以为得到的df为23（24个效应量减1）。我运行的代码如下：

mod_multilevel_pdf = rma.mv(yi = yi, V = vi, mods = ~1, 
                            random=list(~1|study_id,~1|es_id), 
                            data=dat, test="t", dfs = "contain")
summary(mod_multilevel_pdf)

输出结果为：

Multivariate Meta-Analysis Model (k = 24; method: REML)

  logLik  Deviance       AIC       BIC      AICc 
-30.2270   60.4540   66.4540   69.8604   67.7171   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.6783  0.8236     18     no  study_id 
sigma^2.2  0.1416  0.3763     24     no     es_id 

Test for Heterogeneity:
Q(df = 23) = 167.2145, p-val < .0001

Model Results:

estimate      se     tval  df    pval    ci.lb   ci.ub 
 -0.3508  0.2219  -1.5809  17  0.1323  -0.8190  0.1174    

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我对此现象十分困惑，恳请大家提供帮助，非常感谢！

解答 #

你的运行结果完全符合metafor开发版dfs="contain"参数的设计逻辑，不存在异常，困惑来源于混淆了两类不同检验的自由度计算规则：

• 异质性Q检验的自由度固定为总效应量数减1，也就是你输出中看到的23，该计算逻辑不受dfs参数调整。Q检验的本质是检验所有效应量的变异是否显著高于抽样误差，本身就需要基于所有效应量的数量计算自由度，这部分没有问题。
• 模型固定效应（你这里设置的是仅包含截距，即总平均效应）的t检验自由度，才是dfs="contain"参数的调整对象。你输出中模型结果部分的自由度为17，正好对应18篇独立研究减1，完全匹配Nakagawa提到的校正逻辑，目的就是避免单个研究贡献多个效应量时高估自由度，进而控制I类错误的膨胀。
• dfs="contain"的工作逻辑是自动识别多层模型的最高聚类层级，也就是你随机效应结构里的study_id层，取该层级的组数减1作为系数检验的自由度，你输出的方差组件部分已经明确标注study_id的nlvls为18，所以最终自由度为17，完全符合你的预期。
  你最初误以为会得到23，是把异质性检验的自由度和固定效应系数检验的自由度搞混了，属于常见的认知误区。如果需要验证校正效果，可以删掉dfs = "contain"参数重新运行模型，此时默认采用残差法计算自由度，你就会看到固定效应检验的自由度变为23，和你最初误以为的结果一致。

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