Here's a minimal example using the (not-run) English example from the documentation for i2_ml()
:
library(metafor)
#> Loading required package: Matrix
#> Loading required package: metadat
#> Loading required package: numDeriv
#>
#> Loading the 'metafor' package (version 4.0-0). For an
#> introduction to the package please type: help(metafor)
library(orchaRd)
#>
#> Loading the 'orchaRd' package (version 2.0). For an
#> introduction and vignette to the package please see: https://daniel1noble.github.io/orchaRd/
# English example
data(english)
english <- escalc(measure = "SMD", n1i = NStartControl,
sd1i = SD_C, m1i = MeanC, n2i = NStartExpt, sd2i = SD_E,
m2i = MeanE, var.names=c("SMD","vSMD"),data = english)
english_MA <- rma.mv(yi = SMD, V = vSMD,
random = list( ~ 1 | StudyNo, ~ 1 | EffectID), data = english)
I2_eng_1 <- i2_ml(english_MA, data = english, boot = 10)
#> Error: metafor::rma.mv(yi = ysim, V = vi, mods = mods_formula, random = random_formula, :
#> The object/variable ('mods_formula') specified for the 'mods' argument is NULL.
library(metafor)
#> Loading required package: Matrix
#> Loading required package: metadat
#> Loading required package: numDeriv
#>
#> Loading the 'metafor' package (version 4.0-0). For an
#> introduction to the package please type: help(metafor)
library(orchaRd)
#>
#> Loading the 'orchaRd' package (version 2.0). For an
#> introduction and vignette to the package please see: https://daniel1noble.github.io/orchaRd/
# English example
data(english)
english <- escalc(measure = "SMD", n1i = NStartControl,
sd1i = SD_C, m1i = MeanC, n2i = NStartExpt, sd2i = SD_E,
m2i = MeanE, var.names=c("SMD","vSMD"),data = english)
english_MA <- rma.mv(yi = SMD, V = vSMD,
random = list( ~ 1 | StudyNo, ~ 1 | EffectID), data = english)
I2_eng_1 <- i2_ml(english_MA, data = english, boot = 10)
#> Error: metafor::rma.mv(yi = ysim, V = vi, mods = mods_formula, random = random_formula, :
#> The object/variable ('mods_formula') specified for the 'mods' argument is NULL.
# Define new fun:
i2_ml <- function (model, method = c("ratio", "matrix"), data, boot = NULL) {
if (all(class(model) %in% c("robust.rma", "rma.mv", "rma",
"rma.uni")) == FALSE) {
stop("Sorry, you need to fit a metafor model of class robust.rma, rma.mv, rma, rma.uni")
}
if (any(model$tau2 > 0)) {
stop("Sorry. At the moment i2_ml cannot take models with heterogeneous variance.")
}
method <- match.arg(method)
if (method == "matrix") {
I2s <- matrix_i2(model)
}
else {
I2s <- ratio_i2(model)
}
if (!is.null(boot)) {
sim <- metafor::simulate.rma(model, nsim = boot)
mods_formula <- metafor::formula.rma(model, type = "mods")
random_formula <- model$random
vi <- model$vi
pb <- progress::progress_bar$new(total = boot, format = "Bootstrapping [:bar] :percent ETA: :eta",
show_after = 0)
if (is.null(mods_formula)) {
I2_each <- sapply(sim, function(ysim) {
tmp <- tryCatch(metafor::rma.mv(ysim, vi,
random = random_formula, data = data))
pb$tick()
Sys.sleep(1/boot)
if (method == "matrix") {
I2 <- matrix_i2(tmp)
}
else {
I2 <- ratio_i2(tmp)
}
return(I2)
})
} else{
I2_each <- sapply(sim, function(ysim) {
tmp <- tryCatch(metafor::rma.mv(ysim, vi, mods = mods_formula,
random = random_formula, data = data))
pb$tick()
Sys.sleep(1/boot)
if (method == "matrix") {
I2 <- matrix_i2(tmp)
}
else {
I2 <- ratio_i2(tmp)
}
return(I2)
})
}
I2s_each_95 <- data.frame(t(apply(I2_each, 1, stats::quantile,
probs = c(0.5, 0.025, 0.975))))
I2s <- round(I2s_each_95, digits = 3)
colnames(I2s) = c("Est.", "2.5%", "97.5%")
}
return(I2s)
}
# Apply updated fun:
I2_eng_1 <- i2_ml(english_MA, data = english, boot = 10)
I2_eng_1
#> Est. 2.5% 97.5%
#> I2_Total 74.643 69.327 81.686
#> I2_StudyNo 31.608 13.877 51.806
#> I2_EffectID 40.196 19.659 57.785
# Now check that function works when 'mods' is specified, using fish example from
# i2_ml() documentation:
## Fish example
data(fish)
warm_dat <- fish
model <- metafor::rma.mv(yi = lnrr, V = lnrr_vi,
random = list(~1 | group_ID, ~1 | es_ID),
mods = ~ experimental_design + trait.type + deg_dif + treat_end_days,
method = "REML", test = "t", data = warm_dat,
control=list(optimizer="optim", optmethod="Nelder-Mead"))
I2_fish_1 <- i2_ml(model, data = warm_dat, boot = 10)
I2_fish_1
#> Est. 2.5% 97.5%
#> I2_Total 99.943 99.936 99.958
#> I2_group_ID 34.485 27.059 47.810
#> I2_es_ID 65.457 52.147 72.878