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

# R script to go along with the talk:
#
# Viechtbauer, W. (2023). Location-scale models for meta-analysis using the
# metafor package. Evidence Synthesis & Meta-Analysis in R Conference 2023.

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

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

# load the metafor package
library(metafor)

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

# look at the BCG vaccine meta-analysis dataset by Colditz et al. (1994)
dat.bcg

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

# fit random-effects model (using log risk ratios and variances as input)
res <- rma(yi, vi, data=dat)
res

# dichotomize 'alloc' variable
dat$random <- ifelse(dat$alloc == "random", 1, 0)

# fit meta-regression model
res <- rma(yi, vi, mods = ~ random, data=dat)
res

# fit meta-analytic location-scale model
res <- rma(yi, vi, mods = ~ random, scale = ~ random, data=dat)
res

# fit separate random-effects models within subgroups
res <- list(rma(yi, vi, data=dat, subset=random==0),
            rma(yi, vi, data=dat, subset=random==1))
tab <- data.frame(k        = sapply(res, \(x) x$k),
                  estimate = sapply(res, coef),
                  se       = sqrt(sapply(res, vcov)),
                  tau2     = sapply(res, \(x) x$tau2))
rownames(tab) <- c("non-random", "random")
round(tab, digits=4)

# Wald-type test
res1 <- rma(yi, vi, mods = ~ random, scale = ~ random, data=dat)
res1

# likelihood ratio test
res0 <- rma(yi, vi, mods = ~ random, scale = ~ 1, data=dat)
anova(res0, res1)

# permutation test
permutest(res1, seed=1234)

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

# a more elaborate example

# look at the dataset from the meta-analysis by Bangert-Drowns et al. (2004)
dat <- dat.bangertdrowns2004
dat

# collapse subjects down to three main groups
dat$subj <- factor(ifelse(dat$subject %in% c("Algebra", "Calculus", "Math", "Math in Science", "Statistics"), "Math",
                   ifelse(dat$subject %in% c("Chemistry", "Comp Science and Chemistry", "Biology", "Earth Science", "Natural Resources", "Nursing", "Science"), "Sci", "Soc")))

# fit meta-analytic location-scale model
res <- rma(yi, vi, mods = ~ subj + ni, scale = ~ subj + ni, data=dat)
res

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