Added 'decreasing' argument to selmodel().

wviechtb · Mar 27, 2024 · 1c3f7b2 · 1c3f7b2
1 parent f6b7c07
commit 1c3f7b2
Show file tree

Hide file tree

Showing 149 changed files with 624 additions and 334 deletions.
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -4,9 +4,9 @@
 ^_pkgdown\.yml$
 ^docs$
 ^docs_old$
-^vignettes/pkgdown$
 ^data-raw$
 ^.*\.Rproj$
 ^\.Rproj\.user$
-^CITATION.cff$
+^CITATION\.cff$
 ^tests/testthat/images
+^vignettes/pkgdown$
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: metafor
-Version: 4.5-12
-Date: 2024-03-19
+Version: 4.5-13
+Date: 2024-03-27
 Title: Meta-Analysis Package for R
 Authors@R: person(given = "Wolfgang", family = "Viechtbauer", role = c("aut","cre"), email = "wvb@metafor-project.org", comment = c(ORCID = "0000-0003-3463-4063"))
 Depends: R (>= 4.0.0), methods, Matrix, metadat, numDeriv

diff --git a/NEWS.md b/NEWS.md
@@ -1,4 +1,4 @@
-# metafor 4.5-12 (2024-03-19)
+# metafor 4.5-13 (2024-03-27)
 
 - the `steps` argument in the various `profile()` functions can now also be a numeric vector to specify for which parameter values the likelihood should be evaluated
 
@@ -16,7 +16,11 @@
 
 - in `forest()` and `regplot()`, observation limits set via `olim` are now properly applied to all elements
 
-- `selmodel()` no longer automatically stops with an error when one or more intervals defined by the `steps` argument do not contain any observed p-values; instead a warning is issued and model fitting proceeds (but may fail)
+- various internal improvements to `selmodel()`
+
+- `selmodel()` no longer stops with an error when one or more intervals defined by the `steps` argument do not contain any observed p-values (instead a warning is issued and model fitting proceeds, but may fail)
+
+- added `decreasing` argument to `selmodel()` for enforcing that the delta estimates must be a monotonically decreasing function of the p-values in the step function model
 
 - added the undocumented argument `pval` to `selmodel()` for passing p-values directly to the function (doing this is highly experimental)
 

diff --git a/R/confint.rma.ls.r b/R/confint.rma.ls.r
@@ -47,7 +47,7 @@ confint.rma.ls <- function(object, parm, level, fixed=FALSE, alpha, digits, tran
    }
 
    if (x$optbeta)
-      stop(mstyle$stop("CI calculation not yet implemented when model was fit with 'optbeta=TRUE'."))
+      stop(mstyle$stop("CI calculation not yet implemented for models fitted with 'optbeta=TRUE'."))
 
    ### check if user has specified alpha argument
 
