Further improvements to predstyle.

wviechtb · Sep 9, 2024 · 80f55ad · 80f55ad
1 parent 0603761
commit 80f55ad
Show file tree

Hide file tree

Showing 138 changed files with 517 additions and 249 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: metafor
-Version: 4.7-29
-Date: 2024-09-04
+Version: 4.7-30
+Date: 2024-09-09
 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.7-29 (2024-09-04)
+# metafor 4.7-30 (2024-09-09)
 
 - some general changes to the various `forest()` functions: argument `header` is now `TRUE` by default, the y-axis is now created with `yaxs="i"`, and the y-axis limits have been tweaked slightly in accordance
 

diff --git a/R/addpoly.default.r b/R/addpoly.default.r
@@ -1,5 +1,5 @@
 addpoly.default     <- function(x, vi, sei, ci.lb, ci.ub, pi.lb, pi.ub,
-rows=-1, level,         annotate,                predstyle, digits, width, mlab,
+rows=-1, level,         annotate,                predstyle, predlim, digits, width, mlab,
 transf, atransf, targs, efac, col, border, lty, fonts, cex, constarea=FALSE, ...) {
 
    #########################################################################
@@ -514,6 +514,22 @@ transf, atransf, targs, efac, col, border, lty, fonts, cex, constarea=FALSE, ...
 
       if (predstyle %in% c("shade","dist")) {
 
+         if (is.null(pi.se))
+            stop(mstyle$stop("Cannot extract SE of the prediction interval."))
+
+         if (!missing(predlim) && length(predlim) != 2L)
+            stop(mstyle$stop("Argument 'predlim' must be of length 2."))
+
+         if (is.function(transf)) {
+            funlist <- lapply(list("1"=exp, "2"=transf.ztor, "3"=tanh, "4"=transf.ilogit, "5"=plogis, "6"=transf.iarcsin), deparse)
+            funmatch <- sapply(funlist, identical, transf.char)
+            if (!any(funmatch))
+               stop(mstyle$stop("Chosen transformation not (currently) possible with this 'predstyle'."))
+         }
+
+         if (pi.dist != "norm" && pi.ddf <= 1L)
+            stop(mstyle$stop("Cannot shade/draw prediction distribution when df <= 1."))
+
          if (predstyle == "shade") {
             x.len <- 100
             q.lo <- pi.level/2
@@ -524,58 +540,69 @@ transf, atransf, targs, efac, col, border, lty, fonts, cex, constarea=FALSE, ...
             q.hi <- 0.9999
          }
 
-         if (is.null(pi.se))
-            stop(mstyle$stop("Cannot extract SE of the prediction interval."))
-
-         if (pi.dist == "norm") {
-            crits <- qnorm(c(q.lo,q.hi), mean=yi.utransf[i], sd=pi.se)
-            xs <- seq(crits[1], crits[2], length.out=x.len)
-            ys <- dnorm(xs, mean=yi.utransf[i], sd=pi.se)
+         if (missing(predlim) || predstyle == "shade") {
+            if (pi.dist == "norm") {
+               crits <- qnorm(c(q.lo,q.hi), mean=yi.utransf[i], sd=pi.se)
+               xs <- seq(crits[1], crits[2], length.out=x.len)
+               ys <- dnorm(xs, mean=yi.utransf[i], sd=pi.se)
+            } else {
+               crits <- qt(c(q.lo,q.hi), df=pi.ddf) * pi.se + yi.utransf[i]
+               xs <- seq(crits[1], crits[2], length.out=x.len)
+               ys <- dt((xs - yi.utransf[i]) / pi.se, df=pi.ddf) / pi.se
+            }
          } else {
-            if (pi.ddf <= 1L)
-               stop(mstyle$stop("Cannot draw prediction distribution when df <= 1."))
-            crits <- qt(c(q.lo,q.hi), df=pi.ddf) * pi.se + yi.utransf[i]
-            xs <- seq(crits[1], crits[2], length.out=x.len)
-            ys <- dt((xs - yi.utransf[i]) / pi.se, df=pi.ddf) / pi.se
+            xs <- seq(predlim[1], predlim[2], length.out=x.len)
+            if (is.function(transf)) {
+               if (funmatch[1])
+                  xs <- suppressWarnings(log(xs))
+               if (any(funmatch[2:3]))
+                  xs <- suppressWarnings(atanh(xs))
+               if (any(funmatch[4:5]))
+                  xs <- suppressWarnings(qlogis(xs))
+               if (funmatch[6])
+                  xs <- suppressWarnings(transf.arcsin(xs))
+               sel <- is.finite(xs) # FALSE for +-Inf and NA/NaN
+               x.len <- sum(sel)
+               xs <- xs[sel]
+            }
+            if (pi.dist == "norm") {
+               ys <- dnorm(xs, mean=yi.utransf[i], sd=pi.se)
+            } else {
+               ys <- dt((xs - yi.utransf[i]) / pi.se, df=pi.ddf) / pi.se
+            }
          }
 
          sel.l0 <- xs < 0
          sel.g0 <- xs > 0
 
          if (is.function(transf)) {
             xs <- sapply(xs, transf)
-            funlist <- lapply(list("1"=exp, "2"=transf.ztor, "3"=tanh, "4"=transf.ilogit, "5"=plogis, "6"=transf.iarcsin), deparse)
-            funmatch <- sapply(funlist, identical, transf.char)
-            if (any(funmatch)) {
-               if (funmatch[1])
-                  ys <- ys / xs
-               if (any(funmatch[2:3]))
-                  ys <- ys / (1-xs^2)
-               if (any(funmatch[4:5]))
-                  ys <- ys / (xs*(1-xs))
-               if (funmatch[6])
-                  ys <- ys / (2*sqrt(xs*(1-xs)))
-               # can run into numerical problems at the extremes (esp. when pi.dist="t"),
-               # so we need to find the range where ys are first increasing and then decreasing
-               sdys <- sign(diff(ys))
-               y.lo <- which(sdys == 1)
-               if (length(y.lo) == 0L) {
-                  y.lo <- 1
-               } else {
-                  y.lo <- min(y.lo) + 1
-               }
-               y.hi <- which(sdys == -1)
-               if (length(y.hi) == 0L) {
-                  y.hi <- x.len
-               } else {
-                  y.hi <- max(y.hi) + 1
-               }
-               sel.l0 <- sel.l0[y.lo:y.hi]
-               sel.g0 <- sel.g0[y.lo:y.hi]
-               xs <- xs[y.lo:y.hi]
-               ys <- ys[y.lo:y.hi]
-            } else {
-               stop(mstyle$stop("Chosen transformation not (currently) possible with this 'predstyle'."))
+            if (funmatch[1]) {
+               ys <- ys / xs
+               x.lo <- 0.01
+               x.hi <- Inf
+            }
+            if (any(funmatch[2:3])) {
+               ys <- ys / (1-xs^2)
+               x.lo <- -0.99
+               x.hi <-  0.99
+            }
+            if (any(funmatch[4:5])) {
+               ys <- ys / (xs*(1-xs))
+               x.lo <- 0.01
+               x.hi <- 0.99
+            }
+            if (funmatch[6]) {
+               ys <- ys / (2*sqrt(xs*(1-xs)))
+               x.lo <- 0.01
+               x.hi <- 0.99
+            }
+            if (missing(predlim)) {
+               sel <- xs > x.lo & xs < x.hi
+               sel.l0 <- sel.l0[sel]
+               sel.g0 <- sel.g0[sel]
+               ys <- ys[sel]
+               xs <- xs[sel]
             }
          }
 
@@ -587,7 +614,7 @@ transf, atransf, targs, efac, col, border, lty, fonts, cex, constarea=FALSE, ...
 
          colfun <- colorRamp(c(col[2], col[3]))
          rectcol <- colfun(intensity)
-         rectcol <- apply(rectcol, 1, function(x) rgb(x[1], x[2], x[3], maxColorValue=255))
+         rectcol <- apply(rectcol, 1, function(x) if (anyNA(x)) NA else rgb(x[1], x[2], x[3], maxColorValue=255))
 
          lrect(xs[-1], rows[2]-barheight, xs[-length(xs)], rows[2]+barheight, col=rectcol, border=rectcol, ...)
 
@@ -597,7 +624,11 @@ transf, atransf, targs, efac, col, border, lty, fonts, cex, constarea=FALSE, ...
 
          ys <- ys / max(ys) * efac[2]
 
-         sel <- ys > 0.005
+         if (missing(predlim)) {
+            sel <- ys > 0.005
+         } else {
+            sel <- rep(TRUE, length(ys))
+         }
 
          xs.sel.l0 <- xs[sel.l0 & sel]
          xs.sel.g0 <- xs[sel.g0 & sel]

diff --git a/R/forest.cumul.rma.r b/R/forest.cumul.rma.r
@@ -501,6 +501,9 @@ lty, fonts, cex, cex.lab, cex.axis, ...) {
       xlim <- round(xlim, digits[[2]])
       #xlim[1] <- xlim[1]*max(1, digits[[2]]/2)
       #xlim[2] <- xlim[2]*max(1, digits[[2]]/2)
+   } else {
+      if (length(xlim) != 2L)
+         stop(mstyle$stop("Argument 'xlim' must be of length 2."))
    }
 
    xlim <- sort(xlim)

diff --git a/R/forest.default.r b/R/forest.default.r
@@ -645,6 +645,9 @@ lty, fonts, cex, cex.lab, cex.axis, ...) {
       xlim <- round(xlim, digits[[2]])
       #xlim[1] <- xlim[1]*max(1, digits[[2]]/2)
       #xlim[2] <- xlim[2]*max(1, digits[[2]]/2)
+   } else {
+      if (length(xlim) != 2L)
+         stop(mstyle$stop("Argument 'xlim' must be of length 2."))
    }
 
    xlim <- sort(xlim)

diff --git a/R/forest.rma.r b/R/forest.rma.r
@@ -1,6 +1,6 @@
 forest.rma       <- function(x,
 annotate=TRUE, addfit=TRUE, addpred=FALSE, predstyle="line", showweights=FALSE, header=TRUE,
-xlim, alim, olim, ylim, at, steps=5, level=x$level, refline=0, digits=2L, width,
+xlim, alim, olim, ylim, predlim, at, steps=5, level=x$level, refline=0, digits=2L, width,
 xlab, slab, mlab, ilab, ilab.lab, ilab.xpos, ilab.pos, order,
 transf, atransf, targs, rows,
 efac=1, pch, psize, plim=c(0.5,1.5), colout, col, border, shade, colshade,
@@ -160,6 +160,9 @@ lty, fonts, cex, cex.lab, cex.axis, ...) {
       if (missing(border))
          border <- .coladj(par("bg","fg"), dark=0.3, light=-0.3) # border color of the fitted value polygons
 
+      if (predstyle %in% c("bar","shade","dist"))
+         warning(mstyle$warning("Argument 'predstyle' not relevant for meta-regression models."), call.=FALSE)
+
    }
 
    ### set default line types if user has not specified 'lty' argument
@@ -767,6 +770,9 @@ lty, fonts, cex, cex.lab, cex.axis, ...) {
       xlim <- round(xlim, digits[[2]])
       #xlim[1] <- xlim[1]*max(1, digits[[2]]/2)
       #xlim[2] <- xlim[2]*max(1, digits[[2]]/2)
+   } else {
+      if (length(xlim) != 2L)
+         stop(mstyle$stop("Argument 'xlim' must be of length 2."))
    }
 
    xlim <- sort(xlim)
@@ -1063,6 +1069,19 @@ lty, fonts, cex, cex.lab, cex.axis, ...) {
 
          if (predstyle %in% c("shade","dist")) {
 
+            if (!missing(predlim) && length(predlim) != 2L)
+               stop(mstyle$stop("Argument 'predlim' must be of length 2."))
+
+            if (is.function(transf)) {
+               funlist <- lapply(list("1"=exp, "2"=transf.ztor, "3"=tanh, "4"=transf.ilogit, "5"=plogis, "6"=transf.iarcsin), deparse)
+               funmatch <- sapply(funlist, identical, transf.char)
+               if (!any(funmatch))
+                  stop(mstyle$stop("Chosen transformation not (currently) possible with this 'predstyle'."))
+            }
+
+            if (predres$pi.dist != "norm" && predres$pi.ddf <= 1L)
+               stop(mstyle$stop("Cannot shade/draw prediction distribution when df <= 1."))
+
             if (predstyle == "shade") {
                x.len <- 100
                q.lo <- level/2
@@ -1073,55 +1092,69 @@ lty, fonts, cex, cex.lab, cex.axis, ...) {
                q.hi <- 0.9999
             }
 
-            if (predres$pi.dist == "norm") {
-               crits <- qnorm(c(q.lo,q.hi), mean=predres$pred, sd=predres$pi.se)
-               xs <- seq(crits[1], crits[2], length.out=x.len)
-               ys <- dnorm(xs, mean=predres$pred, sd=predres$pi.se)
+            if (missing(predlim) || predstyle == "shade") {
+               if (predres$pi.dist == "norm") {
+                  crits <- qnorm(c(q.lo,q.hi), mean=predres$pred, sd=predres$pi.se)
+                  xs <- seq(crits[1], crits[2], length.out=x.len)
+                  ys <- dnorm(xs, mean=predres$pred, sd=predres$pi.se)
+               } else {
+                  crits <- qt(c(q.lo,q.hi), df=predres$pi.ddf) * predres$pi.se + predres$pred
+                  xs <- seq(crits[1], crits[2], length.out=x.len)
+                  ys <- dt((xs - predres$pred) / predres$pi.se, df=predres$pi.ddf) / predres$pi.se
+               }
             } else {
-               if (predres$pi.ddf <= 1L)
-                  stop(mstyle$stop("Cannot draw prediction distribution when df <= 1."))
-               crits <- qt(c(q.lo,q.hi), df=predres$pi.ddf) * predres$pi.se + predres$pred
-               xs <- seq(crits[1], crits[2], length.out=x.len)
-               ys <- dt((xs - predres$pred) / predres$pi.se, df=predres$pi.ddf) / predres$pi.se
+               xs <- seq(predlim[1], predlim[2], length.out=x.len)
+               if (is.function(transf)) {
+                  if (funmatch[1])
+                     xs <- suppressWarnings(log(xs))
+                  if (any(funmatch[2:3]))
+                     xs <- suppressWarnings(atanh(xs))
+                  if (any(funmatch[4:5]))
+                     xs <- suppressWarnings(qlogis(xs))
+                  if (funmatch[6])
+                     xs <- suppressWarnings(transf.arcsin(xs))
+                  sel <- is.finite(xs) # FALSE for +-Inf and NA/NaN
+                  x.len <- sum(sel)
+                  xs <- xs[sel]
+               }
+               if (predres$pi.dist == "norm") {
+                  ys <- dnorm(xs, mean=predres$pred, sd=predres$pi.se)
+               } else {
+                  ys <- dt((xs - predres$pred) / predres$pi.se, df=predres$pi.ddf) / predres$pi.se
+               }
             }
 
             sel.l0 <- xs < 0
             sel.g0 <- xs > 0
 
             if (is.function(transf)) {
                xs <- sapply(xs, transf)
-               funlist <- lapply(list("1"=exp, "2"=transf.ztor, "3"=tanh, "4"=transf.ilogit, "5"=plogis, "6"=transf.iarcsin), deparse)
-               funmatch <- sapply(funlist, identical, transf.char)
-               if (any(funmatch)) {
-                  if (funmatch[1])
-                     ys <- ys / xs
-                  if (any(funmatch[2:3]))
-                     ys <- ys / (1-xs^2)
-                  if (any(funmatch[4:5]))
-                     ys <- ys / (xs*(1-xs))
-                  if (funmatch[6])
-                     ys <- ys / (2*sqrt(xs*(1-xs)))
-                  # can run into numerical problems at the extremes (esp. when pi.dist="t"),
-                  # so we need to find the range where ys are first increasing and then decreasing
-                  sdys <- sign(diff(ys))
-                  y.lo <- which(sdys == 1)
-                  if (length(y.lo) == 0L) {
-                     y.lo <- 1
-                  } else {
-                     y.lo <- min(y.lo) + 1
-                  }
-                  y.hi <- which(sdys == -1)
-                  if (length(y.hi) == 0L) {
-                     y.hi <- x.len
-                  } else {
-                     y.hi <- max(y.hi) + 1
-                  }
-                  sel.l0 <- sel.l0[y.lo:y.hi]
-                  sel.g0 <- sel.g0[y.lo:y.hi]
-                  xs <- xs[y.lo:y.hi]
-                  ys <- ys[y.lo:y.hi]
-               } else {
-                  stop(mstyle$stop("Chosen transformation not (currently) possible with this 'predstyle'."))
+               if (funmatch[1]) {
+                  ys <- ys / xs
+                  x.lo <- 0.01
+                  x.hi <- Inf
+               }
+               if (any(funmatch[2:3])) {
+                  ys <- ys / (1-xs^2)
+                  x.lo <- -0.99
+                  x.hi <-  0.99
+               }
+               if (any(funmatch[4:5])) {
+                  ys <- ys / (xs*(1-xs))
+                  x.lo <- 0.01
+                  x.hi <- 0.99
+               }
+               if (funmatch[6]) {
+                  ys <- ys / (2*sqrt(xs*(1-xs)))
+                  x.lo <- 0.01
+                  x.hi <- 0.99
+               }
+               if (missing(predlim)) {
+                  sel <- xs > x.lo & xs < x.hi
+                  sel.l0 <- sel.l0[sel]
+                  sel.g0 <- sel.g0[sel]
+                  ys <- ys[sel]
+                  xs <- xs[sel]
                }
             }
 
@@ -1142,7 +1175,7 @@ lty, fonts, cex, cex.lab, cex.axis, ...) {
 
             colfun <- colorRamp(c(col[2], col[3]))
             rectcol <- colfun(intensity)
-            rectcol <- apply(rectcol, 1, function(x) rgb(x[1], x[2], x[3], maxColorValue=255))
+            rectcol <- apply(rectcol, 1, function(x) if (anyNA(x)) NA else rgb(x[1], x[2], x[3], maxColorValue=255))
 
             lrect(xs[-1], -2-barheight, xs[-length(xs)], -2+barheight, col=rectcol, border=rectcol, ...)
 
@@ -1152,7 +1185,12 @@ lty, fonts, cex, cex.lab, cex.axis, ...) {
 
             ys <- ys / max(ys) * efac[4]
 
-            sel <- ys > 0.005
+            if (missing(predlim)) {
+               sel <- ys > 0.005
+            } else {
+               sel <- rep(TRUE, length(ys))
+            }
+
             sel <- sel & xs >= alim[1] & xs <= alim[2]
 
             if (!missing(olim))

diff --git a/R/zzz.r b/R/zzz.r
@@ -1,6 +1,6 @@
 .onAttach <- function(libname, pkgname) {
 
-   ver <- "4.7-29"
+   ver <- "4.7-30"
 
    loadmsg <- paste0("\nLoading the 'metafor' package (version ", ver, "). For an\nintroduction to the package please type: help(metafor)\n")
 

diff --git a/README.md b/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--29-brightgreen.svg)](https://www.metafor-project.org/doku.php/installation#development_version)
+[![devel Version](https://img.shields.io/badge/devel-4.7--30-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)