CODEWITHSUNDEEP
rmetafor
corn_bunting
corn_bunting

Reputation: 381

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

In the meta-analysis package metafor in R, there is a blup() function for rma.uni() but not rma.mv(). I have an rma.mv model with a random intercept term. If I use:

predict.rma(model,transf=exp)

then I get an estimate from my fixed effect moderators for each datapoint in my original dataset, and:

ranef(model,transf=exp)

gives me predictions for each level of my random effect.

Is there then a way to combine the information from these two functions to get combined random and fixed effect BLUPs, as would be provided by the blup() function?

(I have tried taking the mean of both the fixed effect estimate and the random effect prediction for each datapoint in my original dataset, and this looks right when I plot it ... but surely it's not that simple?)

Upvotes: 1

Views: 746

Answers (2)

corn_bunting
corn_bunting

Reputation: 381

In the process of following Wolfgang's very helpful instructions in the marked answer, I created some code to do this with an example dataset. Posting in case it helps anyone else:

Update: following Wolfgang's comment below, I have removed the code for the CIs and replaced it with code to calculate the SE for the 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()

Upvotes: 1

Wolfgang
Wolfgang

Reputation: 3395

You can add (don't take the mean) what predict() gives you to what ranef() gives you, but you have to be careful in correctly 'lining up' the BLUPs with the rows from predict(). The row names from ranef() tell you the levels for which the BLUPs are computed, so you can use them to match things up correctly. Also, first add them, then you can exponentiate the values, not the other way around.

Upvotes: 1

Related Questions