@@ -187,11 +187,11 @@ confint.rma.ls <- function(object, parm, level, fixed=FALSE, alpha, digits, tran
 
       vc.lb <- NA_real_
       vc.ub <- NA_real_
-      ci.null <- FALSE ### logical if CI is a null set
-      lb.conv <- FALSE ### logical if search converged for lower bound (LB)
-      ub.conv <- FALSE ### logical if search converged for upper bound (UB)
-      lb.sign <- ""    ### for sign in case LB must be below vc.min ("<") or above vc.max (">")
-      ub.sign <- ""    ### for sign in case UB must be below vc.min ("<") or above vc.max (">")
+      ci.null <- FALSE # logical if CI is a null set
+      lb.conv <- FALSE # logical if search converged for lower bound (LB)
+      ub.conv <- FALSE # logical if search converged for upper bound (UB)
+      lb.sign <- ""    # for sign in case LB must be below vc.min ("<") or above vc.max (">")
+      ub.sign <- ""    # for sign in case UB must be below vc.min ("<") or above vc.max (">")
 
       ######################################################################
       ######################################################################
@@ -238,6 +238,7 @@ confint.rma.ls <- function(object, parm, level, fixed=FALSE, alpha, digits, tran
                   }
 
                   ### check if uniroot method converged
+
                   if (!inherits(res, "try-error")) {
                      vc.lb <- res
                      lb.conv <- TRUE
@@ -286,6 +287,7 @@ confint.rma.ls <- function(object, parm, level, fixed=FALSE, alpha, digits, tran
                   }
 
                   ### check if uniroot method converged
+
                   if (!inherits(res, "try-error")) {
                      vc.ub <- res
                      ub.conv <- TRUE

diff --git a/R/confint.rma.mh.r b/R/confint.rma.mh.r
@@ -42,7 +42,7 @@ confint.rma.mh <- function(object, parm, level, digits, transf, targs, ...) {
 
    ### if requested, apply transformation function
 
-   if (.isTRUE(transf) && is.element(x$measure, c("OR","RR","IRR"))) ### if transf=TRUE, apply exp transformation to ORs, RRs, and IRRs
+   if (.isTRUE(transf) && is.element(x$measure, c("OR","RR","IRR"))) # if transf=TRUE, apply exp transformation to ORs, RRs, and IRRs
       transf <- exp
 
    if (is.function(transf)) {

diff --git a/R/confint.rma.mv.r b/R/confint.rma.mv.r
@@ -304,15 +304,15 @@ confint.rma.mv <- function(object, parm, level, fixed=FALSE, sigma2, tau2, rho,
       }
       if (comp == "rho") {
          if (is.element(x$struct[1], c("CS","HCS")))
-            con$vc.min <- -1                                  ### this will fail most of the time but with retries, this may get closer to actual lower bound
-            #con$vc.min <- min(-1/(x$g.nlevels.f[1] - 1), vc) ### this guarantees that cor matrix is semi-positive definite, but since V gets added, this is actually too strict
+            con$vc.min <- -1                                  # this will fail most of the time but with retries, this may get closer to actual lower bound
+            #con$vc.min <- min(-1/(x$g.nlevels.f[1] - 1), vc) # this guarantees that cor matrix is semi-positive definite, but since V gets added, this is actually too strict
          if (is.element(x$struct[1], c("AR","HAR","CAR")))
-            con$vc.min <- min(0, vc)                          ### negative autocorrelation parameters not considered (not even sensible for CAR)
+            con$vc.min <- min(0, vc)                          # negative autocorrelation parameters not considered (not even sensible for CAR)
          if (is.element(x$struct[1], c("UN","UNR","GEN")))
-            con$vc.min <- -1                                  ### TODO: this will often fail! (but with retries, this should still work)
+            con$vc.min <- -1                                  # TODO: this will often fail! (but with retries, this should still work)
          con$vc.max <- 1
          if (is.element(x$struct[1], c("SPEXP","SPGAU","SPLIN","SPRAT","SPSPH"))) {
-            con$vc.min <- 0                                   ### TODO: 0 basically always fails
+            con$vc.min <- 0                                   # TODO: 0 basically always fails
             con$vc.max <- max(10, vc*10)
          }
          if (is.element(x$struct[1], c("PHYPL","PHYPD"))) {
@@ -322,15 +322,15 @@ confint.rma.mv <- function(object, parm, level, fixed=FALSE, sigma2, tau2, rho,
       }
       if (comp == "phi") {
          if (is.element(x$struct[2], c("CS","HCS")))
-            con$vc.min <- -1                                  ### this will fail most of the time but with retries, this may get closer to actual lower bound
-            #con$vc.min <- min(-1/(x$h.nlevels.f[1] - 1), vc) ### this guarantees that cor matrix is semi-positive definite, but since V gets added, this is actually too strict
+            con$vc.min <- -1                                  # this will fail most of the time but with retries, this may get closer to actual lower bound
+            #con$vc.min <- min(-1/(x$h.nlevels.f[1] - 1), vc) # this guarantees that cor matrix is semi-positive definite, but since V gets added, this is actually too strict
          if (is.element(x$struct[2], c("AR","HAR","CAR")))
-            con$vc.min <- min(0, vc)                          ### negative autocorrelation parameters not considered (not even sensible for CAR)
+            con$vc.min <- min(0, vc)                          # negative autocorrelation parameters not considered (not even sensible for CAR)
          if (is.element(x$struct[2], c("UN","UNR","GEN")))
-            con$vc.min <- -1                                  ### TODO: this will often fail! (but with retries, this should still work)
+            con$vc.min <- -1                                  # TODO: this will often fail! (but with retries, this should still work)
          con$vc.max <- 1
          if (is.element(x$struct[2], c("SPEXP","SPGAU","SPLIN","SPRAT","SPSPH"))) {
-            con$vc.min <- 0                                   ### TODO: 0 basically always fails
+            con$vc.min <- 0                                   # TODO: 0 basically always fails
             con$vc.max <- max(10, vc*10)
          }
          if (is.element(x$struct[2], c("PHYPL","PHYPD"))) {
@@ -351,11 +351,11 @@ confint.rma.mv <- function(object, parm, level, fixed=FALSE, sigma2, tau2, rho,
 
       vc.lb <- NA_real_
       vc.ub <- NA_real_
-      ci.null <- FALSE ### logical if CI is a null set
-      lb.conv <- FALSE ### logical if search converged for lower bound (LB)
-      ub.conv <- FALSE ### logical if search converged for upper bound (UB)
-      lb.sign <- ""    ### for sign in case LB must be below vc.min ("<") or above vc.max (">")
-      ub.sign <- ""    ### for sign in case UB must be below vc.min ("<") or above vc.max (">")
+      ci.null <- FALSE # logical if CI is a null set
+      lb.conv <- FALSE # logical if search converged for lower bound (LB)
+      ub.conv <- FALSE # logical if search converged for upper bound (UB)
+      lb.sign <- ""    # for sign in case LB must be below vc.min ("<") or above vc.max (">")
+      ub.sign <- ""    # for sign in case UB must be below vc.min ("<") or above vc.max (">")
 
       ######################################################################
       ######################################################################
@@ -410,6 +410,7 @@ confint.rma.mv <- function(object, parm, level, fixed=FALSE, sigma2, tau2, rho,
                   }
 
                   ### check if uniroot method converged
+
                   if (!inherits(res, "try-error")) {
                      vc.lb <- res
                      lb.conv <- TRUE
@@ -466,6 +467,7 @@ confint.rma.mv <- function(object, parm, level, fixed=FALSE, sigma2, tau2, rho,
                   }
 
                   ### check if uniroot method converged
+
                   if (!inherits(res, "try-error")) {
                      vc.ub <- res
                      ub.conv <- TRUE

diff --git a/R/confint.rma.peto.r b/R/confint.rma.peto.r
@@ -42,7 +42,7 @@ confint.rma.peto <- function(object, parm, level, digits, transf, targs, ...) {
 
    ### if requested, apply transformation function
 
-   if (.isTRUE(transf)) ### if transf=TRUE, apply exp transformation to ORs
+   if (.isTRUE(transf)) # if transf=TRUE, apply exp transformation to ORs
       transf <- exp
 
    if (is.function(transf)) {

diff --git a/R/confint.rma.uni.r b/R/confint.rma.uni.r
@@ -126,11 +126,11 @@ confint.rma.uni <- function(object, parm, level, fixed=FALSE, random=TRUE, type,
 
       tau2.lb <- NA_real_
       tau2.ub <- NA_real_
-      ci.null <- FALSE ### logical if CI is a null set
-      lb.conv <- FALSE ### logical if search converged for lower bound (LB)
-      ub.conv <- FALSE ### logical if search converged for upper bound (UB)
-      lb.sign <- ""    ### for sign in case LB must be below tau2.min ("<") or above tau2.max (">")
-      ub.sign <- ""    ### for sign in case UB must be below tau2.min ("<") or above tau2.max (">")
+      ci.null <- FALSE # logical if CI is a null set
+      lb.conv <- FALSE # logical if search converged for lower bound (LB)
+      ub.conv <- FALSE # logical if search converged for upper bound (UB)
+      lb.sign <- ""    # for sign in case LB must be below tau2.min ("<") or above tau2.max (">")
+      ub.sign <- ""    # for sign in case UB must be below tau2.min ("<") or above tau2.max (">")
 
       ######################################################################
 

diff --git a/R/confint.rma.uni.selmodel.r b/R/confint.rma.uni.selmodel.r
@@ -9,9 +9,12 @@ confint.rma.uni.selmodel <- function(object, parm, level, fixed=FALSE, tau2, del
 
    x <- object
 
-   if (x$betaspec) ### TODO: consider providing CIs also for this case
+   if (x$betaspec) # TODO: consider providing CIs also for this case
       stop(mstyle$stop("Cannot obtain confidence intervals when one or more beta values were fixed."))
 
+   if (x$decreasing || x$type == "stepcon")
+      stop(mstyle$stop("Method not currently implemented for this type of model."))
+
    k <- x$k
    p <- x$p
 
@@ -210,11 +213,11 @@ confint.rma.uni.selmodel <- function(object, parm, level, fixed=FALSE, tau2, del
 
       vc.lb <- NA_real_
       vc.ub <- NA_real_
-      ci.null <- FALSE ### logical if CI is a null set
-      lb.conv <- FALSE ### logical if search converged for lower bound (LB)
-      ub.conv <- FALSE ### logical if search converged for upper bound (UB)
-      lb.sign <- ""    ### for sign in case LB must be below vc.min ("<") or above vc.max (">")
-      ub.sign <- ""    ### for sign in case UB must be below vc.min ("<") or above vc.max (">")
+      ci.null <- FALSE # logical if CI is a null set
+      lb.conv <- FALSE # logical if search converged for lower bound (LB)
+      ub.conv <- FALSE # logical if search converged for upper bound (UB)
+      lb.sign <- ""    # for sign in case LB must be below vc.min ("<") or above vc.max (">")
+      ub.sign <- ""    # for sign in case UB must be below vc.min ("<") or above vc.max (">")
 
       ######################################################################
       ######################################################################
@@ -268,6 +271,7 @@ confint.rma.uni.selmodel <- function(object, parm, level, fixed=FALSE, tau2, del
                   }
 
                   ### check if uniroot method converged
+
                   if (!inherits(res, "try-error")) {
                      vc.lb <- res
                      lb.conv <- TRUE
@@ -321,6 +325,7 @@ confint.rma.uni.selmodel <- function(object, parm, level, fixed=FALSE, tau2, del
                   }
 
                   ### check if uniroot method converged
+
                   if (!inherits(res, "try-error")) {
                      vc.ub <- res
                      ub.conv <- TRUE

diff --git a/R/misc.func.hidden.profile.r b/R/misc.func.hidden.profile.r
@@ -259,7 +259,7 @@
 
       res <- try(suppressWarnings(
          selmodel(obj, obj$type, alternative=obj$alternative, prec=obj$prec, scaleprec=obj$scaleprec,
-                  tau2=tau2.arg, delta=delta.arg, steps=obj$steps, verbose=FALSE, control=obj$control,
+                  tau2=tau2.arg, delta=delta.arg, steps=obj$steps, decreasing=obj$decreasing, verbose=FALSE, control=obj$control,
                   skiphes=confint, skiphet=TRUE, defmap=obj$defmap, mapfun=obj$mapfun, mapinvfun=obj$mapinvfun)), silent=TRUE)
 
    }

diff --git a/R/misc.func.hidden.r b/R/misc.func.hidden.r
@@ -1585,15 +1585,21 @@
 
 ############################################################################
 
-.chkopt <- function(optimizer, optcontrol) {
+.chkopt <- function(optimizer, optcontrol, ineq=FALSE) {
 
    mstyle <- .get.mstyle()
 
-   ### set NLOPT_LN_BOBYQA as the default algorithm for nloptr optimizer
+   ### set NLOPT_LN_BOBYQA as the default algorithm for the nloptr optimizer when ineq=FALSE
+   ### and otherwise use NLOPT_LN_COBYLA to allow for nonlinear inequality constraints
    ### and by default use a relative convergence criterion of 1e-8 on the function value
 
-   if (optimizer == "nloptr" && !is.element("algorithm", names(optcontrol)))
-      optcontrol$algorithm <- "NLOPT_LN_BOBYQA"
+   if (optimizer == "nloptr" && !is.element("algorithm", names(optcontrol))) {
+      if (ineq) {
+         optcontrol$algorithm <- "NLOPT_LN_COBYLA"
+      } else {
+         optcontrol$algorithm <- "NLOPT_LN_BOBYQA"
+      }
+   }
 
    if (optimizer == "nloptr" && !is.element("ftol_rel", names(optcontrol)))
       optcontrol$ftol_rel <- 1e-8

diff --git a/R/misc.func.hidden.selmodel.r b/R/misc.func.hidden.selmodel.r
@@ -20,7 +20,11 @@
 
 .mapfun <- function(x, lb, ub, fun=NA) {
    if (is.na(fun)) {
-      lb + (ub-lb) / (1 + exp(-x)) # map (-inf,inf) to (lb,ub)
+      if (lb==0 && ub==1) {
+         plogis(x)
+      } else {
+         lb + (ub-lb) / (1 + exp(-x)) # map (-inf,inf) to (lb,ub)
+      }
    } else {
       x <- sapply(x, fun)
       pmin(pmax(x, lb), ub)
@@ -29,7 +33,11 @@
 
 .mapinvfun <- function(x, lb, ub, fun=NA) {
    if (is.na(fun)) {
-      log((x-lb)/(ub-x))
+      if (lb==0 && ub==1) {
+         qlogis(x)
+      } else {
+         log((x-lb)/(ub-x)) # map (lb,ub) to (-inf,inf)
+      }
    } else {
       sapply(x, fun)
    }
@@ -45,7 +53,7 @@
 }
 
 .selmodel.ll.cont <- function(par, yi, vi, X, preci, k, pX, pvals,
-                              deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max,
+                              deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
                               tau2.arg, tau2.transf, tau2.max, beta.arg,
                               wi.fun, steps, pgrp,
                               alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {
@@ -108,7 +116,7 @@
 ############################################################################
 
 .selmodel.ll.stepfun <- function(par, yi, vi, X, preci, k, pX, pvals,
-                                 deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max,
+                                 deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
                                  tau2.arg, tau2.transf, tau2.max, beta.arg,
                                  wi.fun, steps, pgrp,
                                  alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {
@@ -129,11 +137,26 @@
    tau2[tau2 < .Machine$double.eps*10] <- 0
    tau2[tau2 > tau2.max] <- tau2.max
 
-   if (delta.transf)
-      delta <- mapply(.mapfun, delta, delta.min, delta.max, mapfun)
+   if (decreasing) {
+
+      if (delta.transf) {
+         delta <- exp(delta)
+         delta <- cumsum(c(0, -delta[-1]))
+         delta <- exp(delta)
+      }
+
+   } else {
+
+      if (delta.transf)
+         delta <- mapply(.mapfun, delta, delta.min, delta.max, mapfun)
+
+   }
 
    delta <- ifelse(is.na(delta.arg), delta, delta.arg)
 
+   if (decreasing && any(!is.na(delta.arg[-1])))
+      delta <- rev(cummax(rev(delta)))
+
    yhat <- c(X %*% beta)
 
    sei <- sqrt(vi+tau2)
@@ -206,7 +229,7 @@
 ############################################################################
 
 .selmodel.ll.trunc <- function(par, yi, vi, X, preci, k, pX, pvals,
-                               deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max,
+                               deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
                                tau2.arg, tau2.transf, tau2.max, beta.arg,
                                wi.fun, steps, pgrp,
                                alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {
@@ -271,3 +294,41 @@
 }
 
 ############################################################################
+
+.rma.selmodel.ineqfun.pos <- function(par, yi, vi, X, preci, k, pX, pvals,
+                                      deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
+                                      tau2.arg, tau2.transf, tau2.max, beta.arg,
+                                      wi.fun, steps, pgrp, alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {
+
+   delta <- par[-seq_len(pX+1)]
+
+   if (delta.transf)
+      delta <- mapply(.mapfun, delta, delta.min, delta.max, mapfun)
+
+   delta <- ifelse(is.na(delta.arg), delta, delta.arg)
+
+   diffs <- -diff(delta) # -1 * differences (delta1-delta2, delta2-delta3, ...) must be positive
+
+   return(diffs)
+
+}
+
+.rma.selmodel.ineqfun.neg <- function(par, yi, vi, X, preci, k, pX, pvals,
+                                      deltas, delta.arg, delta.transf, mapfun, delta.min, delta.max, decreasing,
+                                      tau2.arg, tau2.transf, tau2.max, beta.arg,
+                                      wi.fun, steps, pgrp, alternative, pval.min, intCtrl, verbose, digits, dofit=FALSE) {
+
+   delta <- par[-seq_len(pX+1)]
+
+   if (delta.transf)
+      delta <- mapply(.mapfun, delta, delta.min, delta.max, mapfun)
+
+   delta <- ifelse(is.na(delta.arg), delta, delta.arg)
+
+   diffs <- diff(delta) # differences (delta1-delta2, delta2-delta3, ...) must be negative
+
+   return(diffs)
+
+}
+
+############################################################################