如何在 metafor 中计算 rma.mv 的组合随机和固定效应 BLUP?

How to calculate combined random and fixed effect BLUPs for rma.mv in metafor?

在 R 的元分析包 metafor 中,有一个 blup() 函数用于 rma.uni() 而不是 rma.mv()。我有一个带有随机截距项的 rma.mv 模型。如果我使用:

predict.rma(model,transf=exp)

然后我从我的固定效应调节器得到我原始数据集中每个数据点的估计值,并且:

ranef(model,transf=exp)

给我随机效应的每个级别的预测。

有没有办法将来自这两个函数的信息结合起来,得到组合的随机和固定效果 BLUP,就像 blup() 函数提供的那样?

(我已经尝试对原始数据集中的每个数据点取固定效应估计和随机效应预测的平均值,而这个 在我绘制它时看起来 是正确的。 .. 但肯定没那么简单?)

您可以将 predict() 给您的内容加到 ranef() 给您的内容上(不要取平均值),但是您必须小心正确 'lining up' BLUPs来自 predict() 的行。 ranef() 中的行名称告诉您计算 BLUP 的级别,因此您可以使用它们来正确匹配。另外,首先添加它们,然后你可以对值取幂,而不是相反。

在遵循标记答案中 Wolfgang 非常有用的说明的过程中,我创建了一些代码来使用示例数据集执行此操作。发帖以防对其他人有帮助:

更新:根据 Wolfgang 在下面的评论,我删除了 CI 的代码并用计算 BLUP 的 SE 的代码替换它。

#load packages and example data    
library(metafor)
library(plyr)
library(ggplot2)
dat<-dat.konstantopoulos2011
dat

#make a model - note I have nested random effects
#(no idea if it actually makes sense with this example dataset!)
mod<-rma.mv(yi,vi,mods=~year,random=~1|district/school,data=dat)

#predict from the model (without exponentiating) and attach to original data
preds<-predict.rma(mod,addx=TRUE)
dat$pred<-preds$pred
dat$ci.ub<-preds$ci.ub
dat$ci.lb<-preds$ci.lb
dat$fe_se<-preds$se
dat$dist.sch<-interaction(dat$district,dat$school) #create a label for each district/school

#get the district random effects and label them
dist_re<-data.frame(ranef(mod)$district)
dist_re$district<-rownames(dist_re)

#get the school random effects and label them
sch_re<-data.frame(ranef(mod)$`district/school`)
sch_re$dist.sch<-rownames(sch_re)
sch_re$dist.sch<-gsub("/",".",sch_re$dist.sch)
colnames(sch_re)<-c("intrcpt2","se2","pi.lb2","pi.ub2","dist.sch") #to avoid duplicate colnames later

#join the district and school random effects to the data by labels
plotdat<-join(dat,dist_re,by="district")
plotdat2<-join(plotdat,sch_re,by="dist.sch")

#create the blups and intervals by adding the fixed effect estimates and random effect predictions,
#and exponentiating:
plotdat2$blup<-exp(plotdat2$pred+plotdat2$intrcpt+plotdat2$intrcpt2)
plotdat2$blup.se<-sqrt((plotdat2$fe_se^2)+(plotdat2$se^2)+(plotdat2$se2^2))

#forest plot of BLUPs and their SEs just to check they make sense:
ggplot(plotdat2, aes(y=dist.sch, x=blup, xmin=blup-blup.se, xmax=blup+blup.se))+
  geom_point()+
  geom_errorbarh(height=.2)+
  ylab('District and school')+
  geom_vline(xintercept=1,linetype='dashed')+
  xlim(0,5)+
  theme_bw()