powerReplicationSuccess.R

.powerReplicationSuccess_ <- function(
    zo,
    c = 1,
    level = 0.025,
    designPrior = c("conditional", "predictive", "EB"),
    alternative = c("one.sided", "two.sided"),
    type = c("golden", "nominal", "controlled"),
    shrinkage = 0,
    h = 0,
    strict = FALSE) {

    stopifnot(is.numeric(zo),
              length(zo) == 1,
              is.finite(zo),

              is.numeric(c),
              length(c) == 1,
              is.finite(c),
              0 <= c,

              is.numeric(level),
              length(level) == 1,
              is.finite(level),
              0 < level, level < 1,

              !is.null(designPrior))

    designPrior <- match.arg(designPrior)

    stopifnot(!is.null(alternative))
    alternative <- match.arg(alternative)

    stopifnot(!is.null(type))
    type <- match.arg(type)

    stopifnot(is.numeric(shrinkage),
              length(shrinkage) == 1,
              is.finite(shrinkage),
              0 <= shrinkage, shrinkage < 1,

              is.numeric(h),
              length(h) == 1,
              is.finite(h),
              0 <= h,

              is.logical(strict),
              length(strict) == 1)

    ## take the absolute value of zo for easier computations
    zoabs <- abs(zo)

    ## if zoas < zalphaS, power is zero
    alphaS <- levelSceptical(level = level, alternative = alternative, type = type,
                             c = c)
    zalphaS <- p2z(p = alphaS, alternative = alternative)
    if (zoabs < zalphaS) {
        power <- 0
    } else {
        ## computing minimum zr to achieve replication success
        dmin <- effectSizeReplicationSuccess(zo = zoabs, c = c, level = level,
                                             alternative = alternative, type = type)
        zrmin <- dmin * zoabs * sqrt(c)

        if (designPrior == "conditional") {
            power <- stats::pnorm(
                q = zrmin,
                mean = zoabs * (1 - shrinkage) * sqrt(c),
                sd = 1,
                lower.tail = FALSE
            )

            if (strict && alternative == "two.sided") {
                power2 <- stats::pnorm(
                    q = -zrmin,
                    mean = zoabs * (1 - shrinkage) * sqrt(c),
                    sd = 1,
                    lower.tail = TRUE
                )
                power <- power + power2
            }
        } else if (designPrior == "predictive") {
            power <- stats::pnorm(
                q = zrmin,
                mean = zoabs * (1 - shrinkage) * sqrt(c),
                sd = sqrt(c * (1 + 2 * h) + 1),
                lower.tail = FALSE
            )

            if (strict && alternative == "two.sided") {
                power2 <- stats::pnorm(
                    q = -zrmin,
                    mean = zoabs * (1 - shrinkage) * sqrt(c),
                    sd = sqrt(c * (1 + 2 * h) + 1), lower.tail = TRUE
                )
                power <- power + power2
            }
        } else { ## designPrior == "EB"
            EBshrinkage <- pmin((1 + h) / zoabs^2, 1)
            power <- stats::pnorm(
                q = zrmin,
                mean = zoabs * (1 - EBshrinkage) * sqrt(c),
                sd = sqrt((1 - EBshrinkage) * c * (1 + h) + 1 + c * h),
                lower.tail = FALSE
            )

            if (strict && alternative == "two.sided") {
                power2 <- stats::pnorm(
                    q = -zrmin,
                    mean = zoabs * (1 - EBshrinkage) * sqrt(c),
                    sd = sqrt((1 - EBshrinkage) * c * (1 + h) + 1 + c * h),
                    lower.tail = TRUE
                )
                power <- power + power2
            }
        }
    }
    return(power)
}


