Skip to content

Instantly share code, notes, and snippets.

@wviechtb
Last active February 10, 2023 14:06
Show Gist options
  • Save wviechtb/6fbfca40483cb9744384ab4572639169 to your computer and use it in GitHub Desktop.
Save wviechtb/6fbfca40483cb9744384ab4572639169 to your computer and use it in GitHub Desktop.
Bootstrapping rma.mv() model to get CI for pseudo R^2 statistic
############################################################################
library(metafor)
library(boot)
dat <- dat.crede2010
dat <- escalc(measure="ZCOR", ri=ri, ni=ni, data=dat, subset=criterion=="grade")
############################################################################
# fit multilevel model
res0 <- rma.mv(yi, vi, random = ~ 1 | studyid/sampleid, data=dat)
res0
# fit multilevel meta-regression model
res1 <- rma.mv(yi, vi, mods = ~ year, random = ~ 1 | studyid/sampleid, data=dat)
res1
# compute pseudo R^2 (proportional reduction in the sum of the two variance components)
max(0, 100 * (sum(res0$sigma2) - sum(res1$sigma2)) / sum(res0$sigma2))
############################################################################
boot.func <- function(dat, indices) {
res0 <- try(suppressWarnings(rma.mv(yi, vi, random = ~ 1 | studyid/sampleid, data=dat[indices,])), silent=TRUE)
res1 <- try(suppressWarnings(rma.mv(yi, vi, mods = ~ year, random = ~ 1 | studyid/sampleid, data=dat[indices,])), silent=TRUE)
if (inherits(res0, "try-error") || inherits(res1, "try-error")) {
NA
} else {
if (sum(res0$sigma2) > 0) {
max(0, 100 * (sum(res0$sigma2) - sum(res1$sigma2)) / sum(res0$sigma2))
} else {
# catch cases where sum(res0$sigma2) is 0, to avoid -Inf (if sum(res1$sigma2)
# is positive) or NaN (if sum(res1$sigma2) is also 0) for R^2; for such cases,
# define R^2 to be 0
0
}
}
}
set.seed(1234)
# this takes around 210 seconds (on an Intel Xeon CPU E5-2630 v4 @ 2.20GHz)
system.time(res.boot <- boot(dat, boot.func, R=1000))
boot.ci(res.boot, type=c("norm", "basic", "perc", "bca"))
library(parallel)
clust <- makeCluster(20, type="PSOCK")
clusterEvalQ(clust, library(metafor))
# if other objects are needed to fit the model, export them (not needed here)
#clusterExport(clust, c("obj1","obj2"))
set.seed(1234)
# this takes around 14 seconds (workstation has 20 cores)
system.time(res.boot <- boot(dat, boot.func, R=1000, parallel="snow", cl=clust, ncpus=length(clust)))
boot.ci(res.boot, type=c("norm", "basic", "perc", "bca"))
stopCluster(clust)
set.seed(1234)
# this also takes around 14 seconds (multicore only available on Unix-alike platforms)
system.time(res.boot <- boot(dat, boot.func, R=1000, parallel="multicore", ncpus=20))
boot.ci(res.boot, type=c("norm", "basic", "perc", "bca"))
############################################################################
# note: if the lower bound of a CI is negative, we can treat it as if it was 0;
# similarly, if the upper bound is above 100, we can treat it as if it was 100
############################################################################
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment