R/Rfun_itwcHochPower.R

In itrimhoch: Improved Trimmed Weighted Hochberg Procedures and Sample Size Optimization

#' Power for rejecting H1 using various types of the Hochberg Procedure
#'
#' @param n the sample size
#' @param alpha1 the weighted significance levels assigned to H1
#' @param alpha the significance level
#' @param deltavec a numeric vector of two values representing the effect sizes for the two hypotheses
#' @param rho the correlation coefficient between two test statistics
#' @param proctype the improved trimmed weighted Hochberg procedure is denoted by `i`, the trimmed weighted Hochberg procedure is denoted by `t` , the weighted Hochberg procedure is denoted by `w`, and the conservative weighted Hochberg procedure is denoted by `c`
#' @param k a pre-specified constant in the improved trimmed weighted Hochberg procedure
#'
#' @returns the power for rejecting H1 is denoted by `pwr1`, the power for rejecting H2 is denoted by `pwr2`, and the power for rejecting both H1 and H2 is denoted by `pwr12`
#' @export
#' @importFrom mvtnorm pmvnorm
#' @importFrom mvtnorm Miwa
#' @importFrom stats uniroot
#' @author Jiangtao Gou
#' @author Fengqing Zhang
#' @references
#' Gou, J., Chang, Y., Li, T., and Zhang, F. (2025). Improved trimmed weighted Hochberg procedures with two endpoints and sample size optimization. Technical Report.
#'
#' @examples
#' itwcHochPower(n = 100,
#' alpha1 = 0.0125, alpha = 0.025,
#' deltavec = c(0.2, 0.25), rho = 0.2,
#' proctype = "i", k = 0)
#' itwcHochPower(n = 100,
#' alpha1 = 0.0125, alpha = 0.025,
#' deltavec = c(0, 0), rho = 0,
#' proctype = "w", k = 0)
itwcHochPower <- function (n, alpha1, alpha, deltavec, rho, proctype = "i", k = 0) {
  # 2024-12-19
  # Originally in R20241219a.R
  # Applied in R20241219b-e
  pwr1 <- 0
  pwr2 <- 0
  pwr12 <- 0
  #
  sigma <- matrix(c(1, rho, rho, 1), nrow = 2)
  #
  lower <- c( qnorm(1 - alpha) - deltavec[1]*sqrt(n), qnorm(1 - alpha) - deltavec[2]*sqrt(n) )
  upper <- c(1000, 1000)
  pwr12Rslt <- mvtnorm::pmvnorm(lower = lower, upper = upper, mean = rep(0, 2), sigma = sigma, algorithm = mvtnorm::Miwa(steps = 128))
  pwr12 <- as.numeric(pwr12Rslt)
  #
  if (abs(alpha1) < 1e-15) {
    pwr1 <- 0
  } else {
    if (proctype == "i") {
      lower <- c( qnorm(1 - alpha1*(1+k*(alpha - alpha1))) - deltavec[1]*sqrt(n), qnorm(alpha - alpha1) - deltavec[2]*sqrt(n) )
    } else if (proctype == "t") {
      lower <- c( qnorm(1 - alpha1) - deltavec[1]*sqrt(n), qnorm(alpha - alpha1) - deltavec[2]*sqrt(n) )
    } else if (proctype == "w") {
      lower <- c( qnorm(1 - alpha1) - deltavec[1]*sqrt(n), -1000 )
    } else if (proctype == "c") {
      lower <- c( qnorm(1 - alpha1*(1-alpha)) - deltavec[1]*sqrt(n), -1000 )
    }
    upper <- c(1000, qnorm(1 - alpha) - deltavec[2]*sqrt(n) )
    #
    pwr1Rslt <- mvtnorm::pmvnorm(lower = lower, upper = upper, mean = rep(0, 2), sigma = sigma, algorithm = mvtnorm::Miwa(steps = 128))
    pwr1 <- as.numeric(pwr1Rslt)
  }
  #
  #
  if (abs(alpha - alpha1) < 1e-15) {
    pwr2 <- 0
  } else {
    if (proctype == "i") {
      lower <- c( qnorm(alpha1) - deltavec[1]*sqrt(n), qnorm(1 - (alpha - alpha1)*(1+k*alpha1)) - deltavec[2]*sqrt(n) )
    } else if (proctype == "t") {
      lower <- c( qnorm(alpha1) - deltavec[1]*sqrt(n), qnorm(1 - alpha + alpha1) - deltavec[2]*sqrt(n) )
    } else if (proctype == "w") {
      lower <- c( -1000, qnorm(1 - alpha + alpha1) - deltavec[2]*sqrt(n) )
    } else if (proctype == "c") {
      lower <- c( -1000, qnorm(1 - (alpha - alpha1) * (1 - alpha) ) - deltavec[2]*sqrt(n))
    }
    upper <- c(qnorm(1 - alpha) - deltavec[1]*sqrt(n) , 1000 )
    #
    pwr2Rslt <- mvtnorm::pmvnorm(lower = lower, upper = upper, mean = rep(0, 2), sigma = sigma, algorithm = mvtnorm::Miwa(steps = 128))
    pwr2 <- as.numeric(pwr2Rslt)
  }
  #
  return (c(pwr1, pwr2, pwr12))
}