#' Computes the power for replication success with the sceptical p-value
#'
#' Computes the power for replication success with the sceptical
#' p-value  based on the result of the
#' original study, the corresponding variance ratio, and the design prior.
#' @param zo Numeric vector of z-values from original studies.
#' @param c Numeric vector of variance ratios of the original and replication
#'     effect estimates. This is usually the ratio of the sample size of the
#'     replication study to the sample size of the original study.
#' @param level Threshold for the calibrated sceptical p-value.
#'  Default is 0.025.
#' @param designPrior Either "conditional" (default), "predictive", or "EB". If
#'     "EB", the power is computed under a predictive distribution, where
#'     the contribution of the original study is shrunken towards zero based on
#'     the evidence in the original study (with an empirical Bayes shrinkage
#'     estimator).
#' @param alternative Specifies if \code{level} is "one.sided" (default) or "two.sided".
#'     If "one.sided" then power calculations are based
#'     on a one-sided assessment of replication success in the direction of the
#'     original effect estimates.
#' @param type Type of recalibration. Can be either "golden" (default), "nominal" (no recalibration),
#'  or "controlled". "golden" ensures that for an original study just significant at
#' the specified \code{level}, replication success is only possible for
#' replication effect estimates larger than the original one.
#' "controlled" ensures exact overall Type-I error control at level \code{level}^2.
#' @param shrinkage Numeric vector with values in [0,1). Defaults to 0.
#'     Specifies the shrinkage of the original effect estimate towards zero,
#'     e.g., the effect is shrunken by a factor of 25\% for
#'     \code{shrinkage = 0.25}. Is only taken into account if the
#'     \code{designPrior} is "conditional" or "predictive".
#' @param h Numeric vector of relative heterogeneity variances i.e., the ratios
#'     of the heterogeneity variance to the variance of the original effect
#'     estimate. Default is 0 (no heterogeneity). Is only taken into account
#'     when \code{designPrior} = "predictive" or \code{designPrior} = "EB".
#' @param strict Logical vector indicating whether the probability for
#'     replication success in the opposite direction of the original effect
#'     estimate should also be taken into account. Default is \code{FALSE}.
#'     Only taken into account when \code{alternative} = "two.sided".
#' @return The power for replication success with the sceptical p-value
#' @details \code{powerReplicationSuccess} is the vectorized version of
#' the internal function \code{.powerReplicationSuccess_}.
#' \code{\link[base]{Vectorize}} is used to vectorize the function.
#' @author Leonhard Held, Charlotte Micheloud, Samuel Pawel
#' @references
#' Held, L. (2020). A new standard for the analysis and design of replication
#' studies (with discussion). \emph{Journal of the Royal Statistical Society:
#' Series A (Statistics in Society)}, \bold{183}, 431-448.
#' \doi{10.1111/rssa.12493}
#'
#' Held, L., Micheloud, C., Pawel, S. (2022). The assessment of replication
#' success based on relative effect size. \emph{The Annals of Applied
#' Statistics}. 16:706-720. \doi{10.1214/21-AOAS1502}
#'
#' Micheloud, C., Balabdaoui, F., Held, L. (2023). Assessing replicability
#' with the sceptical p-value: Type-I error control and
#' sample size planning. \emph{Statistica Neerlandica}. \doi{10.1111/stan.12312}
#'
#' @seealso \code{\link{sampleSizeReplicationSuccess}}, \code{\link{pSceptical}},
#' \code{\link{levelSceptical}}
#' @examples
#' ## larger sample size in replication (c > 1)
#' powerReplicationSuccess(zo = p2z(0.005), c = 2, level = 0.025, designPrior = "conditional")
#' powerReplicationSuccess(zo = p2z(0.005), c = 2, level = 0.025, designPrior = "predictive")
#'
#' ## smaller sample size in replication (c < 1)
#' powerReplicationSuccess(zo = p2z(0.005), c = 1/2, level = 0.025, designPrior = "conditional")
#' powerReplicationSuccess(zo = p2z(0.005), c = 1/2, level = 0.025, designPrior = "predictive")
#'
#' powerReplicationSuccess(zo = p2z(0.00005), c = 2, level = 0.05,
#'                         alternative = "two.sided",  strict = TRUE, shrinkage = 0.9)
#' powerReplicationSuccess(zo = p2z(0.00005), c = 2, level = 0.05,
#'                         alternative = "two.sided", strict = FALSE, shrinkage = 0.9)
#'
#' @export
powerReplicationSuccess <- Vectorize(.powerReplicationSuccess_)


