R/event_fun.R

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

#' Return expected event number at interim analysis and final analysis.
#'
#' Based on the overall median times (median.1, median.2), the most likely separation timepoint, accrual rate and follow-up duration, this function will return the expected number of events in total at each interim analysis and final analysis based on the simulations.
#' @param median.1 Numeric. The overall median survival time for SOC.
#' @param median.2 Numeric. The overall median survival time for the experimental arm.
#' @param S_likely Numeric. The most likely separation time. Defaults to the midpoint of \code{L} and \code{U}.
#' @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 n.sim Integer. Number of simulations to generate. Default is 1000.
#' @param seed Optional integer. If provided, sets the seed for reproducibility.
#' @return A data frame with two columns:
#'   the first column indicates the stage label,
#'   and the second column contains the expected number of events for each stage.
#' @examples
#' \donttest{
#' result <- event_fun(
#' median.1 = 2.8,
#' median.2 = 3.5,
#' S_likely = 2.1,
#' n.interim = c(28,40),
#' rate = 6,
#' FUP = 6
#' )
#' }
#' @export
event_fun <- function(median.1, median.2, S_likely=(L+U)/2, n.interim, rate, FUP,n.sim=1000,seed=NULL) {

  nmax <- max(n.interim)

  K <- length(n.interim)

  median_inuse <- function(median_0, median_1, S) {
    if (median_0 < S) {
      median.1 <- median_0
    } else {
      median.1 <- (median_1 - S) / (1 - S / median_0)
    }
    return(c(median_0, median.1))
  }

  event_H0 <- matrix(0, nrow = n.sim, ncol = K)
  event_H1 <- matrix(0, nrow = n.sim, ncol = K)
  if (!is.null(seed)) {
    set.seed(seed)
  }
  for (i in 1:n.sim) {
    medians <- median_inuse(median.1, median.2, S_likely)
    median_0 <- medians[1]
    median_1_adj <- medians[2]

    lambda_1 <- log(2) / median_0

    ## --- H0 ---
    wait.t <- rexp(nmax, rate = rate)
    arrival.t <- cumsum(wait.t)
    nobs <- n.interim + 1
    nobs[length(nobs)] <- nmax
    tobs <- arrival.t[nobs]
    tobs[length(tobs)] <- tobs[length(tobs)] + FUP

    event.t.E <- generate_pe(nmax, t = S_likely, lambda1 = lambda_1, lambda2 = lambda_1)
    event.t.C <- rexp(nmax, rate = lambda_1)

    for (k in 1:K) {
      inds <- 1:n.interim[k]
      Ind_event_E <- ifelse(arrival.t[inds] + event.t.E[inds] <= tobs[k], 1, 0)
      Ind_event_C <- ifelse(arrival.t[inds] + event.t.C[inds] <= tobs[k], 1, 0)
      event_H0[i, k] <- sum(Ind_event_E) + sum(Ind_event_C)
    }

    ## --- H1 ---
    lambda_2 <- log(2) / median_1_adj
    wait.t <- rexp(nmax, rate = rate)
    arrival.t <- cumsum(wait.t)
    tobs <- arrival.t[nobs]
    tobs[length(tobs)] <- tobs[length(tobs)] + FUP
    event.t.E <- generate_pe(nmax, t = S_likely, lambda1 = lambda_1, lambda2 = lambda_2)
    event.t.C <- rexp(nmax, rate = lambda_1)

    for (k in 1:K) {
      inds <- 1:n.interim[k]
      Ind_event_E <- ifelse(arrival.t[inds] + event.t.E[inds] <= tobs[k], 1, 0)
      Ind_event_C <- ifelse(arrival.t[inds] + event.t.C[inds] <= tobs[k], 1, 0)
      event_H1[i, k] <- sum(Ind_event_E) + sum(Ind_event_C)
    }
  }

  # Average of both conditions
  total_event <- 0.5 * event_H0 + 0.5 * event_H1
  total_expected <- ceiling(apply(total_event, 2, mean))

  # Label stages
  label_stages <- function(K) {
    suffix <- function(k) {
      if (k == 1) return("1st")
      if (k == 2) return("2nd")
      if (k == 3) return("3rd")
      return(paste0(k, "th"))
    }
    sapply(1:K, function(k) {
      if (k < K) paste0(suffix(k), " interim") else "Final stage"
    })
  }

  result <- data.frame(Events = as.vector(total_expected))
  rownames(result) <- label_stages(length(total_expected))
  return(result)
}