Getting R squared from a mixed effects multilevel model in metafor

Question

I am doing a meta-analysis in R of a specific treatment on forests. For this model I need to fit random effects to account for between study differences in method and variation in age of sites, since both of these are confounding variables and I am not explicitly interested in investigating the variation caused by them.

However, as far as I can tell the package [metfor] does not permit you to calculate a R squared type statistic when you have a multilevel model.

Anyway, to describe my problem more clearly here is a mock dataset

Log<-data.frame(Method=rep(c("RIL","Conv"),each=10),
     RU=runif(n=20,min=10,max=50),SDU=runif(n=20,5,20),
     NU=round(runif(n=20,10,20),0))
Log$Study<-rep(1:4,each=5)
Log$Age<-rep(c(0,10,15,10),times=5)
RIL<-(Log$RU-(Log$RU*(abs(rnorm(n=20,mean=.6,sd=0.1)))))+(0.5*Log$Age)
Conv<-(Log$RU-(Log$RU*(abs(rnorm(n=20,mean=.2,sd=0.1)))))+(0.2*Log$Age)
Log$RL<-ifelse(Log$Method=="RIL",RIL,Conv)
Log$SDL<-Log$SDU
Log$NL<-Log$NU

#now we perform a meta-analysis using metafor
require(metafor)
ROM<-escalc(data=Log,measure="ROM",m2i=RU,
sd2i=SDU,n2i=NU,m1i=RL,sd1i=SDL,n1i=NL,append=T)
Model1<-rma.mv(yi,vi,random=~(1|Study)+(1|Age),method="ML",data=ROM)
summary(Model1)
forest(Model1)

The above model is a null model looking at whether the intercept is statistically significantly different from zero. In our case it is. However, I also want to see whether differences in treatments describe the differences in effect sizes I see on the forest plot that you can see here

So I run this model:

Model2<-rma.mv(yi,vi,mods=~Method,random=~(1|Study)+(1|Age),method="ML",data=ROM)
summary(Model2)

Which looks good.

    Multivariate Meta-Analysis Model (k = 20; method: ML)

  logLik  Deviance       AIC       BIC      AICc  
  0.4725   19.8422    7.0550   11.0380    9.7217  

Variance Components: 

outer factor: Age   (nlvls = 3)
inner factor: Study (nlvls = 4)

            estim    sqrt  fixed
tau^2      0.0184  0.1357     no
rho        1.0000             no

Test for Residual Heterogeneity: 
QE(df = 18) = 23.3217, p-val = 0.1785

Test of Moderators (coefficient(s) 2): 
QM(df = 1) = 19.6388, p-val < .0001

Model Results:

           estimate      se     zval    pval    ci.lb    ci.ub     
intrcpt     -0.1975  0.1007  -1.9622  0.0497  -0.3948  -0.0002    *
MethodRIL   -0.4000  0.0903  -4.4316  <.0001  -0.5768  -0.2231  ***

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

However, I would like to get a goodness of fit measure out of this model equivalent to an R squared. People have had these problems in the past with GLMMs but there are now ways of doing this. I was wondering if anyone knew of a good way of doing something similar with meta-analyses? I have reviewers asking for this and I'm not sure whether I should just tell them that it can't be done or not.

Thanks in advance for the help!

Answer 1

First of all, you are not quite using the right syntax for the rma.mv() function. For the two models, I assume you actually meant to use:

Model1 <- rma.mv(yi, vi, random = list(~ 1 | Study, ~ 1 | Age), method="ML", data=ROM)
Model2 <- rma.mv(yi, vi, mods = ~ Method, random = list(~ 1 | Study, ~ 1 | Age), method="ML", data=ROM)

Now, as for R-squared, you can compute the proportional reduction in the variance components as a sort of pseudo R-squared value. This is just the logical extension of what is typically done in regular meta-regression. So, based on the models above:

(Model1$sigma2[1] - Model2$sigma2[1]) / Model1$sigma2[1]
(Model1$sigma2[2] - Model2$sigma2[2]) / Model1$sigma2[2]

If a value should turn out to be negative, it is typically set to zero.

If you want a single value, you could also compute the proportional reduction in the total variance with:

(sum(Model1$sigma2) - sum(Model2$sigma2)) / sum(Model1$sigma2)