Model Selection using the glmulti and MuMIn Packages with a rma.mv() Model

rcode_rma_mv_example.r

############################################################################

library(metafor)
library(ape)

############################################################################

# read the documentation for this dataset
help(dat.moura2021)

# get the data and the tree
dat  <- dat.moura2021$dat
tree <- dat.moura2021$tree

# calculate r-to-z transformed correlations and corresponding sampling variances
dat <- escalc(measure="ZCOR", ri=ri, ni=ni, data=dat)

# make the tree ultrametric and compute the phylogenetic correlation matrix
tree <- compute.brlen(tree)
A <- vcv(tree, corr=TRUE)

# make a copy of the species.id variable
dat$species.id.phy <- dat$species.id

# fit the full model (multilevel phylogenetic meta-analytic model)
full <- rma.mv(yi, vi, mods = ~ spatially.pooled * temporally.pooled,
   random = list(~ 1 | study.id, ~ 1 | effect.size.id,
                 ~ 1 | species.id, ~ 1 | species.id.phy),
   R=list(species.id.phy=A), data=dat, method="ML")
full

############################################################################

# model selection using glmulti

library(glmulti)

rma.glmulti <- function(formula, data, ...) {
   rma.mv(formula, vi,
          random = list(~ 1 | study.id, ~ 1 | effect.size.id,
                        ~ 1 | species.id, ~ 1 | species.id.phy),
          R=list(species.id.phy=A), data=data, method="ML", ...)
}

# fit all possible models; since level=2, the two-way interaction between the
# two predictors is also considered; and by setting marginality=TRUE the model
# with the interaction must include the two main effects; this leads to a
# total of 5 possible models
system.time(res1 <- glmulti(yi ~ spatially.pooled + temporally.pooled, data=dat,
                level=2, marginality=TRUE, fitfunction=rma.glmulti,
                crit="aicc", confsetsize=5, plotty=FALSE))

# short output
print(res1)

# table with the information criteria for each model
weightable(res1)

# multimodel inference
eval(metafor:::.glmulti)
round(coef(res1, varweighting="Johnson"), 4)

# process the output into a more familiar form
mmi <- as.data.frame(coef(res1, varweighting="Johnson"))
mmi <- data.frame(Estimate=mmi$Est, SE=sqrt(mmi$Uncond),
                  Importance=mmi$Importance, row.names=row.names(mmi))
mmi$z <- mmi$Estimate / mmi$SE
mmi$p <- 2*pnorm(abs(mmi$z), lower.tail=FALSE)
names(mmi) <- c("Estimate", "Std. Error", "Importance", "z value", "Pr(>|z|)")
mmi$ci.lb <- mmi[[1]] - qnorm(.975) * mmi[[2]]
mmi$ci.ub <- mmi[[1]] + qnorm(.975) * mmi[[2]]
mmi <- mmi[order(mmi$Importance, decreasing=TRUE), c(1,2,4:7,3)]
round(mmi, 4)

############################################################################

# model selection using MuMIn

library(MuMIn)
eval(metafor:::.MuMIn)

# fit all possible models
system.time(res2 <- dredge(full, trace=2))
res2

# multimodel inference
summary(model.avg(res2))

# for easier comparison with the results from glmulti
round(mmi[colnames(model.avg(res2)$coefficients),], 4)

############################################################################

# MuMIn with parallel processing

library(parallel)
clust <- makeCluster(2, type="PSOCK")

clusterExport(clust, c("dat","A"))
clusterEvalQ(clust, library(metafor))

system.time(res3 <- dredge(full, trace=2, cluster=clust))
res3

stopCluster(clust)

############################################################################