Nothing
#' harmonic mean chi-squared test
#'
#' p-values and confidence intervals from the harmonic mean chi-squared test.
#' @rdname hMeanChiSq
#' @param z Numeric vector of z-values.
#' @param w Numeric vector of weights.
#' @param alternative Either "greater" (default), "less", "two.sided", or "none".
#' Specifies the alternative to be considered in the computation of the p-value.
#' @param bound If \code{FALSE} (default), p-values that cannot be computed are reported as \code{NaN}.
#' If \code{TRUE}, they are reported as "> bound".
#' @return \code{hMeanChiSq}: returns the p-values from the harmonic mean chi-squared test
#' based on the study-specific z-values.
#' @references Held, L. (2020). The harmonic mean chi-squared test to substantiate scientific findings.
#' \emph{Journal of the Royal Statistical Society: Series C (Applied Statistics)}, \bold{69}, 697-708.
#' \doi{10.1111/rssc.12410}
#' @author Leonhard Held, Florian Gerber
#' @examples
#' ## Example from Fisher (1999) as discussed in Held (2020)
#' pvalues <- c(0.0245, 0.1305, 0.00025, 0.2575, 0.128)
#' lower <- c(0.04, 0.21, 0.12, 0.07, 0.41)
#' upper <- c(1.14, 1.54, 0.60, 3.75, 1.27)
#' se <- ci2se(lower = lower, upper = upper, ratio = TRUE)
#' thetahat <- ci2estimate(lower = lower, upper = upper, ratio = TRUE)
#'
#' ## hMeanChiSq() --------
#' hMeanChiSq(z = p2z(p = pvalues, alternative = "less"),
#' alternative = "less")
#' hMeanChiSq(z = p2z(p = pvalues, alternative = "less"),
#' alternative = "two.sided")
#' hMeanChiSq(z = p2z(p = pvalues, alternative = "less"),
#' alternative = "none")
#'
#' hMeanChiSq(z = p2z(p = pvalues, alternative = "less"),
#' w = 1 / se^2, alternative = "less")
#' hMeanChiSq(z = p2z(p = pvalues, alternative = "less"),
#' w = 1 / se^2, alternative = "two.sided")
#' hMeanChiSq(z = p2z(p = pvalues, alternative = "less"),
#' w = 1 / se^2, alternative = "none")
#' @export
hMeanChiSq <- function(z,
w = rep(1, length(z)),
alternative = c("greater", "less", "two.sided", "none"),
bound = FALSE) {
stopifnot(is.numeric(z),
length(z) > 0,
is.finite(z),
is.numeric(w),
length(w) == length(z),
is.finite(w),
min(w) > 0,
!is.null(alternative))
alternative <- match.arg(alternative)
stopifnot(is.logical(bound),
length(bound) == 1,
is.finite(bound))
n <- length(z)
zH2 <- sum(sqrt(w))^2 / sum(w / z^2)
res <- stats::pchisq(zH2, df = 1, lower.tail = FALSE)
check_greater <- min(z) >= 0
check_less <- max(z) <= 0
break_p <- 1 / (2^n)
if (alternative == "greater") {
if (bound)
res <- if (check_greater) res / (2^n) else paste(">", format(break_p, scientific = FALSE))
else
res <- if (check_greater) res / (2^n) else NaN
}
if (alternative == "less") {
if (bound)
res <- if (check_less) res / (2^n) else paste(">", format(break_p, scientific = FALSE))
else
res <- if (check_less) res / (2^n) else NaN
}
if (alternative == "two.sided") {
if (bound)
res <- if (check_greater || check_less) res / (2^(n - 1)) else
paste(">", format(2 * break_p, scientific = FALSE))
else
res <- if (check_greater || check_less) res / (2^(n - 1)) else NaN
}
## no chnage to res for alternative == "none"
return(res)
}
#' @rdname hMeanChiSq
#' @param thetahat Numeric vector of parameter estimates.
#' @param se Numeric vector of standard errors.
#' @param mu The null hypothesis value. Defaults to 0.
#' @return \code{hMeanChiSqMu}: returns the p-value from the harmonic mean chi-squared test
#' based on study-specific estimates and standard errors.
#' @examples
#'
#'
#' ## hMeanChiSqMu() --------
#' hMeanChiSqMu(thetahat = thetahat, se = se, alternative = "two.sided")
#' hMeanChiSqMu(thetahat = thetahat, se = se, w = 1 / se^2,
#' alternative = "two.sided")
#' hMeanChiSqMu(thetahat = thetahat, se = se, alternative = "two.sided",
#' mu = -0.1)
#' @export
hMeanChiSqMu <- function(
thetahat, se, w = rep(1, length(thetahat)), mu = 0,
alternative = c("greater", "less", "two.sided", "none"),
bound = FALSE) {
stopifnot(is.numeric(thetahat),
length(thetahat) > 0,
is.finite(thetahat),
is.numeric(se),
length(se) == 1 || length(se) == length(thetahat),
is.finite(se),
min(se) > 0,
is.numeric(w),
length(w) == length(thetahat),
is.finite(w),
min(w) > 0,
is.numeric(mu),
length(mu) > 0,
is.finite(mu),
!is.null(alternative))
alternative <- match.arg(alternative)
stopifnot(is.logical(bound),
length(bound) == 1,
is.finite(bound))
n <- length(thetahat)
m <- length(mu)
if (alternative != "none") {
z <- (thetahat - mu) / se
zH2 <- sum(sqrt(w))^2 / sum(w / z^2)
res <- stats::pchisq(zH2, df = 1, lower.tail = FALSE)
check_greater <- min(z) >= 0
check_less <- max(z) <= 0
break_p <- 1 / (2^n)
if (alternative == "greater") {
if (bound)
res <- if (check_greater) res / (2^n) else paste(">", format(break_p, scientific = FALSE))
else
res <- if (check_greater || check_less) res/(2^n) else NaN
}
if (alternative == "less") {
if (bound)
res <- if (check_less) res / (2^n) else paste(">", format(break_p, scientific = FALSE))
else
res <- if (check_greater || check_less) res / (2^n) else NaN
}
if (alternative == "two.sided") {
if (bound)
res <- if (check_greater || check_less) res / (2^(n - 1)) else
paste(">", format(2*break_p, scientific = FALSE))
else
res <- if (check_greater || check_less) res / (2^(n - 1)) else NaN
}
}
if (alternative == "none") {
zH2 <- numeric(m)
sw <- sum(sqrt(w))^2
for (i in seq_along(m)) {
z <- (thetahat - mu[i]) / se
zH2[i] <- sw / sum(w / z^2)
}
res <- stats::pchisq(q = zH2, df = 1, lower.tail = FALSE)
}
return(res)
}
#' @rdname hMeanChiSq
#' @param conf.level Numeric vector specifying the conf.level of the confidence interval. Defaults to 0.95.
#' summarize the gamma values, i.e.,
#' the local minima of the p-value function between the thetahats. Defaults is a vector of 1s.
#' @return \code{hMeanChiSqCI}: returns a list containing confidence interval(s)
#' obtained by inverting the harmonic mean chi-squared test based on study-specific
#' estimates and standard errors. The list contains:
#' \item{CI}{Confidence interval(s).}\cr\cr
#' If the \code{alternative} is "none", the list also contains:
#' \item{gamma}{Local minima of the p-value function between the thetahats.}
#' @examples
#'
#' ## hMeanChiSqCI() --------
#' ## two-sided
#' CI1 <- hMeanChiSqCI(thetahat = thetahat, se = se, w = 1 / se^2,
#' alternative = "two.sided")
#' CI2 <- hMeanChiSqCI(thetahat = thetahat, se = se, w = 1 / se^2,
#' alternative = "two.sided", conf.level = 0.99875)
#' ## one-sided
#' CI1b <- hMeanChiSqCI(thetahat = thetahat, se = se, w = 1 / se^2,
#' alternative = "less", conf.level = 0.975)
#' CI2b <- hMeanChiSqCI(thetahat = thetahat, se = se, w = 1 / se^2,
#' alternative = "less", conf.level = 1 - 0.025^2)
#'
#' ## confidence intervals on hazard ratio scale
#' print(exp(CI1$CI), digits = 2)
#' print(exp(CI2$CI), digits = 2)
#' print(exp(CI1b$CI), digits = 2)
#' print(exp(CI2b$CI), digits = 2)
#'
#'
#' ## example with confidence region consisting of disjunct intervals
#' thetahat2 <- c(-3.7, 2.1, 2.5)
#' se2 <- c(1.5, 2.2, 3.1)
#' conf.level <- 0.95; alpha <- 1 - conf.level
#' muSeq <- seq(-7, 6, length.out = 1000)
#' pValueSeq <- hMeanChiSqMu(thetahat = thetahat2, se = se2,
#' alternative = "none", mu = muSeq)
#' (hm <- hMeanChiSqCI(thetahat = thetahat2, se = se2, alternative = "none"))
#'
#' plot(x = muSeq, y = pValueSeq, type = "l", panel.first = grid(lty = 1),
#' xlab = expression(mu), ylab = "p-value")
#' abline(v = thetahat2, h = alpha, lty = 2)
#' arrows(x0 = hm$CI[, 1], x1 = hm$CI[, 2], y0 = alpha,
#' y1 = alpha, col = "darkgreen", lwd = 3, angle = 90, code = 3)
#' points(hm$gamma, col = "red", pch = 19, cex = 2)
#'
#' @export
hMeanChiSqCI <- function(
thetahat, se, w = rep(1, length(thetahat)),
alternative = c("two.sided", "greater", "less", "none"),
conf.level = 0.95) {
stopifnot(is.numeric(thetahat),
length(thetahat) > 0,
is.finite(thetahat))
stopifnot(is.numeric(se),
length(se) == 1 || length(se) == length(thetahat),
is.finite(se),
min(se) > 0,
is.numeric(w),
length(w) == length(thetahat),
is.finite(w),
min(w) > 0,
!is.null(alternative))
alternative <- match.arg(alternative)
stopifnot(is.numeric(conf.level),
length(conf.level) == 1,
is.finite(conf.level),
0 < conf.level, conf.level < 1,
is.finite(w),
min(w) > 0)
## target function to compute the limits of the CI
target <- function(limit) {
hMeanChiSqMu(thetahat = thetahat, se = se, w = w, mu = limit,
alternative = alternative, bound = FALSE) - alpha
}
## sort 'thetahat', 'se', 'w'
indOrd <- order(thetahat)
thetahat <- thetahat[indOrd]
se <- se[indOrd]
w <- w[indOrd]
## minima are only search between distinct thetahat elements
thetahatUnique <- unique(thetahat)
nThetahatUnique <- length(thetahatUnique)
mini <- which.min(thetahat)
maxi <- which.max(thetahat)
mint <- thetahat[mini]
maxt <- thetahat[maxi]
minse <- se[mini]
maxse <- se[maxi]
alpha <- 1 - conf.level
z1 <- max(-stats::qnorm(alpha), 1)
eps <- 1e-6
factor <- 5
if (alternative == "none") {
## ----------------------------
## find lower bound such that: lower < thetahat[1] AND target(lower) < 0
lower <- mint - z1 * minse
while (target(lower) > 0) lower <- lower - minse
## find root between 'lower' and 'thetahat[1]'
CIlower <- stats::uniroot(
f = target, lower = lower, upper = thetahat[1]
)$root
## -------------------------
## check between thetahats whether 'target' goes below 'alpha'
## if so, search CI limits
CImiddle <- matrix(NA, nrow = 2, ncol = nThetahatUnique - 1)
gam <- matrix(NA, nrow = nThetahatUnique - 1, ncol = 2)
colnames(gam) <- c("minimum", "pvalue_fun/gamma")
for (i in 1:(nThetahatUnique - 1)) {
opt <- stats::optimize(
f = target, lower = thetahatUnique[i],
upper = thetahatUnique[i + 1]
)
gam[i, ] <- c(opt$minimum, opt$objective + alpha)
if (opt$objective <= 0) {
CImiddle[1, i] <- stats::uniroot(
f = target, lower = thetahatUnique[i],
upper = opt$minimum
)$root
CImiddle[2, i] <- stats::uniroot(
f = target, lower = opt$minimum,
upper = thetahatUnique[i + 1]
)$root
}
}
CImiddle <- CImiddle[!is.na(CImiddle)]
## -------------------------
## find upper bound such that:
## upper > thetahat[length(thetahat)] AND target(upper) < 0
upper <- maxt + maxse
while (target(upper) > 0) upper <- upper + z1 * maxse
## find root between 'lower' and 'thetahat[1]'
CIupper <- stats::uniroot(
f = target, lower = thetahat[length(thetahat)],
upper = upper
)$root
CI <- matrix(c(CIlower, CImiddle, CIupper), ncol = 2, byrow = TRUE)
colnames(CI) <- c("lower", "upper")
return(list(CI = CI, gamma = gam))
}
if (alternative == "two.sided") {
lower <- stats::uniroot(f = target,
lower = mint - factor * z1 * minse,
upper = mint - eps * minse)$root
upper <- stats::uniroot(f = target, lower = maxt + eps * maxse,
upper = maxt + factor * z1 * maxse)$root
return(list(CI = cbind(lower, upper)))
}
if (alternative == "greater") {
lower <- stats::uniroot(f = target,
lower = mint - factor * z1 * minse,
upper = mint - eps * minse)$root
upper <- Inf
return(list(CI = cbind(lower, upper)))
}
if (alternative == "less") {
lower <- -Inf
upper <- stats::uniroot(f = target,
lower = maxt + eps * maxse,
upper = maxt + factor * z1 * maxse)$root
return(list(CI = cbind(lower, upper)))
}
stop("function not get here.")
}
#' Find multiple roots in interval
#'
#' Searches the interval from lower to upper for several roots
#' (i.e., zero's) of a univariate function \code{f}.
#' @param f the function for which the root is sought.
#' \code{f} should be vectorized in the first argument.
#' @param interval A vector containing the end-points of the interval to be
#' searched for the root.
#' @param n Number of subintervals on which \code{link[stats]{uniroot}} is called.
#' Default is 1000.
#' @param lower The lower end point of the interval to be searched.
#' @param upper The upper end point of the interval to be searched.
#' @param tol See help of \code{link[stats]{uniroot}}.
#' @param maxiter See help of \code{link[stats]{uniroot}}.
#' @param trace See help of \code{link[stats]{uniroot}}.
#' @param ... Additional named or unnamed arguments to be passed to \code{f}.
#' @return A numeric vector of the roots found in the interval.
#' @references This function is inspired by rootSolve::uniroot.all(),
#' package version 1.8.2.2.
#' @author Florian Gerber
#' @noRd
#' @seealso \code{\link[base]{Vectorize}}
#' @examples
#' f <- function (x) cos(2*x)^3
#' (roots <- unirootAll(f = f, interval = c(0, 10)))
#' f(roots)
unirootAll <- function(
f, interval, lower = min(interval), upper = max(interval),
n = 1000,
tol = .Machine$double.eps^0.2,
maxiter = 1000, trace = 0,
...) {
stopifnot(is.function(f), length(formals(f)) >= 1,
is.numeric(lower), length(lower) == 1, is.finite(lower),
is.numeric(upper), length(upper) == 1, is.finite(upper),
lower < upper,
is.numeric(n), length(n) == 1, is.finite(n), n >= 2)
## 'tol', 'maxiter', 'trace' are checked in uniroot()
n <- round(n)
xseq <- seq(lower, upper, length.out = n + 1)
mod <- f(xseq, ...)
index <- which(mod[1:n] * mod[2:(n + 1)] < 0)
rootsOffGrid <- numeric(length = length(index))
for (i in seq_along(index)) {
rootsOffGrid[i] <- stats::uniroot(
f = f, lower = xseq[index[i]],
upper = xseq[index[i] + 1], maxiter = maxiter,
tol = tol, trace = trace, ...
)$root
}
rootsGrid <- xseq[which(mod == 0)]
if (length(rootsGrid) > 0) return(sort(c(rootsGrid, rootsOffGrid)))
rootsOffGrid
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.