R/get.optimal_2arm_piecewise.R

In DTEBOP2: Bayesian Optimal Phase II Randomized Clinical Trial Design with Delayed Outcomes

#' Obtain Optimal Design Parameters for DTE-BOP2 Design
#'
#' Computes the optimal tuning parameters \eqn{\lambda} and \eqn{\gamma} for use in posterior decision rules of two-arm, multi-stage survival trials under a piecewise exponential model with delayed treatment effect. These parameters are selected to control the type I error and maximize power, given user-specified interim sample sizes, prior distribution for the delayed separation timepoint, and hypothesized median survival times.
#' @importFrom stats pbeta rexp runif
#' @importFrom truncdist rtrunc
#' @param median.1 Numeric. The overall median survival time for the standard-of-care (SOC) arm under the null hypothesis.
#' @param median.2 Numeric. The overall median survival time for the experimental arm under the alternative hypothesis.
#' @param gprior.E_1 Optional. A numeric vector of length two specifying the shape and scale parameters of the inverse-gamma prior for the pre-separation mean survival time (i.e., 1/hazard rate). If NULL, the default is  \code{c(4,3/log(2)*median.1)}
#' @param gprior.E_2 Optional. A numeric vector of length two specifying the shape and scale parameters of the inverse-gamma prior for the post-separation mean survival time (i.e., 1/hazard rate). If NULL, the default is  \code{c(4,6/log(2)*median.1)}
#' @param L Numeric. The lower bound of the delayed treatment effect (DTE) separation timepoint.
#' @param U Numeric. The upper bound of the delayed treatment effect (DTE) separation timepoint.
#' @param S_likely Numeric. The most likely value of the DTE separation timepoint. Defaults to the midpoint of \code{L} and \code{U}.
#' @param lambda.seq A vector. The sequence of \eqn{\lambda} in \eqn{1-\lambda(n/N)^{\gamma}}
#' @param gamma.seq A vector. The sequence of \eqn{\gamma} in \eqn{1-\lambda(n/N)^{\gamma}}
#' @param Uniform Logical value.
#'  \itemize{
#' \item \code{Default} FALSE. The truncated gamma distribution for the separation time will be utilized.
#' \item \code{If TRUE} the average type I error and power are calculated based on 20 evenly divided points in the interval \eqn{[L,U]}.
#' }
#' @param trunc.para Vector value with two elements. The first element is the shape parameter for the truncated gamma prior and the second one is the scale parameter.
#' @param err1 Type I error rate.
#' @param nsim Integer. Number of simulations to generate. Default is 10000.
#' @param n.interim A vector of sample sizes per arm at each interim analysis.
#' \itemize{
#' \item Each element except the last represents an interim sample size per arm.
#' \item The final element is the total sample size N per arm.
#' }
#' @param rate Numeric value. Patient accrual rate (e.g., patients per month).
#' @param FUP Numeric value. Duration of follow-up. Default is 6 month/year in the context.
#' @param track Logical value. If true, it will report the current lambda value in the running simulations.
#' @param control Logical value. If TRUE, the type I error will be additionally controlled at control point. Default is FALSE.
#' @param control.point Vector. It specifies the points where the type I error is controlled.
#' @return A list with the following components:
#' \describe{
#'   \item{optimal}{A numeric vector of length four. Contains the optimal values of \eqn{\lambda} and \eqn{\gamma},
#'     along with the average type I error and average power achieved under these tuning parameters.}
#'
#'   \item{oc.mat}{A numeric matrix containing all evaluated combinations of \eqn{\lambda} and \eqn{\gamma}
#'     that satisfy the type I error constraint. Each row includes values for \eqn{\lambda}, \eqn{\gamma},
#'     average type I error, and average power.}
#' }

