Skip to content

Instantly share code, notes, and snippets.

@wviechtb
Last active March 8, 2024 18:45
Show Gist options
  • Save wviechtb/891483eea79da21d057e60fd1e28856b to your computer and use it in GitHub Desktop.
Save wviechtb/891483eea79da21d057e60fd1e28856b to your computer and use it in GitHub Desktop.
Model Selection using the glmulti and MuMIn Packages with a rma.mv() Model
############################################################################
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)
############################################################################
@Yefeng0920
Copy link

Hi Wolfgang,
Thank you very much for providing this useful and excellent example. I really learn a lot from your open science practices (e.g., sharing code, showing examples). I just have a request - If your time permits, I would be grateful if you would like to show an example of using model selection to select random-effects structure in the context of multivariate/multilevel models (definitely by rma.mv()). The random-effects structure is really worth noting when conducting a meta-analysis. Many papers have touched upon this topic, like

Barr D J, Levy R, Scheepers C, et al. Random effects structure for confirmatory hypothesis testing: Keep it maximal[J]. Journal of memory and language, 2013, 68(3): 255-278.
Bates D, Kliegl R, Vasishth S, et al. Parsimonious mixed models[J]. arXiv preprint arXiv:1506.04967, 2015.
Matuschek H, Kliegl R, Vasishth S, et al. Balancing Type I error and power in linear mixed models[J]. Journal of memory and language, 2017, 94: 305-315.
But we know a little in the context of the meta-analytic models. I am very curious about your answers and your practices to deal with random-effects structure when conducting a (complex) meta-analysis.

Best,
Yefeng Yang PhD
Research Associate
University of New South Wales - Australia
Biological, Earth & Environmental Science

@wviechtb
Copy link
Author

Not sure if/when I will get to this. In most cases, the random effects structure is something that is determined by the data structure and the dependencies that need to be handled/accounted for. In most cases, I would not recommend making any changes to that anyway once the structure has been chosen.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment