Running 'testthat.R' [10s]
Running the tests in 'tests/testthat.R' failed.
Complete output:
  > # This file is part of the standard setup for testthat.
  > # It is recommended that you do not modify it.
  > #
  > # Where should you do additional test configuration?
  > # Learn more about the roles of various files in:
  > # * https://r-pkgs.org/tests.html
  > # * https://testthat.r-lib.org/reference/test_package.html#special-files
  > 
  > library(testthat)
  > library(confMeta)
  > 
  > test_check("confMeta")
  Saving _problems/test-confMeta.weights-19.R
  [ FAIL 1 | WARN 0 | SKIP 0 | PASS 79 ]
  
  ══ Failed tests ════════════════════════════════════════════════════════════════
  ── Error ('test-confMeta.weights.R:13:5'): weights argument is validated correctly ──
  Error in `(function (yi, vi, sei, weights, ai, bi, ci, di, n1i, n2i, x1i, x2i, t1i, t2i, m1i, m2i, sd1i, sd2i, xi, mi, ri, ti, fi, pi, sdi, r2i, ni, mods, scale, measure = "GEN", data, slab, subset, add = 1/2, to = "only0", drop00 = FALSE, intercept = TRUE, method = "REML", weighted = TRUE, test = "z", level = 95, btt, att, tau2, verbose = FALSE, digits, control, ...)  {     mstyle <- .get.mstyle()     if (!is.element(measure, c("RR", "OR", "PETO", "RD", "AS", "PHI", "ZPHI", "YUQ", "YUY", "RTET", "ZTET", "PBIT", "OR2D", "OR2DN", "OR2DL", "MPRD", "MPRR", "MPOR", "MPORC", "MPPETO", "MPORM", "IRR", "IRD", "IRSD", "MD", "SMD", "SMDH", "SMD1", "SMD1H", "ROM", "VR", "CVR", "RPB", "ZPB", "RBIS", "ZBIS", "D2OR", "D2ORN", "D2ORL", "COR", "UCOR", "ZCOR", "PCOR", "ZPCOR", "SPCOR", "ZSPCOR", "R2", "ZR2", "R2F", "ZR2F", "PR", "PLN", "PLO", "PRZ", "PAS", "PFT", "IR", "IRLN", "IRS", "IRFT", "MN", "SMN", "MNLN",          "SDLN", "CVLN", "MC", "SMCC", "SMCR", "SMCRH", "SMCRP", "SMCRPH", "CLESCN", "AUCCN", "ROMC", "VRC", "CVRC", "ARAW", "AHW", "ABT", "REH", "CLES", "CLESN", "AUC", "AUCN", "HR", "HD", "GEN")))          stop(mstyle$stop("Unknown 'measure' specified."))     if (!is.element(method[1], c("FE", "EE", "CE", "HS", "HSk", "HE", "CO", "VC", "DL", "DLIT", "GENQ", "GENQM", "SJ", "SJIT", "PM", "MP", "PMM", "ML", "REML", "EB")))          stop(mstyle$stop("Unknown 'method' specified."))     if (length(add) > 1L)          add <- add[1]     if (length(to) > 1L)          to <- to[1]     na.act <- getOption("na.action")     if (!is.element(na.act, c("na.omit", "na.exclude", "na.fail", "na.pass")))          stop(mstyle$stop("Unknown 'na.action' specified under options()."))     if (missing(tau2)) {         tau2 <- NULL     }     else {         if (!is.null(tau2) && length(tau2) != 1L)              stop(mstyle$stop("Argument 'tau2' must be a scalar."))     }     if (missing(control))          control <- list()     time.start <- proc.time()     ddd <- list(...)     .chkdots(ddd, c("vtype", "knha", "onlyo1", "addyi", "addvi", "correct", "i2def", "r2def", "skipr2", "abbrev", "dfs", "time", "outlist", "link", "optbeta", "alpha", "beta", "skiphes", "retopt", "randhet", "omega2", "pleasedonotreportI2thankyouverymuch"))     if (is.null(ddd$vtype)) {         vtype <- "LS"     }     else {         vtype <- ddd$vtype     }     if (isFALSE(ddd$knha))          test <- "z"     if (isTRUE(ddd$knha))          test <- "knha"     test <- tolower(test)     if (!is.element(test, c("z", "t", "knha", "hksj", "adhoc")))          stop(mstyle$stop("Unknown option specified for the 'test' argument."))     if (test == "hksj")          test <- "knha"     if (missing(scale)) {         model <- "rma.uni"     }     else {         model <- "rma.ls"     }     onlyo1 <- .chkddd(ddd$onlyo1, FALSE)     addyi <- .chkddd(ddd$addyi, TRUE)     addvi <- .chkddd(ddd$addvi, TRUE)     correct <- .chkddd(ddd$correct, TRUE)     i2def <- .chkddd(ddd$i2def, "1")     r2def <- .chkddd(ddd$r2def, "1")     link <- .chkddd(ddd$link, "log", match.arg(ddd$link, c("log", "identity")))     optbeta <- .chkddd(ddd$optbeta, FALSE, isTRUE(ddd$optbeta))     if (model == "rma.uni" && optbeta) {         warning(mstyle$warning("Argument 'optbeta' only relevant for location-scale models and hence ignored."), call. = FALSE)         optbeta <- FALSE     }     if (optbeta && !weighted)          stop(mstyle$stop("Must use 'weighted=TRUE' when 'optbeta=TRUE'."))     alpha <- .chkddd(ddd$alpha, NA_real_)     beta <- .chkddd(ddd$beta, NA_real_)     omega2 <- .chkddd(ddd$omega2, NA_real_)     if (model == "rma.uni" && !missing(att))          warning(mstyle$warning("Argument 'att' only relevant for location-scale models and hence ignored."), call. = FALSE)     if (missing(digits)) {         digits <- .set.digits(dmiss = TRUE)     }     else {         digits <- .set.digits(digits, dmiss = FALSE)     }     formula.yi <- NULL     formula.mods <- NULL     formula.scale <- NULL     if (verbose > 2) {         opwarn <- options(warn = 1)         on.exit(options(warn = opwarn$warn), add = TRUE)     }     randhet <- isTRUE(ddd$randhet)     if (model == "rma.ls") {         if (randhet && link == "identity")              stop(mstyle$stop("Cannot use identity link when allowing for residual heteroscedasticity."))         if (randhet) {             method <- "ML"             optbeta <- TRUE         }     }     else {         if (randhet)              warning(mstyle$warning("Argument 'randhet' only relevant for location-scale models and hence ignored."), call. = FALSE)     }     if (verbose)          .space()     if (verbose > 1)          message(mstyle$message("Extracting/computing the yi/vi values ..."))     if (missing(data))          data <- NULL     if (is.null(data)) {         data <- sys.frame(sys.parent())     }     else {         if (!is.data.frame(data))              data <- data.frame(data)     }     mf <- match.call()     addval <- mf[[match("add", names(mf))]]     if (is.element(measure, c("AS", "PHI", "ZPHI", "RTET", "ZTET", "IRSD", "PAS", "PFT", "IRS", "IRFT")) && is.null(addval))          add <- 0     yi <- .getx("yi", mf = mf, data = data)     if (!is.null(yi) && inherits(yi, "escalc"))          data <- yi     weights <- .getx("weights", mf = mf, data = data, checknumeric = TRUE)     slab <- .getx("slab", mf = mf, data = data)     subset <- .getx("subset", mf = mf, data = data)     mods <- .getx("mods", mf = mf, data = data)     scale <- .getx("scale", mf = mf, data = data)     ai <- bi <- ci <- di <- x1i <- x2i <- t1i <- t2i <- NA_real_     if (!is.null(yi)) {         if (inherits(yi, "formula")) {             formula.yi <- yi             formula.mods <- formula.yi[-2]             options(na.action = "na.pass")             mods <- model.matrix(yi, data = data)             attr(mods, "assign") <- NULL             attr(mods, "contrasts") <- NULL             yi <- model.response(model.frame(yi, data = data))             options(na.action = na.act)             names(yi) <- NULL             intercept <- FALSE         }         if (inherits(yi, "escalc")) {             if (!is.null(attr(yi, "yi.names"))) {                 yi.name <- attr(yi, "yi.names")[1]             }             else {                 if (!is.element("yi", names(yi)))                    stop(mstyle$stop("Cannot determine the name of the 'yi' variable."))                 yi.name <- "yi"             }             if (!is.null(attr(yi, "vi.names"))) {                 vi.name <- attr(yi, "vi.names")[1]             }             else {                 if (!is.element("vi", names(yi)))                    stop(mstyle$stop("Cannot determine the name of the 'vi' variable."))                 vi.name <- "vi"             }             vi <- yi[[vi.name]]             yi <- yi[[yi.name]]             if (is.null(yi))                  stop(mstyle$stop(paste0("Cannot find variable '", yi.name, "' in the object.")))             if (is.null(vi))                  stop(mstyle$stop(paste0("Cannot find variable '", vi.name, "' in the object.")))             yi.escalc <- TRUE         }         else {             yi.escalc <- FALSE         }         if (is.data.frame(yi)) {             if (ncol(yi) == 1L) {                 yi <- yi[[1]]             }             else {                 stop(mstyle$stop("The object/variable specified for the 'yi' argument is a data frame with multiple columns."))             }         }         if (.is.matrix(yi)) {             if (nrow(yi) == 1L || ncol(yi) == 1L) {                 yi <- as.vector(yi)             }             else {                 stop(mstyle$stop("The object/variable specified for the 'yi' argument is a matrix with multiple rows/columns."))             }         }         if (inherits(yi, "array"))              stop(mstyle$stop("The object/variable specified for the 'yi' argument is an array."))         if (!is.numeric(yi))              stop(mstyle$stop("The object/variable specified for the 'yi' argument is not numeric."))         k <- length(yi)         k.all <- k         if (measure == "GEN" && !is.null(attr(yi, "measure")))              measure <- attr(yi, "measure")         attr(yi, "measure") <- measure         if (!yi.escalc) {             vi <- .getx("vi", mf = mf, data = data, checknumeric = TRUE)             sei <- .getx("sei", mf = mf, data = data, checknumeric = TRUE)         }         ni <- .getx("ni", mf = mf, data = data, checknumeric = TRUE)         if (is.null(vi)) {             if (is.null(sei)) {                 stop(mstyle$stop("Must specify the 'vi' or 'sei' argument."))             }             else {                 vi <- sei^2             }         }         .chkviarg(mf$vi)         if (.is.matrix(vi)) {             if (nrow(vi) == 1L || ncol(vi) == 1L) {                 vi <- as.vector(vi)             }             else {                 if (.is.square(vi) && isSymmetric(unname(vi))) {                   vi <- as.matrix(vi)                   if (any(vi[!diag(nrow(vi))] != 0))                      warning(mstyle$warning("Using only the diagonal elements from the 'vi' argument as the sampling variances."), call. = FALSE)                   vi <- diag(vi)                 }                 else {                   stop(mstyle$stop("The object/variable specified for the 'vi' argument is a matrix with multiple rows/columns."))                 }             }         }         if (inherits(vi, "array"))              stop(mstyle$stop("The object/variable specified for the 'vi' argument is an array."))         if ((length(vi) == 1L && vi == 0) || (length(vi) == k && all(vi == 0, na.rm = TRUE))) {             vi0 <- TRUE         }         else {             vi0 <- FALSE         }         vi <- .expand1(vi, k)         if (length(vi) != k)              stop(mstyle$stop("Length of 'yi' and 'vi' (or 'sei') are not the same."))         if (is.null(ni))              ni <- attr(yi, "ni")         if (!is.null(ni) && length(ni) != k)              ni <- NULL         if (!is.null(ni))              attr(yi, "ni") <- ni         if (is.null(slab)) {             slab <- attr(yi, "slab")             if (!is.null(slab) && length(slab) != k)                  slab <- NULL         }         if (!is.null(subset)) {             subset <- .chksubset(subset, k)             yi <- .getsubset(yi, subset)             vi <- .getsubset(vi, subset)             ni <- .getsubset(ni, subset)             attr(yi, "measure") <- measure             attr(yi, "ni") <- ni         }     }     else {         if (is.element(measure, c("RR", "OR", "PETO", "RD", "AS", "PHI", "ZPHI", "YUQ", "YUY", "RTET", "ZTET", "PBIT", "OR2D", "OR2DN", "OR2DL", "MPRD", "MPRR", "MPOR", "MPORC", "MPPETO", "MPORM"))) {             ai <- .getx("ai", mf = mf, data = data, checknumeric = TRUE)             bi <- .getx("bi", mf = mf, data = data, checknumeric = TRUE)             ci <- .getx("ci", mf = mf, data = data, checknumeric = TRUE)             di <- .getx("di", mf = mf, data = data, checknumeric = TRUE)             n1i <- .getx("n1i", mf = mf, data = data, checknumeric = TRUE)             n2i <- .getx("n2i", mf = mf, data = data, checknumeric = TRUE)             ri <- .getx("ri", mf = mf, data = data, checknumeric = TRUE)             pi <- .getx("pi", mf = mf, data = data, checknumeric = TRUE)             if (is.null(bi))                  bi <- n1i - ai             if (is.null(di))                  di <- n2i - ci             k <- length(ai)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 ai <- .getsubset(ai, subset)                 bi <- .getsubset(bi, subset)                 ci <- .getsubset(ci, subset)                 di <- .getsubset(di, subset)                 ri <- .getsubset(ri, subset)                 pi <- .getsubset(pi, subset)             }             args <- list(ai = ai, bi = bi, ci = ci, di = di, ri = ri, pi = pi, add = add, to = to, drop00 = drop00, onlyo1 = onlyo1, addyi = addyi, addvi = addvi)         }         if (is.element(measure, c("IRR", "IRD", "IRSD"))) {             x1i <- .getx("x1i", mf = mf, data = data, checknumeric = TRUE)             x2i <- .getx("x2i", mf = mf, data = data, checknumeric = TRUE)             t1i <- .getx("t1i", mf = mf, data = data, checknumeric = TRUE)             t2i <- .getx("t2i", mf = mf, data = data, checknumeric = TRUE)             k <- length(x1i)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 x1i <- .getsubset(x1i, subset)                 x2i <- .getsubset(x2i, subset)                 t1i <- .getsubset(t1i, subset)                 t2i <- .getsubset(t2i, subset)             }             args <- list(x1i = x1i, x2i = x2i, t1i = t1i, t2i = t2i, add = add, to = to, drop00 = drop00, addyi = addyi, addvi = addvi)         }         if (is.element(measure, c("MD", "SMD", "SMDH", "SMD1", "SMD1H", "ROM", "RPB", "ZPB", "RBIS", "ZBIS", "D2OR", "D2ORN", "D2ORL", "VR", "CVR"))) {             m1i <- .getx("m1i", mf = mf, data = data, checknumeric = TRUE)             m2i <- .getx("m2i", mf = mf, data = data, checknumeric = TRUE)             sd1i <- .getx("sd1i", mf = mf, data = data, checknumeric = TRUE)             sd2i <- .getx("sd2i", mf = mf, data = data, checknumeric = TRUE)             n1i <- .getx("n1i", mf = mf, data = data, checknumeric = TRUE)             n2i <- .getx("n2i", mf = mf, data = data, checknumeric = TRUE)             di <- .getx("di", mf = mf, data = data, checknumeric = TRUE)             ti <- .getx("ti", mf = mf, data = data, checknumeric = TRUE)             pi <- .getx("pi", mf = mf, data = data, checknumeric = TRUE)             ri <- .getx("ri", mf = mf, data = data, checknumeric = TRUE)             if (is.element(measure, c("SMD", "RPB", "ZPB", "RBIS", "ZBIS", "D2OR", "D2ORN", "D2ORL"))) {                 if (!.equal.length(m1i, m2i, sd1i, sd2i, n1i, n2i, di, ti, pi, ri))                    stop(mstyle$stop("Supplied data vectors are not all of the same length."))                 ti <- replmiss(ti, .convp2t(pi, df = n1i + n2i - 2))                 di <- replmiss(di, ti * sqrt(1/n1i + 1/n2i))                 mi <- n1i + n2i - 2                 hi <- mi/n1i + mi/n2i                 di <- replmiss(di, sqrt(hi) * ri/sqrt(1 - ri^2))                 m1i[!is.na(di)] <- di[!is.na(di)]                 m2i[!is.na(di)] <- 0                 sd1i[!is.na(di)] <- 1                 sd2i[!is.na(di)] <- 1             }             k <- length(n1i)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 m1i <- .getsubset(m1i, subset)                 m2i <- .getsubset(m2i, subset)                 sd1i <- .getsubset(sd1i, subset)                 sd2i <- .getsubset(sd2i, subset)                 n1i <- .getsubset(n1i, subset)                 n2i <- .getsubset(n2i, subset)             }             args <- list(m1i = m1i, m2i = m2i, sd1i = sd1i, sd2i = sd2i, n1i = n1i, n2i = n2i)         }         if (is.element(measure, c("COR", "UCOR", "ZCOR"))) {             ri <- .getx("ri", mf = mf, data = data, checknumeric = TRUE)             ni <- .getx("ni", mf = mf, data = data, checknumeric = TRUE)             ti <- .getx("ti", mf = mf, data = data, checknumeric = TRUE)             pi <- .getx("pi", mf = mf, data = data, checknumeric = TRUE)             if (!.equal.length(ri, ni, ti, pi))                  stop(mstyle$stop("Supplied data vectors are not all of the same length."))             ti <- replmiss(ti, .convp2t(pi, df = ni - 2))             ri <- replmiss(ri, ti/sqrt(ti^2 + ni - 2))             k <- length(ri)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 ri <- .getsubset(ri, subset)                 ni <- .getsubset(ni, subset)             }             args <- list(ri = ri, ni = ni)         }         if (is.element(measure, c("PCOR", "ZPCOR", "SPCOR", "ZSPCOR"))) {             ri <- .getx("ri", mf = mf, data = data, checknumeric = TRUE)             ti <- .getx("ti", mf = mf, data = data, checknumeric = TRUE)             mi <- .getx("mi", mf = mf, data = data, checknumeric = TRUE)             ni <- .getx("ni", mf = mf, data = data, checknumeric = TRUE)             pi <- .getx("pi", mf = mf, data = data, checknumeric = TRUE)             r2i <- .getx("r2i", mf = mf, data = data, checknumeric = TRUE)             if (!.equal.length(ri, ti, mi, ni, pi, r2i))                  stop(mstyle$stop("Supplied data vectors are not all of the same length."))             ti <- replmiss(ti, .convp2t(pi, df = ni - mi - 1))             if (is.element(measure, c("PCOR", "ZPCOR")))                  ri <- replmiss(ri, ti/sqrt(ti^2 + ni - mi - 1))             if (is.element(measure, c("SPCOR", "ZSPCOR")))                  ri <- replmiss(ri, ti * sqrt(1 - r2i)/sqrt(ni - mi - 1))             k <- length(ri)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 ri <- .getsubset(ri, subset)                 mi <- .getsubset(mi, subset)                 ni <- .getsubset(ni, subset)                 r2i <- .getsubset(r2i, subset)             }             args <- list(ri = ri, mi = mi, ni = ni, r2i = r2i)         }         if (is.element(measure, c("R2", "ZR2", "R2F", "ZR2F"))) {             r2i <- .getx("r2i", mf = mf, data = data, checknumeric = TRUE)             mi <- .getx("mi", mf = mf, data = data, checknumeric = TRUE)             ni <- .getx("ni", mf = mf, data = data, checknumeric = TRUE)             fi <- .getx("fi", mf = mf, data = data, checknumeric = TRUE)             pi <- .getx("pi", mf = mf, data = data, checknumeric = TRUE)             if (!.equal.length(r2i, mi, ni))                  stop(mstyle$stop("Supplied data vectors are not all of the same length."))             fi <- replmiss(fi, .convp2f(pi, df1 = mi, df2 = ni - mi - 1))             r2i <- replmiss(r2i, mi * fi/(mi * fi + (ni - mi - 1)))             k <- length(r2i)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 r2i <- .getsubset(r2i, subset)                 mi <- .getsubset(mi, subset)                 ni <- .getsubset(ni, subset)             }             args <- list(r2i = r2i, mi = mi, ni = ni)         }         if (is.element(measure, c("PR", "PLN", "PLO", "PRZ", "PAS", "PFT"))) {             xi <- .getx("xi", mf = mf, data = data, checknumeric = TRUE)             mi <- .getx("mi", mf = mf, data = data, checknumeric = TRUE)             ni <- .getx("ni", mf = mf, data = data, checknumeric = TRUE)             if (is.null(mi))                  mi <- ni - xi             k <- length(xi)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 xi <- .getsubset(xi, subset)                 mi <- .getsubset(mi, subset)             }             args <- list(xi = xi, mi = mi, add = add, to = to, addyi = addyi, addvi = addvi)         }         if (is.element(measure, c("IR", "IRLN", "IRS", "IRFT"))) {             xi <- .getx("xi", mf = mf, data = data, checknumeric = TRUE)             ti <- .getx("ti", mf = mf, data = data, checknumeric = TRUE)             k <- length(xi)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 xi <- .getsubset(xi, subset)                 ti <- .getsubset(ti, subset)             }             args <- list(xi = xi, ti = ti, add = add, to = to, addyi = addyi, addvi = addvi)         }         if (is.element(measure, c("MN", "SMN", "MNLN", "SDLN", "CVLN"))) {             mi <- .getx("mi", mf = mf, data = data, checknumeric = TRUE)             sdi <- .getx("sdi", mf = mf, data = data, checknumeric = TRUE)             ni <- .getx("ni", mf = mf, data = data, checknumeric = TRUE)             k <- length(ni)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 mi <- .getsubset(mi, subset)                 sdi <- .getsubset(sdi, subset)                 ni <- .getsubset(ni, subset)             }             args <- list(mi = mi, sdi = sdi, ni = ni)         }         if (is.element(measure, c("MC", "SMCC", "SMCR", "SMCRH", "SMCRP", "SMCRPH", "CLESCN", "AUCCN", "ROMC", "VRC", "CVRC"))) {             m1i <- .getx("m1i", mf = mf, data = data, checknumeric = TRUE)             m2i <- .getx("m2i", mf = mf, data = data, checknumeric = TRUE)             sd1i <- .getx("sd1i", mf = mf, data = data, checknumeric = TRUE)             sd2i <- .getx("sd2i", mf = mf, data = data, checknumeric = TRUE)             ri <- .getx("ri", mf = mf, data = data, checknumeric = TRUE)             ni <- .getx("ni", mf = mf, data = data, checknumeric = TRUE)             di <- .getx("di", mf = mf, data = data, checknumeric = TRUE)             ti <- .getx("ti", mf = mf, data = data, checknumeric = TRUE)             pi <- .getx("pi", mf = mf, data = data, checknumeric = TRUE)             ri <- .expand1(ri, list(m1i, m2i, sd1i, sd2i, ni, di, ti, pi))             if (measure == "SMCC") {                 if (!.equal.length(m1i, m2i, sd1i, sd2i, ri, ni, di, ti, pi))                    stop(mstyle$stop("Supplied data vectors are not all of the same length."))                 ti <- replmiss(ti, .convp2t(pi, df = ni - 1))                 di <- replmiss(di, ti * sqrt(1/ni))                 m1i[!is.na(di)] <- di[!is.na(di)]                 m2i[!is.na(di)] <- 0                 sd1i[!is.na(di)] <- 1                 sd2i[!is.na(di)] <- 1                 ri[!is.na(di)] <- 0.5             }             k <- length(m1i)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 m1i <- .getsubset(m1i, subset)                 m2i <- .getsubset(m2i, subset)                 sd1i <- .getsubset(sd1i, subset)                 sd2i <- .getsubset(sd2i, subset)                 ni <- .getsubset(ni, subset)                 ri <- .getsubset(ri, subset)             }             args <- list(m1i = m1i, m2i = m2i, sd1i = sd1i, sd2i = sd2i, ri = ri, ni = ni)         }         if (is.element(measure, c("ARAW", "AHW", "ABT"))) {             ai <- .getx("ai", mf = mf, data = data, checknumeric = TRUE)             mi <- .getx("mi", mf = mf, data = data, checknumeric = TRUE)             ni <- .getx("ni", mf = mf, data = data, checknumeric = TRUE)             k <- length(ai)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 ai <- .getsubset(ai, subset)                 mi <- .getsubset(mi, subset)                 ni <- .getsubset(ni, subset)             }             args <- list(ai = ai, mi = mi, ni = ni)         }         if (measure == "REH") {             ai <- .getx("ai", mf = mf, data = data, checknumeric = TRUE)             bi <- .getx("bi", mf = mf, data = data, checknumeric = TRUE)             ci <- .getx("ci", mf = mf, data = data, checknumeric = TRUE)             k <- length(ai)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 ai <- .getsubset(ai, subset)                 bi <- .getsubset(bi, subset)                 ci <- .getsubset(ci, subset)             }             args <- list(ai = ai, bi = bi, ci = ci)         }         if (is.element(measure, c("CLES", "AUC"))) {             ai <- .getx("ai", mf = mf, data = data, checknumeric = TRUE)             n1i <- .getx("n1i", mf = mf, data = data, checknumeric = TRUE)             n2i <- .getx("n2i", mf = mf, data = data, checknumeric = TRUE)             mi <- .getx("mi", mf = mf, data = data, checknumeric = TRUE)             if (is.null(mi))                  mi <- rep(0, length(ai))             mi[is.na(mi)] <- 0             k <- length(ai)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 ai <- .getsubset(ai, subset)                 n1i <- .getsubset(n1i, subset)                 n2i <- .getsubset(n2i, subset)                 mi <- .getsubset(mi, subset)             }             args <- list(ai = ai, n1i = n1i, n2i = n2i, mi = mi)         }         if (is.element(measure, c("CLESN", "AUCN"))) {             m1i <- .getx("m1i", mf = mf, data = data, checknumeric = TRUE)             m2i <- .getx("m2i", mf = mf, data = data, checknumeric = TRUE)             sd1i <- .getx("sd1i", mf = mf, data = data, checknumeric = TRUE)             sd2i <- .getx("sd2i", mf = mf, data = data, checknumeric = TRUE)             n1i <- .getx("n1i", mf = mf, data = data, checknumeric = TRUE)             n2i <- .getx("n2i", mf = mf, data = data, checknumeric = TRUE)             di <- .getx("di", mf = mf, data = data, checknumeric = TRUE)             ti <- .getx("ti", mf = mf, data = data, checknumeric = TRUE)             pi <- .getx("pi", mf = mf, data = data, checknumeric = TRUE)             ai <- .getx("ai", mf = mf, data = data, checknumeric = TRUE)             if (!.equal.length(m1i, m2i, sd1i, sd2i, n1i, n2i, di, ti, pi, ai))                  stop(mstyle$stop("Supplied data vectors are not all of the same length."))             if (!.all.specified(n1i, n2i))                  stop(mstyle$stop("Cannot compute outcomes. Check that all of the required information is specified\n  via the appropriate arguments."))             k.all <- .maxlength(m1i, m2i, sd1i, sd2i, n1i, n2i, di, ti, pi, ai)             vtype <- .expand1(vtype, k.all)             if (is.null(sd1i) || is.null(sd2i)) {                 sd1i <- .expand1(NA_real_, k.all)                 sd2i <- .expand1(NA_real_, k.all)             }             ti <- replmiss(ti, .convp2t(pi, df = n1i + n2i - 2))             di <- replmiss(di, ti * sqrt(1/n1i + 1/n2i))             if (!is.null(di))                  vtype[!is.na(di)] <- "HO"             sdpi <- ifelse(vtype == "HO", sqrt(((n1i - 1) * sd1i^2 + (n2i - 1) * sd2i^2)/(n1i + n2i - 2)), sqrt((sd1i^2 + sd2i^2)/2))             di <- replmiss(di, (m1i - m2i)/sdpi)             ai <- replmiss(ai, pnorm(di/sqrt(2)))             di <- replmiss(di, qnorm(ai) * sqrt(2))             k.all <- length(ai)             sdsmiss <- is.na(sd1i) | is.na(sd2i)             sd1i <- ifelse(sdsmiss, 1, sd1i)             sd2i <- ifelse(sdsmiss, 1, sd2i)             vtype[sdsmiss] <- "HO"             k <- length(ai)             k.all <- k             if (!is.null(subset)) {                 subset <- .chksubset(subset, k)                 vtype <- .getsubset(vtype, subset)                 ai <- .getsubset(ai, subset)                 sd1i <- .getsubset(sd1i, subset)                 sd2i <- .getsubset(sd2i, subset)                 n1i <- .getsubset(n1i, subset)                 n2i <- .getsubset(n2i, subset)             }             args <- list(ai = ai, sd1i = sd1i, sd2i = sd2i, n1i = n1i, n2i = n2i)         }         args <- c(args, list(measure = measure, vtype = vtype, correct = correct))         dat <- .do.call(escalc, args)         if (is.element(measure, "GEN"))              stop(mstyle$stop("Specify the desired outcome measure via the 'measure' argument."))         yi <- dat$yi         vi <- dat$vi         ni <- attr(yi, "ni")     }     if (!is.null(weights)) {         if (verbose > 1)              message(mstyle$message("Extracting custom weights ..."))         if (optbeta)              stop(mstyle$stop("Cannot use custom weights when 'optbeta=TRUE'."))         weights <- .expand1(weights, k)         if (length(weights) != k)              stop(mstyle$stop("Length of 'yi' and 'weights' are not the same."))         if (!is.null(subset))              weights <- .getsubset(weights, subset)         if (any(weights < 0, na.rm = TRUE))              stop(mstyle$stop("Negative weights are not allowed."))         if (any(is.infinite(weights)))              stop(mstyle$stop("Infinite weights are not allowed."))     }     if (verbose > 1)          message(mstyle$message("Creating the model matrix ..."))     if (inherits(mods, "formula")) {         formula.mods <- mods         if (.is.tilde1(formula.mods)) {             mods <- matrix(1, nrow = k, ncol = 1)             intercept <- FALSE         }         else {             options(na.action = "na.pass")             mods <- model.matrix(mods, data = data)             attr(mods, "assign") <- NULL             attr(mods, "contrasts") <- NULL             options(na.action = na.act)             intercept <- FALSE         }     }     if (.is.vector(mods))          mods <- cbind(mods)     if (is.data.frame(mods))          mods <- as.matrix(mods)     if (is.character(mods))          stop(mstyle$stop("The model matrix contains character variables."))     if (!is.null(mods) && nrow(mods) != k)          stop(mstyle$stop(paste0("Number of rows in the model matrix (", nrow(mods), ") does not match the length of the outcome vector (", k, ").")))     if (model == "rma.ls") {         if (inherits(scale, "formula")) {             formula.scale <- scale             if (.is.tilde1(formula.scale)) {                 Z <- matrix(1, nrow = k, ncol = 1)                 colnames(Z) <- "intrcpt"             }             else {                 options(na.action = "na.pass")                 Z <- model.matrix(scale, data = data)                 colnames(Z)[grep("(Intercept)", colnames(Z), fixed = TRUE)] <- "intrcpt"                 attr(Z, "assign") <- NULL                 attr(Z, "contrasts") <- NULL                 options(na.action = na.act)             }         }         else {             Z <- scale             if (.is.vector(Z))                  Z <- cbind(Z)             if (is.data.frame(Z))                  Z <- as.matrix(Z)             if (is.character(Z))                  stop(mstyle$stop("Scale model matrix contains character variables."))         }         if (nrow(Z) != k)              stop(mstyle$stop(paste0("Number of rows in the model matrix specified via the 'scale' argument (", nrow(Z), ") does not match the length of the outcome vector (", k, ").")))     }     else {         Z <- NULL     }     if (!is.null(subset)) {         if (verbose > 1)              message(mstyle$message("Subsetting ..."))         mods <- .getsubset(mods, subset)         Z <- .getsubset(Z, subset)     }     if (verbose > 1)          message(mstyle$message("Generating/extracting the study labels ..."))     ids <- seq_len(k)     if (is.null(slab)) {         slab.null <- TRUE         slab <- ids     }     else {         if (anyNA(slab))              stop(mstyle$stop("NAs in study labels."))         if (length(slab) != k)              stop(mstyle$stop(paste0("Length of the 'slab' argument (", length(slab), ") does not correspond to the size of the dataset (", k, ").")))         slab.null <- FALSE     }     if (!is.null(subset)) {         if (verbose > 1)              message(mstyle$message("Subsetting ..."))         slab <- .getsubset(slab, subset)         ids <- .getsubset(ids, subset)     }     if (anyDuplicated(slab))          slab <- .make.unique(slab)     attr(yi, "slab") <- slab     k <- length(yi)     outdat.f <- list(ai = ai, bi = bi, ci = ci, di = di, x1i = x1i, x2i = x2i, t1i = t1i, t2i = t2i)     yi.f <- yi     vi.f <- vi     weights.f <- weights     ni.f <- ni     mods.f <- mods     Z.f <- Z     k.f <- k     has.na <- is.na(yi) | is.na(vi) | (if (is.null(mods))          FALSE     else apply(is.na(mods), 1, any)) | (if (is.null(Z))          FALSE     else apply(is.na(Z), 1, any)) | (if (is.null(weights))          FALSE     else is.na(weights))     not.na <- !has.na     if (any(has.na)) {         if (verbose > 1)              message(mstyle$message("Handling NAs ..."))         if (na.act == "na.omit" || na.act == "na.exclude" || na.act == "na.pass") {             yi <- yi[not.na]             vi <- vi[not.na]             weights <- weights[not.na]             ni <- ni[not.na]             mods <- mods[not.na, , drop = FALSE]             Z <- Z[not.na, , drop = FALSE]             k <- length(yi)             warning(mstyle$warning(paste(sum(has.na), ifelse(sum(has.na) > 1, "studies", "study"), "with NAs omitted from model fitting.")), call. = FALSE)             attr(yi, "measure") <- measure             attr(yi, "ni") <- ni         }         if (na.act == "na.fail")              stop(mstyle$stop("Missing values in data."))     }     if (k < 1L)          stop(mstyle$stop("Processing terminated since k = 0."))     if (any(vi <= 0)) {         allvipos <- FALSE         if (!vi0)              warning(mstyle$warning("There are outcomes with non-positive sampling variances."), call. = FALSE)         vi.neg <- vi < 0         if (any(vi.neg)) {             vi[vi.neg] <- 0             warning(mstyle$warning("Negative sampling variances constrained to 0."), call. = FALSE)         }     }     else {         allvipos <- TRUE     }     if (k == 1L && test != "z") {         warning(mstyle$warning("Setting argument test=\"z\" since k=1."), call. = FALSE)         test <- "z"     }     if (is.null(mods) && !intercept) {         warning(mstyle$warning("Must either include an intercept and/or moderators in the model.\nCoerced an intercept into the model."), call. = FALSE)         intercept <- TRUE     }     if (!is.null(mods) && ncol(mods) == 0L) {         warning(mstyle$warning("Cannot fit model with an empty model matrix. Coerced an intercept into the model."), call. = FALSE)         intercept <- TRUE     }     if (intercept) {         X <- cbind(intrcpt = rep(1, k), mods)         X.f <- cbind(intrcpt = rep(1, k.f), mods.f)     }     else {         X <- mods         X.f <- mods.f     }     tmp <- try(lm(yi ~ 0 + X), silent = TRUE)     if (inherits(tmp, "lm")) {         coef.na <- is.na(coef(tmp))     }     else {         coef.na <- rep(FALSE, NCOL(X))     }     if (any(coef.na)) {         warning(mstyle$warning("Redundant predictors dropped from the model."), call. = FALSE)         X <- X[, !coef.na, drop = FALSE]         X.f <- X.f[, !coef.na, drop = FALSE]     }     is.int <- apply(X, 2, .is.intercept)     if (any(is.int)) {         int.incl <- TRUE         int.indx <- which(is.int, arr.ind = TRUE)         X <- cbind(intrcpt = 1, X[, -int.indx, drop = FALSE])         X.f <- cbind(intrcpt = 1, X.f[, -int.indx, drop = FALSE])         intercept <- TRUE     }     else {         int.incl <- FALSE     }     p <- NCOL(X)     colnames(X) <- colnames(X.f) <- .make.unique(colnames(X))     colnames(Z) <- colnames(Z.f) <- .make.unique(colnames(Z))     if ((p == 1L) && .is.intercept(X)) {         int.only <- TRUE     }     else {         int.only <- FALSE     }     if (!(int.only && k == 1L)) {         if (is.element(method[1], c("FE", "EE", "CE"))) {             if (p > k)                  stop(mstyle$stop("Number of parameters to be estimated is larger than the number of observations."))         }         else {             if (!is.null(tau2) && !is.na(tau2)) {                 if (p > k)                    stop(mstyle$stop("Number of parameters to be estimated is larger than the number of observations."))             }             else {                 if ((p + 1) > k)                    stop(mstyle$stop("Number of parameters to be estimated is larger than the number of observations."))             }         }     }     btt <- .set.btt(btt, p, int.incl, colnames(X))     m <- length(btt)     con <- list(verbose = FALSE, evtol = 1e-07, REMLf = TRUE)     if (model == "rma.uni") {         con <- c(con, list(tau2.init = NULL, tau2.min = 0, tau2.max = 100, threshold = 10^-5, tol = .Machine$double.eps^0.25, ll0check = TRUE, maxiter = 100, stepadj = 1))         con$tau2.max <- max(con$tau2.max, 10 * mad(yi)^2)     }     if (model == "rma.ls") {         con <- c(con, list(beta.init = NULL, hesspack = "numDeriv", optimizer = "nlminb", optmethod = "BFGS", parallel = list(), cl = NULL, ncpus = 1L, tau2.min = 0, tau2.max = Inf, alpha.init = NULL, alpha.min = -Inf, alpha.max = Inf, omega2.init = NULL, hessianCtrl = NULL, htransf = FALSE, hes.beta.fix = FALSE, hes.alpha.fix = FALSE, hes.omega2.fix = FALSE, omega2tol = 1e-05, scaleZ = TRUE, mfmaxit = 10^5))     }     con.pos <- pmatch(names(control), names(con))     con[c(na.omit(con.pos))] <- control[!is.na(con.pos)]     if (verbose)          con$verbose <- verbose     verbose <- con$verbose     if (model == "rma.ls") {         con$hesspack <- match.arg(con$hesspack, c("numDeriv", "pracma", "calculus"))         if (!isTRUE(ddd$skiphes) && !requireNamespace(con$hesspack, quietly = TRUE))              stop(mstyle$stop(paste0("Please install the '", con$hesspack, "' package to compute the Hessian.")))         if (con$hesspack == "numDeriv") {             if (is.null(con$hessianCtrl$r))                  con$hessianCtrl$r <- 8         }         if (con$hesspack == "pracma") {             if (is.null(con$hessianCtrl$h))                  con$hessianCtrl$h <- .Machine$double.eps^(1/4)         }         if (con$hesspack == "calculus") {             if (is.null(con$hessianCtrl$accuracy))                  con$hessianCtrl$accuracy <- 4         }     }     if (model == "rma.uni") {         if (con$tau2.min < 0 && (con$tau2.min < -min(vi))) {             con$tau2.min <- -min(vi)             warning(mstyle$warning(paste0("Value of 'tau2.min' constrained to -min(vi) = ", fmtx(-min(vi), digits[["est"]]), ".")), call. = FALSE)         }     }     else {         if (is.element("tau2.min", names(control)))              con$tau2.min[con$tau2.min < 0] <- 0     }     if (!.chkpd(crossprod(X), tol = con$evtol))          stop(mstyle$stop("Model matrix not of full rank. Cannot fit model."))     vimaxmin <- max(vi)/min(vi)     if (is.finite(vimaxmin) && vimaxmin >= 1/con$evtol)          warning(mstyle$warning("Ratio of largest to smallest sampling variance extremely large. May not be able to obtain stable results."), call. = FALSE)     se.tau2 <- I2 <- H2 <- QE <- QEp <- NA_real_     s2w <- 1     level <- .level(level)     Y <- as.matrix(yi)     ymci <- scale(yi, center = TRUE, scale = FALSE)     Ymc <- as.matrix(ymci)     tau2.inf <- FALSE     if (model == "rma.uni") {         if (!is.null(tau2) && !is.na(tau2) && !is.element(method[1], c("FE", "EE", "CE"))) {             tau2.fix <- TRUE             tau2.arg <- tau2             tau2.inf <- identical(tau2, Inf)         }         else {             tau2.fix <- FALSE             tau2.arg <- NA_real_         }         if (verbose > 1 && !tau2.fix && !is.element(method[1], c("FE", "EE", "CE")))              message(mstyle$message("Estimating the tau^2 value ...\n"))         if (k == 1L) {             method.sav <- method[1]             method <- "k1"             if (!tau2.fix)                  tau2 <- 0         }         conv <- FALSE         while (!conv && !tau2.inf) {             conv <- TRUE             change <- con$threshold + 1             iter <- 0             if (is.element(method[1], c("HS", "HSk"))) {                 if (!allvipos)                    stop(mstyle$stop(paste0(method[1], " estimator cannot be used when there are non-positive sampling variances in the data.")))                 wi <- 1/vi                 W <- .diag(wi)                 stXWX <- .invcalc(X = X, W = W, k = k)                 P <- W - W %*% X %*% stXWX %*% crossprod(X, W)                 RSS <- crossprod(Ymc, P) %*% Ymc                 if (method[1] == "HS") {                   tau2 <- ifelse(tau2.fix, tau2.arg, (RSS - k)/sum(wi))                 }                 else {                   tau2 <- ifelse(tau2.fix, tau2.arg, (k/(k - p) * RSS - k)/sum(wi))                 }             }             if (is.element(method[1], c("HE", "CO", "VC", "ML", "REML", "EB"))) {                 stXX <- .invcalc(X = X, W = diag(k), k = k)                 P <- diag(k) - X %*% tcrossprod(stXX, X)                 RSS <- crossprod(Ymc, P) %*% Ymc                 V <- .diag(vi)                 PV <- P %*% V                 trPV <- .tr(PV)                 tau2 <- ifelse(tau2.fix, tau2.arg, (RSS - trPV)/(k - p))             }             if (method[1] == "DL") {                 if (!allvipos)                    stop(mstyle$stop("DL estimator cannot be used when there are non-positive sampling variances in the data."))                 wi <- 1/vi                 W <- .diag(wi)                 stXWX <- .invcalc(X = X, W = W, k = k)                 P <- W - W %*% X %*% stXWX %*% crossprod(X, W)                 RSS <- crossprod(Ymc, P) %*% Ymc                 trP <- .tr(P)                 tau2 <- ifelse(tau2.fix, tau2.arg, (RSS - (k - p))/trP)             }             if (method[1] == "DLIT") {                 if (is.null(con$tau2.init)) {                   tau2 <- 0                 }                 else {                   tau2 <- con$tau2.init                 }                 while (change > con$threshold) {                   if (verbose)                      cat(mstyle$verbose(paste("Iteration", formatC(iter, width = 5, flag = "-", format = "f", digits = 0), "tau^2 =", fmtx(tau2, digits[["var"]]), "\n")))                   iter <- iter + 1                   old2 <- tau2                   wi <- 1/(vi + tau2)                   if (any(tau2 + vi < 0))                      stop(mstyle$stop("Some marginal variances are negative."))                   if (any(is.infinite(wi)))                      stop(mstyle$stop("Division by zero when computing the inverse variance weights."))                   W <- .diag(wi)                   stXWX <- .invcalc(X = X, W = W, k = k)                   P <- W - W %*% X %*% stXWX %*% crossprod(X, W)                   V <- .diag(vi)                   trP <- .tr(P)                   trPV <- .tr(P %*% V)                   RSS <- crossprod(Ymc, P) %*% Ymc                   tau2 <- ifelse(tau2.fix, tau2.arg, (RSS - trPV)/trP)                   tau2[tau2 < con$tau2.min] <- con$tau2.min                   change <- abs(old2 - tau2)                   if (iter > con$maxiter) {                     conv <- FALSE                     break                   }                 }                 if (!conv) {                   if (length(method) == 1L) {                     stop(mstyle$stop("Iterative DL estimator did not converge."))                   }                   else {                     if (verbose)                        warning(mstyle$warning("Iterative DL estimator did not converge."), call. = FALSE)                   }                 }             }             if (method[1] == "GENQ") {                 if (is.null(weights))                    stop(mstyle$stop("Must specify the 'weights' argument when method='GENQ'."))                 A <- .diag(weights)                 stXAX <- .invcalc(X = X, W = A, k = k)                 P <- A - A %*% X %*% stXAX %*% crossprod(X, A)                 V <- .diag(vi)                 PV <- P %*% V                 trP <- .tr(P)                 trPV <- .tr(PV)                 RSS <- crossprod(Ymc, P) %*% Ymc                 tau2 <- ifelse(tau2.fix, tau2.arg, (RSS - trPV)/trP)             }             if (method[1] == "GENQM") {                 if (is.null(weights))                    stop(mstyle$stop("Must specify the 'weights' argument when method='GENQM'."))                 A <- .diag(weights)                 stXAX <- .invcalc(X = X, W = A, k = k)                 P <- A - A %*% X %*% stXAX %*% crossprod(X, A)                 V <- .diag(vi)                 PV <- P %*% V                 trP <- .tr(P)                 if (!tau2.fix) {                   RSS <- crossprod(Ymc, P) %*% Ymc                   if (.GENQ.func(con$tau2.min, P = P, vi = vi, Q = RSS, level = 0, k = k, p = p, getlower = TRUE) > 0.5) {                     tau2 <- con$tau2.min                   }                   else {                     if (.GENQ.func(con$tau2.max, P = P, vi = vi, Q = RSS, level = 0, k = k, p = p, getlower = TRUE) < 0.5) {                       conv <- FALSE                       if (length(method) == 1L) {                         stop(mstyle$stop("Value of 'tau2.max' too low. Try increasing 'tau2.max' or switch to another 'method'."))                       }                       else {                         if (verbose)                            warning(mstyle$warning("Value of 'tau2.max' too low. Try increasing 'tau2.max' or switch to another 'method'."), call. = FALSE)                       }                     }                     else {                       tau2 <- try(uniroot(.GENQ.func, interval = c(con$tau2.min, con$tau2.max), tol = con$tol, maxiter = con$maxiter, P = P, vi = vi, Q = RSS, level = 0.5, k = k, p = p, getlower = FALSE, verbose = verbose, digits = digits, extendInt = "no")$root, silent = TRUE)                       if (inherits(tau2, "try-error")) {                         conv <- FALSE                         if (length(method) == 1L) {                           stop(mstyle$stop("Error in iterative search for tau^2 using uniroot()."))                         }                         else {                           if (verbose)                              warning(mstyle$warning("Error in iterative search for tau^2 using uniroot()."), call. = FALSE)                         }                       }                     }                   }                 }                 else {                   tau2 <- tau2.arg                 }             }             if (method[1] == "SJ") {                 if (is.null(con$tau2.init)) {                   tau2.0 <- c(var(ymci) * (k - 1)/k)                 }                 else {                   tau2.0 <- con$tau2.init                 }                 wi <- 1/(vi + tau2.0)                 W <- .diag(wi)                 stXWX <- .invcalc(X = X, W = W, k = k)                 P <- W - W %*% X %*% stXWX %*% crossprod(X, W)                 RSS <- crossprod(Ymc, P) %*% Ymc                 V <- .diag(vi)                 PV <- P %*% V                 tau2 <- ifelse(tau2.fix, tau2.arg, tau2.0 * RSS/(k - p))             }             if (method[1] == "SJIT") {                 if (is.null(con$tau2.init)) {                   tau2 <- c(var(ymci) * (k - 1)/k)                 }                 else {                   tau2 <- con$tau2.init                 }                 tau2.0 <- tau2                 while (change > con$threshold) {                   if (verbose)                      cat(mstyle$verbose(paste("Iteration", formatC(iter, width = 5, flag = "-", format = "f", digits = 0), "tau^2 =", fmtx(tau2, digits[["var"]]), "\n")))                   iter <- iter + 1                   old2 <- tau2                   wi <- 1/(vi + tau2)                   W <- .diag(wi)                   stXWX <- .invcalc(X = X, W = W, k = k)                   P <- W - W %*% X %*% stXWX %*% crossprod(X, W)                   RSS <- crossprod(Ymc, P) %*% Ymc                   V <- .diag(vi)                   PV <- P %*% V                   tau2 <- ifelse(tau2.fix, tau2.arg, tau2 * RSS/(k - p))                   change <- abs(old2 - tau2)                   if (iter > con$maxiter) {                     conv <- FALSE                     break                   }                 }                 if (!conv) {                   if (length(method) == 1L) {                     stop(mstyle$stop("Iterative SJ estimator did not converge."))                   }                   else {                     if (verbose)                        warning(mstyle$warning("Iterative SJ estimator did not converge."), call. = FALSE)                   }                 }             }             if (is.element(method[1], c("PM", "MP", "PMM"))) {                 if (!allvipos)                    stop(mstyle$stop(method[1], " estimator cannot be used when there are non-positive sampling variances in the data."))                 if (method[1] == "PMM") {                   target <- qchisq(0.5, df = k - p)                 }                 else {                   target <- k - p                 }                 if (!tau2.fix) {                   if (.QE.func(con$tau2.min, Y = Ymc, vi = vi, X = X, k = k, objective = 0) < target) {                     tau2 <- con$tau2.min                   }                   else {                     if (.QE.func(con$tau2.max, Y = Ymc, vi = vi, X = X, k = k, objective = 0) > target) {                       conv <- FALSE                       if (length(method) == 1L) {                         stop(mstyle$stop("Value of 'tau2.max' too low. Try increasing 'tau2.max' or switch to another 'method'."))                       }                       else {                         if (verbose)                            warning(mstyle$warning("Value of 'tau2.max' too low. Try increasing 'tau2.max' or switch to another 'method'."), call. = FALSE)                       }                     }                     else {                       tau2 <- try(uniroot(.QE.func, interval = c(con$tau2.min, con$tau2.max), tol = con$tol, maxiter = con$maxiter, Y = Ymc, vi = vi, X = X, k = k, objective = target, verbose = verbose, digits = digits, extendInt = "no")$root, silent = TRUE)                       if (inherits(tau2, "try-error")) {                         conv <- FALSE                         if (length(method) == 1L) {                           stop(mstyle$stop("Error in iterative search for tau^2 using uniroot()."))                         }                         else {                           if (verbose)                              warning(mstyle$warning("Error in iterative search for tau^2 using uniroot()."), call. = FALSE)                         }                       }                     }                   }                 }                 else {                   tau2 <- tau2.arg                 }             }             if (is.element(method[1], c("ML", "REML", "EB"))) {                 if (is.null(con$tau2.init)) {                   tau2 <- max(0, tau2, con$tau2.min)                 }                 else {                   tau2 <- con$tau2.init                 }                 while (change > con$threshold) {                   if (verbose)                      cat(mstyle$verbose(paste(mstyle$verbose(paste("Iteration", formatC(iter, width = 5, flag = "-", format = "f", digits = 0), "tau^2 =", fmtx(tau2, digits[["var"]]), "\n")))))                   iter <- iter + 1                   old2 <- tau2                   wi <- 1/(vi + tau2)                   if (any(tau2 + vi < 0))                      stop(mstyle$stop("Some marginal variances are negative."))                   if (any(is.infinite(wi)))                      stop(mstyle$stop("Division by zero when computing the inverse variance weights."))                   W <- .diag(wi)                   stXWX <- .invcalc(X = X, W = W, k = k)                   P <- W - W %*% X %*% stXWX %*% crossprod(X, W)                   if (method[1] == "ML") {                     PP <- P %*% P                     adj <- c(crossprod(Ymc, PP) %*% Ymc - sum(wi))/sum(wi^2)                   }                   if (method[1] == "REML") {                     PP <- P %*% P                     adj <- c(crossprod(Ymc, PP) %*% Ymc - .tr(P))/.tr(PP)                   }                   if (method[1] == "EB") {                     adj <- c(crossprod(Ymc, P) %*% Ymc * k/(k - p) - k)/sum(wi)                   }                   adj <- c(adj) * con$stepadj                   if (is.na(adj))                      adj <- 0                   while (tau2 + adj < con$tau2.min) adj <- adj/2                   tau2 <- ifelse(tau2.fix, tau2.arg, tau2 + adj)                   change <- abs(old2 - tau2)                   if (iter > con$maxiter) {                     conv <- FALSE                     break                   }                 }                 if (!conv) {                   if (length(method) == 1L) {                     stop(mstyle$stop("Fisher scoring algorithm did not converge. See 'help(rma)' for possible remedies."))                   }                   else {                     if (verbose)                        warning(mstyle$warning("Fisher scoring algorithm did not converge. See 'help(rma)' for possible remedies."), call. = FALSE)                   }                 }                 if (conv && is.element(method[1], c("ML", "REML")) && con$ll0check && allvipos && !tau2.fix) {                   wi <- 1/vi                   W <- .diag(wi)                   stXWX <- .invcalc(X = X, W = W, k = k)                   beta <- stXWX %*% crossprod(X, W) %*% Ymc                   RSS <- sum(wi * (ymci - X %*% beta)^2)                   if (method[1] == "ML")                      ll0 <- -1/2 * (k) * log(2 * base::pi) - 1/2 * sum(log(vi)) - 1/2 * RSS                   if (method[1] == "REML")                      ll0 <- -1/2 * (k - p) * log(2 * base::pi) - 1/2 * sum(log(vi)) - 1/2 * determinant(crossprod(X, W) %*% X, logarithm = TRUE)$modulus - 1/2 * RSS                   wi <- 1/(vi + tau2)                   if (any(tau2 + vi < 0))                      stop(mstyle$stop("Some marginal variances are negative."))                   if (any(is.infinite(wi)))                      stop(mstyle$stop("Division by zero when computing the inverse variance weights."))                   W <- .diag(wi)                   stXWX <- .invcalc(X = X, W = W, k = k)                   beta <- stXWX %*% crossprod(X, W) %*% Ymc                   RSS <- sum(wi * (ymci - X %*% beta)^2)                   if (method[1] == "ML")                      ll <- -1/2 * (k) * log(2 * base::pi) - 1/2 * sum(log(vi + tau2)) - 1/2 * RSS                   if (method[1] == "REML")                      ll <- -1/2 * (k - p) * log(2 * base::pi) - 1/2 * sum(log(vi + tau2)) - 1/2 * determinant(crossprod(X, W) %*% X, logarithm = TRUE)$modulus - 1/2 * RSS                   if (ll0 - ll > con$tol && tau2 > con$threshold) {                     warning(mstyle$warning("Fisher scoring algorithm may have gotten stuck at a local maximum.\nSetting tau^2 = 0. Check the profile likelihood plot with profile()."), call. = FALSE)                     tau2 <- 0                   }                 }                 if (conv) {                   wi <- 1/(vi + tau2)                   if (any(tau2 + vi < 0))                      stop(mstyle$stop("Some marginal variances are negative."))                   if (any(is.infinite(wi)))                      stop(mstyle$stop("Division by zero when computing the inverse variance weights."))                   W <- .diag(wi)                   stXWX <- .invcalc(X = X, W = W, k = k)                   P <- W - W %*% X %*% stXWX %*% crossprod(X, W)                 }             }             if (conv) {                 tau2 <- max(con$tau2.min, c(tau2))                 if (!is.na(tau2) && any(tau2 + vi < 0))                    stop(mstyle$stop("Some marginal variances are negative."))                 if (verbose && is.element(method[1], c("ML", "REML", "EB"))) {                   cat(mstyle$verbose(paste("Iteration", formatC(iter, width = 5, flag = "-", format = "f", digits = 0), "tau^2 =", fmtx(tau2, digits[["var"]]), "\n")))                   cat(mstyle$verbose(paste("Fisher scoring algorithm converged after", iter, "iterations.\n")))                 }                 if (method[1] == "HS")                    se.tau2 <- sqrt(1/sum(wi)^2 * (2 * (k - p) + 4 * max(tau2, 0) * .tr(P) + 2 * max(tau2, 0)^2 * sum(P * P)))                 if (method[1] == "HSk")                    se.tau2 <- k/(k - p) * sqrt(1/sum(wi)^2 * (2 * (k - p) + 4 * max(tau2, 0) * .tr(P) + 2 * max(tau2, 0)^2 * sum(P * P)))                 if (is.element(method[1], c("HE", "CO", "VC")))                    se.tau2 <- sqrt(1/(k - p)^2 * (2 * sum(PV * t(PV)) + 4 * max(tau2, 0) * trPV + 2 * max(tau2, 0)^2 * (k - p)))                 if (method[1] == "DL")                    se.tau2 <- sqrt(1/trP^2 * (2 * (k - p) + 4 * max(tau2, 0) * trP + 2 * max(tau2, 0)^2 * sum(P * P)))                 if (is.element(method[1], c("GENQ", "GENQM")))                    se.tau2 <- sqrt(1/trP^2 * (2 * sum(PV * t(PV)) + 4 * max(tau2, 0) * sum(PV * P) + 2 * max(tau2, 0)^2 * sum(P * P)))                 if (method[1] == "SJ")                    se.tau2 <- sqrt(tau2.0^2/(k - p)^2 * (2 * sum(PV * t(PV)) + 4 * max(tau2, 0) * sum(PV * P) + 2 * max(tau2, 0)^2 * sum(P * P)))                 if (method[1] == "ML")                    se.tau2 <- sqrt(2/sum(wi^2))                 if (method[1] == "REML")                    se.tau2 <- sqrt(2/sum(P * P))                 if (is.element(method[1], c("EB", "PM", "MP", "PMM", "DLIT", "SJIT"))) {                   wi <- 1/(vi + tau2)                   se.tau2 <- sqrt(2 * k^2/(k - p)/sum(wi)^2)                 }             }             else {                 method <- method[-1]             }         }         if (k == 1L)              method <- method.sav     }     if (model == "rma.ls") {         if (!is.element(method[1], c("ML", "REML")))              stop(mstyle$stop("Location-scale models can only be fitted with ML or REML estimation."))         tau2.fix <- FALSE         if (!is.null(tau2) && !is.na(tau2))              warning(mstyle$warning("Argument 'tau2' ignored for location-scale models."), call. = FALSE)         optimizer <- match.arg(con$optimizer, c("optim", "nlminb", "uobyqa", "newuoa", "bobyqa", "nloptr", "nlm", "hjk", "nmk", "mads", "ucminf", "lbfgsb3c", "subplex", "BBoptim", "optimParallel", "constrOptim", "solnp", "alabama", "constrOptim.nl", "auglag", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent", "Rcgmin", "Rvmmin"))         optmethod <- match.arg(con$optmethod, c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"))         if (optimizer %in% c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent")) {             optmethod <- optimizer             optimizer <- "optim"         }         parallel <- con$parallel         cl <- con$cl         ncpus <- con$ncpus         optcontrol <- control[is.na(con.pos)]         if (length(optcontrol) == 0L)              optcontrol <- list()         if (ncpus > 1L)              optimizer <- "optimParallel"         if (optimizer == "alabama")              optimizer <- "constrOptim.nl"         if (link == "identity") {             if (optimizer == "nlminb") {                 optimizer <- "constrOptim"             }             else {                 if (!is.element(optimizer, c("constrOptim", "solnp", "nloptr", "constrOptim.nl", "auglag"))) {                   optimizer <- "constrOptim"                   warning(mstyle$warning(paste0("Can only use optimizers 'constrOptim', 'solnp', 'nloptr', 'constrOptim.nl', or 'auglag' when link='identity' (resetting to '", optimizer, "').")), call. = FALSE)                 }             }         }         if (link == "log" && is.element(optimizer, c("constrOptim", "constrOptim.nl", "auglag")))              stop(mstyle$stop(paste0("Cannot use '", optimizer, "' optimizer when using a log link.")))         reml <- ifelse(method[1] == "REML", TRUE, FALSE)         tmp <- try(lm(yi ~ 0 + Z), silent = TRUE)         if (inherits(tmp, "lm")) {             coef.na.Z <- is.na(coef(tmp))         }         else {             coef.na.Z <- rep(FALSE, NCOL(Z))         }         if (any(coef.na.Z)) {             warning(mstyle$warning("Redundant predictors dropped from the scale model."), call. = FALSE)             Z <- Z[, !coef.na.Z, drop = FALSE]             Z.f <- Z.f[, !coef.na.Z, drop = FALSE]         }         is.int <- apply(Z, 2, .is.intercept)         if (any(is.int)) {             Z.int.incl <- TRUE             int.indx <- which(is.int, arr.ind = TRUE)             Z <- cbind(intrcpt = 1, Z[, -int.indx, drop = FALSE])             Z.f <- cbind(intrcpt = 1, Z.f[, -int.indx, drop = FALSE])             Z.intercept <- TRUE         }         else {             Z.int.incl <- FALSE         }         q <- NCOL(Z)         if (!.chkpd(crossprod(Z), tol = con$evtol))              stop(mstyle$stop("Model matrix for scale part of the model not of full rank. Cannot fit model."))         is.int <- apply(Z, 2, .is.intercept)         if (q == 1L && is.int) {             Z.int.only <- TRUE         }         else {             Z.int.only <- FALSE         }         if (missing(alpha) || is.null(alpha) || all(is.na(alpha))) {             alpha <- rep(NA_real_, q)         }         else {             if (any(is.infinite(alpha)))                  stop(mstyle$stop("Infinite values in 'alpha' argument not allowed."))             alpha <- .expand1(alpha, q)             if (length(alpha) != q)                  stop(mstyle$stop(paste0("Length of the 'alpha' argument (", length(alpha), ") does not match the actual number of parameters (", q, ").")))         }         if (optbeta) {             if (missing(beta) || is.null(beta) || all(is.na(beta))) {                 beta <- rep(NA_real_, p)             }             else {                 beta <- .expand1(beta, p)                 if (length(beta) != p)                    stop(mstyle$stop(paste0("Length of the 'beta' argument (", length(beta), ") does not match the actual number of parameters (", p, ").")))             }             X0 <- X             X0[] <- 0         }         else {             X0 <- NULL         }         if (missing(omega2) || is.null(omega2) || all(is.na(omega2))) {             omega2 <- NA_real_         }         else {             if (length(omega2) != 1)                  stop(mstyle$stop("Length of the 'omega2' argument must be 1."))         }         if (!Z.int.only && Z.int.incl && con$scaleZ && is.na(alpha[1]) && all(is.infinite(con$alpha.min)) && all(is.infinite(con$alpha.max)) && !optbeta) {             Zsave <- Z             meanZ <- colMeans(Z[, 2:q, drop = FALSE])             sdZ <- apply(Z[, 2:q, drop = FALSE], 2, sd)             is.d <- apply(Z, 2, .is.dummy)             mZ <- rbind(c(intrcpt = 1, -1 * ifelse(is.d[-1], 0, meanZ/sdZ)), cbind(0, diag(ifelse(is.d[-1], 1, 1/sdZ), nrow = length(is.d) - 1, ncol = length(is.d) - 1)))             imZ <- try(suppressWarnings(solve(mZ)), silent = TRUE)             Z[, !is.d] <- apply(Z[, !is.d, drop = FALSE], 2, scale)             if (any(!is.na(alpha))) {                 if (inherits(imZ, "try-error"))                    stop(mstyle$stop("Unable to rescale starting values for the scale parameters."))                 alpha <- diag(imZ) * alpha             }         }         else {             mZ <- NULL         }         if (k == 1L && Z.int.only) {             if (link == "log")                  con$alpha.init <- -10000             if (link == "identity")                  con$alpha.init <- 1e-05         }         if (verbose > 1)              message(mstyle$message("Extracting/computing the initial values ..."))         if (is.null(con$alpha.init)) {             fit <- suppressWarnings(rma.uni(yi, vi, mods = X, intercept = FALSE, method = "HE", skipr2 = TRUE))             tmp <- rstandard(fit)             if (link == "log") {                 tmp <- suppressWarnings(rma.uni(log(tmp$resid^2), 4/tmp$resid^2 * tmp$se^2, mods = Z, intercept = FALSE, method = "FE"))                 alpha.init <- coef(tmp)             }             if (link == "identity") {                 tmp <- suppressWarnings(rma.uni(tmp$resid^2, tmp$se^2, mods = Z, intercept = FALSE, method = "FE"))                 alpha.init <- coef(tmp)                 if (any(Z %*% alpha.init < 0))                    alpha.init <- ifelse(is.int, fit$tau2 + 0.01, 0)                 if (any(Z %*% alpha.init < 0))                    stop(mstyle$stop("Unable to find suitable starting values for the scale parameters."))             }         }         else {             alpha.init <- con$alpha.init             if (!is.null(mZ)) {                 if (inherits(imZ, "try-error"))                    stop(mstyle$stop("Unable to rescale the starting values for the scale parameters."))                 alpha.init <- c(imZ %*% cbind(alpha.init))             }             if (link == "identity" && any(Z %*% alpha.init < 0))                  stop(mstyle$stop("Starting values for the scale parameters lead to one or more negative tau^2 values."))             if (optbeta)                  fit <- suppressWarnings(rma.uni(yi, vi, mods = X, intercept = FALSE, method = "HE", skipr2 = TRUE))         }         if (length(alpha.init) != q)              stop(mstyle$stop(paste0("Length of the 'alpha.init' argument (", length(alpha.init), ") does not match the actual number of parameters (", q, ").")))         if (anyNA(alpha.init))              stop(mstyle$stop("No missing values allowed in 'alpha.init'."))         if (optbeta) {             if (is.null(con$beta.init)) {                 beta.init <- c(fit$beta)             }             else {                 beta.init <- con$beta.init                 if (length(beta.init) != p)                    stop(mstyle$stop(paste0("Length of the 'beta.init' argument (", length(beta.init), ") does not match the actual number of parameters (", p, ").")))                 if (anyNA(beta.init))                    stop(mstyle$stop("No missing values allowed in 'beta.init'."))             }         }         else {             beta.init <- NULL         }         if (is.null(con$omega2.init)) {             omega2.init <- log(0.1)         }         else {             if (length(con$omega2.init) != 1L)                  stop(mstyle$stop("Argument 'omega2.init' should specify a single value."))             if (is.na(con$omega2.init))                  stop(mstyle$stop("No missing value allowed in 'omega2.init'."))             if (con$omega2.init <= 0)                  stop(mstyle$stop("Value of 'omega2.init' must be > 0."))             omega2.init <- log(con$omega2.init)         }         con$alpha.min <- .expand1(con$alpha.min, q)         con$alpha.max <- .expand1(con$alpha.max, q)         if (length(con$alpha.min) != q)              stop(mstyle$stop(paste0("Length of the 'alpha.min' argument (", length(alpha.min), ") does not match the actual number of parameters (", q, ").")))         if (length(con$alpha.max) != q)              stop(mstyle$stop(paste0("Length of the 'alpha.max' argument (", length(alpha.max), ") does not match the actual number of parameters (", q, ").")))         if (any(xor(is.infinite(con$alpha.min), is.infinite(con$alpha.max))))              stop(mstyle$stop("Constraints on the scale coefficients must be placed on both the lower and upper bound."))         alpha.min <- con$alpha.min         alpha.max <- con$alpha.max         if (link == "identity" && (any(alpha.min != -Inf) || any(alpha.max != Inf)))              stop(mstyle$stop("Cannot use constraints on scale coefficients when using an identity link."))         alpha.init <- pmax(alpha.init, alpha.min)         alpha.init <- pmin(alpha.init, alpha.max)         alpha.init <- mapply(.mapinvfun.alpha, alpha.init, alpha.min, alpha.max)         if (verbose > 1)              message(mstyle$message("Estimating the scale parameters ...\n"))         tmp <- .chkopt(optimizer, optcontrol, ineq = link == "identity")         optimizer <- tmp$optimizer         optcontrol <- tmp$optcontrol         par.arg <- tmp$par.arg         ctrl.arg <- tmp$ctrl.arg         if (optimizer == "optimParallel::optimParallel") {             parallel$cl <- NULL             if (is.null(cl)) {                 ncpus <- as.integer(ncpus)                 if (ncpus < 1L)                    stop(mstyle$stop("Control argument 'ncpus' must be >= 1."))                 cl <- parallel::makePSOCKcluster(ncpus)                 on.exit(parallel::stopCluster(cl), add = TRUE)             }             else {                 if (!inherits(cl, "SOCKcluster"))                    stop(mstyle$stop("Specified cluster is not of class 'SOCKcluster'."))             }             parallel$cl <- cl             if (is.null(parallel$forward))                  parallel$forward <- FALSE             if (is.null(parallel$loginfo)) {                 if (verbose) {                   parallel$loginfo <- TRUE                 }                 else {                   parallel$loginfo <- FALSE                 }             }         }         if (link == "log") {             if (randhet) {                 init.str <- "=c(beta.init, alpha.init, omega2.init), "             }             else {                 init.str <- "=c(beta.init, alpha.init), "             }             optcall <- paste0(optimizer, "(", par.arg, init.str, " .ll.rma.ls, ", ifelse(optimizer == "optim", "method=optmethod, ", ""), "yi=yi, vi=vi, X=X, Z=Z, reml=reml, alpha.arg=alpha, beta.arg=beta, omega2.arg=omega2, verbose=verbose, digits=digits,\n                          REMLf=con$REMLf, link=link, mZ=mZ, alpha.min=alpha.min, alpha.max=alpha.max, alpha.transf=TRUE, omega2.transf=TRUE,\n                          tau2.min=con$tau2.min, tau2.max=con$tau2.max, optbeta=optbeta, randhet=randhet, mfmaxit=con$mfmaxit",                  ctrl.arg, ")\n")         }         if (link == "identity") {             if (optimizer == "constrOptim")                  optcall <- paste0("constrOptim(theta=c(beta.init, alpha.init), f=.ll.rma.ls, grad=NULL, ui=cbind(X0,Z), ci=rep(0,k),\n                              yi=yi, vi=vi, X=X, Z=Z, reml=reml, alpha.arg=alpha, beta.arg=beta, omega2.arg=omega2, verbose=verbose, digits=digits,\n                              REMLf=con$REMLf, link=link, mZ=mZ, alpha.min=alpha.min, alpha.max=alpha.max, alpha.transf=TRUE, omega2.transf=TRUE,\n                              tau2.min=con$tau2.min, tau2.max=con$tau2.max, optbeta=optbeta, randhet=randhet, mfmaxit=con$mfmaxit",                    ctrl.arg, ")\n")             if (optimizer == "Rsolnp::solnp")                  optcall <- paste0("Rsolnp::solnp(pars=c(beta.init, alpha.init), fun=.ll.rma.ls, ineqfun=.rma.ls.ineqfun.pos, ineqLB=rep(0,k), ineqUB=rep(Inf,k),\n                              yi=yi, vi=vi, X=X, Z=Z, reml=reml, alpha.arg=alpha, beta.arg=beta, omega2.arg=omega2, verbose=verbose, digits=digits,\n                              REMLf=con$REMLf, link=link, mZ=mZ, alpha.min=alpha.min, alpha.max=alpha.max, alpha.transf=TRUE, omega2.transf=TRUE,\n                              tau2.min=con$tau2.min, tau2.max=con$tau2.max, optbeta=optbeta, randhet=randhet, mfmaxit=con$mfmaxit",                    ctrl.arg, ")\n")             if (optimizer == "nloptr::nloptr")                  optcall <- paste0("nloptr::nloptr(x0=c(beta.init, alpha.init), eval_f=.ll.rma.ls, eval_g_ineq=.rma.ls.ineqfun.neg,\n                              yi=yi, vi=vi, X=X, Z=Z, reml=reml, alpha.arg=alpha, beta.arg=beta, omega2.arg=omega2, verbose=verbose, digits=digits,\n                              REMLf=con$REMLf, link=link, mZ=mZ, alpha.min=alpha.min, alpha.max=alpha.max, alpha.transf=TRUE, omega2.transf=TRUE,\n                              tau2.min=con$tau2.min, tau2.max=con$tau2.max, optbeta=optbeta, randhet=randhet, mfmaxit=con$mfmaxit",                    ctrl.arg, ")\n")             if (is.element(optimizer, c("alabama::constrOptim.nl", "alabama::auglag")))                  optcall <- paste0(optimizer, "(par=c(beta.init, alpha.init), fn=.ll.rma.ls, hin=.rma.ls.ineqfun.pos,\n                              yi=yi, vi=vi, X=X, Z=Z, reml=reml, alpha.arg=alpha, beta.arg=beta, omega2.arg=omega2, verbose=verbose, digits=digits,\n                              REMLf=con$REMLf, link=link, mZ=mZ, alpha.min=alpha.min, alpha.max=alpha.max, alpha.transf=TRUE, omega2.transf=TRUE,\n                              tau2.min=con$tau2.min, tau2.max=con$tau2.max, optbeta=optbeta, randhet=randhet, mfmaxit=con$mfmaxit",                    ctrl.arg, ")\n")         }         iteration <- 0         try(assign("iteration", iteration, envir = .metafor), silent = TRUE)         if (verbose) {             opt.res <- try(eval(str2lang(optcall)), silent = !verbose)         }         else {             opt.res <- try(suppressWarnings(eval(str2lang(optcall))), silent = !verbose)         }         if (isTRUE(ddd$retopt))              return(opt.res)         opt.res <- .chkconv(optimizer = optimizer, opt.res = opt.res, optcontrol = optcontrol, fun = "rma", verbose = verbose, paronly = FALSE)         if (optbeta) {             pos.beta <- 1:p             pos.alpha <- (p + 1):(p + q)             if (randhet) {                 pos.omega2 <- p + q + 1             }             else {                 pos.omega2 <- integer(0)             }         }         else {             pos.alpha <- 1:q             pos.beta <- integer(0)             pos.omega2 <- integer(0)         }         opt.res$par[pos.alpha] <- mapply(.mapfun.alpha, opt.res$par[pos.alpha], alpha.min, alpha.max)         opt.res$par[pos.beta] <- ifelse(is.na(beta), opt.res$par[pos.beta], beta)         opt.res$par[pos.alpha] <- ifelse(is.na(alpha), opt.res$par[pos.alpha], alpha)         beta.arg <- beta         alpha.arg <- alpha         omega2.arg <- omega2         beta <- cbind(opt.res$par[pos.beta])         alpha <- cbind(opt.res$par[pos.alpha])         lnomega2 <- unname(opt.res$par[pos.omega2])         omega2 <- exp(lnomega2)         H <- NA_real_         if (con$hes.beta.fix) {             beta.hes <- beta         }         else {             beta.hes <- beta.arg         }         if (con$hes.alpha.fix) {             alpha.hes <- alpha         }         else {             alpha.hes <- alpha.arg         }         if (con$hes.omega2.fix || isTRUE(omega2 < con$omega2tol)) {             omega2.hes <- omega2         }         else {             omega2.hes <- omega2.arg         }         if (optbeta) {             if (randhet) {                 vcovmat <- matrix(NA_real_, nrow = p + q + 1, ncol = p + q + 1)                 hest <- c(is.na(beta.hes), is.na(alpha.hes), is.na(omega2.hes))                 hesspars <- c(beta, alpha, omega2)             }             else {                 vcovmat <- matrix(NA_real_, nrow = p + q, ncol = p + q)                 hest <- c(is.na(beta.hes), is.na(alpha.hes))                 hesspars <- c(beta, alpha)             }         }         else {             vcovmat <- matrix(NA_real_, nrow = q, ncol = q)             hest <- is.na(alpha.hes)             hesspars <- c(alpha)         }         if (any(hest) && !isTRUE(ddd$skiphes)) {             if (verbose > 1)                  message(mstyle$message("\nComputing the Hessian ..."))             if (con$htransf) {                 if (con$hesspack == "numDeriv")                    H <- try(numDeriv::hessian(func = .ll.rma.ls, x = opt.res$par, method.args = con$hessianCtrl, yi = yi, vi = vi, X = X, Z = Z, reml = reml, alpha.arg = alpha.hes, beta.arg = beta.hes, omega2.arg = omega2.hes, verbose = FALSE, digits = digits, REMLf = con$REMLf, link = link, mZ = mZ, alpha.min = alpha.min, alpha.max = alpha.max, alpha.transf = FALSE, omega2.transf = TRUE, tau2.min = con$tau2.min, tau2.max = con$tau2.max, optbeta = optbeta, randhet = randhet, mfmaxit = con$mfmaxit),                      silent = TRUE)                 if (con$hesspack == "pracma")                    H <- try(pracma::hessian(f = .ll.rma.ls, x0 = opt.res$par, h = con$hessianCtrl$h, yi = yi, vi = vi, X = X, Z = Z, reml = reml, alpha.arg = alpha.hes, beta.arg = beta.hes, omega2.arg = omega2.hes, verbose = FALSE, digits = digits, REMLf = con$REMLf, link = link, mZ = mZ, alpha.min = alpha.min, alpha.max = alpha.max, alpha.transf = FALSE, omega2.transf = TRUE, tau2.min = con$tau2.min, tau2.max = con$tau2.max, optbeta = optbeta, randhet = randhet, mfmaxit = con$mfmaxit), silent = TRUE)                 if (con$hesspack == "calculus")                    H <- try(calculus::hessian(f = .ll.rma.ls, var = opt.res$par, accuracy = con$hessianCtrl$accuracy, stepsize = con$hessianCtrl$stepsize, params = list(yi = yi, vi = vi, X = X, Z = Z, reml = reml, alpha.arg = alpha.hes, beta.arg = beta.hes, omega2.arg = omega2.hes, verbose = FALSE, digits = digits, REMLf = con$REMLf, link = link, mZ = mZ, alpha.min = alpha.min, alpha.max = alpha.max, alpha.transf = FALSE, omega2.transf = TRUE, tau2.min = con$tau2.min, tau2.max = con$tau2.max, optbeta = optbeta,                      randhet = randhet, mfmaxit = con$mfmaxit)), silent = TRUE)             }             else {                 if (con$hesspack == "numDeriv")                    H <- try(numDeriv::hessian(func = .ll.rma.ls, x = hesspars, method.args = con$hessianCtrl, yi = yi, vi = vi, X = X, Z = Z, reml = reml, alpha.arg = alpha.hes, beta.arg = beta.hes, omega2.arg = omega2.hes, verbose = FALSE, digits = digits, REMLf = con$REMLf, link = link, mZ = mZ, alpha.min = alpha.min, alpha.max = alpha.max, alpha.transf = FALSE, omega2.transf = FALSE, tau2.min = con$tau2.min, tau2.max = con$tau2.max, optbeta = optbeta, randhet = randhet, mfmaxit = con$mfmaxit),                      silent = TRUE)                 if (con$hesspack == "pracma")                    H <- try(pracma::hessian(f = .ll.rma.ls, x0 = hesspars, h = con$hessianCtrl$h, yi = yi, vi = vi, X = X, Z = Z, reml = reml, alpha.arg = alpha.hes, beta.arg = beta.hes, omega2.arg = omega2.hes, verbose = FALSE, digits = digits, REMLf = con$REMLf, link = link, mZ = mZ, alpha.min = alpha.min, alpha.max = alpha.max, alpha.transf = FALSE, omega2.transf = FALSE, tau2.min = con$tau2.min, tau2.max = con$tau2.max, optbeta = optbeta, randhet = randhet, mfmaxit = con$mfmaxit), silent = TRUE)                 if (con$hesspack == "calculus")                    H <- try(calculus::hessian(f = .ll.rma.ls, var = hesspars, accuracy = con$hessianCtrl$accuracy, stepsize = con$hessianCtrl$stepsize, params = list(yi = yi, vi = vi, X = X, Z = Z, reml = reml, alpha.arg = alpha.hes, beta.arg = beta.hes, omega2.arg = omega2.hes, verbose = FALSE, digits = digits, REMLf = con$REMLf, link = link, mZ = mZ, alpha.min = alpha.min, alpha.max = alpha.max, alpha.transf = FALSE, omega2.transf = FALSE, tau2.min = con$tau2.min, tau2.max = con$tau2.max, optbeta = optbeta,                      randhet = randhet, mfmaxit = con$mfmaxit)), silent = TRUE)             }             if (inherits(H, "try-error")) {                 warning(mstyle$warning("Error when trying to compute the Hessian."), call. = FALSE)             }             else {                 if (optbeta) {                   if (randhet) {                     rownames(H) <- colnames(H) <- c(paste0("X", colnames(X)), paste0("Z", colnames(Z)), "omega2")                   }                   else {                     rownames(H) <- colnames(H) <- c(paste0("X", colnames(X)), paste0("Z", colnames(Z)))                   }                 }                 else {                   rownames(H) <- colnames(H) <- colnames(Z)                 }                 H.hest <- H[hest, hest, drop = FALSE]                 iH.hest <- try(suppressWarnings(chol2inv(chol(H.hest))), silent = TRUE)                 if (inherits(iH.hest, "try-error") || anyNA(iH.hest) || any(is.infinite(iH.hest))) {                   warning(mstyle$warning("Error when trying to invert the Hessian."), call. = FALSE)                 }                 else {                   vcovmat[hest, hest] <- iH.hest                 }             }         }         vb <- vcovmat[pos.beta, pos.beta, drop = FALSE]         va <- vcovmat[pos.alpha, pos.alpha, drop = FALSE]         vo <- vcovmat[pos.omega2, pos.omega2]         vo <- max(0, vo)         crit.omega2 <- qnorm(level/2, lower.tail = FALSE)         if (con$htransf) {             se.omega2 <- sqrt(vo) * omega2             ci.lb.omega2 <- exp(lnomega2 - crit.omega2 * sqrt(vo))             ci.ub.omega2 <- exp(lnomega2 + crit.omega2 * sqrt(vo))         }         else {             se.omega2 <- sqrt(vo)             ci.lb.omega2 <- omega2 - crit.omega2 * se.omega2             ci.ub.omega2 <- omega2 + crit.omega2 * se.omega2         }         ci.lb.omega2 <- max(0, ci.lb.omega2)         if (any(alpha <= alpha.min + 10 * .Machine$double.eps^0.25) || any(alpha >= alpha.max - 10 * .Machine$double.eps^0.25))              warning(mstyle$warning("One or more 'alpha' estimates are (almost) equal to their lower or upper bound.\nTreat results with caution (or consider adjusting 'alpha.min' and/or 'alpha.max')."), call. = FALSE)         if (!is.null(mZ)) {             alpha <- mZ %*% alpha             va[!hest, ] <- 0             va[, !hest] <- 0             va <- mZ %*% va %*% t(mZ)             va[!hest, ] <- NA_real_             va[, !hest] <- NA_real_             Z <- Zsave         }         att <- .set.btt(att, q, Z.int.incl, colnames(Z))         m.alpha <- length(att)         if (is.element(test, c("knha", "adhoc", "t"))) {             ddf.alpha <- k - q         }         else {             ddf.alpha <- NA_integer_         }         QS <- try(as.vector(t(alpha)[att] %*% chol2inv(chol(va[att, att])) %*% alpha[att]), silent = TRUE)         if (inherits(QS, "try-error"))              QS <- NA_real_         se.alpha <- sqrt(diag(va))         if (isTRUE(ddd$abbrev)) {             tmp <- colnames(Z)             tmp <- gsub("relevel(factor(", "", tmp, fixed = TRUE)             tmp <- gsub("\\), ref = \"[[:alnum:]]*\")", "", tmp)             tmp <- gsub("poly(", "", tmp, fixed = TRUE)             tmp <- gsub(", degree = [[:digit:]], raw = TRUE)", "^", tmp)             tmp <- gsub(", degree = [[:digit:]], raw = TRUE)", "^", tmp)             tmp <- gsub(", degree = [[:digit:]])", "^", tmp)             tmp <- gsub("rcs\\([[:alnum:]]*, [[:digit:]]\\)", "", tmp)             tmp <- gsub("factor(", "", tmp, fixed = TRUE)             tmp <- gsub("I(", "", tmp, fixed = TRUE)             tmp <- gsub(")", "", tmp, fixed = TRUE)             colnames(Z) <- tmp         }         rownames(alpha) <- rownames(va) <- colnames(va) <- colnames(Z)         names(se.alpha) <- NULL         zval.alpha <- c(alpha/se.alpha)         if (is.element(test, c("knha", "adhoc", "t"))) {             QS <- QS/m.alpha             QSdf <- c(m.alpha, ddf.alpha)             QSp <- if (QSdf[2] > 0)                  pf(QS, df1 = QSdf[1], df2 = QSdf[2], lower.tail = FALSE)             else NA_real_             pval.alpha <- if (ddf.alpha > 0)                  2 * pt(abs(zval.alpha), df = ddf.alpha, lower.tail = FALSE)             else rep(NA_real_, q)             crit.alpha <- if (ddf.alpha > 0)                  qt(level/2, df = ddf.alpha, lower.tail = FALSE)             else NA_real_         }         else {             QSdf <- c(m.alpha, ddf.alpha)             QSp <- pchisq(QS, df = QSdf[1], lower.tail = FALSE)             pval.alpha <- 2 * pnorm(abs(zval.alpha), lower.tail = FALSE)             crit.alpha <- qnorm(level/2, lower.tail = FALSE)         }         ci.lb.alpha <- c(alpha - crit.alpha * se.alpha)         ci.ub.alpha <- c(alpha + crit.alpha * se.alpha)         if (link == "log")              tau2 <- exp(as.vector(Z %*% alpha))         if (link == "identity")              tau2 <- as.vector(Z %*% alpha)     }     if (is.element(method[1], c("FE", "EE", "CE")))          tau2 <- 0     if (verbose > 1)          message(mstyle$message("\nModel fitting ..."))     wi <- 1/(vi + tau2)     W <- .diag(wi)     M <- .diag(vi + tau2)     if (weighted) {         if (is.null(weights) || is.element(test, c("knha", "adhoc"))) {             if (any(is.infinite(wi)))                  stop(mstyle$stop("Division by zero when computing the inverse variance weights."))             if (!optbeta) {                 if (tau2.inf) {                   beta <- cbind(coef(lm(yi ~ 0 + X)))                   vb <- .diag(rep(Inf, p))                 }                 else {                   stXWX <- .invcalc(X = X, W = W, k = k)                   beta <- stXWX %*% crossprod(X, W) %*% Y                   vb <- stXWX                 }             }             RSS.f <- sum(wi * c(yi - X %*% beta)^2)             RSS.knha <- RSS.f         }         if (!is.null(weights)) {             A <- .diag(weights)             stXAX <- .invcalc(X = X, W = A, k = k)             beta <- stXAX %*% crossprod(X, A) %*% Y             vb <- stXAX %*% t(X) %*% A %*% M %*% A %*% X %*% stXAX             RSS.f <- sum(wi * c(yi - X %*% beta)^2)         }         if (is.element(test, c("knha", "adhoc"))) {             if (RSS.knha <= .Machine$double.eps) {                 s2w <- 0             }             else {                 s2w <- RSS.knha/(k - p)             }         }     }     else {         stXX <- .invcalc(X = X, W = diag(k), k = k)         beta <- stXX %*% crossprod(X, Y)         vb <- tcrossprod(stXX, X) %*% M %*% X %*% stXX         RSS.f <- sum(wi * (yi - X %*% beta)^2)         if (is.element(test, c("knha", "adhoc"))) {             if (any(is.infinite(wi)))                  stop(mstyle$stop("Division by zero when computing the inverse variance weights."))             stXWX <- .invcalc(X = X, W = W, k = k)             beta.knha <- stXWX %*% crossprod(X, W) %*% Y             RSS.knha <- sum(wi * (yi - X %*% beta.knha)^2)             if (RSS.knha <= .Machine$double.eps) {                 s2w <- 0             }             else {                 s2w <- RSS.knha/(k - p)             }         }     }     if (verbose > 1)          message(mstyle$message("Conducting the tests of the fixed effects ..."))     if (is.element(method[1], c("FE", "EE", "CE")) && is.element(test, c("knha", "adhoc")))          warning(mstyle$warning(paste0("Knapp and Hartung method is not meant to be used in the context of '", method[1], "' models.")), call. = FALSE)     if (test == "adhoc")          s2w[s2w < 1] <- 1     vb <- s2w * vb     if (tau2.inf)          vb <- diag(rep(Inf, p))     if (is.element(test, c("knha", "adhoc", "t"))) {         ddf <- .chkddd(ddd$dfs, k - p, ddd$dfs[[1]])     }     else {         ddf <- NA_integer_     }     QM <- try(as.vector(t(beta)[btt] %*% chol2inv(chol(vb[btt, btt])) %*% beta[btt]), silent = TRUE)     if (inherits(QM, "try-error"))          QM <- NA_real_     if (isTRUE(ddd$abbrev)) {         tmp <- colnames(X)         tmp <- gsub("relevel(factor(", "", tmp, fixed = TRUE)         tmp <- gsub("\\), ref = \"[[:alnum:]]*\")", "", tmp)         tmp <- gsub("poly(", "", tmp, fixed = TRUE)         tmp <- gsub(", degree = [[:digit:]], raw = TRUE)", "^", tmp)         tmp <- gsub(", degree = [[:digit:]], raw = TRUE)", "^", tmp)         tmp <- gsub(", degree = [[:digit:]])", "^", tmp)         tmp <- gsub("rcs\\([[:alnum:]]*, [[:digit:]]\\)", "", tmp)         tmp <- gsub("factor(", "", tmp, fixed = TRUE)         tmp <- gsub("I(", "", tmp, fixed = TRUE)         tmp <- gsub(")", "", tmp, fixed = TRUE)         colnames(X) <- tmp     }     rownames(beta) <- rownames(vb) <- colnames(vb) <- colnames(X.f) <- colnames(X)     se <- sqrt(diag(vb))     names(se) <- NULL     zval <- c(beta/se)     if (is.element(test, c("knha", "adhoc", "t"))) {         QM <- QM/m         QMdf <- c(m, ddf)         QMp <- if (QMdf[2] > 0)              pf(QM, df1 = QMdf[1], df2 = QMdf[2], lower.tail = FALSE)         else NA_real_         pval <- if (ddf > 0)              2 * pt(abs(zval), df = ddf, lower.tail = FALSE)         else rep(NA_real_, p)         crit <- if (ddf > 0)              qt(level/2, df = ddf, lower.tail = FALSE)         else NA_real_     }     else {         QMdf <- c(m, ddf)         QMp <- pchisq(QM, df = QMdf[1], lower.tail = FALSE)         pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)         crit <- qnorm(level/2, lower.tail = FALSE)     }     ci.lb <- c(beta - crit * se)     ci.ub <- c(beta + crit * se)     if (verbose > 1)          message(mstyle$message("Conducting the heterogeneity test ..."))     if (allvipos) {         if (k > p) {             wi <- 1/vi             W.FE <- .diag(wi)             stXWX <- .invcalc(X = X, W = W.FE, k = k)             P <- W.FE - W.FE %*% X %*% stXWX %*% crossprod(X, W.FE)             QE <- max(0, c(crossprod(Ymc, P) %*% Ymc))             QEp <- pchisq(QE, df = k - p, lower.tail = FALSE)             if (i2def == "1")                  vt <- (k - p)/.tr(P)             if (i2def == "2")                  vt <- 1/mean(wi)             if (is.element(method[1], c("FE", "EE", "CE"))) {                 I2 <- max(0, 100 * (QE - (k - p))/QE)                 H2 <- QE/(k - p)             }             else {                 I2 <- 100 * tau2/(vt + tau2)                 H2 <- tau2/vt + 1             }         }         else {             QE <- 0             QEp <- 1             I2 <- 0             H2 <- 1             vt <- 0         }     }     else {         if (!vi0)              warning(mstyle$warning(paste0("Cannot compute ", ifelse(int.only, "Q", "QE"), "-test, I^2, or H^2 when there are non-positive sampling variances in the data.")), call. = FALSE)         vt <- NA_real_     }     if (verbose > 1)          message(mstyle$message("Computing fit statistics and log-likelihood ..."))     q.est <- ifelse(model == "rma.uni", 0, sum(is.na(alpha.arg)) + (randhet && is.na(omega2.arg)))     parms <- ifelse(optbeta, sum(is.na(beta.arg)), p) + ifelse(model == "rma.uni", ifelse(is.element(method[1], c("FE", "EE", "CE")) || tau2.fix, 0, 1), q.est)     if (randhet) {         ll.ML <- -opt.res$value         ll.REML <- NA_real_     }     else {         ll.ML <- -1/2 * (k) * log(2 * base::pi) - 1/2 * sum(log(vi + tau2)) - 1/2 * RSS.f         ll.REML <- -1/2 * (k - p) * log(2 * base::pi) + ifelse(con$REMLf, 1/2 * determinant(crossprod(X), logarithm = TRUE)$modulus, 0) + -1/2 * sum(log(vi + tau2)) - 1/2 * determinant(crossprod(X, W) %*% X, logarithm = TRUE)$modulus - 1/2 * RSS.f     }     if (k > p) {         if (allvipos) {             dev.ML <- -2 * (ll.ML - sum(dnorm(yi, mean = yi, sd = sqrt(vi), log = TRUE)))         }         else {             dev.ML <- -2 * ll.ML         }     }     else {         dev.ML <- 0     }     AIC.ML <- -2 * ll.ML + 2 * parms     BIC.ML <- -2 * ll.ML + parms * log(k)     AICc.ML <- -2 * ll.ML + 2 * parms * max(k, parms + 2)/(max(k, parms + 2) - parms - 1)     dev.REML <- -2 * (ll.REML - 0)     AIC.REML <- -2 * ll.REML + 2 * parms     BIC.REML <- -2 * ll.REML + parms * log(k - p)     AICc.REML <- -2 * ll.REML + 2 * parms * max(k - p, parms + 2)/(max(k - p, parms + 2) - parms - 1)     fit.stats <- matrix(c(ll.ML, dev.ML, AIC.ML, BIC.ML, AICc.ML, ll.REML, dev.REML, AIC.REML, BIC.REML, AICc.REML), ncol = 2, byrow = FALSE)     dimnames(fit.stats) <- list(c("ll", "dev", "AIC", "BIC", "AICc"), c("ML", "REML"))     fit.stats <- data.frame(fit.stats)     is.nested.in.int.only <- .is.nested(X, cbind(rep(1, k)))     if (model == "rma.uni" && !int.only && is.nested.in.int.only && !isTRUE(ddd$skipr2)) {         if (verbose > 1)              message(mstyle$message("Computing R^2 ..."))         if (is.element(method[1], c("FE", "EE", "CE"))) {             if (identical(var(yi), 0)) {                 R2 <- 0             }             else {                 if (weighted) {                   if (is.null(weights)) {                     R2 <- max(0, 100 * summary(lm(yi ~ X, weights = wi))$adj.r.squared)                   }                   else {                     R2 <- max(0, 100 * summary(lm(yi ~ X, weights = weights))$adj.r.squared)                   }                 }                 else {                   R2 <- max(0, 100 * summary(lm(yi ~ X))$adj.r.squared)                 }             }         }         else {             if (r2def %in% c("1", "1v", "3", "3v", "5", "6", "7", "8")) {                 args <- list(yi = yi, vi = vi, weights = weights, method = method, weighted = weighted, test = test, verbose = ifelse(verbose, TRUE, FALSE), control = con, digits = digits, outlist = "minimal")                 if (verbose > 1) {                   res0 <- try(.do.call(rma.uni, args), silent = FALSE)                 }                 else {                   res0 <- try(suppressWarnings(.do.call(rma.uni, args)), silent = TRUE)                 }                 if (!inherits(res0, "try-error")) {                   tau2.RE <- res0$tau2                   if (identical(tau2.RE, 0) && r2def %in% c("1", "3")) {                     R2 <- 0                   }                   else {                     ll0 <- logLik(res0)                     ll1 <- ifelse(method[1] == "ML", ll.ML, ll.REML)                     lls <- (ifelse(method[1] == "ML", dev.ML, dev.REML) + 2 * ll1)/2                     if (r2def == "1")                        R2 <- (tau2.RE - tau2)/tau2.RE                     if (r2def == "1v")                        R2 <- (tau2.RE - tau2)/(tau2.RE + 1/mean(1/vi))                     if (r2def == "3")                        R2 <- var(c(X %*% beta))/tau2.RE                     if (r2def == "3v")                        R2 <- var(c(X %*% beta))/(tau2.RE + 1/mean(1/vi))                     if (r2def == "5")                        R2 <- 1 - ll1/ll0                     if (r2def == "6")                        R2 <- 1 - (exp(ll0)/exp(ll1))^(2/k)                     if (r2def == "7")                        R2 <- (1 - (exp(ll0)/exp(ll1))^(2/k))/(1 - exp(ll0)^(2/k))                     if (r2def == "8")                        R2 <- (ll1 - ll0)/(lls - ll0)                   }                 }                 else {                   R2 <- NA_real_                 }             }             else {                 if (r2def == "2")                    R2 <- var(c(X %*% beta))/(var(c(X %*% beta)) + tau2)                 if (r2def == "2v")                    R2 <- var(c(X %*% beta))/(var(c(X %*% beta)) + tau2 + 1/mean(1/vi))                 if (r2def == "4")                    R2 <- cor(yi, c(X %*% beta))^2                 if (r2def == "4w") {                   if (is.null(weights)) {                     R2 <- cov.wt(cbind(yi, c(X %*% beta)), cor = TRUE, wt = 1/(vi + tau2))$cor[1, 2]^2                   }                   else {                     R2 <- cov.wt(cbind(yi, c(X %*% beta)), cor = TRUE, wt = weights)$cor[1, 2]^2                   }                 }             }             R2 <- max(0, 100 * R2)         }     }     else {         R2 <- NULL     }     if (isTRUE(ddd$pleasedonotreportI2thankyouverymuch)) {         I2 <- NA         H2 <- NA     }     if (verbose > 1)          message(mstyle$message("Preparing the output ..."))     p.eff <- p     k.eff <- k     if (is.null(ddd$outlist) || ddd$outlist == "nodata") {         res <- list(b = beta, beta = beta, se = se, zval = zval, pval = pval, ci.lb = ci.lb, ci.ub = ci.ub, vb = vb, tau2 = tau2, se.tau2 = se.tau2, tau2.fix = tau2.fix, tau2.f = tau2, I2 = I2, H2 = H2, R2 = R2, vt = vt, QE = QE, QEp = QEp, QM = QM, QMdf = QMdf, QMp = QMp, k = k, k.f = k.f, k.eff = k.eff, k.all = k.all, p = p, p.eff = p.eff, parms = parms, int.only = int.only, int.incl = int.incl, intercept = intercept, allvipos = allvipos, coef.na = coef.na, yi = yi, vi = vi, X = X, weights = weights,              yi.f = yi.f, vi.f = vi.f, X.f = X.f, weights.f = weights.f, M = M, chksumyi = digest::digest(as.vector(yi)), chksumvi = digest::digest(as.vector(vi)), chksumX = digest::digest(X), outdat.f = outdat.f, ni = ni, ni.f = ni.f, ids = ids, not.na = not.na, subset = subset, slab = slab, slab.null = slab.null, measure = measure, method = method[1], model = model, weighted = weighted, test = test, dfs = ddf, ddf = ddf, s2w = s2w, btt = btt, m = m, digits = digits, level = level, control = control,              verbose = verbose, add = add, to = to, drop00 = drop00, fit.stats = fit.stats, formula.yi = formula.yi, formula.mods = formula.mods, version = packageVersion("metafor"), call = mf)         if (is.null(ddd$outlist))              res <- append(res, list(data = data), which(names(res) == "fit.stats"))     }     else {         if (ddd$outlist == "minimal") {             res <- list(b = beta, beta = beta, se = se, zval = zval, pval = pval, ci.lb = ci.lb, ci.ub = ci.ub, vb = vb, tau2 = tau2, se.tau2 = se.tau2, tau2.fix = tau2.fix, I2 = I2, H2 = H2, R2 = R2, QE = QE, QEp = QEp, QM = QM, QMdf = QMdf, QMp = QMp, k = k, k.f = k.f, k.eff = k.eff, p = p, p.eff = p.eff, parms = parms, int.only = int.only, int.incl = int.incl, intercept = intercept, chksumyi = digest::digest(as.vector(yi)), chksumvi = digest::digest(as.matrix(vi)), chksumX = digest::digest(X),                  measure = measure, method = method[1], model = model, weighted = weighted, test = test, dfs = ddf, ddf = ddf, btt = btt, m = m, digits = digits, level = level, fit.stats = fit.stats)         }         else {             res <- eval(str2lang(paste0("list(", ddd$outlist, ")")))         }     }     if (model == "rma.ls") {         res$alpha <- alpha         res$va <- va         res$se.alpha <- se.alpha         res$zval.alpha <- zval.alpha         res$pval.alpha <- pval.alpha         res$ci.lb.alpha <- ci.lb.alpha         res$ci.ub.alpha <- ci.ub.alpha         res$alpha.fix <- !is.na(alpha.arg)         res$optbeta <- optbeta         if (optbeta) {             res$vcovmat <- vcovmat             res$beta.fix <- !is.na(beta.arg)         }         res$randhet <- randhet         if (randhet) {             res$omega2 <- omega2             res$omega2.fix <- !is.na(omega2.arg)             res$se.omega2 <- se.omega2             res$ci.lb.omega2 <- ci.lb.omega2             res$ci.ub.omega2 <- ci.ub.omega2         }         res$q <- q         res$alphas <- q         res$link <- link         res$Z <- Z         res$Z.f <- Z.f         res$tau2.f <- rep(NA_real_, k.f)         res$tau2.f[not.na] <- tau2         res$att <- att         res$m.alpha <- m.alpha         res$ddf.alpha <- ddf.alpha         res$QS <- QS         res$QSdf <- QSdf         res$QSp <- QSp         res$formula.scale <- formula.scale         res$Z.int.incl <- Z.int.incl         res$Z.intercept <- Z.int.incl         res$Z.int.only <- Z.int.only         res$coef.na.Z <- coef.na.Z         res$H <- H     }     time.end <- proc.time()     res$time <- unname(time.end - time.start)[3]     if (isTRUE(ddd$time))          .print.time(res$time)     if (verbose || isTRUE(ddd$time))          cat("\n")     if (model == "rma.ls") {         class(res) <- c("rma.ls", "rma.uni", "rma")     }     else {         class(res) <- c("rma.uni", "rma")     }     return(res) })(yi = c(0.608633627839699, 0.829063699821454, 0.758777742442889, -0.632681644068139, 0.583780204394991), sei = c(1.65045067235645, 0.926849181557613, 1.39949966536894, 0.39530076550589, 1.37058246455349), method = "REML")`: Fisher scoring algorithm did not converge. See 'help(rma)' for possible remedies.
  Backtrace:
       ▆
    1. └─confMeta::confMeta(...) at test-confMeta.weights.R:13:5
    2.   └─confMeta:::new_confMeta(...)
    3.     └─confMeta:::get_stats_others(...)
    4.       └─base::lapply(...)
    5.         └─confMeta (local) FUN(X[[i]], ...)
    6.           └─confMeta:::get_obj_fe(estimates = estimates, SEs = SEs, conf_level = conf_level)
    7.             └─meta::metagen(...)
    8.               └─meta:::hetcalc(...)
    9.                 └─meta:::runNN(...)
   10.                   ├─base::do.call(func, args)
   11.                   └─metafor (local) `<fn>`(yi = `<dbl>`, sei = `<dbl>`, method = "REML")
  
  [ FAIL 1 | WARN 0 | SKIP 0 | PASS 79 ]
  Error:
  ! Test failures.
  Execution halted