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:
library(metafor)
dat<- dat.konstantopoulos2011
dat$year <-rnorm(nrow(dat))
dat$yi <- dat$yi + dat$year/3 + dat$district/200
dat$study <- factor(dat$study)
dat$district<-factor(dat$district)
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):
#libraries
library(RCurl)
library(glmmTMB)
library(nlme)
#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)]
studies<-studies[order(studies$order),]
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
vi2<-1/vi
vi2<-vi2/sum(vi2)
#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.
#Plotting:
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)
#models
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))
#coefficients
c(coef(M)[2], summary(M2)$coefficients$cond[2,1], M3$coefficients$fixed[2])
#randoms:
sqrt(M$sigma2)
VarCorr(M2) #order reversed in comparison to the other 2
VarCorr(M3)
# 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
