R/getoc_2arm_piecewise.R

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

#' Compute Operating Characteristics for Two-Arm Piecewise Exponential Designs
#'
#' Computes the operating characteristics (e.g., type I error, power, and other performance metrics) of a two-arm survival trial design under the DTE-BOP2 framework with delayed treatment effect, based on a piecewise exponential model. This function is typically used after selecting design parameters (e.g., lambda, gamma) to assess the statistical properties of the proposed design.
#' @importFrom stats pbeta rexp runif
#' @importFrom truncdist rtrunc
#' @param median.true A numeric vector of length two.
#' \itemize{
#' \item The first element is the overall median survival time (in months) for the standard-of-care (SOC) arm under the null hypothesis.
#' \item The second element is the overall median survival time for the experimental arm.
#' }
#' @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.true[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.true[1])}
#' @param lambda Numeric. The tuning parameter \eqn{\lambda} in the decision boundary function \eqn{1 - \lambda \cdot (n/N)^\gamma}.
#' @param gamma Numeric. The tuning parameter \eqn{\gamma} in the decision boundary function \eqn{1 - \lambda \cdot (n/N)^\gamma}.
#' @param n.interim A numeric vector. Specifies the sample sizes per arm at each analysis. The last element represents the final sample size per arm.
#' @param L Numeric. Lower bound of the delayed treatment effect (DTE) separation time.
#' @param U Numeric. Upper bound of the DTE separation time.
#' @param S_likely Numeric. The most likely value of the separation time. Defaults to the midpoint of \code{L} and \code{U}.
#'   Default is the midpoint of \code{L} and \code{U}.
#' @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 rate Numeric. 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 nsim Integer. Number of simulations to generate. Default is 10000.
#' @param track Logical value.
#' \itemize{
#' \item \code{Default} FALSE
#' \item \code{If TRUE}, prints progress updates every 1000 simulations.
#' }
#' @param seed Optional integer. If provided, sets the seed for reproducibility.
#' @return A list with the following components:
#' \describe{
#'   \item{earlystop}{Numeric. The proportion of simulated trials that stopped early
#'     due to futility or predefined stopping rules.}
#'
#'   \item{reject}{Numeric. The proportion of trials in which the null hypothesis ($H_0$) was rejected.}
#'
#'   \item{average.patients}{Numeric. The average number of patients enrolled across all simulated trials.}
#'
#'   \item{trial.duration}{Numeric. The average total duration of the trial (in months), including follow-up.}
#' }
#' @examples
#' # Define trial parameters
#' median.1 <- 6
#' median.2 <- 9
#' trunc.para <- c(1, 1)
#' rate <- 3
#' FUP <- 9
#' n.interim <- c(30, 50)  # Each arm: 30 pts at interim, 50 pts at final
#' # Run operating characteristics computation
#' getoc_2arm_piecewise(median.true = c(median.1, median.2),
#'Uniform = FALSE,lambda = 0.9,gamma = 1,n.interim = n.interim,
#'   L = 3,U = 3,S_likely = 3,FUP = FUP,trunc.para = trunc.para,
#'   rate = rate,nsim = 1)
#' @export
getoc_2arm_piecewise <- function(median.true, gprior.E_1=NULL, gprior.E_2=NULL, lambda, gamma, n.interim, L, U,S_likely=(L+U)/2, Uniform=FALSE,trunc.para, rate, FUP = 6, nsim = 1000, track = FALSE,seed=NULL) {
  median_inuse <- function(median_0, median_1, S) {
    median.0 <- median_0


    if (median_0 < S) {
      median.1 <- median_0
    } else {
      median.1 <- (median_1-S)/(1-S/median_0)
    }

    return(c(median.0,median.1))
  }
  median_true=median_inuse(median.true[1],median.true[2],S=S_likely)
  if(is.null(gprior.E_1)){
    gprior.E_1=c(4,3*median_true[1]/log(2))
  }
  if(is.null(gprior.E_2)){
    gprior.E_2=c(4,6*median_true[1]/log(2))
  }
  lambda_true=log(2)/median_true
  nsim.T=1
  nmax = max(n.interim)
  nobs = n.interim+1
  nobs[length(nobs)] = nmax
  earlystop=rej=0
  t.trial = pts = rep(0,nsim)

  if (!is.null(seed)) {
    set.seed(seed)
  }

  for (sim in 1:nsim)
  {
    if(track) {if(sim%%1000==0) cat(sim,"\n")}

    wait.t = rexp(nmax,rate = rate)
    arrival.t = cumsum(wait.t)
    #for (i in 1:length(nobs)) {arrival.t = c(arrival.t,cumsum(wait.t)[(nobs.n[i]+1):(nobs.n[i+1])]+(i-1)*layover)}

    #   event.t.C = rexp(nmax,rate=log(2)/mu_true[2])
    tobs = arrival.t[nobs]
    tobs[length(tobs)] = tobs[length(tobs)] + FUP  ## For the last patient, we allow FUP of follow-up duration

    abort = 0
    k = 1

    data_E=data_C=matrix(0,nrow=nmax,ncol=nsim.T)

    for(j in 1:nsim.T){
      if(Uniform==FALSE){
        T.sep=rtrunc_gamma(n = 1,  L, U, shape = trunc.para[1], scale = trunc.para[2])   ##the event should be defined before the while loop
        event.t.E = generate_pe(nmax,T.sep,lambda1=lambda_true[1],lambda2=lambda_true[2])
        event.t.C=rexp(nmax,rate=lambda_true[1])
        data_E[,j]=event.t.E
        data_C[,j]=event.t.C
      }
      if(Uniform==TRUE){
        n.seq=20
        T.sep=sample(L+(U-L)/n.seq*seq(1,n.seq-1),1,replace = T)
        event.t.E = generate_pe(nmax,T.sep,lambda1=lambda_true[1],lambda2=lambda_true[2])
        event.t.C=rexp(nmax,rate=lambda_true[1])
        data_E[,j]=event.t.E
        data_C[,j]=event.t.C
      }
    }
    while(abort == 0)
    {temp=rep(0,nsim.T)

    for(j in 1:nsim.T){

      event.t.E=data_E[,j]
      event.t.C=data_C[,j]
      Ind_event_E = ifelse(arrival.t[1:n.interim[k]] + event.t.E[1:n.interim[k]] <= tobs[k],1,0)   ##event indicator for th treatment, 1: event.
      Ind_event_C = ifelse(arrival.t[1:n.interim[k]] + event.t.C[1:n.interim[k]] <= tobs[k],1,0)
      t.event.E  = rep(0,n.interim[k])
      t.event.C=rep(0,n.interim[k])
      for(l in 1:length(t.event.E)){
        t.event.E[l] = ifelse(arrival.t[l]+event.t.E[l]<=tobs[k],event.t.E[l],tobs[k]-arrival.t[l])
        t.event.C[l] = ifelse(arrival.t[l]+event.t.C[l]<=tobs[k],event.t.C[l],tobs[k]-arrival.t[l])

      }
      time_data=data.frame(Time=t.event.E,Status=Ind_event_E)
      TTOT=TTOT_fun(time_data,t=T.sep)
      para1=gprior.E_1[1]+TTOT$event_count_1+sum(Ind_event_C)
      para2=gprior.E_2[1]+TTOT$event_count_2
      cutoff=(gprior.E_1[2]+TTOT$TTOT_1+sum(t.event.C))/(gprior.E_1[2]+gprior.E_2[2]+TTOT$TTOT_1+TTOT$TTOT_2+sum(t.event.C))
      temp[j]=pbeta(cutoff,para1,para2)
    }

    if(mean(temp) > 1 - lambda * (n.interim[k]/nmax)^gamma) abort = 1
    else k = k+1

    if(k > length(n.interim)) break
    }

    if(abort==0) rej = rej + 1
    else {
      if(k < length(n.interim)) {earlystop = earlystop + 1}
    }

    t.trial[sim] = ifelse(k<=length(n.interim),tobs[k],tobs[length(tobs)])
    pts[sim] = ifelse(k<=length(n.interim),n.interim[k],nmax)

  }
  list(earlystop=earlystop/nsim,reject=rej/nsim,average.patients=mean(pts),trial.duration=mean(t.trial))
}