Bootstrapping rma.mv() model to get CI for pseudo R^2 statistic

boot_rma_mv_parallel.r

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

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

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