Added pairwise() and allow p-dimensional 'newmods' matrix.

wviechtb · Jun 10, 2024 · 0b9744c · 0b9744c
1 parent feacbe6
commit 0b9744c
Show file tree

Hide file tree

Showing 149 changed files with 1,113 additions and 324 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: metafor
-Version: 4.7-14
-Date: 2024-06-07
+Version: 4.7-15
+Date: 2024-06-10
 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

diff --git a/NEWS.md b/NEWS.md
@@ -1,12 +1,16 @@
-# metafor 4.7-14 (2024-06-07)
+# metafor 4.7-15 (2024-06-10)
+
+- the `predict.rma()` and `predict.rma.ls()` functions now also accept a matrix as input that includes a column for the intercept term (in which case the `intercept` argument is ignored)
+
+- added `pairwise()` function to construct a matrix of pairwise contrasts
 
 - `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
 
-- argument `shade` in `funnel()` now automatically uses a color gradient for the regions when multiple `level` values are specified
+- when unspecified, argument `shade` in `funnel()` now automatically uses a color gradient for the regions when multiple `level` values are specified
 
-- added extractor for function `se()` for extracting standard errors from model objects
+- added extractor function `se()` for extracting standard errors from model objects
 
 - added `lim`, `ci`, `pi`, `legend`, and `flip` arguments to `labbe()`
 

diff --git a/R/anova.rma.r b/R/anova.rma.r
@@ -12,7 +12,7 @@ anova.rma <- function(object, object2, btt, X, att, Z, rhs, digits, refit=FALSE,
 
    ddd <- list(...)
 
-   .chkdots(ddd, c("test", "L", "verbose", "fixed", "df"))
+   .chkdots(ddd, c("test", "L", "verbose", "fixed", "df", "abbrev"))
 
    if (!is.null(ddd$L))
       X <- ddd$L
@@ -25,9 +25,12 @@ anova.rma <- function(object, object2, btt, X, att, Z, rhs, digits, refit=FALSE,
    if (!missing(Z) && !inherits(object, "rma.ls"))
       stop(mstyle$stop("Can only specify 'Z' for location-scale models."))
 
-   #mf <- match.call()
-   #if (any(grepl("pairwise(", as.character(mf), fixed=TRUE)))
-   #   try(assign("pairwise", object, envir=.metafor), silent=TRUE)
+   mf <- match.call()
+
+   if (any(grepl("pairwise(", as.character(mf), fixed=TRUE))) {
+      try(assign("pairwise", object, envir=.metafor), silent=TRUE)
+      on.exit(suppressWarnings(rm("pairwise", envir=.metafor)))
+   }
 
    if (missing(object2)) {
 
@@ -256,6 +259,13 @@ anova.rma <- function(object, object2, btt, X, att, Z, rhs, digits, refit=FALSE,
             colnames(hyp) <- ""
             rownames(hyp) <- paste0(seq_len(m), ":") # add '1:', '2:', ... as row names
 
+            ### abbreviate some hyp elements
+
+            if (.isTRUE(ddd$abbrev)) {
+               hyp[,1] <- gsub("factor(", "", hyp[,1], fixed=TRUE)
+               hyp[,1] <- gsub(")", "", hyp[,1], fixed=TRUE)
+            }
+
             res <- list(QS=QS, QSdf=QSdf, QSp=QSp, hyp=hyp, Za=Za, se=se, zval=zval, pval=pval, k=x$k, q=x$q, m=m, test=x$test, ddf=x$ddf.alpha, digits=digits, type="Wald.Za")
 
          } else {
@@ -415,6 +425,13 @@ anova.rma <- function(object, object2, btt, X, att, Z, rhs, digits, refit=FALSE,
             colnames(hyp) <- ""
             rownames(hyp) <- paste0(seq_len(m), ":") # add '1:', '2:', ... as row names
 
+            ### abbreviate some hyp elements
+
+            if (.isTRUE(ddd$abbrev)) {
+               hyp[,1] <- gsub("factor(", "", hyp[,1], fixed=TRUE)
+               hyp[,1] <- gsub(")", "", hyp[,1], fixed=TRUE)
+            }
+
             res <- list(QM=QM, QMdf=QMdf, QMp=QMp, hyp=hyp, Xb=Xb, se=se, zval=zval, pval=pval, k=x$k, p=x$p, m=m, test=x$test, ddf=ddf, digits=digits, type="Wald.Xb")
 
          }

diff --git a/R/misc.func.hidden.r b/R/misc.func.hidden.r
@@ -1920,7 +1920,7 @@
 
    }
 
-   if (themeopt != "default")
+   if (themeopt != "default" && isFALSE(par("new")))
       par(fg=fg, bg=bg, col=fg, col.axis=fg, col.lab=fg, col.main=fg, col.sub=fg)
 
    invisible()

diff --git a/R/pairwise.r b/R/pairwise.r
@@ -0,0 +1,85 @@
+pairwise <- function(x, btt, btt2, ...) {
+
+   mstyle <- .get.mstyle()
+
+   if (missing(x)) {
+      x <- .getfromenv("pairwise", envir=.metafor)
+   }# else {
+   #   if (is.atomic(x)) {
+   #      btt <- x
+   #      x <- .getfromenv("pairwise", envir=.metafor)
+   #   }
+   #}
+
+   if (is.null(x))
+      stop(mstyle$stop("Need to specify 'x' argument."), call.=FALSE)
+
+   .chkclass(class(x), must="rma")
+
+   if (x$int.only)
+      stop(mstyle$stop("Cannot construct contrast matrices for intercept-only models."))
+
+   if (missing(btt) || is.null(btt))
+      stop(mstyle$stop("Need to specify 'btt' argument."), call.=FALSE)
+
+   ddd <- list(...)
+
+   .chkdots(ddd, c("fixed"))
+
+   fixed <- .chkddd(ddd$fixed, FALSE, .isTRUE(ddd$fixed))
+
+   #########################################################################
+
+   btt <- .set.btt(btt, x$p, x$int.incl, colnames(x$X), fixed=fixed)
+
+   p <- length(btt)
+
+   if (p == 1L)
+      stop(mstyle$stop("Need to specify multiple coefficients via argument 'btt' for pairwise comparisons."), call.=FALSE)
+
+   names <- rownames(x$beta)
+   connames <- rep("", p*(p-1)/2)
+
+   X <- matrix(0, nrow=p*(p-1)/2, ncol=x$p)
+   row <- 0
+
+   for (i in 1:(p-1)) {
+      btti <- btt[i]
+      for (j in (i+1):p) {
+         bttj <- btt[j]
+         row <- row + 1
+         X[row,btti] <- -1
+         X[row,bttj] <- +1
+         connames[row] <- paste0(names[btti], "-", names[bttj])
+      }
+   }
+
+   rownames(X) <- connames
+
+   #########################################################################
+
+   ### in case btt2 is specified, add these coefficients to X
+
+   if (!missing(btt2)) {
+
+      btt <- .set.btt(btt2, x$p, x$int.incl, colnames(x$X), fixed=fixed)
+
+      p <- length(btt)
+
+      Xadd <- matrix(0, nrow=p, ncol=x$p)
+
+      for (i in 1:p) {
+         Xadd[i,btt[i]] <- 1
+      }
+
+      rownames(Xadd) <- names[btt]
+
+      X <- rbind(Xadd, X)
+
+   }
+
+   #########################################################################
+
+   return(X)
+
+}
diff --git a/R/predict.rma.ls.r b/R/predict.rma.ls.r
@@ -14,11 +14,22 @@ level, digits, transf, targs, vcov=FALSE, ...) {
 
    x <- object
 
+   mf <- match.call()
+
+   if (any(grepl("pairwise(", as.character(mf), fixed=TRUE))) {
+      try(assign("pairwise", object, envir=.metafor), silent=TRUE)
+      on.exit(suppressWarnings(rm("pairwise", envir=.metafor)))
+   }
+
    if (missing(newmods))
       newmods <- NULL
 
-   if (missing(intercept))
+   if (missing(intercept)) {
       intercept <- x$intercept
+      int.spec <- FALSE
+   } else {
+      int.spec <- TRUE
+   }
 
    if (missing(newscale))
       newscale <- NULL
@@ -63,8 +74,10 @@ level, digits, transf, targs, vcov=FALSE, ...) {
       if (!(.is.vector(newmods) || inherits(newmods, "matrix")))
          stop(mstyle$stop(paste0("Argument 'newmods' should be a vector or matrix, but is of class '", class(newmods), "'.")))
 
-      if ((!x$int.incl && x$p == 1L) || (x$int.incl && x$p == 2L)) {
-         k.new <- length(newmods)                                # if single moderator (multiple k.new possible) (either without or with intercept in the model)
+      singlemod <- (NCOL(newmods) == 1L) && ((!x$int.incl && x$p == 1L) || (x$int.incl && x$p == 2L))
+
+      if (singlemod) {                                           # if single moderator (multiple k.new possible) (either without or with intercept in the model)
+         k.new <- length(newmods)                                # (but when specifying a matrix, it must be a column vector for this work)
          X.new <- cbind(c(newmods))                              #
          if (.is.vector(newmods)) {                              #
             rnames <- names(newmods)                             #
@@ -121,12 +134,23 @@ level, digits, transf, targs, vcov=FALSE, ...) {
       ### but user can also decide to remove the intercept from the predictions with intercept=FALSE
       ### one special case: when the location model is an intercept-only model, one can set newmods=1 to obtain the predicted intercept
 
-      if (x$int.incl && !(x$int.only && ncol(X.new) == 1L && nrow(X.new) == 1L && X.new[1,1] == 1)) {
-         if (intercept) {
-            X.new <- cbind(intrcpt=1, X.new)
-         } else {
-            X.new <- cbind(intrcpt=0, X.new)
+      if (inherits(newmods, "matrix") && ncol(newmods) == x$p) {
+
+         if (int.spec)
+            warning(mstyle$warning("Arguments 'intercept' ignored when 'newmods' includes 'p' columns."), call.=FALSE)
+
+      } else {
+
+         if (x$int.incl && !(x$int.only && ncol(X.new) == 1L && nrow(X.new) == 1L && X.new[1,1] == 1)) {
+
+            if (intercept) {
+               X.new <- cbind(intrcpt=1, X.new)
+            } else {
+               X.new <- cbind(intrcpt=0, X.new)
+            }
+
          }
+
       }
 
       if (ncol(X.new) != x$p)
@@ -139,8 +163,10 @@ level, digits, transf, targs, vcov=FALSE, ...) {
       if (!(.is.vector(newscale) || inherits(newscale, "matrix")))
          stop(mstyle$stop(paste0("Argument 'newscale' should be a vector or matrix, but is of class '", class(newscale), "'.")))
 
-      if ((!x$Z.int.incl && x$q == 1L) || (x$Z.int.incl && x$q == 2L)) {
-         Z.k.new <- length(newscale)                             # if single moderator (multiple k.new possible) (either without or with intercept in the model)
+      singlescale <- (NCOL(newscale) == 1L) && ((!x$Z.int.incl && x$q == 1L) || (x$Z.int.incl && x$q == 2L))
+
+      if (singlescale) {                                         # if single moderator (multiple k.new possible) (either without or with intercept in the model)
+         Z.k.new <- length(newscale)                             #
          Z.new <- cbind(c(newscale))                             #
          if (is.null(rnames)) {                                  #
             if (.is.vector(newscale)) {                          #
@@ -201,16 +227,25 @@ level, digits, transf, targs, vcov=FALSE, ...) {
       ### one special case: when the scale model is an intercept-only model, one can set newscale=1 to obtain the predicted intercept
       ### (which can be converted to tau^2 with transf=exp when using a log link)
 
-      if (x$Z.int.incl && !(x$Z.int.only && ncol(Z.new) == 1L && nrow(Z.new) == 1L && Z.new[1,1] == 1)) {
-         if (is.null(newmods)) {
-            if (intercept) {
-               Z.new <- cbind(intrcpt=1, Z.new)
+      if (inherits(newscale, "matrix") && ncol(newscale) == x$q) {
+
+         if (int.spec)
+            warning(mstyle$warning("Arguments 'intercept' ignored when 'newscale' includes 'q' columns."), call.=FALSE)
+
+      } else {
+
+         if (x$Z.int.incl && !(x$Z.int.only && ncol(Z.new) == 1L && nrow(Z.new) == 1L && Z.new[1,1] == 1)) {
+            if (is.null(newmods)) {
+               if (intercept) {
+                  Z.new <- cbind(intrcpt=1, Z.new)
+               } else {
+                  Z.new <- cbind(intrcpt=0, Z.new)
+               }
             } else {
-               Z.new <- cbind(intrcpt=0, Z.new)
+               Z.new <- cbind(intrcpt=1, Z.new)
             }
-         } else {
-            Z.new <- cbind(intrcpt=1, Z.new)
          }
+
       }
 
       if (ncol(Z.new) != x$q)
@@ -271,7 +306,12 @@ level, digits, transf, targs, vcov=FALSE, ...) {
             }
          }
 
-         if (length(tau2.f) == 1L) {
+         if (k.new == 1L && Z.k.new > 1L) {
+            X.new <- X.new[rep(1,Z.k.new),,drop=FALSE]
+            k.new <- Z.k.new
+         }
+
+         if (length(tau2.f) == 1L && k.new > 1L) {
             Z.new  <- Z.new[rep(1,k.new),,drop=FALSE]
             tau2.f <- rep(tau2.f, k.new)
          }