Skip to content

Commit

Permalink
Add 'cluster' argument to leave1out() functions.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wviechtb
wviechtb committed Nov 24, 2024
1 parent 00d51a1 commit d0c00ed
Show file tree
Hide file tree
Showing 130 changed files with 534 additions and 356 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: metafor
Version: 4.7-55
Date: 2024-11-22
Version: 4.7-56
Date: 2024-11-24
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
Expand Down
4 changes: 3 additions & 1 deletion NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# metafor 4.7-55 (2024-11-22)
# metafor 4.7-56 (2024-11-24)

- 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

Expand Down Expand Up @@ -26,6 +26,8 @@

- added a `collapse` argument to the various `cumul()` functions (to specify whether studies with the same value of the `order` variable should be added simultaneously)

- the various `leave1out()` functions gain a `cluster` argument

- `rma.mv()` now counts the number of levels of a random effect more appropriately; this may trigger more often the check that 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
Expand Down
116 changes: 79 additions & 37 deletions R/leave1out.rma.mh.r
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
leave1out.rma.mh <- function(x, digits, transf, targs, progbar=FALSE, ...) {
leave1out.rma.mh <- function(x, cluster, digits, transf, targs, progbar=FALSE, ...) {

mstyle <- .get.mstyle()

Expand Down Expand Up @@ -39,43 +39,80 @@ leave1out.rma.mh <- function(x, digits, transf, targs, progbar=FALSE, ...) {

#########################################################################

beta <- rep(NA_real_, x$k.f)
se <- rep(NA_real_, x$k.f)
zval <- rep(NA_real_, x$k.f)
pval <- rep(NA_real_, x$k.f)
ci.lb <- rep(NA_real_, x$k.f)
ci.ub <- rep(NA_real_, x$k.f)
QE <- rep(NA_real_, x$k.f)
QEp <- rep(NA_real_, x$k.f)
#tau2 <- rep(NA_real_, x$k.f)
I2 <- rep(NA_real_, x$k.f)
H2 <- rep(NA_real_, x$k.f)
### process cluster variable

misscluster <- ifelse(missing(cluster), TRUE, FALSE)

if (misscluster) {
cluster <- seq_len(x$k.all)
} else {
mf <- match.call()
cluster <- .getx("cluster", mf=mf, data=x$data)
}

### note: cluster variable must be of the same length as the original dataset
### so we have to apply the same subsetting (if necessary) and removing
### of NAs as was done during model fitting

if (length(cluster) != x$k.all)
stop(mstyle$stop(paste0("Length of the variable specified via 'cluster' (", length(cluster), ") does not match the length of the data (", x$k.all, ").")))

cluster <- .getsubset(cluster, x$subset)

cluster.f <- cluster

cluster <- cluster[x$not.na]

### checks on cluster variable

if (anyNA(cluster.f))
stop(mstyle$stop("No missing values allowed in 'cluster' variable."))

if (length(cluster.f) == 0L)
stop(mstyle$stop(paste0("Cannot find 'cluster' variable (or it has zero length).")))

### cluster ids and number of clusters

ids <- unique(cluster)
n <- length(ids)

if (!misscluster)
ids <- sort(ids)

#########################################################################

beta <- rep(NA_real_, n)
se <- rep(NA_real_, n)
zval <- rep(NA_real_, n)
pval <- rep(NA_real_, n)
ci.lb <- rep(NA_real_, n)
ci.ub <- rep(NA_real_, n)
QE <- rep(NA_real_, n)
QEp <- rep(NA_real_, n)
#tau2 <- rep(NA_real_, n)
I2 <- rep(NA_real_, n)
H2 <- rep(NA_real_, n)

### elements that need to be returned

outlist <- "beta=beta, se=se, zval=zval, pval=pval, ci.lb=ci.lb, ci.ub=ci.ub, QE=QE, QEp=QEp, tau2=tau2, I2=I2, H2=H2"

### note: skipping NA cases

if (progbar)
pbar <- pbapply::startpb(min=0, max=x$k.f)
pbar <- pbapply::startpb(min=0, max=n)

for (i in seq_len(x$k.f)) {
for (i in seq_len(n)) {

if (progbar)
pbapply::setpb(pbar, i)

if (!x$not.na[i])
next

if (is.element(x$measure, c("RR","OR","RD"))) {
args <- list(ai=x$outdat.f$ai, bi=x$outdat.f$bi, ci=x$outdat.f$ci, di=x$outdat.f$di,
measure=x$measure, add=x$add, to=x$to, drop00=x$drop00,
correct=x$correct, level=x$level, subset=-i, outlist=outlist)
args <- list(ai=x$outdat$ai, bi=x$outdat$bi, ci=x$outdat$ci, di=x$outdat$di,
measure=x$measure, add=x$add, to=x$to, drop00=x$drop00, correct=x$correct, level=x$level,
subset=ids[i]!=cluster, outlist=outlist)
} else {
args <- list(x1i=x$outdat.f$x1i, x2i=x$outdat.f$x2i, t1i=x$outdat.f$t1i, t2i=x$outdat.f$t2i,
measure=x$measure, add=x$add, to=x$to, drop00=x$drop00,
correct=x$correct, level=x$level, subset=-i, outlist=outlist)
args <- list(x1i=x$outdat$x1i, x2i=x$outdat$x2i, t1i=x$outdat$t1i, t2i=x$outdat$t2i,
measure=x$measure, add=x$add, to=x$to, drop00=x$drop00, correct=x$correct, level=x$level,
subset=ids[i]!=cluster, outlist=outlist)
}
res <- try(suppressWarnings(.do.call(rma.mh, args)), silent=TRUE)

Expand All @@ -90,9 +127,8 @@ leave1out.rma.mh <- function(x, digits, transf, targs, progbar=FALSE, ...) {
ci.ub[i] <- res$ci.ub
QE[i] <- res$QE
QEp[i] <- res$QEp
#tau2[i] <- res$tau2
I2[i] <- res$I2
H2[i] <- res$H2
I2[i] <- res$I2
H2[i] <- res$H2

}

Expand All @@ -109,14 +145,14 @@ leave1out.rma.mh <- function(x, digits, transf, targs, progbar=FALSE, ...) {
if (is.function(transf)) {
if (is.null(targs)) {
beta <- sapply(beta, transf)
se <- rep(NA_real_, x$k.f)
se <- rep(NA_real_, n)
ci.lb <- sapply(ci.lb, transf)
ci.ub <- sapply(ci.ub, transf)
} else {
if (!is.primitive(transf) && !is.null(targs) && length(formals(transf)) == 1L)
stop(mstyle$stop("Function specified via 'transf' does not appear to have an argument for 'targs'."))
beta <- sapply(beta, transf, targs)
se <- rep(NA_real_, x$k.f)
se <- rep(NA_real_, n)
ci.lb <- sapply(ci.lb, transf, targs)
ci.ub <- sapply(ci.ub, transf, targs)
}
Expand All @@ -131,19 +167,25 @@ leave1out.rma.mh <- function(x, digits, transf, targs, progbar=FALSE, ...) {

#########################################################################

out <- list(estimate=beta, se=se, zval=zval, pval=pval, ci.lb=ci.lb, ci.ub=ci.ub, Q=QE, Qp=QEp, I2=I2, H2=H2)

if (na.act == "na.omit") {
out <- list(estimate=beta[x$not.na], se=se[x$not.na], zval=zval[x$not.na], pval=pval[x$not.na], ci.lb=ci.lb[x$not.na], ci.ub=ci.ub[x$not.na], Q=QE[x$not.na], Qp=QEp[x$not.na], I2=I2[x$not.na], H2=H2[x$not.na])
out$slab <- x$slab[x$not.na]
if (misscluster) {
out$slab <- paste0("-", x$slab[x$not.na])
} else {
out$slab <- paste0("-", ids)
}
}

if (na.act == "na.exclude" || na.act == "na.pass") {
out <- list(estimate=beta, se=se, zval=zval, pval=pval, ci.lb=ci.lb, ci.ub=ci.ub, Q=QE, Qp=QEp, I2=I2, H2=H2)
out$slab <- x$slab
if (misscluster) {
out <- .expandna(out, x$not.na)
out$slab <- paste0("-", x$slab)
} else {
out$slab <- paste0("-", ids)
}
}

if (na.act == "na.fail" && any(!x$not.na))
stop(mstyle$stop("Missing values in results."))

out$digits <- digits
out$transf <- transf

Expand Down
108 changes: 76 additions & 32 deletions R/leave1out.rma.peto.r
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
leave1out.rma.peto <- function(x, digits, transf, targs, progbar=FALSE, ...) {
leave1out.rma.peto <- function(x, cluster, digits, transf, targs, progbar=FALSE, ...) {

mstyle <- .get.mstyle()

Expand Down Expand Up @@ -39,36 +39,75 @@ leave1out.rma.peto <- function(x, digits, transf, targs, progbar=FALSE, ...) {

#########################################################################

beta <- rep(NA_real_, x$k.f)
se <- rep(NA_real_, x$k.f)
zval <- rep(NA_real_, x$k.f)
pval <- rep(NA_real_, x$k.f)
ci.lb <- rep(NA_real_, x$k.f)
ci.ub <- rep(NA_real_, x$k.f)
QE <- rep(NA_real_, x$k.f)
QEp <- rep(NA_real_, x$k.f)
#tau2 <- rep(NA_real_, x$k.f)
I2 <- rep(NA_real_, x$k.f)
H2 <- rep(NA_real_, x$k.f)
### process cluster variable

misscluster <- ifelse(missing(cluster), TRUE, FALSE)

if (misscluster) {
cluster <- seq_len(x$k.all)
} else {
mf <- match.call()
cluster <- .getx("cluster", mf=mf, data=x$data)
}

### note: cluster variable must be of the same length as the original dataset
### so we have to apply the same subsetting (if necessary) and removing
### of NAs as was done during model fitting

if (length(cluster) != x$k.all)
stop(mstyle$stop(paste0("Length of the variable specified via 'cluster' (", length(cluster), ") does not match the length of the data (", x$k.all, ").")))

cluster <- .getsubset(cluster, x$subset)

cluster.f <- cluster

cluster <- cluster[x$not.na]

### checks on cluster variable

if (anyNA(cluster.f))
stop(mstyle$stop("No missing values allowed in 'cluster' variable."))

if (length(cluster.f) == 0L)
stop(mstyle$stop(paste0("Cannot find 'cluster' variable (or it has zero length).")))

### cluster ids and number of clusters

ids <- unique(cluster)
n <- length(ids)

if (!misscluster)
ids <- sort(ids)

#########################################################################

beta <- rep(NA_real_, n)
se <- rep(NA_real_, n)
zval <- rep(NA_real_, n)
pval <- rep(NA_real_, n)
ci.lb <- rep(NA_real_, n)
ci.ub <- rep(NA_real_, n)
QE <- rep(NA_real_, n)
QEp <- rep(NA_real_, n)
#tau2 <- rep(NA_real_, n)
I2 <- rep(NA_real_, n)
H2 <- rep(NA_real_, n)

### elements that need to be returned

outlist <- "beta=beta, se=se, zval=zval, pval=pval, ci.lb=ci.lb, ci.ub=ci.ub, QE=QE, QEp=QEp, tau2=tau2, I2=I2, H2=H2"

### note: skipping NA cases

if (progbar)
pbar <- pbapply::startpb(min=0, max=x$k.f)
pbar <- pbapply::startpb(min=0, max=n)

for (i in seq_len(x$k.f)) {
for (i in seq_len(n)) {

if (progbar)
pbapply::setpb(pbar, i)

if (!x$not.na[i])
next

args <- list(ai=x$outdat.f$ai, bi=x$outdat.f$bi, ci=x$outdat.f$ci, di=x$outdat.f$di, add=x$add, to=x$to, drop00=x$drop00, level=x$level, subset=-i, outlist=outlist)
args <- list(ai=x$outdat$ai, bi=x$outdat$bi, ci=x$outdat$ci, di=x$outdat$di,
add=x$add, to=x$to, drop00=x$drop00, level=x$level,
subset=ids[i]!=cluster, outlist=outlist)
res <- try(suppressWarnings(.do.call(rma.peto, args)), silent=TRUE)

if (inherits(res, "try-error"))
Expand All @@ -82,9 +121,8 @@ leave1out.rma.peto <- function(x, digits, transf, targs, progbar=FALSE, ...) {
ci.ub[i] <- res$ci.ub
QE[i] <- res$QE
QEp[i] <- res$QEp
#tau2[i] <- res$tau2
I2[i] <- res$I2
H2[i] <- res$H2
I2[i] <- res$I2
H2[i] <- res$H2

}

Expand All @@ -101,14 +139,14 @@ leave1out.rma.peto <- function(x, digits, transf, targs, progbar=FALSE, ...) {
if (is.function(transf)) {
if (is.null(targs)) {
beta <- sapply(beta, transf)
se <- rep(NA_real_, x$k.f)
se <- rep(NA_real_, n)
ci.lb <- sapply(ci.lb, transf)
ci.ub <- sapply(ci.ub, transf)
} else {
if (!is.primitive(transf) && !is.null(targs) && length(formals(transf)) == 1L)
stop(mstyle$stop("Function specified via 'transf' does not appear to have an argument for 'targs'."))
beta <- sapply(beta, transf, targs)
se <- rep(NA_real_, x$k.f)
se <- rep(NA_real_, n)
ci.lb <- sapply(ci.lb, transf, targs)
ci.ub <- sapply(ci.ub, transf, targs)
}
Expand All @@ -123,19 +161,25 @@ leave1out.rma.peto <- function(x, digits, transf, targs, progbar=FALSE, ...) {

#########################################################################

out <- list(estimate=beta, se=se, zval=zval, pval=pval, ci.lb=ci.lb, ci.ub=ci.ub, Q=QE, Qp=QEp, I2=I2, H2=H2)

if (na.act == "na.omit") {
out <- list(estimate=beta[x$not.na], se=se[x$not.na], zval=zval[x$not.na], pval=pval[x$not.na], ci.lb=ci.lb[x$not.na], ci.ub=ci.ub[x$not.na], Q=QE[x$not.na], Qp=QEp[x$not.na], I2=I2[x$not.na], H2=H2[x$not.na])
out$slab <- x$slab[x$not.na]
if (misscluster) {
out$slab <- paste0("-", x$slab[x$not.na])
} else {
out$slab <- paste0("-", ids)
}
}

if (na.act == "na.exclude" || na.act == "na.pass") {
out <- list(estimate=beta, se=se, zval=zval, pval=pval, ci.lb=ci.lb, ci.ub=ci.ub, Q=QE, Qp=QEp, I2=I2, H2=H2)
out$slab <- x$slab
if (misscluster) {
out <- .expandna(out, x$not.na)
out$slab <- paste0("-", x$slab)
} else {
out$slab <- paste0("-", ids)
}
}

if (na.act == "na.fail" && any(!x$not.na))
stop(mstyle$stop("Missing values in results."))

out$digits <- digits
out$transf <- transf

Expand Down
Loading

0 comments on commit d0c00ed

Please sign in to comment.