1

在 R 的元分析包 metafor 中，rma.uni() 有一个 blup() 函数，但 rma.mv() 没有。我有一个带有随机截距项的 rma.mv 模型。如果我使用：

predict.rma(model,transf=exp)

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

ranef(model,transf=exp)

为我的随机效应的每个级别提供预测。

那么有没有一种方法可以将来自这两个函数的信息组合起来以获得组合的随机和固定效果 BLUP，就像 blup() 函数提供的那样？

（我尝试对原始数据集中的每个数据点取固定效应估计和随机效应预测的平均值，当我绘制它时这看起来是正确的......但肯定不是那么简单吗？）

4

2 回答 2

1

您可以添加（不要吝啬）什么predict()给了你什么ranef()给了你，但你必须小心正确地将 BLUP 与predict(). 行名称ranef()告诉您计算 BLUP 的级别，因此您可以使用它们来正确匹配事物。此外，首先添加它们，然后您可以对这些值取幂，而不是相反。

于 2021-03-15T06:57:27.433 回答
1

在遵循 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()
于 2021-03-17T02:00:20.540 回答