R/SimData_per_group.R

In winratiosim: Simulate Power for Hierarchical Win Ratio Endpoints

#' Simulate Individual-Level Trial Data for One Arm
#'
#' Generates individual-level simulated data for a treatment or control arm in
#' a hierarchical win ratio trial. The simulation includes frailty-adjusted
#' time to death, recurrent event counts, censoring times, and a continuous
#' quality-of-life change score.
#'
#' @param treatment Integer. Treatment group indicator, usually 1 for the
#'   active treatment arm and 0 for the control arm.
#' @param ngroup Integer. Number of subjects to simulate in this arm.
#' @param alpha.JFM Numeric. Alpha parameter for the joint frailty model.
#' @param theta.JFM Numeric. Frailty variance parameter for the joint frailty
#'   model. Must be positive.
#' @param lambda Numeric. Annual mortality probability. Must be in
#'   \code{[0, 1)}.
#' @param ann.icr Numeric. Annual incidence rate of recurrent events.
#' @param censorrate Numeric. Annual censoring probability. Must be in
#'   \code{[0, 1)}.
#' @param xbase Numeric. Baseline value of the continuous outcome.
#' @param xfinal Numeric. Expected final value of the continuous outcome among
#'   subjects followed through 360 days.
#' @param sd.delta.x Numeric. Standard deviation of the change in the
#'   continuous outcome.
#'
#' @return A named list. If \code{treatment = 1}, the list contains
#'   \code{surv_1}; otherwise, it contains \code{surv_0}. The data frame has
#'   one row per subject and includes subject ID, treatment indicator, death
#'   time, censoring time, death indicator, recurrent event count, and
#'   continuous outcome value.
#'
#' @examples
#' set.seed(1)
#' sim <- SimData_per_group(
#'   treatment = 1, ngroup = 5,
#'   alpha.JFM = 0, theta.JFM = 1,
#'   lambda = 0.13, ann.icr = 0.32,
#'   censorrate = 0.2, xbase = 45, xfinal = 52.5,
#'   sd.delta.x = 20
#' )
#' str(sim$surv_1)
#'
#' @export
SimData_per_group <- function(treatment, ngroup, alpha.JFM, theta.JFM,
                              lambda, ann.icr, censorrate, xbase, xfinal,
                              sd.delta.x) {
  if (length(ngroup) != 1L || !is.finite(ngroup) || ngroup < 1) {
    stop("ngroup must be a positive length-one integer.", call. = FALSE)
  }
  if (length(theta.JFM) != 1L || !is.finite(theta.JFM) || theta.JFM <= 0) {
    stop("theta.JFM must be positive.", call. = FALSE)
  }
  if (lambda < 0 || lambda >= 1 || censorrate < 0 || censorrate >= 1) {
    stop("lambda and censorrate must be in [0, 1).", call. = FALSE)
  }
  if (ann.icr < 0 || sd.delta.x < 0) {
    stop("ann.icr and sd.delta.x must be non-negative.", call. = FALSE)
  }

  n <- as.integer(ngroup)
  subject_ids <- seq_len(n)
  followup <- rep(12, n)
  followdays_subject <- followup * 30

  frailty <- stats::rgamma(n, shape = 1 / theta.JFM, rate = 1 / theta.JFM)

  beta_hfh <- ann.icr / 360
  repeated_ids <- rep(subject_ids, each = 100)
  repeated_frailty <- rep(frailty, each = 100)
  hfh_rate <- repeated_frailty * beta_hfh
  gap_times <- stats::rexp(length(repeated_ids), rate = hfh_rate)

  gap_by_subject <- split(gap_times, repeated_ids)
  start_times <- unlist(
    lapply(gap_by_subject, function(x) cumsum(c(0, x[-length(x)]))),
    use.names = FALSE
  )
  stop_times <- unlist(lapply(gap_by_subject, cumsum), use.names = FALSE)

  beta_mortality <- -log(1 - lambda) / 360
  mortality_rate <- (frailty ^ alpha.JFM) * beta_mortality
  mortd <- stats::rexp(n, rate = mortality_rate)

  censoring_rate <- -log(1 - censorrate) / 360
  ctime <- stats::rexp(n, rate = censoring_rate)

  df_event <- data.frame(
    subjid = repeated_ids,
    treatment = treatment,
    mortd = rep(mortd, each = 100),
    followup = rep(followup, each = 100),
    group = treatment,
    ftime = 360,
    ctime = rep(ctime, each = 100),
    followdays = rep(followdays_subject, each = 100),
    startm = start_times,
    stopm = stop_times
  )

  df_event$fudays <- pmin(df_event$followdays, df_event$ctime)
  df_event$censortime <- pmin(df_event$fudays, df_event$mortd)
  df_event$finalstopm <- pmin(df_event$stopm, df_event$censortime)
  df_event$censor_new <- ifelse(df_event$stopm > df_event$censortime, 0, 1)
  df_event$deathdays <- ifelse(
    df_event$mortd > df_event$fudays,
    floor(df_event$fudays),
    floor(df_event$mortd)
  )
  df_event$death <- ifelse(df_event$mortd > df_event$fudays, 0, 1)

  df_event <- df_event[!(df_event$startm > df_event$finalstopm), , drop = FALSE]
  df_event <- df_event[order(df_event$subjid, df_event$startm), , drop = FALSE]
  df_event$HFH <- as.numeric(stats::ave(
    df_event$censor_new,
    df_event$subjid,
    FUN = sum
  ))

  df_final <- df_event[!duplicated(df_event$subjid, fromLast = TRUE), ,
                       drop = FALSE]
  df_final <- df_final[order(df_final$subjid), , drop = FALSE]
  df_final <- df_final[, c("subjid", "treatment", "mortd", "followup",
                           "group", "ftime", "ctime", "followdays",
                           "fudays", "censortime", "deathdays", "death",
                           "HFH")]

  mean_x_change <- xfinal - xbase
  kccq_change <- pmax(
    pmin(stats::rnorm(n, mean = mean_x_change, sd = sd.delta.x), 100),
    -100
  )
  df_final$kccq <- ifelse(df_final$deathdays < 360, NA_real_, kccq_change)

  if (treatment == 1) {
    list(surv_1 = df_final)
  } else {
    list(surv_0 = df_final)
  }
}