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

# R script to go along with the talk:
#
# Viechtbauer, W. (2024). Model selection in complex meta-analyses. Symposium
# on 'Recent Advances in Meta-Analysis', Department of Medical Statistics,
# University Medical Center Göttingen, Germany.

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

# install (if not already installed) the metafor package
#install.packages("metafor")

# load the metafor package
library(metafor)

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

# Example: Relationship between Class Attendance and Grade

# copy the dataset to dat
dat <- dat.crede2010

# select rows from the dataset that report the correlation between class
# attendance and performance (grade) in the class
dat <- subset(dat, criterion=="grade")

# compute r-to-z transformed correlations and corresponding sampling variances
dat <- escalc(measure="ZCOR", ri=ri, ni=ni, data=dat,
              slab=paste("Study:", studyid, " / Sample:", sampleid))

# forest plot
forest(dat$yi, dat$vi, xlim=c(-2.2,3.3), cex=0.7, psize=1, efac=0,
       header="Study Id / Sample Id", digits=c(2,1),
       shade=dat$studyid %in% names(which(table(dat$studyid) > 1)))

# fit multilevel model
res <- rma.mv(yi, vi, random = ~ 1 | studyid/sampleid, data=dat)
res

# back-transform results
predict(res, transf=transf.ztor, digits=2)

# estimated the intraclass correlation coefficient
round(res$sigma2[1] / sum(res$sigma2), digits=2)

# test H0: sigma^2_1 = 0 (same as testing H0: rho = 0)
res0 <- update(res, sigma2=c(0,NA))
anova(res0, res)

# test H0: sigma^2_2 = 0 (same as testing H0: rho = 1)
res0 <- update(res, sigma2=c(NA,0))
anova(res0, res)

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

# Example: Handedness and Eye-Dominance

# copy the dataset to dat
dat <- dat.bourassa1996

# restructure the dataset to keep only the male/female data when it is reported
# separately and the combined data when this is the only data reported
dat <- escalc(measure="OR", ai=rh.re, bi=rh.le,
                            ci=lh.re, di=lh.le,
              data=dat, add=1/2, to="all")
dat <- lapply(split(dat, dat$sample), function(x) {
   if (nrow(x) == 3L) {
      x[-which(x$sex == "combined"),]
   } else {
      x
   }
})
dat <- do.call(rbind, dat)
rownames(dat) <- NULL
dat

# forest plot
forest(dat$yi, dat$vi, xlim=c(-17,14), at=seq(-4,8,by=2), cex=0.6, psize=1, efac=0, header="Study",
       shade=as.numeric(factor(dat$study)) %% 2 != 0, slab=paste("Study:", dat$study),
       ilab=cbind(paste("Sample:", dat$sample), dat$sex), ilab.xpos=c(-13,-8.5), ilab.pos=4,
       ilab.lab=c("Sample","Sex"))

# fit multilevel model
res <- rma.mv(yi, vi, random = ~ 1 | study/sample/sex, data=dat)
res

# back-transform results
predict(res, transf=exp, digits=2)

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

## Example: BCG Vaccine against Tuberculosis

# copy the dataset to dat
dat <- dat.bcg

# examine the dataset
dat

# compute log odds ratios and corresponding sampling variances
dat <- escalc(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

# fit a standard random-effects model
rma(yi, vi, data=dat)

# restructure into long format
dat2 <- to.long(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat)
dat2$group <- factor(dat2$group, levels=c(2,1), labels=c("con","exp"))
dat2$id <- 1:nrow(dat2)

# compute log odds for each treatment arm
dat2 <- escalc(measure="PLO", xi=out1, mi=out2, data=dat2)

# fit model with fixed study effects
rma.mv(yi, vi, mods = ~ 0 + factor(study) + group, random = ~ 1 | id, data=dat2)

# fit model with random study effects
rma.mv(yi, vi, mods = ~ group, random = ~ 1 | study/group, data=dat2)

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

## Example: Treatment for Periodontal Disease

# copy the dataset to dat
dat <- dat.berkey1998

# examine the dataset
dat

# construct the variance-covariance matrix of the observed outcomes
V <- vcalc(vi=1, cluster=author, rvars=c(v1i, v2i), data=dat)
V

# fit the multivariate model
res <- rma.mv(yi, V, mods = ~ outcome - 1, data = dat,
              random = ~ outcome | trial, struct = "UN")
res

# estimate the difference between the two pooled effects
predict(res, newmods=c(1,-1))

# test the difference between the two pooled effects
anova(res, L=c(1,-1))

# test the correlation among the true effects
res0 <- update(res, rho=0)
anova(res0, res)

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

## Example: Association between Recidivism and Mental Health

# copy the dataset to dat
dat <- dat.assink2016

# examine the dataset
dat

# assume a correlation of 0.7 for effect sizes corresponding to the same
# type of delinquent behavior and a correlation of 0.5 for effect sizes
# corresponding to different types of delinquent behavior
V <- vcalc(vi, cluster=study, type=deltype, obs=esid,
           data=dat, rho=c(0.7, 0.5))

# fit multilevel model using this approximate V matrix
res <- rma.mv(yi, V, random = ~ 1 | study/esid, data=dat)
res

# apply cluster-robust inference methods
robust(res, cluster=study, clubSandwich=TRUE)

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

## Example: Effects of Mycorrhizal Fungi Inoculation

# this example is skipped here because the data preparation takes quite a bit
# more effort and the data are not available via the metadat package; but for
# those interested, the data are freely available; for details, see:
#
# Chaudhary, V. B., Rúa, M. A., Antoninka, A., Bever, J. D., Cannon, J.,
# Craig, A., Duchicela, J., Frame, A., Gardes, M., Gehring, C., Ha, M., Hart,
# M., Hopkins, J., Ji, B., Collins Johnson, N., Kaonongbua, W., Karst, J.,
# Koide, R. T., Lamit, L. J., Meadow, J., Milligan, B. G., Moore, J. C.,
# Pendergast IV, T. H., Piculell, B., Ramsby, B., Simard, S., Shrestha, S.,
# Umbanhowar, J., Viechtbauer, W., Walters, L., Wilson, G. W. T., Zee, P. C.,
# & Hoeksema, J. D. (2016). MycoDB, a global database of plant response to
# mycorrhizal fungi. Scientific Data, 3, 160028.
# https://doi.org/10.1038/sdata.2016.28

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

## Example: Generation Effect

# copy the dataset to dat
dat <- dat.mccurdy2020

# fit multilevel model
res <- rma.mv(yi, vi, mods = ~ condition,
              random = list(~ 1 | article/experiment/sample/id, ~ 1 | pairing),
              data=dat, sparse=TRUE)
res

# apply cluster-robust inference methods
robust(res, cluster=article, clubSandwich=TRUE)

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