R/TwoSample.generate.sequential.R

In gsMeanFreq: Group Sequential Clinical Trial Designs for Composite Endpoints

#' @title Function to simulate two-sample composite endpoint data under staggered entry.
#' @description Simulate two-sample composite endpoints data with recurrent events and a terminal event
#' under two time scales: event time \code{t} and calendar time \code{s}. A uniform recruitment period is assumed,
#' and the function returns all observed data available at a specified calendar time. Recurrent event occurrences
#' are generated from an underlying Poisson process with subject-specific Gamma frailty.
#'
#' @param lambda_0vec Numeric vector of length 2 giving the baseline recurrent event rate parameters for two arms, default is \code{c(1.15, 1.15)}.
#' @param sizevec Integer vector of size 2 giving the group sizes.
#' @param beta.trt Numeric value giving the treatment effect coefficient applied to the treatment indicator (\code{Z1}) in the proportional mean model (arm 2 vs arm 1).
#' @param calendar Calendar time of the end of the trial (in years), default is \code{5}.
#' @param recruitment Length of the recruitment period (in years), default is \code{3}.
#' @param random.censor.rate Rate parameter for independent random right censoring.
#' @param seed Seed for reproducibility.
#' @returns A data frame in long format containing simulated composite endpoint data.
#' Each subject may contributing multiple rows corresponding to recurrent events, a terminal event (death), or censoring. The data include:
#' \itemize{
#' \item \code{group}: Arm indicator (\code{1} = control, \code{2}= treatment).
#' \item \code{id}: Subject identifier (unique across both arms).
#' \item \code{e}: Enrollment time on the calendar scale.
#' \item \code{event_time_cal}: Cumulative event time on the calendar scale.
#' \item \code{status}: Event indicator with values
#' \code{2}=recurrent event, \code{1}=death, and \code{0}=censoring.
#' \item \code{Z1}, \code{Z2}: Simulated covariates, where \code{Z1} is the treatment indicator (\code{0} for arm 1, \code{1} for arm 2).
#' \item \code{tau_star}: Subject-specific stopping time, the last event observed in \code{[0, tau_star]} is classified as death.
#' \item \code{death}: Binary indicator for death.
#' \item \code{recurrent}: Binary indicator for recurrent events.
#' \item \code{event}: Binary event indicator, \code{event = death + recurrent}.
#' \item \code{calendar}: Calendar time cutoff used to generate the returned data.
#' \item \code{lambda_0}: Baseline Poisson process rate parameter.
#' \item \code{lambda_star}: Rate parameter of an exponential distribution in generating \code{tau_star}.
#' \item \code{gamma_scale}, \code{gamma_shape}: Parameters of the Gamma distribution used to generate subject-specific frailty terms.
#' }
#' @export
#' @importFrom dplyr mutate
#' @importFrom stats rgamma rexp rnorm rbinom runif
#' @importFrom pracma newtonRaphson
#' @references Mao L, Lin DY. Semiparametric regression for the weighted composite endpoint of recurrent and terminal events. \emph{Biostatistics}. 2016 Apr; \strong{17(2)} :390-403.
#'
#' @examples
#' # Generate two-sample composite endpoint data
#' df <- TwoSample.generate.sequential(sizevec = c(200, 200),
#' beta.trt = 0.8, calendar = 5, recruitment = 3,
#' random.censor.rate = 0.05, seed = 2026)
#'
TwoSample.generate.sequential <- function(lambda_0vec = c(1.15, 1.15), sizevec, beta.trt,
                                          calendar = 5, recruitment = 3,
                                          random.censor.rate, seed){
  set.seed(seed)
  df <- NULL

  lambda_starvec = c(0.1, 0.1)
  gamma_shape = 2
  gamma_scale = 0.5

  study_end <- 5 # Administrative censoring time, assuming the trial ends at year 5
  admin_censoring <- min(study_end, calendar)

  # Two arms are simulated separately
  for (a in 1:2){
    # df.group <- NULL
    n <- sizevec[a]

    # Step1: simulate a recruitment time 'e'
    enroll <- runif(n, 0, recruitment)

    #Qinghua 04/07/2025 Update: Assume certain percentage of random censoring per year
    if (random.censor.rate > 0){
      random.censor.lambda <- -log(1- random.censor.rate)
      random.censor.cal <- enroll + rexp(n, random.censor.lambda)
    }

    lambda_star <- lambda_starvec[a]
    lambda_0 <- lambda_0vec[a]

    pre <- c(0, sizevec[2])
    for (j in 1:n){

      # Step 2: Simulate recurrent events from a homogeneous Poisson process

      # gamma distribution: mean = shape*scale, variance = shape*scale^2
      xi <- rgamma(1, shape = gamma_shape, scale = gamma_scale)
      # beta <- c(0,0)

      # Qinghua 04/23/2025 Update: set group 1 as the placebo group and group 2 as the treatment group
      # set Z1, the binary covariate as the treatment indicator and use beta_1 to differentiate the two groups
      if (a == 1){
        beta = c(0, 0)
      } else {
        beta <- c(beta.trt,0)
      }

      # Qinghua 5/01/2025 Update: changed Z1 to binary and Z1 = 0 for group 1 (placabo), Z1 = 1 for group 2 (treatment).
      # Z1 <- rbinom(1, 1, 0.5)
      Z1 = ifelse(a == 2, 1, 0)
      Z2 <- rnorm(1)
      Z <- c(Z1,Z2)
      times <- rexp(100, rate = lambda_0*as.vector(exp(Z%*%beta))*xi)

      # Step 3: Simulate tau_star, the stopping time (unique for each patient), the last event in [0, tau_star] is labeled as death.
      # Step 3.1: simulate tau_tilde

      F_tau_tilde2 <- function(t){
        (1 - exp(-lambda_star*xi*t))/(1 - exp(-lambda_0*xi*as.vector(exp(Z%*%beta))*t))
      }

      u <- runif(1, 0, 1)
      if (u <= lambda_star/(lambda_0*exp(as.vector(Z%*%beta)))){
        tau_tilde <- 0
      }else{
        tau_tilde <- newtonRaphson(function(t) {F_tau_tilde2(t)-u}, 0.05)$root
      }

      # Step 3.2: find tau_1 (the smallest event time)
      tau_1 <- times[1]
      tau_star <- max(tau_1, tau_tilde)
      event_times <- cumsum(times)
      event_times <- event_times[event_times <= tau_star] # a vector of all event times (recurrent and death)

      # Events starts happening after the patient is enrolled
      # 'event_times_calendar' is a vector of event on the calendar scale
      event_times_calendar = enroll[j] + event_times

      # Step 4: Apply the censoring at calendar time
      if (random.censor.rate > 0){
        censoring_time = min(admin_censoring, random.censor.cal[j])
      } else {
        censoring_time = admin_censoring
      }


      obs_time <- event_times_calendar[event_times_calendar <= censoring_time] # Observed times

      if (enroll[j] > censoring_time) {
        # patient has not been enrolled yet, should not contribute to the total sample size
        obs_time <- NA
        status <- NA
      } else if (length(obs_time) == 0) {
        # patient was enrolled and no events happened before given calendar time
        obs_time <- censoring_time
        status <- 0
      } else if (length(obs_time) < length(event_times_calendar)) {
        # patient was enrolled and experienced at least one event and censored at this given calendar time
        obs_time <- c(obs_time, censoring_time)
        status <- c(rep(2, length(obs_time) - 1),0)
      } else {
        # patient is enrolled and died before this given calendar time
        status <- c(rep(2, length(obs_time) - 1),1)
      }
      id <- rep((j+pre[a]), length(obs_time))
      e <- rep(enroll[j], length(obs_time))
      Z1 <- rep(Z1, length(obs_time))
      Z2 <- rep(Z2, length(obs_time))
      temp <- data.frame(group = a, id = id, e = e, event_time_cal = obs_time, status = status,
                         Z1 = Z1,Z2 = Z2,
                         tau_star = tau_star,
                         lambda_0 = lambda_0, lambda_star = lambda_star,
                         gamma_scale = gamma_scale, gamma_shape = gamma_shape)
      df <- rbind(df, temp)

      # Qinghua 05/07/2025 Update: Add clean up functions to save memory
      rm(times, xi, Z, event_times, event_times_calendar, temp, obs_time, status, id, e, Z1, Z2)
      # gc(verbose = FALSE) # Running the 'gc' function free up the memory but seems to be time consuming

    } # End of the "j" loop

  } # End of the "a" loop

  # df <- df %>% mutate(death = ifelse(status == 1,1,0)) %>%
  #   mutate(recurrent = ifelse(status == 2,1,0)) %>%
  #   mutate(event = death + recurrent) %>%
  #   mutate(calendar = calendar)

  df <- df %>%
    dplyr::mutate(
      death     = ifelse(.data$status == 1, 1, 0),
      recurrent = ifelse(.data$status == 2, 1, 0),
      event     = ifelse(.data$status %in% c(1, 2), 1, 0),
      calendar  = calendar
    )


  return(df)
}