#' @examples
#' # Define design and simulation parameters
#' median.1 <- 6
#' median.2 <- 10
#' L <- 2.6
#' U <- 3.4
#' S_likely <- 2.7
#' trunc.para <- c(1, 1)
#' rate <- 3
#' FUP <- 9
#' err1 <- 0.15
#' n.interim <- c(30, 50)  # Each arm has 30 patients at interim and 50 at final
#' \donttest{
#' get.optimal_2arm_piecewise(
#'   median.1 = median.1,
#'   median.2 = median.2,
#'   L = L,
#'   U = U,
#'   Uniform = FALSE,
#'   S_likely = S_likely,
#'   trunc.para = trunc.para,
#'   err1 = err1,
#'   n.interim = n.interim,
#'   rate = rate,
#'   FUP = FUP,
#'   track = TRUE,
#'   nsim = 10000
#' )
#' }
#' @export
get.optimal_2arm_piecewise <- function(median.1,median.2, gprior.E_1=NULL, gprior.E_2=NULL,L,U,S_likely=(L+U)/2,lambda.seq = seq(0.5,0.975,by=0.025),gamma.seq = seq(0,1,by=0.1),trunc.para,Uniform=FALSE, err1,nsim,n.interim, rate, FUP,track=TRUE,control=FALSE,control.point=NULL)
{
  #cut.seq = c(0.05*1/2^(c(4,3,2,1)),seq(0.05,0.95,by=0.025))
  # lambda.seq = seq(0.5,0.975,by=0.025)
  # gamma.seq = seq(0,1,by=0.1)
  oc.mat = c()
  cut.start = 1

  for (j in 1:length(lambda.seq)){

    if(track) cat("Cutoff:",lambda.seq[j],"\n")
    for (k in cut.start:length(gamma.seq)){ #cat("cutoff & power:",c(cut.seq[j],power.seq[k]),"\n")
      temp1=getoc_2arm_piecewise(gprior.E_1=gprior.E_1,gprior.E_2=gprior.E_2,median.true = c(median.1,median.1),lambda=lambda.seq[j],gamma=gamma.seq[k],
                                 n.interim=n.interim, L=L,U=U,S_likely=S_likely,trunc.para=trunc.para,rate=rate, FUP=FUP,Uniform=Uniform, nsim=nsim, track = FALSE)
      if(control){
        temp1.temp_1=getoc_2arm_piecewise(gprior.E_1=gprior.E_1,gprior.E_2=gprior.E_2,median.true = c(median.1,median.1),lambda=lambda.seq[j],gamma=gamma.seq[k],
                                          n.interim=n.interim, L=min(control.point),U=min(control.point),S_likely=S_likely,trunc.para=trunc.para,rate=rate, FUP=FUP,Uniform=Uniform, nsim=nsim, track = FALSE)
        temp1.temp_2=getoc_2arm_piecewise(gprior.E_1=gprior.E_1,gprior.E_2=gprior.E_2,median.true = c(median.1,median.1),lambda=lambda.seq[j],gamma=gamma.seq[k],
                                          n.interim=n.interim, L=max(control.point),U=max(control.point),S_likely=S_likely,trunc.para=trunc.para,rate=rate, FUP=FUP,Uniform=Uniform, nsim=nsim, track = FALSE)
        temp.error=max(temp1.temp_1$reject,temp1.temp_2$reject)

        if ((temp1$reject>err1)|(temp.error>err1)) break;}
      if (temp1$reject>err1) break;
      temp2=getoc_2arm_piecewise(gprior.E_1=gprior.E_1,gprior.E_2=gprior.E_2,median.true = c(median.1,median.2),lambda=lambda.seq[j],gamma=gamma.seq[k],
                                 n.interim=n.interim, L=L,U=U,S_likely=S_likely,trunc.para=trunc.para,rate=rate, FUP=FUP,Uniform=Uniform, nsim=nsim, track = FALSE)
      #if (temp1$new_accept<=0.1&temp2$new_accept>=0.9)
      oc.mat <- rbind(oc.mat,c(lambda.seq[j],gamma.seq[k],temp1$reject,temp2$reject))
    }
    cut.start = k
    if(cut.start == length(gamma.seq)) break
  }
  # Find the maximum value in the 4th column
  max_col4 <- max(oc.mat[, 4])

  # Filter rows where the 4th column equals the maximum value
  rows_with_max_col4 <- oc.mat[oc.mat[, 4] == max_col4, ]

  # Find the row with the maximum value in the 2nd column among the filtered rows
  if (is.matrix(rows_with_max_col4)) {
    # If the result is a matrix (multiple rows), find the row with the max value in the 2nd column
    optimal <- rows_with_max_col4[which.max(rows_with_max_col4[, 2]), ]
  } else {
    # If there's only one row, return it directly
    optimal <- rows_with_max_col4
  }

  # Print the optimal row
  return(list(optimal=optimal,oc.mat=oc.mat))

}