Блог

В пакете метаанализа metafor в R есть функция blup() для rma.uni(), но нет rma.mv (). У меня есть rma.mv модель со случайным термином перехвата. Если я использую:

 predict.rma(model,transf=exp)
 

затем я получаю оценку от моих модераторов с фиксированным эффектом для каждой точки данных в моем исходном наборе данных и:

 ranef(model,transf=exp)
 

дает мне прогнозы для каждого уровня моего случайного эффекта.

Есть ли тогда способ объединить информацию из этих двух функций, чтобы получить комбинированные случайные и фиксированные эффекты, как это было бы предусмотрено функцией blup ()?

(Я попытался получить среднее значение как для оценки фиксированного эффекта, так и для прогнозирования случайного эффекта для каждой точки данных в моем исходном наборе данных, и это выглядит правильно, когда я его строю… но, конечно, это не так просто?)

Вы можете добавить (не принимайте среднее значение) то, что predict() дает вам, к тому, что ranef() дает вам, но вы должны быть осторожны в правильном «выстраивании» блефов со строками из predict() . Имена строк из ranef() указывают вам уровни, для которых вычисляются провалы, так что вы можете использовать их для правильного сопоставления. Кроме того, сначала добавьте их, затем вы можете возвести значения в степень, а не наоборот.

В процессе следования очень полезным инструкциям Вольфганга в отмеченном ответе я создал некоторый код, чтобы сделать это с помощью примера набора данных. Публикация на случай, если это поможет кому-то еще:

Обновление: следуя комментарию Вольфганга ниже, я удалил код для CIs и заменил его кодом для вычисления SE для BLUP.

 #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()