Better counting of levels in rma.mv() and added 'check.k.gtr.1' contr…

…ol argument.

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: metafor
-Version: 4.7-9
-Date: 2024-04-19
+Version: 4.7-10
+Date: 2024-04-29
 Title: Meta-Analysis Package for R
 Authors@R: person(given = "Wolfgang", family = "Viechtbauer", role = c("aut","cre"), email = "[email protected]", comment = c(ORCID = "0000-0003-3463-4063"))
 Depends: R (>= 4.0.0), methods, Matrix, metadat, numDeriv

NEWS.md

@@ -1,4 +1,6 @@
-# metafor 4.7-9 (2024-04-19)
+# metafor 4.7-10 (2024-04-29)
+
+- `rma.mv()` now counts the number of levels of a random effect more appropriately; this may trigger more often the check whether the number of levels is equal to 1, in which case the corresponding variance component is automatically fixed to 0; this check can be omitted with `control=list(check.k.gtr.1=FALSE)`
 
 - made optimizers `Rcgmin` and `Rvmmin` available again via the `optimx` package
 

R/confint.rma.ls.r

@@ -146,8 +146,8 @@ confint.rma.ls <- function(object, parm, level, fixed=FALSE, alpha, digits, tran
 
       ######################################################################
 
-      ### set control parameters for uniroot() and possibly replace with user-defined values
-      ### set vc.min and vc.max and possibly replace with user-defined values
+      ### set defaults for control parameters for uniroot() and replace with any user-defined values
+      ### set vc.min and vc.max and possibly replace with any user-defined values
 
       con <- list(tol=.Machine$double.eps^0.25, maxiter=1000, verbose=FALSE, eptries=10)
 

R/confint.rma.mv.r

@@ -293,8 +293,8 @@ confint.rma.mv <- function(object, parm, level, fixed=FALSE, sigma2, tau2, rho,
 
       ######################################################################
 
-      ### set control parameters for uniroot() and possibly replace with user-defined values
-      ### set vc.min and vc.max and possibly replace with user-defined values
+      ### set defaults for control parameters for uniroot() and replace with any user-defined values
+      ### set vc.min and vc.max and possibly replace with any user-defined values
 
       con <- list(tol=.Machine$double.eps^0.25, maxiter=1000, verbose=FALSE, eptries=10)
 

R/confint.rma.uni.r

@@ -97,8 +97,8 @@ confint.rma.uni <- function(object, parm, level, fixed=FALSE, random=TRUE, type,
 
       ######################################################################
 
-      ### set control parameters for uniroot() and possibly replace with user-defined values
-      ### set tau2.min and tau2.max and possibly replace with user-defined values
+      ### set defaults for control parameters for uniroot() and replace with any user-defined values
+      ### set tau2.min and tau2.max and replace with any user-defined values
       ### note: default tau2.min is smaller of 0 or tau2, since tau2 could in principle be negative
       ### note: default tau2.max must be larger than tau2 and tau2.min and really should be much larger (at least 100)
 

R/confint.rma.uni.selmodel.r

@@ -186,8 +186,8 @@ confint.rma.uni.selmodel <- function(object, parm, level, fixed=FALSE, tau2, del
 
       ######################################################################
 
-      ### set control parameters for uniroot() and possibly replace with user-defined values
-      ### set vc.min and vc.max and possibly replace with user-defined values
+      ### set defaults for control parameters for uniroot() and replace with any user-defined values
+      ### set vc.min and vc.max and possibly replace with any user-defined values
 
       con <- list(tol=.Machine$double.eps^0.25, maxiter=1000, verbose=FALSE, eptries=10)
 

R/hc.rma.uni.r

@@ -38,7 +38,8 @@ hc.rma.uni <- function(object, digits, transf, targs, control, ...) {
 
    #########################################################################
 
-   ### set control parameters for uniroot() and possibly replace with user-defined values
+   ### set defaults for control parameters for uniroot() and replace with any user-defined values
+
    con <- list(tol=.Machine$double.eps^0.25, maxiter=1000, verbose=FALSE)
    con.pos <- pmatch(names(control), names(con))
    con[c(na.omit(con.pos))] <- control[!is.na(con.pos)]

R/misc.func.hidden.mv.r

@@ -252,7 +252,7 @@
 
 ############################################################################
 
-.process.G.afterrmna <- function(mf.g, g.nlevels, g.levels, g.values, struct, formula, tau2, rho, Z.G1, Z.G2, isG, sparse, distspec, verbose) {
+.process.G.afterrmna <- function(mf.g, g.nlevels, g.levels, g.values, struct, formula, tau2, rho, Z.G1, Z.G2, isG, sparse, distspec, check.k.gtr.1, verbose) {
 
    mstyle <- .get.mstyle()
 
@@ -327,18 +327,18 @@
    if (is.element(struct, c("SPEXP","SPGAU","SPLIN","SPRAT","SPSPH","PHYBM","PHYPL","PHYPD","GEN","GDIAG"))) {
       g.levels.k <- table(factor(apply(mf.g[-nvars], 1, paste, collapse=" + "), levels=g.levels.f[[1]]))
    } else {
-      g.levels.k <- table(factor(mf.g[[1]], levels=g.levels.f[[1]]))
+      #g.levels.k <- table(factor(mf.g[[1]], levels=g.levels.f[[1]]))
+      g.levels.k <- apply(table(factor(mf.g[[1]], levels=g.levels.f[[1]]), mf.g[[2]]), 1, function(x) sum(x>0L))
    }
 
    ### for "HCS","UN","DIAG","HAR": if a particular level of the inner factor only occurs once, then set corresponding tau2 value to 0 (if not already fixed)
-   ### note: no longer done; variance component should still be (weakly) identifiable
-
-   #if (is.element(struct, c("HCS","UN","DIAG","HAR"))) {
-   #   if (any(is.na(tau2) & g.levels.k == 1)) {
-   #      tau2[is.na(tau2) & g.levels.k == 1] <- 0
-   #      warning(mstyle$warning("Inner factor has k=1 for one or more levels. Corresponding 'tau2' value(s) fixed to 0."), call.=FALSE)
-   #   }
-   #}
+
+   if (is.element(struct, c("HCS","UN","DIAG","HAR")) && check.k.gtr.1) {
+      if (any(is.na(tau2) & g.levels.k == 1)) {
+         tau2[is.na(tau2) & g.levels.k == 1] <- 0
+         warning(mstyle$warning("Inner factor has k=1 for one or more levels. Corresponding 'tau2' value(s) fixed to 0."), call.=FALSE)
+      }
+   }
 
    ### check if each study has only a single arm (could be different arms!)
    ### for "CS","HCS","AR","HAR","CAR" must then fix rho to 0 (if not already fixed)

R/permutest.rma.ls.r

@@ -50,7 +50,7 @@ permutest.rma.ls <- function(x, exact=FALSE, iter=1000, progbar=TRUE, digits, co
    if (skip.beta && skip.alpha)
       stop(mstyle$stop("Must run permutation test for at least one part of the model."))
 
-   ### set control parameters and possibly replace with user-defined values
+   ### set defaults for control parameters and replace with any user-defined values
 
    if (missing(control))
       control <- list()

R/permutest.rma.uni.r

@@ -95,7 +95,7 @@ permutest.rma.uni <- function(x, exact=FALSE, iter=1000, permci=FALSE, progbar=T
    if (!X.exact)
       seed <- as.integer(runif(1)*2e9)
 
-   ### set control parameters and possibly replace with user-defined values
+   ### set defaults for control parameters and replace with any user-defined values
 
    if (missing(control))
       control <- list()

R/rma.glmm.r

@@ -623,7 +623,7 @@ test="z", level=95, btt, nAGQ=7, verbose=FALSE, digits, control, ...) {
 
    #########################################################################
 
-   ### set default control parameters
+   ### set defaults for control parameters
 
    con <- list(verbose = FALSE,            # also passed on to glm/glmer/optim/nlminb/minqa (uobyqa/newuoa/bobyqa)
                package="lme4",             # package for fitting logistic mixed-effects models ("lme4", "GLMMadaptive", "glmmTMB")

R/rma.mv.r

@@ -200,6 +200,32 @@ cvvc=FALSE, sparse=FALSE, verbose=FALSE, digits, control, ...) {
    if (inherits(verbose, "try-error") || is.function(verbose) || length(verbose) != 1L || !(is.logical(verbose) || is.numeric(verbose)))
       stop(mstyle$stop("Argument 'verbose' must be a scalar (logical or numeric/integer)."))
 
+   ### set defaults for control parameters (part 1)
+
+   con <- list(verbose = FALSE,
+               optimizer = "nlminb",      # optimizer to use ("optim","nlminb","uobyqa","newuoa","bobyqa","nloptr","nlm","hjk","nmk","mads","ucminf","lbfgsb3c","subplex","BBoptim","optimParallel","Rcgmin","Rvmmin")
+               optmethod = "BFGS",        # argument 'method' for optim() ("Nelder-Mead" and "BFGS" are sensible options)
+               parallel = list(),         # parallel argument for optimParallel() (note: 'cl' argument in parallel is not passed; this is directly specified via 'cl')
+               cl = NULL,                 # arguments for optimParallel()
+               ncpus = 1L,                # arguments for optimParallel()
+               REMLf = TRUE,              # full REML likelihood (including all constants)
+               evtol = 1e-07,             # lower bound for eigenvalues to determine if model matrix is positive definite
+               nearpd = FALSE,            # to force G and H matrix to become positive definite
+               hessianCtrl = list(r=8),   # arguments passed on to 'method.args' of hessian()
+               hesstol = .Machine$double.eps^0.5, # threshold for detecting fixed elements in Hessian
+               hesspack = "numDeriv",     # package for computing the Hessian (numDeriv or pracma)
+               check.k.gtr.1 = TRUE)      # check that s.nlevels > 1 and g.levels.k > 1
+
+   ### replace defaults with any user-defined values
+
+   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
+
    ### set options(warn=1) if verbose > 2
 
    if (verbose > 2) {
@@ -1239,7 +1265,7 @@ cvvc=FALSE, sparse=FALSE, verbose=FALSE, digits, control, ...) {
          V[vi.neg,] <- 0 # note: entire row set to 0 (so covariances are also 0)
          V[,vi.neg] <- 0 # note: entire col set to 0 (so covariances are also 0)
          vi[vi.neg] <- 0
-         warning(mstyle$warning("Negative sampling variances constrained to zero."), call.=FALSE)
+         warning(mstyle$warning("Negative sampling variances constrained to 0."), call.=FALSE)
       }
    } else {
       allvipos <- TRUE
@@ -1312,6 +1338,11 @@ cvvc=FALSE, sparse=FALSE, verbose=FALSE, digits, control, ...) {
       int.incl <- FALSE
    }
 
+   ### check whether model matrix is of full rank
+
+   if (!.chkpd(crossprod(X), tol=con$evtol))
+      stop(mstyle$stop("Model matrix not of full rank. Cannot fit model."))
+
    ### number of columns in X (including the intercept if it is included)
 
    p <- NCOL(X)
@@ -1369,7 +1400,7 @@ cvvc=FALSE, sparse=FALSE, verbose=FALSE, digits, control, ...) {
 
       ### for any single-level factor with unfixed sigma2, fix the sigma2 value to 0
 
-      if (any(is.na(sigma2) & s.nlevels == 1)) {
+      if (any(is.na(sigma2) & s.nlevels == 1) && con$check.k.gtr.1) {
          sigma2[is.na(sigma2) & s.nlevels == 1] <- 0
          warning(mstyle$warning("Single-level factor(s) found in 'random' argument. Corresponding 'sigma2' value(s) fixed to 0."), call.=FALSE)
       }
@@ -1477,7 +1508,7 @@ cvvc=FALSE, sparse=FALSE, verbose=FALSE, digits, control, ...) {
 
    if (withG) {
 
-      tmp <- .process.G.afterrmna(mf.g, g.nlevels, g.levels, g.values, struct[1], formulas[[1]], tau2, rho, Z.G1, Z.G2, isG=TRUE, sparse, ddd$dist[[1]], verbose)
+      tmp <- .process.G.afterrmna(mf.g, g.nlevels, g.levels, g.values, struct[1], formulas[[1]], tau2, rho, Z.G1, Z.G2, isG=TRUE, sparse, ddd$dist[[1]], con$check.k.gtr.1, verbose)
 
       mf.g <- tmp$mf.g
 
@@ -1515,7 +1546,7 @@ cvvc=FALSE, sparse=FALSE, verbose=FALSE, digits, control, ...) {
 
    if (withH) {
 
-      tmp <- .process.G.afterrmna(mf.h, h.nlevels, h.levels, h.values, struct[2], formulas[[2]], gamma2, phi, Z.H1, Z.H2, isG=FALSE, sparse, ddd$dist[[2]], verbose)
+      tmp <- .process.G.afterrmna(mf.h, h.nlevels, h.levels, h.values, struct[2], formulas[[2]], gamma2, phi, Z.H1, Z.H2, isG=FALSE, sparse, ddd$dist[[2]], con$check.k.gtr.1, verbose)
 
       mf.h <- tmp$mf.g
 
@@ -1628,37 +1659,20 @@ cvvc=FALSE, sparse=FALSE, verbose=FALSE, digits, control, ...) {
 
    #########################################################################
 
-   ### set default control parameters
+   ### set default control parameters (part 2)
 
-   con <- list(verbose = FALSE,
-               optimizer = "nlminb",      # optimizer to use ("optim","nlminb","uobyqa","newuoa","bobyqa","nloptr","nlm","hjk","nmk","mads","ucminf","lbfgsb3c","subplex","BBoptim","optimParallel","Rcgmin","Rvmmin")
-               optmethod = "BFGS",        # argument 'method' for optim() ("Nelder-Mead" and "BFGS" are sensible options)
-               parallel = list(),         # parallel argument for optimParallel() (note: 'cl' argument in parallel is not passed; this is directly specified via 'cl')
-               cl = NULL,                 # arguments for optimParallel()
-               ncpus = 1L,                # arguments for optimParallel()
-               sigma2.init = sigma2.init, # initial value(s) for sigma2
-               tau2.init = tau2.init,     # initial value(s) for tau2
-               rho.init = rho.init,       # initial value(s) for rho
-               gamma2.init = gamma2.init, # initial value(s) for gamma2
-               phi.init = phi.init,       # initial value(s) for phi
-               REMLf = TRUE,              # full REML likelihood (including all constants)
-               evtol = 1e-07,             # lower bound for eigenvalues to determine if model matrix is positive definite
-               cholesky = ifelse(is.element(struct, c("UN","UNR","GEN")), TRUE, FALSE), # by default, use Cholesky factorization for G and H matrix for "UN", "UNR", and "GEN" structures
-               nearpd = FALSE,            # to force G and H matrix to become positive definite
-               hessianCtrl = list(r=8),   # arguments passed on to 'method.args' of hessian()
-               hesstol = .Machine$double.eps^0.5, # threshold for detecting fixed elements in Hessian
-               hesspack = "numDeriv")     # package for computing the Hessian (numDeriv or pracma)
+   con <- c(con, list(sigma2.init = sigma2.init, # initial value(s) for sigma2
+                      tau2.init = tau2.init,     # initial value(s) for tau2
+                      rho.init = rho.init,       # initial value(s) for rho
+                      gamma2.init = gamma2.init, # initial value(s) for gamma2
+                      phi.init = phi.init,       # initial value(s) for phi
+                      cholesky = ifelse(is.element(struct, c("UN","UNR","GEN")), TRUE, FALSE))) # by default, use Cholesky factorization for G and H matrix for "UN", "UNR", and "GEN" structures
 
    ### replace defaults with any user-defined values
 
    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
-
    ### when restart=TRUE, restart at current estimates
 
    if (isTRUE(ddd$restart)) {
@@ -1893,11 +1907,6 @@ cvvc=FALSE, sparse=FALSE, verbose=FALSE, digits, control, ...) {
 
    #########################################################################
 
-   ### check whether model matrix is of full rank
-
-   if (!.chkpd(crossprod(X), tol=con$evtol))
-      stop(mstyle$stop("Model matrix not of full rank. Cannot fit model."))
-
    ### which variance components are fixed? (TRUE/FALSE or NA if not applicable = not included)
 
    if (withS) {

R/rma.uni.r

@@ -950,7 +950,7 @@ test="z", level=95, btt, att, tau2, verbose=FALSE, digits, control, ...) {
       vi.neg <- vi < 0
       if (any(vi.neg)) {
          vi[vi.neg] <- 0
-         warning(mstyle$warning("Negative sampling variances constrained to zero."), call.=FALSE)
+         warning(mstyle$warning("Negative sampling variances constrained to 0."), call.=FALSE)
       }
    } else {
       allvipos <- TRUE
@@ -1056,7 +1056,7 @@ test="z", level=95, btt, att, tau2, verbose=FALSE, digits, control, ...) {
 
    #########################################################################
 
-   ### set default control parameters
+   ### set defaults for control parameters
 
    con <- list(verbose = FALSE,
                evtol = 1e-07,                  # lower bound for eigenvalues to determine if model matrix is positive definite (also for checking if vimaxmin >= 1/con$evtol)

R/selmodel.rma.uni.r

@@ -210,7 +210,7 @@ selmodel.rma.uni <- function(x, type, alternative="greater", prec, delta, steps,
 
    ############################################################################
 
-   ### set default control parameters
+   ### set defaults for control parameters
 
    con <- list(verbose = FALSE,
                delta.init = NULL,     # initial value(s) for selection model parameter(s)

R/zzz.r

@@ -1,6 +1,6 @@
 .onAttach <- function(libname, pkgname) {
 
-   ver <- "4.7-9"
+   ver <- "4.7-10"
 
    loadmsg <- paste0("\nLoading the 'metafor' package (version ", ver, "). For an\nintroduction to the package please type: help(metafor)\n")
 

README.md

@@ -5,7 +5,7 @@ metafor: A Meta-Analysis Package for R
 [![R build status](https://github.com/wviechtb/metafor/workflows/R-CMD-check/badge.svg)](https://github.com/wviechtb/metafor/actions)
 [![Code Coverage](https://codecov.io/gh/wviechtb/metafor/branch/master/graph/badge.svg)](https://app.codecov.io/gh/wviechtb/metafor)
 [![CRAN Version](https://www.r-pkg.org/badges/version/metafor)](https://cran.r-project.org/package=metafor)
-[![devel Version](https://img.shields.io/badge/devel-4.7--9-brightgreen.svg)](https://www.metafor-project.org/doku.php/installation#development_version)
+[![devel Version](https://img.shields.io/badge/devel-4.7--10-brightgreen.svg)](https://www.metafor-project.org/doku.php/installation#development_version)
 [![Monthly Downloads](https://cranlogs.r-pkg.org/badges/metafor)](https://cranlogs.r-pkg.org/badges/metafor)
 [![Total Downloads](https://cranlogs.r-pkg.org/badges/grand-total/metafor)](https://cranlogs.r-pkg.org/badges/grand-total/metafor)
 

_pkgdown.yml

@@ -32,13 +32,13 @@ home:
 
 figures:
   dev: grDevices::png
-  dpi: 96
+  dpi: 192
   dev.args: []
   fig.ext: png
   fig.width: 9
   fig.height: 8
   fig.retina: 2
-  fig.asp: 1
+  fig.asp: ~
 
 #figures:
 #  dev: ragg::agg_png