I am conducting a meta-analysis and I am struggling with the random structure when there is a continuous moderator and a nested random term. Starting with a simple example:

dat<- dat.konstantopoulos2011
dat$year <-rnorm(nrow(dat))
dat$yi <- dat$yi + dat$year/3 + dat$district/200
dat$study <- factor(dat$study)

This dataset has 11 districts, with 3-11 studies per district. There is one effect size per study. dat$yi are the effect sizes, which correlate with the continous moderator dat$year (slope 0.33) and have varying intercepts for the 11 districts. Is the following model correct?

rma.mv(yi ~ year, vi, random = ~ 1 | district/study, data=dat)

I just want to clarify because the usual model in ecology would be:

lme(yi ~ year, random = ~1|district, data=dat)  #and VarFixed (~ vi), lmecontrol (sigma=1)

corresponding in metafor to

rma.mv(yi ~ age, vi, random = ~ 1 | district, data=dat)

I can intuitively understand that we usually do not want to have all variance ending up in the “study” term, but that this is different in a meta-analysis were variances are known exactly. Just want to make sure that the nesting of study in district is correct.

My actual data for the meta-analysis is a bit more complicated. It consists of effect sizes (mean diapause date) from 447 populations, extracted from 57 studies. The 57 studies were conducted on 46 species of 32 genera in 9 orders. There is a single continuous moderator. A full random term would thus be order/genus/species/study/population. I plan to drop the term study, because there is almost always only one study per species, except for a few cases where the same authors conducted several studies with equal methods on the same species. I’m also thinking about dropping the term genus, as most species come from different genera. This would make the random term order/species/population with sample sizes of 9/46/447. Or would it be only order/species ? The model first seems fine, but to calculate an R² value I need to use a null model with the moderator dropped, and in that case the term order suddenly explains zero variance. Here is the script so far (including access the raw data):


#load data
url <- getURL('https://raw.githubusercontent.com/JensJoschi/variability_timing/master/lit_extract/mcmcresults.txt')
studies <- read.table(text=url, header = TRUE)
studies <- studies[,-c(2:4,6:18,23:25,28,30:32)]
r<-studies$upper_e-studies$lower_e  #credible interval range
r[r<(1/6)]<-1/6 #prevents studies from getting infinite weight
vi<-r    #CI should be (r / (2*1.96))^2  but perhaps this is sufficient for demonstration purposes
#Order, genus, spec, ID and popid are the terms for nesting, med_e the effect sizes, 
#vi the variances, and degN is a moderator (latitude). 
#Vi2 is a scaled inverse variance needed for glmmTMB.

plot(studies$med_e ~ studies$degN, pch=21, col=NA, bg = studies$order)
segments(x0=studies$degN,y0 = studies$med_e-vi/2, y1 = studies$med_e+vi/2,col=studies$order)

M<- rma.mv(med_e ~ degN, vi, random= ~ (1|order/spec/popid), data=studies)
M2<-glmmTMB(med_e~degN +  (1|order/spec/popid),weights = vi2, data= studies, dispformula = ~0)
M3<-lme(med_e~degN, random = ~1|order/spec/popid, weights = varFixed(~vi), data= studies, control = lmeControl(sigma=1))

#null models
M_null<- rma.mv(med_e, vi, random= ~ (1|order/spec/popid), data=studies)
M2_null<-glmmTMB(med_e ~1 +  (1|order/spec/popid),weights = vi2, data= studies, dispformula = ~0)
M3_null<-lme(med_e~1, random = ~1|order/spec/popid, weights = varFixed(~vi), data= studies, control = lmeControl(sigma=1))

c(coef(M)[2], summary(M2)$coefficients$cond[2,1], M3$coefficients$fixed[2])
VarCorr(M2) #order reversed in comparison to the other 2

# R² values (metafor only)
(M_null$sigma2-M$sigma2)/M_null$sigma2 # -4435754, 0.67 and 0.79
(sum(M_null$sigma2)-sum(M$sigma2))/sum(M_null$sigma2) #0.54

I now wonder about the 0 variance term of order. Is this because popid should not be part of the random term, or did I do something else fundamentally wrong in my models? Given that the models are correct, can I use the R²-values (reporting as 0, 0.67 and 0.79; and 0.54 overall)?

Lastly, I wonder why glmmTMB always gives different estimates, no matter which random terms I use. Is there something wrong with my use of the function? I will need it later because one of my effect sizes is beta- distributed. I would really appreciate if someone with more expertise could check the models.

Further background on the study is here


1 Answer 1


Thanks for providing a reproducible example. Indeed, at first sight, it seems a bit odd that the order random effect drops to 0 when removing the degN predictor. But order has only 9 levels, so I would not expect very precise estimates of the corresponding variance component. Let's see:

res1 <- rma.mv(med_e ~ degN, vi, random = ~ 1 | order/spec/popid, data=studies)
res0 <- rma.mv(med_e, vi, random = ~ 1 | order/spec/popid, data=studies)

confint(res1, sigma2=1)
confint(res0, sigma2=1)

This yields:

> confint(res1, sigma2=1)

          estimate  ci.lb  ci.ub
sigma^2.1   0.3629 0.0362 1.7155
sigma.1     0.6024 0.1902 1.3098

> confint(res0, sigma2=1)

          estimate  ci.lb  ci.ub
sigma^2.1   0.0000 0.0000 1.6463
sigma.1     0.0003 0.0000 1.2831

So, the 95% profile likelihood CI for this variance component is very wide and quite similar for the two models. So, I would say it is alright to use:

pmax(0, (res0$sigma2 - res1$sigma2) / res0$sigma2)
sum(res0$sigma2 - res1$sigma2) / sum(res0$sigma2)

to get the pseudo-R^2 values.

As for glmmTMB(): First of all, you should use REML=TRUE since REML estimation is also the default for rma.mv(). But the results are still different. Actually, whether you use weights = vi2 or not makes no difference here to what glmmTMB() returns. I am not an expert on glmmTMB(), but I think the meaning of weights is really different here than how the sampling variances are used by rma.mv() (and lme()). My guess is that the weights are applied to the residual variance component, but since you are basically forcing that component to 0, it makes no difference what weights you specify. So, I think you cannot fit the same model with glmmTMB() that you are fitting with rma.mv().

  • $\begingroup$ Thank you very much for your clarification, this helped a lot and I got sensible output. I am still stuck with the beta-distributed effect size, but I opened it as new question, because the topic is slightly different (not about random terms but about choice of functions/approaches): stats.stackexchange.com/questions/405178/… $\endgroup$
    – Jens
    Apr 26, 2019 at 8:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.