R/Rfun_avetrteff2.R

In cobenrich: Using Multiple Continuous Biomarkers for Patient Enrichment in Two-Stage Clinical Designs

#' @name avetrteff2
#' @title Compute the average subpopulation treatment effect and the standardized average subpopulation treatment effect when two biomarkers are involved
#'
#' @param z1z2 a numeric vector of two numbers that are standardized biomarker values
#' @param kappa a number of the correlation coefficient between two biomarkers
#' @param rhovec a numeric vector of two correlation coefficients between the output and two biomarkers
#' @param sigma a number of the standard deviation of outcome
#' @param muminusmu0 a number of the difference between the mean of outcome and the minimal clinically important treatment effect
#'
#' @returns a list of three numbers: \code{delta} is the average subpopulation treatment effect, \code{lambda} is the standardized average subpopulation treatment effect, and \code{cVar} is the variance with respect to the truncated distribution with specified cutoff values
#' @export
#' @importFrom tmvtnorm mtmvnorm
#'
#' @author Jiangtao Gou
#' @references
#' Zhang, F. and Gou, J. (2025). Using multiple biomarkers for patient enrichment in two-stage clinical designs. Technical Report.
#'
#' @examples
#' x1x2 <- c(2, 1)
#' nu1nu2 <- c(0,0)
#' tau1tau2 <- c(1,1)
#' z1z2 <- (x1x2 - nu1nu2)/tau1tau2
#' muminusmu0 <- 1.8
#' kappa <- 0.1
#' sigma <- 1
#' rhovec <- c(0.1, 0.2)
#' avetrteff2(z1z2, kappa, rhovec, sigma, muminusmu0)
avetrteff2 <- function (z1z2, kappa, rhovec, sigma, muminusmu0) {
  #
  corrmat <- matrix(c(1, kappa, kappa, 1), 2, 2)
  #
  hfunrslt <- tmvtnorm::mtmvnorm(mean = c(0, 0),
                       sigma = corrmat,
                       lower = z1z2,
                       upper = c(Inf, Inf),
                       doComputeVariance = TRUE)
  h10 <- hfunrslt$tmean[1]
  h01 <- hfunrslt$tmean[2]
  h20 <- hfunrslt$tvar[1,1] + (hfunrslt$tmean[1])^2
  h02 <- hfunrslt$tvar[2,2] + (hfunrslt$tmean[2])^2
  h11 <- hfunrslt$tvar[1,2] + hfunrslt$tmean[1] * hfunrslt$tmean[2]
  #
  cVar <- sigma^2/(1-kappa^2)^2*( 2*(rhovec[1] - kappa*rhovec[2])*(rhovec[2] - kappa*rhovec[1])*(h11 - h10*h01) + (rhovec[1] - kappa*rhovec[2])^2*(h20 - h10^2) + (rhovec[2] - kappa*rhovec[1])^2*(h02 - h01^2) ) + sigma^2*(1 - (rhovec[1]^2 + rhovec[2]^2 - 2*kappa*rhovec[1]*rhovec[2])/(1-kappa^2))
  #
  delta <- muminusmu0 + sigma/(1-kappa^2)*((rhovec[1] - kappa*rhovec[2])*h10  + (rhovec[2] - kappa*rhovec[1])*h01)
  #
  lambda <- delta/sqrt(cVar)
  #
  return (list(delta, lambda, cVar))
}