https://github.com/cran/metafor
Tip revision: 5e080e28c25477c201792f082375d718b86ed758 authored by Wolfgang Viechtbauer on 28 March 2024, 16:00:06 UTC
version 4.6-0
version 4.6-0
Tip revision: 5e080e2
confint.rma.mv.r
confint.rma.mv <- function(object, parm, level, fixed=FALSE, sigma2, tau2, rho, gamma2, phi, digits, transf, targs, verbose=FALSE, control, ...) {
mstyle <- .get.mstyle()
.chkclass(class(object), must="rma.mv")
if (!missing(parm))
warning(mstyle$warning("Argument 'parm' (currently) ignored."), call.=FALSE)
x <- object
k <- x$k
p <- x$p
if (missing(level))
level <- x$level
if (missing(digits)) {
digits <- .get.digits(xdigits=x$digits, dmiss=TRUE)
} else {
digits <- .get.digits(digits=digits, xdigits=x$digits, dmiss=FALSE)
}
if (missing(transf))
transf <- FALSE
if (missing(targs))
targs <- NULL
if (missing(control))
control <- list()
ddd <- list(...)
.chkdots(ddd, c("time", "xlim", "extint"))
level <- .level(level, stopon100=.isTRUE(ddd$extint))
if (.isTRUE(ddd$time))
time.start <- proc.time()
if (!is.null(ddd$xlim)) {
if (length(ddd$xlim) != 2L)
stop(mstyle$stop("Argument 'xlim' should be a vector of length 2."))
control$vc.min <- ddd$xlim[1]
control$vc.max <- ddd$xlim[2]
}
### check if user has specified one of the sigma2, tau2, rho, gamma2, or phi arguments
random <- !all(missing(sigma2), missing(tau2), missing(rho), missing(gamma2), missing(phi))
if (!fixed && !random) {
### if both 'fixed' and 'random' are FALSE, obtain CIs for all variance/correlation components
cl <- match.call()
### total number of non-fixed components
comps <- ifelse(x$withS, sum(!x$vc.fix$sigma2), 0) +
ifelse(x$withG, sum(!x$vc.fix$tau2) + sum(!x$vc.fix$rho), 0) +
ifelse(x$withH, sum(!x$vc.fix$gamma2) + sum(!x$vc.fix$phi), 0)
if (comps == 0)
stop(mstyle$stop("No components for which a CI can be obtained."))
res.all <- list()
j <- 0
if (x$withS && any(!x$vc.fix$sigma2)) {
for (pos in seq_len(x$sigma2s)[!x$vc.fix$sigma2]) {
j <- j + 1
cl.vc <- cl
cl.vc$sigma2 <- pos
cl.vc$time <- FALSE
#cl.vc$object <- quote(x)
if (verbose)
cat(mstyle$verbose(paste("\nObtaining CI for sigma2 =", pos, "\n")))
res.all[[j]] <- eval(cl.vc, envir=parent.frame())
}
}
if (x$withG) {
if (any(!x$vc.fix$tau2)) {
for (pos in seq_len(x$tau2s)[!x$vc.fix$tau2]) {
j <- j + 1
cl.vc <- cl
cl.vc$tau2 <- pos
cl.vc$time <- FALSE
#cl.vc$object <- quote(x)
if (verbose)
cat(mstyle$verbose(paste("\nObtaining CI for tau2 =", pos, "\n")))
res.all[[j]] <- eval(cl.vc, envir=parent.frame())
}
}
if (any(!x$vc.fix$rho)) {
for (pos in seq_len(x$rhos)[!x$vc.fix$rho]) {
j <- j + 1
cl.vc <- cl
cl.vc$rho <- pos
cl.vc$time <- FALSE
#cl.vc$object <- quote(x)
if (verbose)
cat(mstyle$verbose(paste("\nObtaining CI for rho =", pos, "\n")))
res.all[[j]] <- eval(cl.vc, envir=parent.frame())
}
}
}
if (x$withH) {
if (any(!x$vc.fix$gamma2)) {
for (pos in seq_len(x$gamma2s)[!x$vc.fix$gamma2]) {
j <- j + 1
cl.vc <- cl
cl.vc$gamma2 <- pos
cl.vc$time <- FALSE
#cl.vc$object <- quote(x)
if (verbose)
cat(mstyle$verbose(paste("\nObtaining CI for gamma2 =", pos, "\n")))
res.all[[j]] <- eval(cl.vc, envir=parent.frame())
}
}
if (any(!x$vc.fix$phi)) {
for (pos in seq_len(x$phis)[!x$vc.fix$phi]) {
j <- j + 1
cl.vc <- cl
cl.vc$phi <- pos
cl.vc$time <- FALSE
#cl.vc$object <- quote(x)
if (verbose)
cat(mstyle$verbose(paste("\nObtaining CI for phi =", pos, "\n")))
res.all[[j]] <- eval(cl.vc, envir=parent.frame())
}
}
}
if (.isTRUE(ddd$time)) {
time.end <- proc.time()
.print.time(unname(time.end - time.start)[3])
}
if (length(res.all) == 1L) {
return(res.all[[1]])
} else {
res.all$digits <- digits
class(res.all) <- "list.confint.rma"
return(res.all)
}
}
#########################################################################
#########################################################################
#########################################################################
if (random) {
type <- "pl"
######################################################################
### check if user has specified more than one of these arguments
if (sum(!missing(sigma2), !missing(tau2), !missing(rho), !missing(gamma2), !missing(phi)) > 1L)
stop(mstyle$stop("Must specify only one of the arguments 'sigma2', 'tau2', 'rho', 'gamma2', or 'phi'."))
### check if model actually contains (at least one) such a component and that it was actually estimated
### note: a component that is not in the model is NA; components that are fixed are TRUE
if (!missing(sigma2) && (all(is.na(x$vc.fix$sigma2)) || all(x$vc.fix$sigma2)))
stop(mstyle$stop("Model does not contain any (estimated) 'sigma2' components."))
if (!missing(tau2) && (all(is.na(x$vc.fix$tau2)) || all(x$vc.fix$tau2)))
stop(mstyle$stop("Model does not contain any (estimated) 'tau2' components."))
if (!missing(rho) && c(all(is.na(x$vc.fix$rho)) || all(x$vc.fix$rho)))
stop(mstyle$stop("Model does not contain any (estimated) 'rho' components."))
if (!missing(gamma2) && (all(is.na(x$vc.fix$gamma2)) || all(x$vc.fix$gamma2)))
stop(mstyle$stop("Model does not contain any (estimated) 'gamma2' components."))
if (!missing(phi) && c(all(is.na(x$vc.fix$phi)) || all(x$vc.fix$phi)))
stop(mstyle$stop("Model does not contain any (estimated) 'phi' components."))
### check if user specified more than one sigma2, tau2, rho, gamma2, or rho component
if (!missing(sigma2) && (length(sigma2) > 1L))
stop(mstyle$stop("Can only specify one 'sigma2' component."))
if (!missing(tau2) && (length(tau2) > 1L))
stop(mstyle$stop("Can only specify one 'tau2' component."))
if (!missing(rho) && (length(rho) > 1L))
stop(mstyle$stop("Can only specify one 'rho' component."))
if (!missing(gamma2) && (length(gamma2) > 1L))
stop(mstyle$stop("Can only specify one 'gamma2' component."))
if (!missing(phi) && (length(phi) > 1L))
stop(mstyle$stop("Can only specify one 'phi' component."))
### check if user specified a logical
if (!missing(sigma2) && is.logical(sigma2))
stop(mstyle$stop("Must specify a number for the 'sigma2' component."))
if (!missing(tau2) && is.logical(tau2))
stop(mstyle$stop("Must specify a number for the 'tau2' component."))
if (!missing(rho) && is.logical(rho))
stop(mstyle$stop("Must specify a number for the 'rho' component."))
if (!missing(gamma2) && is.logical(gamma2))
stop(mstyle$stop("Must specify a number for the 'gamma2' component."))
if (!missing(phi) && is.logical(phi))
stop(mstyle$stop("Must specify a number for the 'phi' component."))
### check if user specified a component that does not exist
if (!missing(sigma2) && (sigma2 > length(x$vc.fix$sigma2) || sigma2 <= 0))
stop(mstyle$stop("No such 'sigma2' component in the model."))
if (!missing(tau2) && (tau2 > length(x$vc.fix$tau2) || tau2 <= 0))
stop(mstyle$stop("No such 'tau2' component in the model."))
if (!missing(rho) && (rho > length(x$vc.fix$rho) || rho <= 0))
stop(mstyle$stop("No such 'rho' component in the model."))
if (!missing(gamma2) && (gamma2 > length(x$vc.fix$gamma2) || gamma2 <= 0))
stop(mstyle$stop("No such 'gamma2' component in the model."))
if (!missing(phi) && (phi > length(x$vc.fix$phi) || phi <= 0))
stop(mstyle$stop("No such 'phi' component in the model."))
### check if user specified a component that was fixed
if (!missing(sigma2) && x$vc.fix$sigma2[sigma2])
stop(mstyle$stop("Specified 'sigma2' component was fixed."))
if (!missing(tau2) && x$vc.fix$tau2[tau2])
stop(mstyle$stop("Specified 'tau2' component was fixed."))
if (!missing(rho) && x$vc.fix$rho[rho])
stop(mstyle$stop("Specified 'rho' component was fixed."))
if (!missing(gamma2) && x$vc.fix$gamma2[gamma2])
stop(mstyle$stop("Specified 'gamma2' component was fixed."))
if (!missing(phi) && x$vc.fix$phi[phi])
stop(mstyle$stop("Specified 'phi' component was fixed."))
### if everything is good so far, get value of the variance component and set 'comp'
sigma2.pos <- NA_integer_
tau2.pos <- NA_integer_
rho.pos <- NA_integer_
gamma2.pos <- NA_integer_
phi.pos <- NA_integer_
if (!missing(sigma2)) {
vc <- x$sigma2[sigma2]
comp <- "sigma2"
sigma2.pos <- sigma2
}
if (!missing(tau2)) {
vc <- x$tau2[tau2]
comp <- "tau2"
tau2.pos <- tau2
}
if (!missing(rho)) {
vc <- x$rho[rho]
comp <- "rho"
rho.pos <- rho
}
if (!missing(gamma2)) {
vc <- x$gamma2[gamma2]
comp <- "gamma2"
gamma2.pos <- gamma2
}
if (!missing(phi)) {
vc <- x$phi[phi]
comp <- "phi"
phi.pos <- phi
}
#return(list(comp=comp, vc=vc, sigma2.pos=sigma2.pos, tau2.pos=tau2.pos, rho.pos=rho.pos, gamma2.pos=gamma2.pos, phi.pos=phi.pos))
######################################################################
### 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
con <- list(tol=.Machine$double.eps^0.25, maxiter=1000, verbose=FALSE, eptries=10)
if (is.element(comp, c("sigma2", "tau2", "gamma2"))) {
con$vc.min <- 0
con$vc.max <- max(ifelse(vc <= .Machine$double.eps^0.5, 10, max(10, vc*100)), con$vc.min)
}
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
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)
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.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.max <- max(10, vc*10)
}
if (is.element(x$struct[1], c("PHYPL","PHYPD"))) {
con$vc.min <- 0
con$vc.max <- max(2, vc*2)
}
}
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
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)
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.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.max <- max(10, vc*10)
}
if (is.element(x$struct[2], c("PHYPL","PHYPD"))) {
con$vc.min <- 0
con$vc.max <- max(2, vc*2)
}
}
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
######################################################################
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 (">")
######################################################################
######################################################################
######################################################################
### Profile Likelihood method
if (type == "pl") {
if (con$vc.min > vc)
stop(mstyle$stop("Lower bound of interval to be searched must be <= estimated value of component."))
if (con$vc.max < vc)
stop(mstyle$stop("Upper bound of interval to be searched must be >= estimated value of component."))
objective <- qchisq(1-level, df=1)
###################################################################
### search for lower bound
### get diff value when setting component to vc.min; this value should be positive (i.e., discrepancy must be larger than critical value)
### if it is not, then the lower bound must be below vc.min
epdiff <- abs(con$vc.min - vc) / con$eptries
for (i in seq_len(con$eptries)) {
res <- try(.profile.rma.mv(con$vc.min, obj=x, comp=comp, sigma2.pos=sigma2.pos, tau2.pos=tau2.pos, rho.pos=rho.pos, gamma2.pos=gamma2.pos, phi.pos=phi.pos, confint=TRUE, objective=objective, verbose=verbose), silent=TRUE)
if (!inherits(res, "try-error") && !is.na(res)) {
if (!.isTRUE(ddd$extint) && res < 0) {
vc.lb <- con$vc.min
lb.conv <- TRUE
if (is.element(comp, c("sigma2", "tau2", "gamma2")) && con$vc.min > 0)
lb.sign <- "<"
if (is.element(comp, c("rho", "phi")) && con$vc.min > -1)
lb.sign <- "<"
if (((comp == "rho" && is.element(x$struct[1], c("SPEXP","SPGAU","SPLIN","SPRAT","SPSPH","PHYPL","PHYPD"))) || (comp == "phi" && is.element(x$struct[2], c("SPEXP","SPGAU","SPLIN","SPRAT","SPSPH","PHYPL","PHYPD")))) && con$vc.min > 0)
lb.sign <- "<"
} else {
if (.isTRUE(ddd$extint)) {
res <- try(uniroot(.profile.rma.mv, interval=c(con$vc.min, vc), tol=con$tol, maxiter=con$maxiter, extendInt="downX", obj=x, comp=comp, sigma2.pos=sigma2.pos, tau2.pos=tau2.pos, rho.pos=rho.pos, gamma2.pos=gamma2.pos, phi.pos=phi.pos, confint=TRUE, objective=objective, verbose=verbose, check.conv=TRUE)$root, silent=TRUE)
} else {
res <- try(uniroot(.profile.rma.mv, interval=c(con$vc.min, vc), tol=con$tol, maxiter=con$maxiter, obj=x, comp=comp, sigma2.pos=sigma2.pos, tau2.pos=tau2.pos, rho.pos=rho.pos, gamma2.pos=gamma2.pos, phi.pos=phi.pos, confint=TRUE, objective=objective, verbose=verbose, check.conv=TRUE)$root, silent=TRUE)
}
### check if uniroot method converged
if (!inherits(res, "try-error")) {
vc.lb <- res
lb.conv <- TRUE
}
}
break
}
con$vc.min <- con$vc.min + epdiff
}
if (verbose)
cat("\n")
###################################################################
### search for upper bound
### get diff value when setting component to vc.max; this value should be positive (i.e., discrepancy must be larger than critical value)
### if it is not, then the upper bound must be above vc.max
epdiff <- abs(con$vc.max - vc) / con$eptries
for (i in seq_len(con$eptries)) {
res <- try(.profile.rma.mv(con$vc.max, obj=x, comp=comp, sigma2.pos=sigma2.pos, tau2.pos=tau2.pos, rho.pos=rho.pos, gamma2.pos=gamma2.pos, phi.pos=phi.pos, confint=TRUE, objective=objective, verbose=verbose), silent=TRUE)
if (!inherits(res, "try-error") && !is.na(res)) {
if (!.isTRUE(ddd$extint) && res < 0) {
vc.ub <- con$vc.max
ub.conv <- TRUE
if (is.element(comp, c("sigma2", "tau2", "gamma2")))
ub.sign <- ">"
if (is.element(comp, c("rho", "phi")) && con$vc.max < 1)
ub.sign <- ">"
if ((comp == "rho" && is.element(x$struct[1], c("SPEXP","SPGAU","SPLIN","SPRAT","SPSPH","PHYPL","PHYPD"))) || (comp == "phi" && is.element(x$struct[2], c("SPEXP","SPGAU","SPLIN","SPRAT","SPSPH","PHYPL","PHYPD"))))
ub.sign <- ">"
} else {
if (.isTRUE(ddd$extint)) {
res <- try(uniroot(.profile.rma.mv, interval=c(vc, con$vc.max), tol=con$tol, maxiter=con$maxiter, extendInt="upX", obj=x, comp=comp, sigma2.pos=sigma2.pos, tau2.pos=tau2.pos, rho.pos=rho.pos, gamma2.pos=gamma2.pos, phi.pos=phi.pos, confint=TRUE, objective=objective, verbose=verbose, check.conv=TRUE)$root, silent=TRUE)
} else {
res <- try(uniroot(.profile.rma.mv, interval=c(vc, con$vc.max), tol=con$tol, maxiter=con$maxiter, obj=x, comp=comp, sigma2.pos=sigma2.pos, tau2.pos=tau2.pos, rho.pos=rho.pos, gamma2.pos=gamma2.pos, phi.pos=phi.pos, confint=TRUE, objective=objective, verbose=verbose, check.conv=TRUE)$root, silent=TRUE)
}
### check if uniroot method converged
if (!inherits(res, "try-error")) {
vc.ub <- res
ub.conv <- TRUE
}
}
break
}
con$vc.max <- con$vc.max - epdiff
}
###################################################################
}
######################################################################
######################################################################
######################################################################
if (!lb.conv)
warning(mstyle$warning("Cannot obtain lower bound of profile likelihood CI due to convergence problems."), call.=FALSE)
if (!ub.conv)
warning(mstyle$warning("Cannot obtain upper bound of profile likelihood CI due to convergence problems."), call.=FALSE)
######################################################################
vc <- c(vc, vc.lb, vc.ub)
if (is.element(comp, c("sigma2", "tau2", "gamma2"))) {
vcsqrt <- sqrt(ifelse(vc >= 0, vc, NA_real_))
res.random <- rbind(vc, vcsqrt)
if (comp == "sigma2") {
if (length(x$sigma2) == 1L) {
rownames(res.random) <- c("sigma^2", "sigma")
} else {
rownames(res.random) <- paste0(c("sigma^2", "sigma"), ".", sigma2.pos)
}
}
if (comp == "tau2") {
if (length(x$tau2) == 1L) {
rownames(res.random) <- c("tau^2", "tau")
} else {
rownames(res.random) <- paste0(c("tau^2", "tau"), ".", tau2.pos)
}
}
if (comp == "gamma2") {
if (length(x$gamma2) == 1L) {
rownames(res.random) <- c("gamma^2", "gamma")
} else {
rownames(res.random) <- paste0(c("gamma^2", "gamma"), ".", gamma2.pos)
}
}
} else {
res.random <- rbind(vc)
if (comp == "rho") {
if (length(x$rho) == 1L) {
rownames(res.random) <- "rho"
} else {
rownames(res.random) <- paste0("rho.", rho.pos)
}
}
if (comp == "phi") {
if (length(x$phi) == 1L) {
rownames(res.random) <- "phi"
} else {
rownames(res.random) <- paste0("phi.", rho.pos)
}
}
}
colnames(res.random) <- c("estimate", "ci.lb", "ci.ub")
}
#########################################################################
#########################################################################
#########################################################################
if (fixed) {
if (is.element(x$test, c("knha","adhoc","t"))) {
crit <- sapply(seq_along(x$ddf), function(j) if (x$ddf[j] > 0) qt(level/2, df=x$ddf[j], lower.tail=FALSE) else NA_real_)
} else {
crit <- qnorm(level/2, lower.tail=FALSE)
}
beta <- c(x$beta)
ci.lb <- c(beta - crit * x$se)
ci.ub <- c(beta + crit * x$se)
if (is.function(transf)) {
if (is.null(targs)) {
beta <- sapply(beta, transf)
ci.lb <- sapply(ci.lb, transf)
ci.ub <- sapply(ci.ub, transf)
} else {
beta <- sapply(beta, transf, targs)
ci.lb <- sapply(ci.lb, transf, targs)
ci.ub <- sapply(ci.ub, transf, targs)
}
}
### make sure order of intervals is always increasing
tmp <- .psort(ci.lb, ci.ub)
ci.lb <- tmp[,1]
ci.ub <- tmp[,2]
res.fixed <- cbind(estimate=beta, ci.lb=ci.lb, ci.ub=ci.ub)
rownames(res.fixed) <- rownames(x$beta)
}
#########################################################################
#########################################################################
#########################################################################
res <- list()
if (fixed)
res$fixed <- res.fixed
if (random)
res$random <- res.random
res$digits <- digits
if (random) {
res$ci.null <- ci.null
res$lb.sign <- lb.sign
res$ub.sign <- ub.sign
#res$vc.min <- con$vc.min
}
if (.isTRUE(ddd$time)) {
time.end <- proc.time()
.print.time(unname(time.end - time.start)[3])
}
class(res) <- "confint.rma"
return(res)
}