# function that returns z_r^2 quantile for given z_o, c
zr2quantile <- function(zo,
                        c,
                        p,
                        designPrior,
                        shrinkage) {

    if (designPrior == "predictive") {
        lambda <- (1 - shrinkage)^2 * zo^2 / (1 + 1 / c)
        factor <- c + 1
        # h <- 0 # relative heterogeneity h, implement later
        # lambda <- (1 - shrinkage)^2*zo^2/(1/c + 1 + 2*h)
        # factor <- 1 + c + 2*h*c
    }
    if (designPrior == "conditional") {
        lambda <- (1 - shrinkage)^2 * c * zo^2
        factor <- 1
    }
    if (designPrior == "EB") {
        s <- pmax(1 - 1 / zo^2, 0)
        lambda <- zo^2 * s^2 / (s + 1 / c)
        factor <- s * c + 1
        # h <- 0 # relative heterogeneity h, implement later
        # s <- pmax(1 - (1 + h)/zo^2, 0)
        # lambda <- zo^2*s^2/(s*(1 + h) + 1/c + h)
        # factor <- s*c*(1 + h) + 1 + h*c
    }
    if (lambda < 100)
        res <- stats::qchisq(p = p, df = 1, ncp = lambda)
    else
        res <- stats::qnorm(p = p, mean = sqrt(lambda), sd = 1)^2
    return(factor * res)
}


powerReplicationSuccessTargetPower <- function(
    power,
    zo,
    c,
    level,
    designPrior,
    alternative,
    type,
    shrinkage) {

    zr2 <- zr2quantile(zo = zo, c = c, p = 1 - power,
                       designPrior = designPrior, shrinkage = shrinkage)
    pC <- pSceptical(zo = zo, zr = sqrt(zr2), c = c,
                     alternative = alternative, type = type)
    return(pC - level)
    # pC <- pSceptical(zo = zo, zr = sqrt(zr2), c = c,
    #                  alternative = alternative)
    # return(pC - levelSceptical(level = level, alternative = alternative,
    #                            type = type))
}

.powerReplicationSuccessNum_ <- function(
    zo,
    c = 1,
    level = 0.025,
    designPrior = c("conditional", "predictive", "EB"),
    alternative = c("one.sided", "two.sided"),
    type = c("golden", "nominal", "controlled"),
    shrinkage = 0) {

    stopifnot(is.numeric(zo),
              length(zo) == 1,
              is.finite(zo),

              is.numeric(c),
              length(c) == 1,
              is.finite(c),
              0 <= c,

              is.numeric(level),
              length(level) == 1,
              is.finite(level),
              0 < level, level < 1,

              !is.null(designPrior))

    designPrior <- match.arg(designPrior)

    stopifnot(!is.null(alternative))
    alternative <- match.arg(alternative)

    stopifnot(!is.null(type))
    type <- match.arg(type)

    stopifnot(is.numeric(shrinkage),
              length(shrinkage) == 1,
              is.finite(shrinkage),
              0 <= shrinkage, shrinkage < 1)


    ## check if original study was not significant, then power is zero
    zo <- abs(zo)
    p <- z2p(z = zo, alternative = alternative)
    eps <- 1e-5
    mylower <- eps
    myupper <- 1 - eps

    if (p > levelSceptical(level = level, alternative = alternative,
                           type = type)) {
        res <- 0
    } else {
        targetLower <- powerReplicationSuccessTargetPower(
            power = mylower,
            zo = zo,
            c = c,
            level = level,
            designPrior = designPrior,
            alternative = alternative,
            type = type,
            shrinkage = shrinkage
        )
        targetUpper <- powerReplicationSuccessTargetPower(
            power = myupper,
            zo = zo,
            c = c,
            level = level,
            designPrior = designPrior,
            alternative = alternative,
            type = type,
            shrinkage = shrinkage
        )

        if (sign(targetLower) == sign(targetUpper)) {
            if ((sign(targetLower) >= 0) && (sign(targetUpper) >= 0))
                res <- 0
            if ((sign(targetLower) < 0) && (sign(targetUpper) < 0))
                res <- 1
        } else {
            res <- stats::uniroot(
                f = powerReplicationSuccessTargetPower,
                lower = mylower,
                upper = myupper,
                zo = zo,
                c = c,
                level = level,
                designPrior = designPrior,
                alternative = alternative,
                type = type,
                shrinkage = shrinkage
            )$root
        }
    }
    return(res)
}

powerReplicationSuccessNum <- Vectorize(.powerReplicationSuccessNum_)