R/simHCE.R

In hce: Design and Analysis of Hierarchical Composite Endpoints

#' Simulate `hce` object with given event rates of time-to-event outcomes (Weibull), mean and SD of the continuous outcome (normal or log-normal) by treatment group
#'
#' @param n sample size in the active treatment group.
#' @param n0 sample size in the placebo group.
#' @param TTE_A event rates per year in the active group for the time-to-event outcomes.
#' @param TTE_P event rates per year in the placebo group for the time-to-event outcomes. Should have the same length as `TTE_A`.
#' @param CM_A mean value for the continuous outcome of the active group.
#' @param CM_P mean value for the continuous outcome of the placebo group.
#' @param CSD_A standard deviation for the continuous outcome of the active group.
#' @param CSD_P standard deviation for the continuous outcome of the placebo group.
#' @param fixedfy length of follow-up in years.
#' @param yeardays number of days in a year.
#' @param pat scale of provided event rates (per `pat`-years).
#' @param shape shape of the Weibull distribution for time-to-event outcomes. Default is exponential distribution with `shape = 1`.
#' @param logC logical, whether to use log-normal distribution for the continuous outcome.
#' @param seed for generating random numbers.
#' @param dec decimal places for the continuous outcome used for rounding. The default is `dec = 2`.
#' @return an object of class `hce` containing the following columns:
#' * ID subject identifier.
#' * TRTP planned treatment group - "A" for active, "P" for Placebo.
#' * GROUP type of the outcome, either "TTE" for time-to-event outcomes or "C" for continuous. 
#' Only one continuous outcome is possible, but no restriction on the number of "TTE" outcomes.
#' * GROUPN order of outcomes in `GROUP`, with a higher value signifying a better outcome.
#' * AVALT the timing of the time-to-event outcomes.
#' * AVAL0 numeric values of the continuous outcome and the timing of "TTE" outcomes.
#' * AVAL analysis values derived as `AVAL0 + GROUPN`. For the continuous outcome the values of `AVAL0` are shifted to start always from 0.
#' * seed the seed of the random sample. If not specified in `seed` argument will be selected based on system time.
#' * PADY primary analysis day, the length of fixed follow-up in days calculated as `yeardays` multiplied by `fixedfy`.
#' @export
#' @md
#' @seealso [hce::hce()], [hce::as_hce()] for the helper a coerce function to `hce` objects.
#' @examples
#' # Example 1
#' Rates_A <- c(1.72, 1.74, 0.58, 1.5, 1)
#' Rates_P <- c(2.47, 2.24, 2.9, 4, 6)
#' dat <- simHCE(n = 2500, TTE_A = Rates_A, TTE_P = Rates_P,
#'               CM_A = -3, CM_P = -6, CSD_A = 16, CSD_P = 15, fixedfy = 3)
#' head(dat)
#' 
#'# Example 2
#' Rates_A <- 10
#' Rates_P <- 15
#' dat <- simHCE(n = 1000, n0 = 500, TTE_A = Rates_A, TTE_P = Rates_P, 
#'               CM_A = 0.1, CM_P = 0, seed = 5, shape = 0.2, logC = TRUE, dec = 0)
#' summaryWO(dat)

simHCE <- function(n, n0 = n, TTE_A, TTE_P, CM_A, CM_P, CSD_A = 1, CSD_P = CSD_A, fixedfy = 1, yeardays = 360, pat = 100, shape = 1, logC = FALSE, seed = NULL, dec = 2 ){
  if(base::length(TTE_P) != base::length(TTE_A))
    stop("Event rate vectors for active and placebo groups should have the same length")
  
  n <- n[1]
  n0 <- n0[1]
  CM_A <- CM_A[1]
  CM_P <- CM_P[1]
  CSD_A <-  CSD_A[1]
  CSD_P <-  CSD_P[1]
  fixedfy <- fixedfy[1]
  yeardays <- yeardays[1]
  pat <- pat[1]
  seed <- seed[1]
  logC <- as.logical(logC[1])
  shape <- shape[1]
  dec <- dec[1]

  fixedf <- yeardays*fixedfy
  TTE_P <- TTE_P/yeardays/pat
  TTE_A <- TTE_A/yeardays/pat

  l <- base::length(TTE_A)
  if(is.null(seed)){
    seed <- base::as.numeric(base::Sys.time())
  }
  seed <- seed[1]
  base::set.seed(seed)


  M_P <- base::sapply(1/TTE_P, stats::rweibull, n = n0, shape = shape)
  M_P <- round(M_P)
  M_P <- base::as.data.frame(M_P)
  base::names(M_P) <- base::paste0("TTE", 1:l)
  M_P$TRTP <- "P"
  M_P$ID <- (n + 1):(n0 + n)

  M_A <- base::sapply(1/TTE_A, stats::rweibull, n = n, shape = shape)
  M_A <- round(M_A)
  M_A <- base::as.data.frame(M_A)
  base::names(M_A) <- base::paste0("TTE", 1:l)
  M_A$TRTP <- "A"
  M_A$ID <- 1:n

  dat <- base::rbind(M_A, M_P)
  dat$PADY <- fixedf
  dat$seed <- seed
  
  
  dat$AVALT <- base::apply(dat[, base::startsWith(names(dat), "TTE"), drop = FALSE], 1, function(x) x[x <= fixedf][1])
  dat$AVALT[is.na(dat$AVALT)] <- fixedf
  dat$EVENTN <- base::apply(dat[, base::startsWith(names(dat), "TTE"), drop = FALSE], 1, function(x) which(x <= fixedf)[1])
  dat$EVENTN[is.na(dat$EVENTN)]  <- l + 1

  dat$EVENT <- base::factor(dat$EVENTN, levels = 1:(l+1))
  base::levels(dat$EVENT) <- base::c(paste0("TTE", 1:l), "C")
  dat$GROUP <- base::as.character(dat$EVENT)
  dat$GROUPN <- dat$EVENTN*dat$PADY

  lXA <- round(stats::rlnorm(n, meanlog = CM_A, sdlog = CSD_A), dec)
  lXP <- round(stats::rlnorm(n0, meanlog = CM_P, sdlog = CSD_P), dec)
  XA <- round(stats::rnorm(n, mean = CM_A, sd = CSD_A), dec)
  XP <- round(stats::rnorm(n0, mean = CM_P, sd = CSD_P), dec)
  
  if(logC) {
    dat[dat$TRTP == "A", "AVAL0"] <- lXA
    dat[dat$TRTP == "P", "AVAL0"] <- lXP
  } else {
    dat[dat$TRTP == "A", "AVAL0"] <- XA
    dat[dat$TRTP == "P", "AVAL0"] <- XP
  }
  dat[dat$GROUP != "C", "AVAL0"] <- dat[dat$GROUP != "C", "AVALT"]
  mcont <- min(dat[dat$GROUP == "C", "AVAL0"])
  dat$AVAL <- ifelse(dat$GROUP != "C", dat$AVAL0 + dat$GROUPN, dat$AVAL0 - mcont + 1 + dat$GROUPN)
  dat <- dat[, c("ID", "TRTP", "GROUP", "GROUPN", "AVALT", "AVAL0", "AVAL", "seed")]
  dat$PADY <- fixedf
  dat <- as_hce(dat)
  return(dat)
}