R/Rfun_find_k_given_rho_target_mvtnorm.R

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

#' Find the difference between the error rate when k and rho are both given and the prespecified alpha level
#'
#' @param k a pre-specified constant in the improved trimmed weighted Hochberg procedure
#' @param rho the correlation coefficient between two test statistics
#' @param alpha the significance level
#' @param alphavec a numeric vector of two values representing the weighted significance levels assigned to the two hypotheses
#'
#' @returns the difference between the error rate when k is specified and rho is optimal and the prespecified alpha level
#' @export
#' @importFrom stats qnorm
#' @importFrom mvtnorm pmvnorm
#' @importFrom mvtnorm Miwa
#' @author Jiangtao Gou
#' @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.
find_k_given_rho_target_mvtnorm <- function (k, rho, alpha, alphavec = c(alpha/2, alpha/2)) {
  # Rfun_find_k_given_rho_target_mvtnorm
  adjFct <- k
  sigma <- matrix(c(1, rho, rho, 1), nrow = 2)
  # Middle
  lower <- c(-1000, -1000)
  upper <- c(qnorm(alpha), qnorm(alpha))
  RegionMiddleRslt <- mvtnorm::pmvnorm(lower = lower, upper = upper, mean = rep(0, 2), sigma = sigma, algorithm = mvtnorm::Miwa(steps = 128))
  RegionMiddle <- as.numeric(RegionMiddleRslt)
  # Left
  lower <- c(-1000, qnorm(alpha))
  upper <- c(qnorm(  alphavec[1]*(1+alphavec[2]*adjFct) ), qnorm(1-alphavec[2]) )
  RegionLeftRslt <- mvtnorm::pmvnorm(lower = lower, upper = upper, mean = rep(0, 2), sigma = sigma, algorithm = mvtnorm::Miwa(steps = 128))
  RegionLeft <- as.numeric(RegionLeftRslt)
  # Bottom
  lower <- c(qnorm(alpha), -1000)
  upper <- c(qnorm(1-alphavec[1]),  qnorm(  alphavec[2]*(1+alphavec[1]*adjFct)) )
  RegionBottomRslt <- mvtnorm::pmvnorm(lower = lower, upper = upper, mean = rep(0, 2), sigma = sigma, algorithm = mvtnorm::Miwa(steps = 128))
  RegionBottom <- as.numeric(RegionBottomRslt)
  #
  Region <- RegionMiddle + RegionLeft + RegionBottom
  #
  return (Region - alpha)
}