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

# Code to go along with the talk "Confidence intervals for the amount of
# heterogeneity accounted for in meta-regression models" by Wolfgang
# Viechtbauer (https://www.wvbauer.com) presented at the 2023 European
# Congress of Methodology on 2023-06-13

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

### standard linear regression

# examine the mtcars dataset
mtcars

# linear regression model to predict mpg (miles per gallon) from:
# - cyl (number of cylinders)
# - hp  (gross horsepower)
# - wt  (weight in 1000 lbs)
# - am  (transmission: 0 = automatic, 1 = manual)
res <- lm(mpg ~ cyl + hp + wt + am, data=mtcars)
summary(res)

# extract the R^2 value from the model object
summary(res)$r.squared

# illustrate the computation of R^2
(var(mtcars$mpg) - var(resid(res))) / var(mtcars$mpg)

# CI for R^2 by test inversion
#install.packages("confintr")
library(confintr)
fstat <- summary(res)$fstatistic
ci_rsquared(fstat["value"], df1=fstat["numdf"], df2=fstat["dendf"])

# illustrate the computation of the CI
pf(fstat[["value"]], df1=fstat[["numdf"]], df2=fstat[["dendf"]], ncp=69.054285)
69.054285 / (69.054285 + fstat[["numdf"]] + fstat[["dendf"]] + 1)
pf(fstat[["value"]], df1=fstat[["numdf"]], df2=fstat[["dendf"]], ncp=256.15667)
256.15667 / (256.15667 + fstat[["numdf"]] + fstat[["dendf"]] + 1)

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

### mixed-effects meta-regression model

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

# load the metafor package
library(metafor)

# results from 58 studies on the effectiveness of cognitive-behavioral
# therapy (CBT) for reducing recidivism (Landenberger & Lipsey, 2005)
dat <- dat.landenberger2005

# look at a subset of the data
dat[c(1:6,56:58),c(1,8:11,18)]

# n.ctrl.rec - number of recidivists in the control group
# n.ctrl.non - number of non-recidivists in the control group
# n.cbt.rec  - number of recidivists in the CBT group
# n.cbt.non  - number of non-recidivists in the CBT group
#
# these variables denote the values in 2x2 tables of the form:
#
#           Non-Rec      Rec
#         +------------+------------+
# CBT     | n.cbt.non  | n.cbt.rec  |
#         +------------+------------+
# control | n.ctrl.non | n.ctrl.rec |
#         +------------+------------+
#
# session - number of CBT sessions per week

# calculate log odds ratios (for non-recidivism in CBT vs. control groups)
# and the 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 a random-effects model (log odds ratios and variances as input)
res1 <- rma(yi, vi, data=dat)
res1

# fit a mixed-effects meta-regression model with number of sessions
# as a moderator (log odds ratios and variances as input)
res2 <- rma(yi, vi, mods = ~ sessions, data=dat)
res2

# extract the R^2 value from the model object
summary(res2)$R2

# illustrate the computation of R^2
100 * (res1$tau2 - res2$tau2) / res1$tau2

# fit a mixed-effects meta-regression model with number of sessions
# as a moderator (use the Knapp & Hartung method for inferences)
res3 <- rma(yi, vi, mods = ~ sessions, data=dat, test="knha")
res3

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

### non-parametric bootstrapping to obtain a CI for R^2

# load the boot package
library(boot)

# function to draw bootstrap replicates, fit the model, and return R^2
boot.func.np <- function(dat, indices, formula) {
   res <- suppressWarnings(rma(yi, vi, mods = formula, data=dat[indices,], method=c("REML","DL")))
   return(res$R2)
}

# run the bootstrapping (10000 replicates)
set.seed(1234)
res.boot <- boot(dat, boot.func.np, R=10000, formula=formula(res2))

# obtain the percentile and BCA CIs
boot.ci(res.boot, type=c("perc","bca"))

# histogram of the bootstrap distribution of R^2
hist(c(res.boot$t), breaks=100*seq(0,1,by=0.025), main="",
     xlab=expression(R^2 * " Value (in %)"), ylim=c(0,1000))

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