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

# install the metafor package
#install.packages("metafor")

# install the clubSandwich package
#install.packages("clubSandwich")

# load the metafor package
library(metafor)

# load the clubSandwich package
library(clubSandwich)

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

# copy dataset to 'dat'
dat <- dat.landenberger2005

# compute log odds ratios and corresponding sampling variances
dat <- escalc(measure="OR", ai=n.cbt.non,  bi=n.cbt.rec,
                            ci=n.ctrl.non, di=n.ctrl.rec, data=dat)

# fit random-effects model
res <- rma(yi, vi, data=dat)
res

# compute average odds ratio (and corresponding 95% CI/PI)
predict(res, transf=exp, digits=2)

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

# copy dataset to 'dat'
dat <- dat.bornmann2007

# compute log odds ratios and corresponding sampling variances
dat <- escalc(measure="OR", ai=waward, n1i=wtotal,
                            ci=maward, n2i=mtotal, data=dat)

# within-study averaging assuming independent sampling errors
agg <- aggregate(dat, cluster=dat$study, struct="ID")

# fit random-effects model
res <- rma(yi, vi, data=agg)
res

# compute average odds ratio (and corresponding 95% CI/PI)
predict(res, transf=exp, digits=2)

# fit random-effects model
res <- rma(yi, vi, data=agg)
res

# compute average odds ratio (and corresponding 95% CI/PI)
predict(res, transf=exp, digits=2)

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

# compute average odds ratio (and corresponding 95% CI/PI)
predict(res, transf=exp, digits=2)

# compute the ICC
round(res$sigma2[1] / sum(res$sigma2), 2)

# fit model assuming homogeneity within studies
res <- rma.mv(yi, vi, random = ~ 1 | study, data=dat)
res

# compute average odds ratio (and corresponding 95% CI/PI)
predict(res, transf=exp, digits=2)

# LRT comparing the two models
res0 <- rma.mv(yi, vi, random = ~ 1 | study, data=dat)
res1 <- rma.mv(yi, vi, random = ~ 1 | study/obs, data=dat)
anova(res0, res1)

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

# copy data into 'dat' and examine data
dat <- dat.berkey1998
dat

# construct the variance-covariance matrices of the observed outcomes
V <- lapply(split(dat[c("v1i", "v2i")], dat$trial), as.matrix)
V <- bldiag(V)
V

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

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

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

# fit random-effects model for the two outcomes separately
res.AL <- rma(yi, vi, data=dat, subset=outcome=="AL")
res.AL
res.PD <- rma(yi, vi, data=dat, subset=outcome=="PD")
res.PD

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

# copy data into 'dat' and examine first 10 rows
dat <- dat.assink2016
head(dat, 10)

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

# use cluster-robust inference methods
robust(res, cluster=dat$study)

# or use coef_test() from clubSandwich
coef_test(res, vcov="CR2", cluster=dat$study)

# construct approx V matrix assuming r=0.6 for the sampling errors
V <- impute_covariance_matrix(dat$vi, cluster=dat$study, r=0.6)

# fit multivariate model
res <- rma.mv(yi, V, random = ~ factor(esid) | study, data=dat)
res

# use cluster-robust inference methods
robust(res, cluster=dat$study)

# or use coef_test() from clubSandwich
coef_test(res, vcov="CR2", cluster=dat$study)

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