i2_ml() fails when `boot` arg is not NULL if optional arg 'mods' not specified during rma model fitting

When attempting to calculate $I^2$ with confidence intervals, i.e. with a non-null boot arg in call to i2_ml(), and when model was fitted without specifying the mods argument (which is optional according to rma documentation), i2_ml() fails.

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.

^{Created on 2023-04-17 with reprex v2.0.2}

Below is a suggested re-write that allows i2_ml() to be applied when metafor::formula.rma(model, type = "mods") evaluates to NULL. I've applied it the reprex below as proof-of-concept:

Local .Rprofile detected at /Users/egould/Documents/code/ManyAnalysts/.Rprofile

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

^{Created on 2023-04-17 with reprex v2.0.2}

daniel1noble / orchaRd

i2_ml() fails when `boot` arg is not NULL if optional arg 'mods' not specified during rma model fitting #22