simIDM coverage - 97.08%

Files
Source

#' Event Times Distributed as Sum of Weibull
#'
#' This returns event times with a distribution resulting from the sum of two Weibull distributed random variables
#' using the inversion method.
#'
#' @param U (`numeric`)\cr uniformly distributed random variables.
#' @param haz1 (positive `number`)\cr  first summand (constant hazard).
#' @param haz2 (positive `number`)\cr  second summand (constant hazard).
#' @param p1 (positive `number`)\cr rate parameter of Weibull distribution for `haz1`.
#' @param p2 (positive `number`)\cr rate parameter of Weibull distribution for `haz2`.
#' @param entry (`numeric`)\cr the entry times in the current state.
#'
#' @return This returns a vector with event times.
#' @export
#'
#' @examples
#' getWaitTimeSum(U = c(0.4, 0.5), haz1 = 0.8, haz2 = 1, p1 = 1.1, p2 = 1.5, entry = c(0, 0))
getWaitTimeSum <- function(U, haz1, haz2, p1, p2, entry) {
  assert_numeric(U, lower = 0, upper = 1, any.missing = FALSE, all.missing = FALSE)
  assert_positive_number(haz1, zero_ok = TRUE)
  assert_positive_number(haz2, zero_ok = TRUE)
  assert_positive_number(p1)
  assert_positive_number(p2)
  assert_numeric(entry, lower = 0, len = length(U), any.missing = FALSE, all.missing = FALSE)

  N <- length(U)
  # Exponential distributed survival times.
  if (p1 == 1 && p2 == 1) {
    return(-log(1 - U) / (haz1 + haz2))
  }
  # Weibull distributed survival times.
  temp <- function(x, y, t0) {
    return(haz1 * (x + t0)^p1 - haz1 * t0^p1 - haz2 * t0^p2 + haz2 * (x + t0)^p2 + y)
  }
  stime <- NULL
  i <- 1
  while (length(stime) < N) {
    u <- U[i]
    t0 <- entry[i]
    if (temp(0, log(1 - u), t0) * temp(10000, log(1 - u), t0) < 0) {
      res <- stats::uniroot(temp, c(0, 10000), tol = 10^-16, y = log(1 - u), t0 = t0)
      stime[i] <- res$root
      i <- i + 1
      if (res$root == 0) {
        warning("uniroot: accuracy (tol=10^-16) is not enough.")
      }
    } else {
      warning("uniroot: Values at interval endpoints (interval=c(0,10000)) are not of opposite signs. \n")
    }
  }
  stime
}

#' Time-point by which a specified number of events occurred.
#'
#' This returns the study time-point by which a specified number of events (PFS or OS) occurred.
#'
#' @param data (`data.frame`)\cr illness-death data set in `1rowPatient` format.
#' @param eventNum (`int`)\cr number of events.
#' @param typeEvent (`string`)\cr type of event. Possible values are `PFS` and `OS`.
#' @param byArm  (`logical`)\cr if `TRUE` time-point per treatment arm, else joint evaluation
#'   of treatment arms.
#'
#' @return This returns  the time-point by which `eventNum` of `typeEvent`-events occurred.
#' @export
#'
#' @examples
#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
#'
#' simStudy <- getOneClinicalTrial(
#'   nPat = c(20, 20), transitionByArm = list(transition1, transition2),
#'   dropout = list(rate = 0.3, time = 10),
#'   accrual = list(param = "time", value = 0)
#' )
#' simStudyWide <- getDatasetWideFormat(simStudy)
#' getTimePoint(simStudyWide, eventNum = 10, typeEvent = "OS", byArm = FALSE)
getTimePoint <- function(data, eventNum, typeEvent, byArm = FALSE) {
  assert_data_frame(data, ncols = 11)
  assert_int(eventNum, lower = 1L)
  assert_choice(typeEvent, c("OS", "PFS"))
  assert_logical(byArm)

  if (!byArm) {
    data$trt <- "all"
  }
  byVar <- unique(data$trt)

  if (typeEvent == "OS") {
    data$event <- data$OSevent
    data$time <- data$OStimeCal
  } else if (typeEvent == "PFS") {
    data$event <- data$PFSevent
    data$time <- data$PFStimeCal
  }

  sapply(byVar, function(trt) {
    dataTemp <- data[data$trt == trt, ]
    sortedTimes <- sort(dataTemp$time[dataTemp$event == 1])
    if (eventNum > length(sortedTimes)) {
      warning("Less events occur till EOS than specified \n")
    }
    timePoint <- sortedTimes[eventNum]
  })
}

#' Helper function for `censoringByNumberEvents`
#'
#' @param time (`numeric`) \cr event times.
#' @param event (`numeric`)\cr event indicator.
#' @param data (`data.frame`)\cr data frame including patient id `id`, recruiting time `recruitTime`
#'  and individual censoring time `censTimeInd`.
#'
#' @return This function returns a data frame with columns:
#' event time, censoring indicator, event indicator and event time
#' in calendar time.
#' @export
#'
#' @examples
#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
#'
#' simStudy <- getOneClinicalTrial(
#'   nPat = c(20, 20), transitionByArm = list(transition1, transition2),
#'   dropout = list(rate = 0.3, time = 10),
#'   accrual = list(param = "time", value = 7)
#' )
#' simStudyWide <- getDatasetWideFormat(simStudy)
#' simStudyWide$censTimeInd <- 5 - simStudyWide$recruitTime
#' NotRecruited <- simStudyWide$id[simStudyWide$censTimeInd < 0]
#' censoredData <- simStudyWide[!(simStudyWide$id %in% NotRecruited), ]
#' getCensoredData(time = censoredData$OStime, event = censoredData$OSevent, data = censoredData)
getCensoredData <- function(time, event, data) {
  assert_numeric(time, lower = 0)
  assert_numeric(event, lower = 0, upper = 1)
  assert_data_frame(data)

  # Event time censored?
  epsilon <- 1e-10
  Censored <- data$id[(time - epsilon) > data$censTimeInd]
  # keep minimum of censoring and event time. epsilon)
  Censoredtime <- pmin(time, data$censTimeInd)
  # adjust event and censoring indicators.
  CensoredEvent <- ifelse(data$id %in% Censored, 0, event)
  Censored <- abs(1 - CensoredEvent)
  # calculate corresponding calendar times.
  timeCal <- Censoredtime + data$recruitTime

  data.frame(
    time = Censoredtime,
    Censored = Censored,
    event = CensoredEvent,
    timeCal = timeCal
  )
}

#' Event-driven censoring.
#'
#' This function censors a study after a pre-specified number of events occurred.
#'
#' @param data (`data.frame`)\cr illness-death data set in `1rowPatient` format.
#' @param eventNum  (`int`)\cr number of events.
#' @param typeEvent (`string`)\cr type of event. Possible values are `PFS` and `OS`.
#'
#' @return This function returns a data set that is censored after `eventNum` of
#'   `typeEvent`-events occurred.
#' @export
#'
#' @examples
#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
#'
#' simStudy <- getOneClinicalTrial(
#'   nPat = c(20, 20), transitionByArm = list(transition1, transition2),
#'   dropout = list(rate = 0.3, time = 10),
#'   accrual = list(param = "time", value = 7)
#' )
#' simStudyWide <- getDatasetWideFormat(simStudy)
#' censoringByNumberEvents(data = simStudyWide, eventNum = 20, typeEvent = "PFS")
censoringByNumberEvents <- function(data, eventNum, typeEvent) {
  assert_data_frame(data, ncols = 11)
  assert_int(eventNum, lower = 1L)
  assert_choice(typeEvent, c("OS", "PFS"))

  # first step get timePoint (calendar time) for censoring.
  censTime <- getTimePoint(data, eventNum, typeEvent, byArm = FALSE)

  # censoring time at individual time scale.
  data$censTimeInd <- censTime - data$recruitTime
  # patients that are not yet recruited at censoring time has to be deleted.
  NotRecruited <- data$id[data$censTimeInd < 0]
  censoredData <- data[!(data$id %in% NotRecruited), ]

  OSdata <- getCensoredData(time = censoredData$OStime, event = censoredData$OSevent, data = censoredData)
  PFSdata <- getCensoredData(time = censoredData$PFStime, event = censoredData$PFSevent, data = censoredData)

  data.frame(
    id = censoredData$id, trt = censoredData$trt,
    PFStime = PFSdata$time, PFSevent = PFSdata$event,
    OStime = OSdata$time, CensoredOS = OSdata$Censored, OSevent = OSdata$event,
    recruitTime = censoredData$recruitTime, OStimeCal = OSdata$timeCal, PFStimeCal = PFSdata$timeCal
  )
}

#' Number of recruited/censored/ongoing Patients.
#'
#'
#' @param data (`data.frame`)\cr illness-death data set in `1rowPatient` format.
#' @param t (`numeric`)\cr  study time-point.
#'
#' @return This function returns number of recruited patients,
#' number of censored and number of patients under observations.
#' @export
#'
#' @examples
#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
#'
#' simStudy <- getOneClinicalTrial(
#'   nPat = c(20, 20), transitionByArm = list(transition1, transition2),
#'   dropout = list(rate = 0.6, time = 10),
#'   accrual = list(param = "time", value = 0)
#' )
#' simStudyWide <- getDatasetWideFormat(simStudy)
#' getEventsAll(data = simStudyWide, t = 1.5)
getEventsAll <- function(data, t) {
  assert_data_frame(data, ncols = 11)
  assert_positive_number(t, zero_ok = FALSE)

  allRecruited <- data$id[data$recruitTime <= t]
  allCensored <- data$id[(data$OStime + data$recruitTime) <= t & data$OSevent == 0]
  allDeath <- data$id[(data$OStime + data$recruitTime) <= t & data$OSevent == 1]
  ## recruited, not death, not censored.
  allUnderObs <- data$id[data$id %in% allRecruited & !(data$id %in% allCensored) &
    !(data$id %in% allDeath)]
  numRecruited <- length(unique(allRecruited))
  numCensored <- length(unique(allCensored))
  numUnderObs <- length(unique(allUnderObs))

  c(
    Recruited = numRecruited,
    Censored = numCensored,
    UnderObs = numUnderObs
  )
}

#' Helper Function for `trackEventsPerTrial`
#'
#' @param event (`numeric`)\cr event indicator.
#' @param time (`numeric`) \cr event times.
#' @param t  (`numeric`)\cr  study time-point.
#'
#' @return This function returns the number of events occurred until time t.
#' @export
#'
#' @examples
#' event <- c(0, 1, 1, 1, 0)
#' time <- c(3, 3.4, 5, 6, 5.5)
#' getNumberEvents(event = event, time = time, t = 5)
getNumberEvents <- function(event, time, t) {
  assert_numeric(event, lower = 0, upper = 1)
  assert_numeric(time, lower = 0)
  assert_number(t, lower = 0)

  sum(event[time <= t])
}


#' Event tracking in an oncology trial.
#'
#' @param data (`data.frame`)\cr illness-death data set in `1rowPatient` format.
#' @param timeP (`numeric`)\cr  vector of study time-points.
#' @param byArm  (`logical`)\cr  if `TRUE` time-point per treatment arm, else joint evaluation of treatment arms.
#'
#' @return This function returns a data frame including number of PFS events, number of OS events,
#'  number of recruited patients, number of censored patients and number of ongoing patients at `timeP`.
#'
#' @export
#'
#' @examples
#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
#'
#' simStudy <- getOneClinicalTrial(
#'   nPat = c(20, 20), transitionByArm = list(transition1, transition2),
#'   dropout = list(rate = 0.3, time = 10),
#'   accrual = list(param = "time", value = 0)
#' )
#' simStudyWide <- getDatasetWideFormat(simStudy)
#' trackEventsPerTrial(data = simStudyWide, timeP = 1.5, byArm = FALSE)
trackEventsPerTrial <- function(data, timeP, byArm = FALSE) {
  assert_data_frame(data, ncols = 11)
  assert_numeric(timeP, lower = 0)
  assert_logical(byArm)

  if (!byArm) {
    data$trt <- "all"
  }
  byVar <- unique(data$trt)
  allNumbers <- lapply(byVar, function(j) {
    datTemp <- data[data$trt == j, ]
    eventsPFS <- sapply(timeP, getNumberEvents, event = datTemp$PFSevent, time = datTemp$PFStimeCal)
    eventsOS <- sapply(timeP, getNumberEvents, event = datTemp$OSevent, time = datTemp$OStimeCal)
    eventsTrial <- sapply(timeP, getEventsAll, data = datTemp)

    allNumbers <- rbind(eventsPFS, eventsOS, eventsTrial)
    rownames(allNumbers) <- c("PFS", "OS", "Recruited", "Censored", "Ongoing")
    colnames(allNumbers) <- paste0("Timepoint: ", round(timeP, 2))
    allNumbers
  })
  names(allNumbers) <- byVar
  allNumbers
}

# survPFS ----

#' PFS Survival Function for Different Transition Models
#'
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
#' @param t (`numeric`)\cr time at which the value of the PFS survival function is to be computed.
#'
#' @return The value of the survival function for the specified transition and time.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' survPFS(transition, 0.4)
survPFS <- function(transition, t) {
  UseMethod("survPFS")
}

#' @describeIn survPFS Survival Function for an exponential transition model.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' survPFS(transition, 0.4)
survPFS.ExponentialTransition <- function(transition, t) {
  ExpSurvPFS(t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02)
}

#' @describeIn survPFS Survival Function for a Weibull transition model.
#' @export
#'
#' @examples
#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
#' survPFS(transition, 0.4)
survPFS.WeibullTransition <- function(transition, t) {
  WeibSurvPFS(
    t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02,
    p01 = transition$weibull_rates$p01, p02 = transition$weibull_rates$p02
  )
}

#' @describeIn survPFS Survival Function for a piecewise constant transition model.
#' @export
#'
#' @examples
#' transition <- piecewise_exponential(
#'   h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1),
#'   pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9)
#' )
#' survPFS(transition, 0.4)
survPFS.PWCTransition <- function(transition, t) {
  PWCsurvPFS(t,
    h01 = transition$hazards$h01, h02 = transition$hazards$h02,
    pw01 = transition$intervals$pw01, pw02 = transition$intervals$pw02
  )
}

# survOS ----

#' OS Survival Function for Different Transition Models
#'
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
#' @param t (`numeric`)\cr time at which the value of the OS survival function is to be computed.
#'
#' @return The value of the survival function for the specified transition and time.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' survOS(transition, 0.4)
survOS <- function(transition, t) {
  UseMethod("survOS")
}

#' @describeIn survOS Survival Function for an exponential transition model.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' survOS(transition, 0.4)
survOS.ExponentialTransition <- function(transition, t) {
  ExpSurvOS(
    t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02,
    h12 = transition$hazards$h12
  )
}

#' @describeIn survOS Survival Function for a Weibull transition model.
#' @export
#'
#' @examples
#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
#' survOS(transition, 0.4)
survOS.WeibullTransition <- function(transition, t) {
  WeibSurvOS(
    t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02,
    h12 = transition$hazards$h12, p01 = transition$weibull_rates$p01,
    p02 = transition$weibull_rates$p02, p12 = transition$weibull_rates$p12
  )
}

#' @describeIn survOS Survival Function for a piecewise constant transition model.
#' @export
#'
#' @examples
#' transition <- piecewise_exponential(
#'   h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1),
#'   pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9)
#' )
#' survOS(transition, 0.4)
survOS.PWCTransition <- function(transition, t) {
  PWCsurvOS(
    t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02,
    h12 = transition$hazards$h12, pw01 = transition$intervals$pw01,
    pw02 = transition$intervals$pw02, pw12 = transition$intervals$pw12
  )
}

# expval ----

#' Helper Function for Computing E(PFS^2)
#'
#' @param x (`numeric`)\cr variable of integration.
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
#'
#' @return Numeric results of the integrand used to calculate E(PFS^2).
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' expvalPFSInteg(0.4, transition)
expvalPFSInteg <- function(x, transition) {
  x * survPFS(transition, x)
}

#' Helper Function for Computing E(OS^2)
#'
#' @param x (`numeric`)\cr variable of integration.
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
#'
#' @return Numeric results of the integrand used to calculate E(OS^2).
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' expvalOSInteg(0.4, transition)
expvalOSInteg <- function(x, transition) {
  x * survOS(transition = transition, t = x)
}

# p11 ----

#' Helper Function for `log_p11()`
#'
#' @param x (`numeric`)\cr variable of integration.
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
#'
#' @return Hazard rate at the specified time for the transition from progression to death.
#' @keywords internal
p11Integ <- function(x, transition) {
  haz(transition = transition, t = x, trans = 3)
}

#' Probability of Remaining in Progression Between Two Time Points for Different Transition Models
#'
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
#' @param s (`numeric`)\cr lower time points.
#' @param t (`numeric`)\cr higher time points.
#' @return This returns the natural logarithm of the probability of remaining in progression (state 1)
#' between two time points, conditional on being in state 1 at the lower time point.
#'
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' log_p11(transition, 1, 3)
log_p11 <- function(transition, s, t) {
  assert_numeric(s, finite = TRUE, any.missing = FALSE, lower = 0)
  assert_numeric(t, finite = TRUE, any.missing = FALSE, lower = 0)
  assert_true(identical(length(s), length(t)))
  assert_true(all(t > s))

  intval <- mapply(function(s, t) {
    stats::integrate(p11Integ,
      lower = s,
      upper = t,
      transition
    )$value
  }, s, t)
  -intval
}

# PFSOS ----

#' Helper Function for `survPFSOS()`
#'
#' @param u (`numeric`)\cr variable of integration.
#' @param t (`numeric`)\cr time at which the value of the PFS*OS survival function is to be computed.
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
#'
#' @note Not all vectors `u` and `t` work here due to assertions in [log_p11()].
#'
#' @return Numeric result of the integrand used to calculate the PFS*OS survival function.
#' @keywords internal
PFSOSInteg <- function(u, t, transition) {
  exp(log_p11(transition, u, t / u) + log(survPFS(transition, u)) + log(haz(transition, u, 1)))
}

#' Survival Function of the Product PFS*OS for Different Transition Models
#'
#' @param t (`numeric`)\cr time at which the value of the PFS*OS survival function is to be computed.
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
#'
#' @return This returns the value of PFS*OS survival function at time t.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' survPFSOS(0.4, transition)
survPFSOS <- function(t, transition) {
  sapply(t, function(x) {
    intval <- stats::integrate(PFSOSInteg, lower = 0, upper = sqrt(x), x, transition)$value
    survPFS(transition, sqrt(x)) + intval
  })
}

# correlation ----

#' Correlation of PFS and OS event times for Different Transition Models
#'
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
#'
#' @return The correlation of PFS and OS.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' corTrans(transition)
corTrans <- function(transition) {
  # E(PFS) & E(OS).
  expvalPFS <- stats::integrate(survPFS,
    lower = 0, upper = Inf,
    transition = transition
  )$value

  expvalOS <- stats::integrate(survOS,
    lower = 0, upper = Inf,
    transition = transition
  )$value

  # Var(PFS) & Var(OS).
  expvalPFS2 <- 2 * stats::integrate(expvalPFSInteg,
    lower = 0, upper = Inf,
    transition = transition
  )$value

  expvalOS2 <- 2 * stats::integrate(expvalOSInteg,
    lower = 0, upper = Inf,
    transition = transition
  )$value

  varPFS <- expvalPFS2 - expvalPFS^2

  varOS <- expvalOS2 - expvalOS^2

  # E(PFS*OS).
  expvalPFSOS <- stats::integrate(survPFSOS,
    lower = 0, upper = Inf,
    transition
  )$value

  # Cor(PFS, OS).
  (expvalPFSOS - expvalPFS * expvalOS) / sqrt(varPFS * varOS)
}

#' Correlation of PFS and OS event times for data from the IDM
#'
#' @param data (`data.frame`)\cr in the format produced by [getOneClinicalTrial()].
#' @param transition (`TransitionParameters` object)\cr specifying the assumed distribution of transition hazards.
#'   Initial parameters for optimization can be specified here.
#'   See [exponential_transition()] or [weibull_transition()] for details.
#' @param bootstrap (`flag`)\cr if `TRUE` computes confidence interval via bootstrap.
#' @param bootstrap_n (`count`)\cr number of bootstrap samples.
#' @param conf_level (`proportion`)\cr confidence level for the confidence interval.
#'
#' @return The correlation of PFS and OS.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' data <- getClinicalTrials(
#'   nRep = 1, nPat = c(100), seed = 1234, datType = "1rowTransition",
#'   transitionByArm = list(transition), dropout = list(rate = 0.5, time = 12),
#'   accrual = list(param = "intensity", value = 7)
#' )[[1]]
#' corPFSOS(data, transition = exponential_transition(), bootstrap = FALSE)
#' \dontrun{
#' corPFSOS(data, transition = exponential_transition(), bootstrap = TRUE)
#' }
corPFSOS <- function(data, transition, bootstrap = TRUE, bootstrap_n = 100, conf_level = 0.95) {
  assert_data_frame(data)
  assert_flag(bootstrap)
  assert_count(bootstrap_n)
  assert_number(conf_level, lower = 0.01, upper = 0.999)

  trans <- estimateParams(data, transition)
  res <- list("corPFSOS" = corTrans(trans))
  if (bootstrap) {
    oplan <- future::plan(future::multisession, workers = parallelly::availableCores(omit = 1))
    on.exit(future::plan(oplan), add = TRUE)
    ids <- lapply(1:bootstrap_n, function(x) sample(seq_len(nrow(data)), nrow(data), replace = TRUE))
    corBootstrap <- furrr::future_map_dbl(ids, ~ {
      furrr::furrr_options(
        globals = list(data = data, transition = transition),
        packages = c("simIDM")
      )
      b_sample <- data[.x, , drop = FALSE]
      b_transition <- estimateParams(b_sample, transition)
      corTrans(b_transition)
    })
    lowerQuantile <- (1 - conf_level) / 2
    upperQuantile <- lowerQuantile + conf_level
    c(stats::quantile(corBootstrap, lowerQuantile),
      "corPFSOS" = res,
      stats::quantile(corBootstrap, upperQuantile)
    )
    res$lower <- stats::quantile(corBootstrap, lowerQuantile)
    res$upper <- stats::quantile(corBootstrap, upperQuantile)
  }
  res
}

#' Log-Rank Test for a Single Trial
#'
#' This function evaluates the significance of either PFS or OS endpoints in a trial,
#' based on a pre-specified critical value.
#'
#' @param data (`data.frame`)\cr data frame containing entry and exit times of an
#' illness-death model. See [getSimulatedData()] for details.
#' @param typeEvent (`string`)\cr endpoint to be evaluated, possible values are `PFS` and `OS`.
#' @param critical (positive `number`)\cr critical value of the log-rank test.
#'
#' @return Logical value indicating log-rank test significance.
#' @export
#'
#' @examples
#' transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
#' transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
#' simTrial <- getClinicalTrials(
#'   nRep = 1, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
#'   transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
#'   accrual = list(param = "intensity", value = 7)
#' )[[1]]
#' logRankTest(data = simTrial, typeEvent = "OS", critical = 3.4)
logRankTest <- function(data, typeEvent = c("PFS", "OS"), critical) {
  # Must be have the format of one row per patient (`datType` must be 1rowPatient`
  # in getClinicalTrials()), i.e. have 10 (if censored) or 11 columns.
  assert_data_frame(data, min.cols = 10, max.cols = 11)
  typeEvent <- match.arg(typeEvent)
  assert_positive_number(critical)

  if (typeEvent == "OS") {
    time <- data$OStime
    event <- data$OSevent
  } else if (typeEvent == "PFS") {
    time <- data$PFStime
    event <- data$PFSevent
  }

  logRank <- survival::survdiff(survival::Surv(time, event) ~ trt, data)
  sqrt(logRank$chisq) > critical
}


#' Helper function to conduct log-rank tests for either PFS or OS
#'
#' This function evaluates the significance of either PFS or OS endpoints for each trial
#' in a list of trials, based on a pre-specified critical value.
#'
#' @param simTrials (`list`)\cr simulated trial data sets, see [getClinicalTrials()].
#' @param typeEvent (`string`)\cr endpoint to be evaluated, possible values are `PFS` and `OS`.
#' @param eventNum (`integer`)\cr number of events required to trigger analysis.
#' @param critical (positive `number`)\cr critical value of the log-rank test.
#'
#' @return Logical vector indicating log-rank test significance for each trial.
#' @export
#'
#' @examples
#' transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
#' transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
#' simTrials <- getClinicalTrials(
#'   nRep = 50, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
#'   transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
#'   accrual = list(param = "intensity", value = 7)
#' )
#' passedLogRank(simTrials = simTrials, typeEvent = "PFS", eventNum = 300, critical = 2.4)
#' @keywords internal
passedLogRank <- function(simTrials, typeEvent, eventNum, critical) {
  assert_list(simTrials, null.ok = FALSE)
  assert_positive_number(critical)
  assert_count(eventNum, positive = TRUE)

  # Censor simulated trials at time-point of OS/PFS analysis.
  trialsAna <- lapply(
    X = simTrials,
    FUN = censoringByNumberEvents,
    eventNum = eventNum,
    typeEvent = typeEvent
  )
  # Compute log-rank test for all trials for OS/PFS.
  unlist(lapply(
    X = trialsAna,
    FUN = logRankTest,
    typeEvent = typeEvent,
    critical = critical
  ))
}

#' Empirical Significance for a List of Simulated Trials
#'
#' This function computes four types of empirical significance — PFS, OS, at-least (significant
#' in at least one of PFS/OS), and joint (significant in both PFS and OS) — using
#' the log-rank test. Empirical significance is calculated as the proportion of significant
#' results in simulated trials, each ending when a set number of PFS/OS events occur.
#' Critical values for PFS and OS test significance must be specified. If trials simulate equal
#' transition hazards across groups (H0), empirical significance estimates type I error;
#' if they simulate differing transition hazards (H1), it estimates power.
#'
#' @param simTrials (`list`)\cr simulated trial data sets, see [getClinicalTrials()].
#' @param criticalPFS (positive `number`)\cr critical value of the log-rank test for PFS.
#' @param criticalOS (positive `number`)\cr critical value of the log-rank test for OS.
#' @param eventNumPFS (`integer`)\cr number of PFS events required to trigger PFS analysis.
#' @param eventNumOS (`integer`)\cr number of OS events required to trigger OS analysis.
#'
#' @return This returns values of four measures of empirical significance.
#' @export
#'
#' @examples
#' transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
#' transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
#' simTrials <- getClinicalTrials(
#'   nRep = 50, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
#'   transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
#'   accrual = list(param = "intensity", value = 7)
#' )
#' empSignificant(
#'   simTrials = simTrials, criticalPFS = 2.4, criticalOS = 2.2,
#'   eventNumPFS = 300, eventNumOS = 500
#' )
empSignificant <- function(simTrials, criticalPFS, criticalOS, eventNumPFS, eventNumOS) {
  assert_list(simTrials, null.ok = FALSE)
  assert_positive_number(criticalPFS)
  assert_positive_number(criticalOS)
  assert_count(eventNumPFS, positive = TRUE)
  assert_count(eventNumOS, positive = TRUE)

  nRep <- length(simTrials)
  simTrials <- lapply(
    X = simTrials,
    FUN = function(x) if (ncol(x) == 9) getDatasetWideFormat(x) else x
  )

  # Which trials passed the log-rank test for PFS/OS?
  passedLogRankPFS <- passedLogRank(
    simTrials = simTrials,
    typeEvent = "PFS",
    eventNum = eventNumPFS,
    critical = criticalPFS
  )
  passedLogRankOS <- passedLogRank(
    simTrials = simTrials,
    typeEvent = "OS",
    eventNum = eventNumOS,
    critical = criticalOS
  )

  # Empirical significance is the fraction of trials with significant log-rank test:
  significantPFS <- sum(passedLogRankPFS) / nRep
  significantOS <- sum(passedLogRankOS) / nRep

  # Derived measures of significance:
  sumPassed <- passedLogRankPFS + passedLogRankOS
  # At least one of PFS/OS has significant log-rank test.
  significantAtLeastOne <- sum(sumPassed >= 1) / nRep
  # Both PFS/OS have significant log-rank tests.
  significantBoth <- sum(sumPassed == 2) / nRep

  list(
    "significantPFS" = significantPFS,
    "significantOS" = significantOS,
    "significantAtLeastOne" = significantAtLeastOne,
    "significantBoth" = significantBoth
  )
}

#' Transitions from the Intermediate State to the Absorbing State
#'
#' This function creates transition entry and exit times from the intermediate state to the absorbing state
#' for an existing data frame containing the exit times out of the initial state.
#'
#' @param simDataOne (`data.frame`)\cr a data frame containing all patients with transitions
#'  into the intermediate state. See [getSimulatedData()] for details.
#' @param transition (`TransitionParameters`)\cr transition parameters comprising
#'  `hazards`, corresponding `intervals` and `weibull_rates`, see [exponential_transition()], [piecewise_exponential()]
#'  and [weibull_transition()] for details.
#'
#' @return This returns a data frame with one row per patient for the second transition,
#'  i.e. the transition out of the intermediate
#'  state. This is a helper function of [getSimulatedData()].
#'
#' @export
#'
#' @examples
#' simDataOne <- data.frame(
#'   id = c(1:3), to = c(1, 1, 1), from = c(0, 0, 0), entry = c(0, 0, 0),
#'   exit = c(3, 5.6, 7.2), censTime = c(6.8, 5.9, 9.4)
#' )
#' transition <- exponential_transition(1, 1.6, 0.3)
#' getOneToTwoRows(simDataOne, transition)
getOneToTwoRows <- function(simDataOne, transition) {
  assert_data_frame(simDataOne, ncols = 6)
  assert_class(transition, "TransitionParameters")

  id1 <- simDataOne$id
  N1 <- nrow(simDataOne)
  U1 <- stats::runif(N1)
  to1 <- rep(2, N1)
  from1 <- rep(1, N1)
  entry1 <- simDataOne$exit
  h12 <- transition$hazard$h12
  p12 <- transition$weibull_rates$p12
  pw12 <- transition$intervals$pw12

  # Create waiting time in state 1.
  wait_time1 <- if (transition$family == "exponential") {
    (-log(1 - U1)) / (h12)
  } else if (transition$family == "Weibull") {
    getWaitTimeSum(U = U1, haz1 = h12, haz2 = 0, p1 = p12, p2 = 1, entry = entry1)
  } else if (transition$family == "piecewise exponential") {
    getPCWDistr(U = U1, haz = h12, pw = pw12, t_0 = entry1)
  }
  exit1 <- entry1 + wait_time1

  # Add censoring.
  censTime1 <- simDataOne$censTime
  to1 <- ifelse(censTime1 < exit1, "cens", to1)
  exit1 <- pmin(censTime1, exit1)

  data.frame(
    id = id1,
    from = from1,
    to = to1,
    entry = entry1,
    exit = exit1,
    censTime = censTime1,
    stringsAsFactors = FALSE
  )
}

#' Staggered Study Entry
#'
#' This function adds staggered study entry times to a simulated data set with illness-death model structure.
#'
#' @param simData (`data.frame`)\cr simulated data frame containing entry and exit times
#'   at individual study time scale. See [getSimulatedData()] for details.
#' @param N (`int`)\cr number of patients.
#' @param accrualParam (`string`)\cr possible values are 'time' or 'intensity'.
#' @param accrualValue  (`number`)\cr specifies the accrual intensity. For `accrualParam` equal time,
#'   it describes the length of the accrual period. For `accrualParam` equal intensity, it describes
#'   the number of patients recruited per time unit.  If `accrualValue` is equal to 0,
#'   all patients start at calendar time 0
#'   in the initial state.
#'
#' @return This returns a data set containing a single simulated study containing accrual times,
#'  i.e. staggered study entry.
#'  This is a helper function of [getSimulatedData()].
#'
#' @export
#'
#' @examples
#' simData <- data.frame(
#'   id = c(1, 1, 2, 3), from = c(0, 1, 0, 0), to = c(1, 2, "cens", 2),
#'   entry = c(0, 3, 0, 0),
#'   exit = c(3, 5.3, 5.6, 7.2), censTime = c(6.8, 6.8, 5.6, 9.4)
#' )
#' addStaggeredEntry(simData, 3, accrualParam = "time", accrualValue = 5)
addStaggeredEntry <- function(simData, N, accrualParam, accrualValue) {
  assert_choice(accrualParam, c("time", "intensity"))
  assert_number(accrualValue, lower = 0)

  # Get accrual times in calendar time per individual.
  # If no staggered study entry is present, all individuals have the same entry time 0.
  entry_act <- if (accrualValue != 0) {
    if (accrualParam == "time") {
      stats::runif(N, 0, accrualValue)
    } else if (accrualParam == "intensity") {
      accrualTime <- N / accrualValue
      stats::runif(N, 0, accrualTime)
    }
  } else if (accrualValue == 0) {
    rep(0, N)
  }
  entryAct <- cbind(id = 1:N, entry_act)
  # Combine simulated data with actual entry time.
  # Generate actual times  (actual entry time + individual time).
  simData <- merge(simData, entryAct)
  simData$entryAct <- simData$entry + simData$entry_act
  simData$exitAct <- simData$exit + simData$entry_act
  simData$censAct <- simData$censTime + simData$entry_act
  # Delete temporary helper columns.
  simData$entry_act <- NULL
  simData[, c("censTime")] <- list(NULL)
  simData
}

#' Simulate Data Set from an Illness-Death Model
#'
#' This function creates a single simulated data set for a single treatment arm. It simulates data
#'  from an illness-death model with one row per transition and subject.
#'
#' @param N (`int`)\cr number of patients.
#' @param transition (`TransitionParameters`)\cr transition parameters comprising
#'   `hazards`, corresponding `intervals` and `weibull_rates`, see [exponential_transition()], [piecewise_exponential()]
#'   and [weibull_transition()] for details.
#' @param dropout (`list`)\cr specifies drop-out probability.
#'   Random censoring times are generated using exponential distribution. `dropout$rate` specifies
#'   the drop-out probability per `dropout$time` time units.
#'   If `dropout$rate` is equal to 0, then no censoring is applied.
#' @param accrual (`list`)\cr specifies accrual intensity. See [addStaggeredEntry()] for details.
#'
#' @return This returns a data frame with one row per transition per individual.
#' @details The output data set contains the following columns:
#' - id (`integer`): patient id.
#' - from (`numeric`): starting state of the transition.
#' - to (`character`): final state of the transition.
#' - entry (`numeric`): entry time of the transition on the individual time scale.
#' - exit (`numeric`): exit time of the transition on the individual time scale.
#' - entryAct (`numeric`): entry time of the transition on study time scale.
#' - exitAct (`numeric`):  exit time  of the transition on study time scale.
#' - censAct (`numeric`):  censoring time of the individual on study time scale.
#' @export
#'
#' @examples
#' getSimulatedData(
#'   N = 10,
#'   transition = exponential_transition(h01 = 1, h02 = 1.5, h12 = 1),
#'   dropout = list(rate = 0.3, time = 1),
#'   accrual = list(param = "time", value = 5)
#' )
getSimulatedData <- function(N,
                             transition = exponential_transition(h01 = 1, h02 = 1, h12 = 1),
                             dropout = list(rate = 0, time = 12),
                             accrual = list(param = "time", value = 0)) {
  # Check input parameters.
  assert_int(N, lower = 1L)
  assert_class(transition, "TransitionParameters")
  assert_list(dropout)

  assert(check_number(dropout$rate, lower = 0, upper = 1),
    check_number(dropout$time),
    check_true(dropout$time > 0, ),
    combine = "and", .var.name = "dropout"
  )
  assert_list(accrual)
  # Initialize transition hazards, Weibull rates and intervals.
  h <- transition$hazard
  p <- transition$weibull_rates
  pw <- transition$intervals

  # Get rate parameter for exponential distributed censoring times.
  # if rate=0, no censoring is applied.
  censRate <- if (dropout$rate > 0) {
    -log(1 - dropout$rate) / dropout$time
  } else {
    0
  }
  # Censoring times for all individuals, infinity if no censoring is applied.
  censTime <- if (dropout$rate > 0) {
    stats::rexp(N, censRate)
  } else {
    rep(Inf, N)
  }

  # All individuals start in the initial state.
  entry <- rep(0, N)
  from <- rep(0, N)
  # Waiting time in the initial state 0.
  U <- stats::runif(N)
  # Exponential or Weibull distributed survival times.
  if (transition$family %in% c("exponential", "Weibull")) {
    # For distribution of waiting time in the initial state the all-cause hazard (h01+h02) is needed.
    wait_time <- getWaitTimeSum(U, h$h01, h$h02, p$p01, p$p02, entry)
    # A binomial experiment decides on death or progression.
    numerator <- p$p01 * h$h01 * wait_time^(p$p01 - 1) # nolint
    denumerator <- numerator + p$p02 * h$h02 * wait_time^(p$p02 - 1) # nolint

    # Piecewise exponential distributed survival times.
  } else if (transition$family == "piecewise exponential") {
    Sum_PCW <- getSumPCW(h$h01, h$h02, pw$pw01, pw$pw02)
    wait_time <- getPCWDistr(U, Sum_PCW$hazards, Sum_PCW$intervals, entry)
    # A binomial experiment decides on death or progression.
    numerator <- getPWCHazard(h$h01, pw$pw01, wait_time)
    denumerator <- getPWCHazard(h$h01, pw$pw01, wait_time) +
      getPWCHazard(h$h02, pw$pw02, wait_time)
  }

  to_prob <- stats::rbinom(N, 1, numerator / denumerator)
  to <- ifelse(to_prob == 0, 2, 1)
  exit <- wait_time

  # Add censoring.
  to <- ifelse(censTime < wait_time, "cens", to)
  exit <- pmin(censTime, wait_time)

  simData <- data.frame(
    id = 1:N, from = from, to = to, entry = entry, exit = exit,
    censTime = censTime, stringsAsFactors = FALSE
  )
  is_intermediate_state <- simData$to == 1
  if (any(is_intermediate_state)) {
    # Add 1 -> 2 transition only for patients that are in the intermediate state 1.
    simDataOne <- simData[is_intermediate_state, ]
    newrows <- getOneToTwoRows(simDataOne, transition)
    simData <- rbind(simData, newrows)
    simData <- simData[order(simData$id), ]
  }
  # Add staggered study entry, i.e. study entry at calendar time.
  simData <- addStaggeredEntry(
    simData = simData,
    N = N,
    accrualParam = accrual$param,
    accrualValue = accrual$value
  )
  simData$to <- as.character(simData$to)
  simData
}

#' Preparation of a Data Set to Compute Log-likelihood
#'
#' @param data (`data.frame`)\cr containing entry and exit times of an illness-death model.
#'   See [getOneClinicalTrial()] for details.
#'
#' @return This function returns a data set with one row per patient and transition, when the patient is at risk.
#' @export
#'
#' @details
#' The output data set contains the following columns:
#' - id (`integer`): patient id.
#' - from (`integer`): start event state.
#' - to (`integer`): end event state.
#' - trans (`integer`): transition (1, 2 or 3) identifier
#'   - `1`: Transition from state 0 (stable) to 1 (progression).
#'   - `2`: Transition from state 0 (stable) to 2 (death).
#'   - `3`: Transition from state 1 (progression) to 2 (death).
#' - entry (`numeric`): time at which the patient begins to be at risk for the transition.
#' - exit (`numeric`): time at which the patient ends to be at risk for the transition.
#' - status (`logical`): event indicator for the transition.
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' simData <- getOneClinicalTrial(
#'   nPat = c(30), transitionByArm = list(transition),
#'   dropout = list(rate = 0.8, time = 12),
#'   accrual = list(param = "time", value = 1)
#' )
#' prepareData(simData)
prepareData <- function(data) {
  assert_data_frame(data, min.cols = 9, max.cols = 11)
  colNames <- c(
    "id", "trt", "PFStime", "CensoredPFS", "PFSevent", "OStime",
    "CensoredOS", "OSevent", "recruitTime",
    "OStimeCal", "PFStimeCal"
  )
  if (!all(names(data) %in% colNames)) {
    data <- getDatasetWideFormat(data)
  }

  # Transform simIDM trial data to log-likelihood-compatible format.
  # Suppress warning about how msprep handles 1 -> 3 transitions.
  dataNew <- suppressWarnings(mstate::msprep(
    time = c("recruitTime", "PFStime", "OStime"),
    status = c("trt", "PFSevent", "OSevent"),
    data = data,
    trans = mstate::trans.illdeath(),
    id = data$id
  ))
  cols <- which(names(dataNew) %in% c("Tstart", "Tstop"))
  names(dataNew)[cols] <- c("entry", "exit")
  # Correct msprep results for uncensored PFS=OS events.
  ids <- data$id[data$PFStimeCal == data$OStimeCal & data$CensoredPFS == 0]
  dataNew <- dataNew[!(dataNew$id %in% ids & dataNew$trans == 3), ]
  dataNew$status[dataNew$id %in% ids] <- abs(dataNew$status[dataNew$id %in% ids] - 1)

  as.data.frame(dataNew[, -which(names(dataNew) == "time")], row.names = seq_len(nrow(dataNew)))
}

#' Compute the Negative Log-Likelihood for a Given Data Set and Transition Model
#'
#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr
#'   see [exponential_transition()] or [weibull_transition()] for details.
#' @param data (`data.frame`)\cr in the format created by [prepareData()].
#'
#' @return The value of the negative log-likelihood.
#' @export
#'
#' @details
#' Calculates the negative log-likelihood for a given data set and transition model. It uses the hazard
#' and survival functions specific to the transition model.
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' simData <- getOneClinicalTrial(
#'   nPat = c(30), transitionByArm = list(transition),
#'   dropout = list(rate = 0.8, time = 12),
#'   accrual = list(param = "time", value = 1)
#' )
#' negLogLik(transition, prepareData(simData))
negLogLik <- function(transition, data) {
  with(
    data,
    -sum(log(
      haz(transition, exit, trans)^status * survTrans(transition, exit, trans) /
        survTrans(transition, entry, trans)
    ))
  )
}

#' Hazard Function for Different Transition Models
#'
#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr
#'   see [exponential_transition()] or [weibull_transition()] for details.
#' @param t (`numeric`)\cr time at which hazard is to be computed.
#' @param trans (`integer`)\cr index specifying the transition type.
#'
#' @details
#' The transition types are:
#' - `1`: Transition from state 0 (stable) to 1 (progression).
#' - `2`: Transition from state 0 (stable) to 2 (death).
#' - `3`: Transition from state 1 (progression) to 2 (death).
#'
#' @return The hazard rate for the specified transition and time.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' haz(transition, 0.4, 2)
haz <- function(transition, t, trans) {
  assert_class(transition, "TransitionParameters")
  assert_numeric(t, lower = 0)
  assert_subset(trans, c(1, 2, 3))
  UseMethod("haz")
}

#' @describeIn haz for an exponential transition model.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' haz(transition, 0.4, 2)
haz.ExponentialTransition <- function(transition, t, trans) {
  # params (in this order): h01, h02, h12.
  params <- unlist(transition$hazards, use.names = FALSE)
  rep(params[trans], length(t) / length(trans))
}

#' @describeIn haz for the Weibull transition model.
#' @export
#'
#' @examples
#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
#' haz(transition, 0.4, 2)
haz.WeibullTransition <- function(transition, t, trans) {
  # params (in this order): h01, h02, h12, p01, p02, p12.
  params <- c(unlist(transition$hazards, use.names = FALSE), unlist(transition$weibull_rates, use.names = FALSE))
  params[trans] * params[trans + 3] * t^(params[trans + 3] - 1)
}

#' @describeIn haz for the piecewise constant transition model.
#' @export
#'
#' @examples
#' transition <- piecewise_exponential(
#'   h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1),
#'   pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9)
#' )
#' haz(transition, 6, 2)
haz.PWCTransition <- function(transition, t, trans) {
  getPWCHazard(unlist(transition$hazards[trans], use.names = FALSE),
    unlist(transition$intervals[trans], use.names = FALSE),
    x = t
  )
}

#' Survival Function for Different Transition Models
#'
#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr
#'   see [exponential_transition()] or [weibull_transition()] for details.
#' @param t (`numeric`)\cr time at which survival probability is to be computed.
#' @param trans (`integer`)\cr index specifying the transition type.
#'
#' @return The survival probability for the specified transition and time.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' survTrans(transition, 0.4, 2)
survTrans <- function(transition, t, trans) {
  assert_class(transition, "TransitionParameters")
  assert_numeric(t, lower = 0)
  assert_subset(trans, c(1, 2, 3))
  UseMethod("survTrans")
}

#' @describeIn survTrans for the Exponential Transition Model
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' survTrans(transition, 0.4, 2)
survTrans.ExponentialTransition <- function(transition, t, trans) {
  # params (in this order): h01, h02, h12.
  params <- unlist(transition$hazards, use.names = FALSE)
  exp(-params[trans] * t)
}

#' @describeIn survTrans for the Weibull Transition Model
#' @export
#'
#' @examples
#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
#' survTrans(transition, 0.4, 2)
survTrans.WeibullTransition <- function(transition, t, trans) {
  # params (in this order): h01, h02, h12, p01, p02, p12.
  params <- c(unlist(transition$hazards, use.names = FALSE), unlist(transition$weibull_rates, use.names = FALSE))
  exp(-params[trans] * t^params[trans + 3])
}

#' Retrieve Initial Parameter Vectors for Likelihood Maximization
#'
#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr containing the initial parameters.
#'   See [exponential_transition()] or [weibull_transition()] for details.
#'
#' @return The numeric vector of initial parameters for likelihood maximization.
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' getInit(transition)
getInit <- function(transition) {
  assert_class(transition, "TransitionParameters")
  UseMethod("getInit")
}

#' @describeIn getInit for the Exponential Transition Model
#' @export
#'
#' @examples
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' getInit(transition)
getInit.ExponentialTransition <- function(transition) {
  unlist(transition$hazards, use.names = FALSE)
}

#' @describeIn getInit for the Weibull Transition Model
#' @export
#'
#' @examples
#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
#' getInit(transition)
getInit.WeibullTransition <- function(transition) {
  c(unlist(transition$hazards, use.names = FALSE), unlist(transition$weibull_rates, use.names = FALSE))
}

#' Generate the Target Function for Optimization
#'
#' @param transition (`TransitionParameters`)\cr
#'   specifying the distribution family. See [exponential_transition()] or [weibull_transition()] for details.
#'
#' @return Function that calculates the negative log-likelihood for the given parameters.
#' @export
#'
#' @details
#' This function creates a target function for optimization, computing the negative log-likelihood for given
#' parameters, data, and transition model type.
#'
#' @examples
#' transition <- exponential_transition(2, 1.3, 0.8)
#' simData <- getOneClinicalTrial(
#'   nPat = c(30), transitionByArm = list(transition),
#'   dropout = list(rate = 0.8, time = 12),
#'   accrual = list(param = "time", value = 1)
#' )
#' params <- c(1.2, 1.5, 1.6) # For ExponentialTransition
#' data <- prepareData(simData)
#' transition <- exponential_transition()
#' fun <- getTarget(transition)
#' fun(params, data)
getTarget <- function(transition) {
  assert_class(transition, "TransitionParameters")
  UseMethod("getTarget", transition)
}

#' @describeIn getTarget for the Exponential Transition Model
#' @export
#'
#' @examples
#' transition <- exponential_transition(2, 1.3, 0.8)
#' simData <- getOneClinicalTrial(
#'   nPat = c(30), transitionByArm = list(transition),
#'   dropout = list(rate = 0.8, time = 12),
#'   accrual = list(param = "time", value = 1)
#' )
#' params <- c(1.2, 1.5, 1.6)
#' data <- prepareData(simData)
#' transition <- exponential_transition()
#' target <- getTarget(transition)
#' target(params, data)
getTarget.ExponentialTransition <- function(transition) {
  function(params, data) {
    negLogLik(transition = exponential_transition(h01 = params[1], h02 = params[2], h12 = params[3]), data = data)
  }
}

#' @describeIn getTarget for the Weibull Transition Model
#' @export
#'
#' @examples
#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
#' simData <- getOneClinicalTrial(
#'   nPat = c(30), transitionByArm = list(transition),
#'   dropout = list(rate = 0.8, time = 12),
#'   accrual = list(param = "time", value = 1)
#' )
#' params <- c(1.2, 1.5, 1.6, 0.8, 1.3, 1.1)
#' data <- prepareData(simData)
#' transition <- weibull_transition()
#' target <- getTarget(transition)
#' target(params, data)
getTarget.WeibullTransition <- function(transition) {
  function(params, data) {
    negLogLik(transition = weibull_transition(
      h01 = params[1], h02 = params[2], h12 = params[3],
      p01 = params[4], p02 = params[5], p12 = params[6]
    ), data = data)
  }
}

#' Format Results of Parameter Estimation for Different Transition Models
#'
#' @param transition (`TransitionParameters`)\cr
#'   see [exponential_transition()] or [weibull_transition()] for details.
#' @param res (`numeric` vector)\cr vector of parameter estimates from the likelihood maximization procedure.
#'
#' @return Returns a `TransitionParameters` object with parameter estimates.
#' @export
#'
#' @examples
#' results <- c(1.2, 1.5, 1.6)
#' getResults(exponential_transition(), results)
getResults <- function(transition, res) {
  UseMethod("getResults")
}

#' @describeIn getResults for the Exponential Transition Model
#' @export
#'
#' @examples
#' results <- c(1.2, 1.5, 1.6)
#' getResults(exponential_transition(), results)
getResults.ExponentialTransition <- function(transition, res) {
  exponential_transition(h01 = res[1], h02 = res[2], h12 = res[3])
}

#' @describeIn getResults for the Weibull Transition Model
#' @export
#'
#' @examples
#' results <- c(1.2, 1.5, 1.6, 2, 2.5, 1)
#' getResults(weibull_transition(), results)
getResults.WeibullTransition <- function(transition, res) {
  weibull_transition(
    h01 = res[1], h02 = res[2], h12 = res[3],
    p01 = res[4], p02 = res[5], p12 = res[6]
  )
}

#' Estimate Parameters of the Multistate Model Using Clinical Trial Data
#'
#' @param data (`data.frame`)\cr in the format produced by [getOneClinicalTrial()].
#' @param transition (`TransitionParameters` object)\cr specifying the assumed distribution of transition hazards.
#'   Initial parameters for optimization can be specified here.
#'   See [exponential_transition()] or [weibull_transition()] for details.
#'
#' @return Returns a `TransitionParameters` object with the estimated parameters.
#' @export
#'
#' @details
#' This function estimates parameters for transition models using clinical trial data.
#' The `transition` object can be initialized with starting values for parameter estimation.
#' It uses [stats::optim()] to optimize the parameters.
#'
#' @examples
#' transition <- exponential_transition(h01 = 2, h02 = 1.4, h12 = 1.6)
#' simData <- getOneClinicalTrial(
#'   nPat = c(30), transitionByArm = list(transition),
#'   dropout = list(rate = 0.3, time = 12),
#'   accrual = list(param = "time", value = 1)
#' )
#' # Initialize transition with desired starting values for optimization:
#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' estimate <- estimateParams(simData, transition)
estimateParams <- function(data, transition) {
  data <- prepareData(data)
  par <- getInit(transition)
  target <- getTarget(transition)

  res <- stats::optim(
    par = par,
    fn = target,
    method = "L-BFGS-B",
    lower = 1e-3,
    data = data
  )$par

  getResults(transition, res)
}

#' Simulation of a Single Oncology Clinical Trial
#'
#' This function creates a data set with a single simulated oncology clinical trial with one row per transition
#' based on an illness-death model. Studies with an arbitrary number of treatment arms are possible.
#'
#' @param nPat (`integer`)\cr numbers of patients per treatment arm.
#' @param transitionByArm (`list`) \cr transition parameters for each treatment group.
#'   See [exponential_transition()], [piecewise_exponential()] and [weibull_transition()] for details.
#' @param dropout  dropout (`list`)\cr specifies drop-out probability. See [getSimulatedData()] for details.
#'  Can be specified either as one list that should be applied to all treatment groups or a separate list
#'  for each treatment group.
#' @param accrual  accrual (`list`)\cr specifies accrual intensity. See [addStaggeredEntry()] for details.
#'  Can be specified either as one list that should be applied to all treatment groups or a separate list
#'  for each treatment group.
#'
#' @return This returns a data frame with one simulated clinical trial and multiple treatment arms.
#'   See [getSimulatedData()] for the explanation of the columns. The column `trt` contains the treatment indicator.
#'   This is a helper function of [getClinicalTrials()].
#' @export
#'
#' @examples
#' transition1 <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' transition2 <- exponential_transition(h01 = 1, h02 = 1.3, h12 = 1.7)
#' transition3 <- exponential_transition(h01 = 1.1, h02 = 1, h12 = 1.5)
#' getOneClinicalTrial(
#'   nPat = c(30, 20, 30), transitionByArm = list(transition1, transition2, transition3),
#'   dropout = list(rate = 0, time = 12),
#'   accrual = list(param = "time", value = 0)
#' )
getOneClinicalTrial <- function(nPat, transitionByArm,
                                dropout = list(rate = 0, time = 12),
                                accrual = list(param = "time", value = 0)) {
  assert_list(transitionByArm)
  nPat <- as.integer(nPat)
  nArm <- length(transitionByArm)
  assert_integer(nPat,
    lower = 1,
    any.missing = FALSE,
    all.missing = FALSE,
    len = nArm,
  )
  assert_list(dropout)
  assert_list(accrual)

  # Same accrual and dropout parameters for each group?
  if (is.list(dropout[[1]])) {
    assert_list(dropout, len = nArm, types = "list")
  } else {
    dropout <- rep(list(dropout), nArm)
  }

  if (is.list(accrual[[1]])) {
    assert_list(accrual, len = nArm, types = "list")
  } else {
    accrual <- rep(list(accrual), nArm)
  }

  # Each loop simulates a single trial arm.
  # Starting values for the loop.
  simdata <- NULL
  previousPts <- 0
  for (i in seq_len(nArm)) {
    group <- getSimulatedData(nPat[i], transitionByArm[[i]], dropout[[i]], accrual[[i]])
    group$trt <- i
    group$id <- group$id + previousPts
    simdata <- rbind(simdata, group)
    previousPts <- previousPts + nPat[i]
  }
  if (!any(simdata$from == 1)) {
    warning("no progression to death transitions included in the simulated data")
  }
  simdata
}

#' Conversion of a Data Set from One Row per Transition to One Row per Patient
#'
#' @param data (`data.frame`)\cr data frame containing entry and exit times of an illness-death model.
#'   See [getSimulatedData()] for details.
#'
#' @return This function returns a data set with one row per patient and endpoints PFS and OS.
#' @export
#'
#' @details
#' The output data set contains the following columns:
#' - id (`integer`): patient id.
#' - trt `integer`): treatment id.
#' - PFStime (`numeric`): event time of PFS event.
#' - CensoredPFS (`logical`): censoring indicator for PFS event.
#' - PFSevent (`logical`): event indicator for PFS event.
#' - OStime (`numeric`): event time of OS event.
#' - CensoredOS (`logical`): censoring indicator for OS event.
#' - OSevent (`logical`): event indicator for OS event.
#' - recruitTime (`numeric`): time of recruitment.
#' - OStimeCal (`numeric`): OS event time at calendar time scale.
#' - PFStimeCal (`numeric`): PFS event time at calendar time scale.
#'
#' @examples
#' transition1 <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' transition2 <- exponential_transition(h01 = 1, h02 = 1.3, h12 = 1.7)
#' transition3 <- exponential_transition(h01 = 1.1, h02 = 1, h12 = 1.5)
#' simData <- getOneClinicalTrial(
#'   nPat = c(30, 20, 30), transitionByArm = list(transition1, transition2, transition3),
#'   dropout = list(rate = 0, time = 12),
#'   accrual = list(param = "time", value = 0)
#' )
#' getDatasetWideFormat(simData)
getDatasetWideFormat <- function(data) {
  assert_data_frame(data, ncols = 9)
  assert_subset(c("id", "from", "to", "entry", "exit", "entryAct", "exitAct", "censAct", "trt"), names(data))

  # Recruitment time is the actual entry time of the initial state.
  recruitTime <- subset(data[, c("id", "entryAct")], data$from == 0)
  names(recruitTime)[names(recruitTime) == "entryAct"] <- "recruitTime"

  # The OS time is the entry time into state 2 or the censoring time whatever occurs first.
  OStime <- subset(data[, c("id", "exit")], data$to == 2 | data$to == "cens")
  names(OStime)[names(OStime) == "exit"] <- "OStime"

  # The PFS time is the entry time into state 1 or state 2 or the censoring time whatever occurs first.
  PFStime <- subset(
    data[, c("id", "exit")],
    (data$to == 2 & data$from == 0) | (data$to == 1) | (data$to == "cens" & data$from == 0)
  )
  names(PFStime)[names(PFStime) == "exit"] <- "PFStime"

  # Add censoring indicators.
  censoredIdsPFS <- data$id[data$to == "cens" & data$from == 0]
  censoredIdsOS <- data$id[data$to == "cens"]
  id <- unique(data$id)
  CensoredOS <- cbind(id = id, CensoredOS = as.integer(id %in% censoredIdsOS))
  CensoredPFS <- cbind(id = id, CensoredPFS = as.integer(id %in% censoredIdsPFS))

  # Merge all data sets to one.
  newdata <- unique(data[, c("id", "trt")])
  newdata <- merge(x = newdata, y = PFStime, by = "id")
  newdata <- merge(x = newdata, y = CensoredPFS, by = "id")

  # Do we have an observed PFS event?
  newdata$PFSevent <- abs(1 - newdata$CensoredPFS)
  newdata <- merge(x = newdata, y = OStime, by = "id")
  newdata <- merge(x = newdata, y = CensoredOS, by = "id")

  # Do we have an observed OS event?
  newdata$OSevent <- abs(1 - newdata$CensoredOS)
  newdata <- merge(x = newdata, y = recruitTime, by = "id")

  # Add variables with event times in calendar time.
  newdata$OStimeCal <- newdata$OStime + newdata$recruitTime
  newdata$PFStimeCal <- newdata$PFStime + newdata$recruitTime

  newdata
}

#' Helper Function for Adding Progress Bar to Trial Simulation
#'
#' @param x (`int`)\cr iteration index within lapply.
#' @param ... parameters transferred to [getOneClinicalTrial()], see [getOneClinicalTrial()] for details.
#'
#' @return This returns the same as [getOneClinicalTrial()], but updates the progress bar.
#' @keywords internal
runTrial <- function(x, pb, ...) {
  utils::setTxtProgressBar(pb, x)
  getOneClinicalTrial(...)
}

#' Simulation of a Large Number of Oncology Clinical Trials
#'
#' @param nRep (`int`)\cr number of simulated trials.
#' @param ... parameters transferred to [getOneClinicalTrial()], see [getOneClinicalTrial()] for details.
#' @param seed (`int`)\cr random seed used for this simulation.
#' @param datType (`string`)\cr possible values are `1rowTransition` and `1rowPatient`.
#'
#' @return This function returns a list with `nRep` simulated data sets in the format specified by `datType`.
#'  See [getDatasetWideFormat()] [getOneClinicalTrial()] for details.
#' @export
#'
#' @examples
#' transition1 <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
#' transition2 <- exponential_transition(h01 = 1, h02 = 1.3, h12 = 1.7)
#' getClinicalTrials(
#'   nRep = 10, nPat = c(20, 20), seed = 1234, datType = "1rowTransition",
#'   transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
#'   accrual = list(param = "intensity", value = 7)
#' )
getClinicalTrials <- function(nRep, ..., seed = 1234, datType = "1rowTransition") {
  assert_number(nRep, lower = 1)
  assert_choice(datType, c("1rowTransition", "1rowPatient"))

  set.seed(seed)
  cat("Simulating", nRep, "trials:\n")
  pb <- utils::txtProgressBar(min = 0, max = nRep, style = 3)
  # getOneClinicalTrial generates a single clinical trial with multiple arms. Generate nRep simulated trials:
  simulatedTrials <- lapply(
    seq_len(nRep),
    FUN = function(x, ...) runTrial(x, pb, ...),
    ...
  )
  close(pb)
  # Final data set format: one row per patient or one row per transition?
  if (datType == "1rowPatient") {
    simulatedTrials <- lapply(simulatedTrials, getDatasetWideFormat)
  }
  simulatedTrials
}

#' PFS Survival Function from Constant Transition Hazards
#'
#' @inheritParams WeibSurvOS
#' @return This returns the value of PFS survival function at time t.
#' @export
#'
#' @examples
#' ExpSurvPFS(c(1:5), 0.2, 0.4)
ExpSurvPFS <- function(t, h01, h02) {
  assert_numeric(t, lower = 0, any.missing = FALSE)
  assert_positive_number(h01, zero_ok = TRUE)
  assert_positive_number(h02, zero_ok = TRUE)

  exp(-(h01 + h02) * t)
}

#' OS Survival Function from Constant Transition Hazards
#'
#' @inheritParams WeibSurvOS
#'
#' @return This returns the value of OS survival function at time t.
#' @export
#'
#' @examples
#' ExpSurvOS(c(1:5), 0.2, 0.4, 0.1)
ExpSurvOS <- function(t, h01, h02, h12) {
  assert_numeric(t, lower = 0, any.missing = FALSE)
  assert_positive_number(h01, zero_ok = TRUE)
  assert_positive_number(h02, zero_ok = TRUE)
  assert_positive_number(h12, zero_ok = TRUE)

  h012 <- h12 - h01 - h02
  ExpSurvPFS(t, h01, h02) + h01 * h012^-1 * (ExpSurvPFS(t, h01, h02) - exp(-h12 * t))
}

#' PFS Survival Function from Weibull Transition Hazards
#'
#' @inheritParams WeibSurvOS
#' @return This returns the value of PFS survival function at time t.
#' @export
#' @export
#'
#' @examples
#' WeibSurvPFS(c(1:5), 0.2, 0.5, 1.2, 0.9)
WeibSurvPFS <- function(t, h01, h02, p01, p02) {
  assert_numeric(t, lower = 0, any.missing = FALSE)
  assert_positive_number(h01, zero_ok = TRUE)
  assert_positive_number(h02, zero_ok = TRUE)
  assert_positive_number(p01)
  assert_positive_number(p02)

  exp(-h01 * t^p01 - h02 * t^p02)
}

#' Helper for Efficient Integration
#'
#' @param integrand (`function`)\cr  to be integrated.
#' @param upper (`numeric`)\cr upper limits of integration.
#' @param ...  additional arguments to be passed to `integrand`.
#'
#' @return This function returns for each upper limit the estimates of the integral.
#' @keywords internal
integrateVector <- function(integrand, upper, ...) {
  assert_true(all(upper >= 0))
  boundaries <- sort(unique(upper))
  nIntervals <- length(boundaries)
  intervals <- mapply(
    FUN = function(f, lower, upper, ...) stats::integrate(f, ..., lower, upper)$value,
    MoreArgs = list(f = integrand, ...),
    lower = c(0, boundaries[-nIntervals]),
    upper = boundaries
  )
  cumIntervals <- cumsum(intervals)
  cumIntervals[match(upper, boundaries)]
}

#' Helper Function for `WeibSurvOS()`
#'
#' @param x (`numeric`)\cr  variable of integration.
#' @inheritParams WeibSurvOS
#'
#' @return Numeric results of the integrand used to calculate
#' the OS survival function for Weibull transition hazards, see  `WeibSurvOS()`.
#'
#' @keywords internal
WeibOSInteg <- function(x, t, h01, h02, h12, p01, p02, p12) {
  assert_numeric(x, finite = TRUE, any.missing = FALSE)
  assert_numeric(t, finite = TRUE, any.missing = FALSE)
  assert_true(test_scalar(x) || identical(length(x), length(t)) || test_scalar(t))
  assert_positive_number(h01, zero_ok = TRUE)
  assert_positive_number(h02, zero_ok = TRUE)
  assert_positive_number(h12, zero_ok = TRUE)
  assert_positive_number(p01)
  assert_positive_number(p02)
  assert_positive_number(p12)

  x^(p01 - 1) * exp(-h01 * x^p01 - h02 * x^p02 - h12 * (t^p12 - x^p12))
}

#' OS Survival Function from Weibull Transition Hazards
#'
#' @param t  (`numeric`)\cr  study time-points.
#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition.
#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition.
#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition.
#' @param p01 (positive `number`)\cr rate parameter of Weibull distribution for `h01`.
#' @param p02 (positive `number`)\cr rate parameter of Weibull distribution for `h02`.
#' @param p12 (positive `number`)\cr rate parameter of Weibull distribution for `h12`.
#'
#' @return This returns the value of OS survival function at time t.
#' @export
#'
#' @examples WeibSurvOS(c(1:5), 0.2, 0.5, 2.1, 1.2, 0.9, 1)
WeibSurvOS <- function(t, h01, h02, h12, p01, p02, p12) {
  assert_numeric(t, lower = 0, any.missing = FALSE)
  assert_positive_number(h01, zero_ok = TRUE)
  assert_positive_number(p01)

  WeibSurvPFS(t, h01, h02, p01, p02) +
    h01 * p01 *
      sapply(t, function(t) {
        integrateVector(WeibOSInteg,
          upper = t,
          t = t,
          h01 = h01,
          h02 = h02,
          h12 = h12,
          p01 = p01,
          p02 = p02,
          p12 = p12
        )
      })
}

#' Cumulative Hazard for Piecewise Constant Hazards
#'
#' @param t (`numeric`)\cr  study time-points.
#' @param haz (`numeric vector`)\cr constant transition hazards.
#' @param pw (`numeric vector`)\cr time intervals for the piecewise constant hazards.
#'
#' @return This returns the value of cumulative hazard at time t.
#' @export
#'
#' @examples
#' pwA(1:5, c(0.5, 0.9), c(0, 4))
pwA <- function(t, haz, pw) {
  assert_numeric(t, lower = 0, any.missing = FALSE)
  assert_numeric(haz, lower = 0, any.missing = FALSE, all.missing = FALSE)
  assert_intervals(pw, length(haz))

  # Time intervals: time-points + cut points.
  pw_new <- c(pw[pw < max(t)])
  pw_time <- sort(unique(c(pw_new, t)))
  haz_all <- getPWCHazard(haz, pw, pw_time)

  # Cumulative hazard.
  i <- 1:(length(pw_time) - 1)
  dA_pw <- haz_all[i] * (pw_time[i + 1] - pw_time[i])
  A_pw <- c(0, cumsum(dA_pw))

  # Match result with input.
  A_pw[match(t, pw_time)]
}

#' PFS Survival Function from Piecewise Constant Hazards
#'
#' @inheritParams PWCsurvOS
#'
#' @return This returns the value of PFS survival function at time t.
#' @export
#'
#' @examples
#' PWCsurvPFS(1:5, c(0.3, 0.5), c(0.5, 0.8), c(0, 4), c(0, 8))
PWCsurvPFS <- function(t, h01, h02, pw01, pw02) {
  assert_numeric(t, lower = 0, any.missing = FALSE)
  assert_numeric(h01, lower = 0, any.missing = FALSE, all.missing = FALSE)
  assert_numeric(h02, lower = 0, any.missing = FALSE, all.missing = FALSE)

  assert_intervals(pw01, length(h01))
  assert_intervals(pw02, length(h02))

  A01 <- pwA(t, h01, pw01)
  A02 <- pwA(t, h02, pw02)

  exp(-A01 - A02)
}

#' Helper Function of `PWCsurvOS()`
#'
#' @param x (`numeric`)\cr  variable of integration.
#' @inheritParams PWCsurvOS
#'
#' @return Numeric results of the integrand used to calculate
#' the OS survival function for piecewise constant transition hazards, see  `PWCsurvOS`.
#'
#' @keywords internal
PwcOSInt <- function(x, t, h01, h02, h12, pw01, pw02, pw12) {
  PWCsurvPFS(x, h01, h02, pw01, pw02) *
    getPWCHazard(h01, pw01, x) *
    exp(pwA(x, h12, pw12) - pwA(t, h12, pw12))
}

#' OS Survival Function from Piecewise Constant Hazards
#'
#' @param t (`numeric`)\cr study time-points.
#' @param h01 (`numeric vector`)\cr constant transition hazards for 0 to 1 transition.
#' @param h02 (`numeric vector`)\cr constant transition hazards for 0 to 2 transition.
#' @param h12 (`numeric vector`)\cr constant transition hazards for 1 to 2 transition.
#' @param pw01 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h01`.
#' @param pw02 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h02`.
#' @param pw12 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h12`.
#'
#' @return This returns the value of OS survival function at time t.
#' @export
#'
#' @examples
#' PWCsurvOS(1:5, c(0.3, 0.5), c(0.5, 0.8), c(0.7, 1), c(0, 4), c(0, 8), c(0, 3))
PWCsurvOS <- function(t, h01, h02, h12, pw01, pw02, pw12) {
  assert_numeric(t, lower = 0, any.missing = FALSE)
  assert_numeric(h01, lower = 0, any.missing = FALSE, all.missing = FALSE)
  assert_numeric(h02, lower = 0, any.missing = FALSE, all.missing = FALSE)
  assert_numeric(h12, lower = 0, any.missing = FALSE, all.missing = FALSE)
  assert_intervals(pw01, length(h01))
  assert_intervals(pw02, length(h02))
  assert_intervals(pw12, length(h12))

  # Assemble unique time points and corresponding hazard values once.
  unique_pw_times <- sort(unique(c(pw01, pw02, pw12)))
  h01_at_times <- getPWCHazard(h01, pw01, unique_pw_times)
  h02_at_times <- getPWCHazard(h02, pw02, unique_pw_times)
  h12_at_times <- getPWCHazard(h12, pw12, unique_pw_times)

  # We work for the integral parts with sorted and unique time points.
  t_sorted <- sort(unique(t))
  t_sorted_zero <- c(0, t_sorted)
  n_t <- length(t_sorted)
  int_sums <- numeric(n_t)
  cum_haz_12 <- pwA(t_sorted, h12, pw12)
  for (i in seq_len(n_t)) {
    t_start <- t_sorted_zero[i]
    t_end <- t_sorted_zero[i + 1]
    # Determine the indices of the time intervals we are working with here.
    start_index <- findInterval(t_start, unique_pw_times)
    # We use here left.open = TRUE, so that when t_end is identical to a change point,
    # then we don't need an additional part.
    end_index <- max(findInterval(t_end, unique_pw_times, left.open = TRUE), 1L)
    index_bounds <- start_index:end_index
    # Determine corresponding time intervals.
    time_bounds <- if (length(index_bounds) == 1L) {
      rbind(c(t_start, t_end))
    } else if (length(index_bounds) == 2L) {
      rbind(
        c(t_start, unique_pw_times[end_index]),
        c(unique_pw_times[end_index], t_end)
      )
    } else {
      n_ind_bounds <- length(index_bounds)
      rbind(
        c(t_start, unique_pw_times[start_index + 1L]),
        cbind(
          unique_pw_times[index_bounds[-c(1, n_ind_bounds)]],
          unique_pw_times[index_bounds[-c(1, 2)]]
        ),
        c(unique_pw_times[end_index], t_end)
      )
    }
    for (j in seq_along(index_bounds)) {
      time_j <- time_bounds[j, 1]
      time_jp1 <- time_bounds[j, 2]
      intercept_a <- pwA(time_j, h12, pw12) - pwA(time_j, h01, pw01) - pwA(time_j, h02, pw02)
      slope_b <- h12_at_times[index_bounds[j]] - h01_at_times[index_bounds[j]] - h02_at_times[index_bounds[j]]
      h01_j <- h01_at_times[index_bounds[j]]
      int_js <- if (slope_b != 0) {
        h01_j / slope_b * exp(intercept_a) *
          (exp(slope_b * (time_jp1 - time_j) - cum_haz_12[i:n_t]) - exp(-cum_haz_12[i:n_t]))
      } else {
        h01_j * exp(intercept_a - cum_haz_12[i:n_t]) * (time_jp1 - time_j)
      }
      int_sums[i:n_t] <- int_sums[i:n_t] + int_js
    }
  }

  # Match back to all time points,
  # which might not be unique or ordered.
  int_sums <- int_sums[match(t, t_sorted)]
  result <- PWCsurvPFS(t, h01, h02, pw01, pw02) + int_sums

  # Cap at 1 in a safe way - i.e. first check.
  assert_numeric(result, finite = TRUE, any.missing = FALSE)
  above_one <- result > 1
  if (any(above_one)) {
    assert_true(all((result[above_one] - 1) < sqrt(.Machine$double.eps)))
    result[above_one] <- 1
  }
  result
}

#' Helper Function for Single Quantile for OS Survival Function
#'
#' @param q (`number`)\cr single quantile at which to compute event time.
#' @inheritParams ExpQuantOS
#'
#' @return Single time t such that the OS survival function at t equals q.
#' @keywords internal
singleExpQuantOS <- function(q, h01, h02, h12) {
  assert_number(q, finite = TRUE, lower = 0, upper = 1)
  toroot <- function(x) {
    return(q - ExpSurvOS(t = x, h01, h02, h12))
  }
  stats::uniroot(
    toroot,
    interval = c(0, 10^3),
    extendInt = "yes"
  )$root
}

#' Quantile function for OS survival function induced by an illness-death model
#'
#' @param q (`numeric`)\cr quantile(s) at which to compute event time (q = 1 / 2 corresponds to median).
#' @param h01 (`numeric vector`)\cr constant transition hazards for 0 to 1 transition.
#' @param h02 (`numeric vector`)\cr constant transition hazards for 0 to 2 transition.
#' @param h12 (`numeric vector`)\cr constant transition hazards for 1 to 2 transition.
#'
#' @return This returns the time(s) t such that the OS survival function at t equals q.
#' @export
#'
#' @examples ExpQuantOS(1 / 2, 0.2, 0.5, 2.1)
ExpQuantOS <- function(q = 1 / 2, h01, h02, h12) {
  assert_numeric(q, lower = 0, upper = 1, any.missing = FALSE)
  assert_numeric(h01, lower = 0, any.missing = FALSE)
  assert_numeric(h02, lower = 0, any.missing = FALSE)
  assert_numeric(h12, lower = 0, any.missing = FALSE)

  vapply(
    q,
    FUN = singleExpQuantOS,
    FUN.VALUE = numeric(1),
    h01 = h01,
    h02 = h02,
    h12 = h12
  )
}

#' Transition Hazards for Exponential Event Times
#'
#' This creates a list with class `TransitionParameters` containing
#' hazards, time intervals and Weibull rates for exponential event times
#' in an illness-death model.
#'
#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition.
#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition.
#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition.
#'
#' @return List with elements `hazards`, `intervals`, `weibull_rates` and `family`
#'   (exponential).
#' @export
#'
#' @examples
#' exponential_transition(1, 1.6, 0.3)
exponential_transition <- function(h01 = 1, h02 = 1, h12 = 1) {
  assert_positive_number(h01, zero_ok = TRUE)
  assert_positive_number(h02, zero_ok = TRUE)
  assert_positive_number(h12, zero_ok = TRUE)
  structure(
    list(
      hazards = list(h01 = h01, h02 = h02, h12 = h12),
      intervals = list(pw01 = 0, pw02 = 0, pw12 = 0),
      weibull_rates = list(p01 = 1, p02 = 1, p12 = 1),
      family = "exponential"
    ),
    class = c("ExponentialTransition", "TransitionParameters")
  )
}

#' Transition Hazards for Piecewise Exponential Event Times
#'
#' This creates a list with class `TransitionParameters` containing
#' hazards, time intervals and Weibull rates for piecewise exponential event times
#' in an illness-death model.
#'
#' @param h01 (`numeric vector`)\cr constant transition hazards for 0 to 1 transition
#' @param h02 (`numeric vector`)\cr  constant transition hazards for 0 to 2 transition
#' @param h12 (`numeric vector`)\cr  constant transition hazards for 1 to 2 transition
#' @param pw01 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h01`
#' @param pw02 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h02`
#' @param pw12 (`numeric vector`)\cr  time intervals for the piecewise constant hazards `h12`
#'
#' @return List with elements `hazards`, `intervals`, `weibull_rates` and `family`
#'   (piecewise exponential).
#' @export
#'
#' @examples
#' piecewise_exponential(
#'   h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1),
#'   pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9)
#' )
piecewise_exponential <- function(h01, h02, h12, pw01, pw02, pw12) {
  assert_numeric(h01, lower = 0, any.missing = FALSE, all.missing = FALSE)
  assert_numeric(h02, lower = 0, any.missing = FALSE, all.missing = FALSE)
  assert_numeric(h12, lower = 0, any.missing = FALSE, all.missing = FALSE)

  assert_intervals(pw01, length(h01))
  assert_intervals(pw02, length(h02))
  assert_intervals(pw12, length(h12))

  structure(
    list(
      hazards = list(h01 = h01, h02 = h02, h12 = h12),
      intervals = list(pw01 = pw01, pw02 = pw02, pw12 = pw12),
      weibull_rates = list(p01 = 1, p02 = 1, p12 = 1),
      family = "piecewise exponential"
    ),
    class = c("PWCTransition", "TransitionParameters")
  )
}

#' Transition Hazards for Weibull Distributed Event Times
#'
#' This creates a list with class `TransitionParameters` containing
#' hazards, time intervals and Weibull rates for Weibull distributed event times
#' in an illness-death model.
#'
#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition
#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition
#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition
#' @param p01 (positive `number`)\cr rate parameter of Weibull distribution for `h01`
#' @param p02 (positive `number`)\cr rate parameter of Weibull distribution for `h02`
#' @param p12 (positive `number`)\cr rate parameter of Weibull distribution for `h12`
#'
#' @return List with elements `hazards`, `intervals`, `weibull_rates` and `family`
#'   (Weibull).
#' @export
#'
#' @examples
#' weibull_transition(h01 = 1, h02 = 1.3, h12 = 0.5, p01 = 1.2, p02 = 1.3, p12 = 0.5)
weibull_transition <- function(h01 = 1, h02 = 1, h12 = 1, p01 = 1, p02 = 1, p12 = 1) {
  assert_positive_number(h01, zero_ok = TRUE)
  assert_positive_number(h02, zero_ok = TRUE)
  assert_positive_number(h12, zero_ok = TRUE)
  assert_positive_number(p01)
  assert_positive_number(p02)
  assert_positive_number(p12)
  structure(
    list(
      hazards = list(h01 = h01, h02 = h02, h12 = h12),
      weibull_rates = list(p01 = p01, p02 = p02, p12 = p12),
      intervals = list(pw01 = 0, pw02 = 0, pw12 = 0),
      family = "Weibull"
    ),
    class = c("WeibullTransition", "TransitionParameters")
  )
}

#' OS Hazard Function from Constant Transition Hazards
#'
#' @param t (`numeric`)\cr study time-points.
#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition.
#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition.
#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition.
#'
#' @return This returns the value of the OS hazard function at time t.
#' @export
#'
#' @examples
#' ExpHazOS(c(1:5), 0.2, 1.1, 0.8)
ExpHazOS <- function(t, h01, h02, h12) {
  assert_numeric(t, lower = 0, any.missing = FALSE)
  assert_positive_number(h01)
  assert_positive_number(h02)
  assert_positive_number(h12)

  h012 <- h12 - h01 - h02
  ((h12 - h02) * (h01 + h02) - h01 * h12 * exp(-h012 * t)) / ((h12 - h02) - h01 * exp(-h012 * t))
}

#' Helper Function for `avgHRExpOS()`
#'
#' It is an integrand of the form OS hazard function with intensities h01, h02, h12
#' at time point t multiplied with a weighted product of the two OS Survival functions
#' at t (one for intensities h0 and one for h1).
#'
#' @param x (`numeric`)\cr variable of integration.
#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition.
#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition.
#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition.
#' @param h0 (`list`)\cr transition parameters for the first treatment group.
#' @param h1 (`list`)\cr transition parameters for the second treatment group.
#' @param alpha (`number`)\cr weight parameter, see [avgHRExpOS()].
#' @return This returns the value of the integrand used to calculate
#' the average hazard ratio for constant transition hazards, see [avgHRExpOS()].
#' @export
#'
#' @examples
#' h0 <- list(h01 = 0.18, h02 = 0.06, h12 = 0.17)
#' h1 <- list(h01 = 0.23, h02 = 0.07, h12 = 0.19)
#' avgHRIntegExpOS(x = 5, h01 = 0.2, h02 = 0.5, h12 = 0.7, h0 = h0, h1 = h1, alpha = 0.5)
avgHRIntegExpOS <- function(x, h01, h02, h12, h0, h1, alpha) {
  assert_positive_number(alpha)

  weightedSurv <- (ExpSurvOS(t = x, h01 = h0$h01, h02 = h0$h02, h12 = h0$h12) *
    ExpSurvOS(t = x, h01 = h1$h01, h02 = h1$h02, h12 = h1$h12))^alpha
  ExpHazOS(x, h01, h02, h12) * weightedSurv
}

#' Average OS Hazard Ratio from Constant Transition Hazards
#'
#' @param transitionByArm (`list`)\cr transition parameters for each treatment group.
#'   See [exponential_transition()], [piecewise_exponential()] and [weibull_transition()] for details.
#' @param alpha (`number`)\cr assigns weights to time points, where values higher than 0.5 assign greater weight
#' to earlier times and values lower than 0.5 assign greater weight to later times.
#' @param upper (`number`)\cr upper (time) limit of integration.
#'
#' @return This returns the value of the average hazard ratio.
#' @export
#'
#' @examples
#' transitionTrt <- exponential_transition(h01 = 0.18, h02 = 0.06, h12 = 0.17)
#' transitionCtl <- exponential_transition(h01 = 0.23, h02 = 0.07, h12 = 0.19)
#' transitionList <- list(transitionCtl, transitionTrt)
#' avgHRExpOS(transitionByArm = transitionList, alpha = 0.5, upper = 100)
avgHRExpOS <- function(transitionByArm, alpha = 0.5, upper = Inf) {
  assert_list(transitionByArm, len = 2)
  assert_class(transitionByArm[[1]], "TransitionParameters")
  assert_class(transitionByArm[[2]], "TransitionParameters")
  assert_positive_number(alpha)
  assert_positive_number(upper)

  h0 <- transitionByArm[[1]]$hazard
  h1 <- transitionByArm[[2]]$hazard

  num <- stats::integrate(avgHRIntegExpOS,
    lower = 0, upper = upper, h01 = h1$h01, h02 = h1$h02, h12 = h1$h12,
    h0 = h0, h1 = h1, alpha = alpha
  )$value
  denom <- stats::integrate(avgHRIntegExpOS,
    lower = 0, upper = upper, h01 = h0$h01, h02 = h0$h02, h12 = h0$h12,
    h0 = h0, h1 = h1, alpha = alpha
  )$value
  num / denom
}

#' Piecewise Exponentially Distributed Event Times
#'
#' This returns event times with a distribution resulting from piece-wise constant hazards
#' using the inversion method.
#'
#' @param U (`numeric`)\cr uniformly distributed random variables.
#' @param haz (`numeric`)\cr piecewise constant hazard.
#' @param pw (`numeric`)\cr time intervals for the piecewise constant hazard.
#' @param t_0 (`numeric`)\cr the starting times.
#'
#' @return This returns a vector with event times.
#' @export
#'
#' @examples
#' getPCWDistr(U = runif(3), haz = c(1.1, 0.5, 0.4), pw = c(0, 7, 10), t_0 = c(0, 1, 4.2))
getPCWDistr <- function(U, haz, pw, t_0) {
  assert_numeric(U, lower = 0, upper = 1, any.missing = FALSE, all.missing = FALSE)
  assert_numeric(haz, lower = 0, any.missing = FALSE, all.missing = FALSE)
  assert_intervals(pw, length(haz))
  assert_numeric(t_0, lower = 0, len = length(U), any.missing = FALSE, all.missing = FALSE)

  N <- length(U)
  t1 <- rep(NA, N) # Initialize event times.
  n2 <- length(haz)
  for (kk in seq_len(N)) {
    # Shift if t_0 !=0.
    cuts_temp <- c(t_0[kk], pw[pw > t_0[kk]])
    cuts_temp <- cuts_temp - t_0[kk]
    # Number of cut-points after shift.
    n <- length(cuts_temp)
    # Number of hazards to remove for shift (t_0 > times).
    remov <- n2 - n
    haz_temp <- haz[(remov + 1):n2]
    LogU <- log(1 - U[kk])

    t1[kk] <- if (n != 1) {
      PCWInversionMethod(haz = haz_temp, pw = cuts_temp, LogU = LogU)
    } else if (n == 1) { # I.e. exponential distribution.
      -LogU / haz_temp
    }
  }
  t1
}

#' Single Piecewise Exponentially Distributed Event Time
#'
#' This returns an event time with a distribution resulting from piece-wise constant hazards
#' using the inversion method.
#'
#' @param pw (`numeric`)\cr time intervals for the piecewise constant hazard.
#' @param haz (`numeric`)\cr piecewise constant hazard.
#' @param LogU (`numeric`)\cr transformed uniformly distributed random variables (log(1-U)).
#'
#' @return This returns one single event time.
#' @export
#'
#' @examples
#' PCWInversionMethod(haz = c(1.1, 0.5, 0.4), pw = c(0, 7, 10), LogU = log(1 - runif(1)))
PCWInversionMethod <- function(haz, pw, LogU) {
  assert_numeric(haz, lower = 0, any.missing = FALSE, all.missing = FALSE)
  assert_intervals(pw, length(haz))
  assert_number(LogU)

  n <- length(pw)
  t1 <- 0
  # Determine sum of alpha*time-interval for all i.
  dt <- pw[2:n] - pw[1:(n - 1)]

  # Helping matrix.
  tempMatrix <- matrix(0, nrow = n, ncol = n - 1)
  tempMatrix[lower.tri(tempMatrix)] <- 1

  sumA <- -as.vector(tempMatrix %*% (haz[1:(n - 1)] * dt))
  # Find the appropriate time interval.
  rev_index <- findInterval(LogU, rev(sumA), rightmost.closed = FALSE, left.open = TRUE, all.inside = FALSE)
  # Use reverse index.
  index <- n - rev_index
  # Calculate event time.
  pw[index] + (sumA[index] - LogU) / haz[index]
}

#' Piecewise Constant Hazard Values
#'
#' This returns piecewise constant hazard values at specified time points.
#'
#' @param haz (`numeric`)\cr piecewise constant input hazard.
#' @param pw (`numeric`)\cr time intervals for the piecewise constant hazard.
#' @param x (`numeric`)\cr time-points.
#'
#' @return Hazard values at input time-points.
#' @export
#'
#' @examples
#' getPWCHazard(c(1, 1.2, 1.4), c(0, 2, 3), c(1, 4, 6))
getPWCHazard <- function(haz, pw, x) { # nolint
  sapply(x, function(jj) {
    y <- NULL
    # Find interval and corresponding hazard value for time x[jj].
    for (ii in seq_along(haz)) {
      if (jj >= pw[ii]) {
        y <- haz[ii]
      }
    }
    y
  })
}

#' Sum of Two Piecewise Constant Hazards
#'
#' This returns the sum of two piecewise constant hazards per interval.
#'
#' @param haz1 (`numeric`)\cr first summand (piecewise constant hazard).
#' @param haz2 (`numeric`)\cr second summand (piecewise constant hazard).
#' @param pw1 (`numeric`)\cr time intervals of first summand.
#' @param pw2  (`numeric`)\cr time intervals of second summand.
#'
#' @return List with elements `hazards` and `intervals` for the sum of two piecewise constant hazards.
#' @export
#'
#' @examples
#' getSumPCW(c(1.2, 0.3, 0.6), c(1.2, 0.7, 1), c(0, 8, 9), c(0, 1, 4))
getSumPCW <- function(haz1, haz2, pw1, pw2) {
  # Get all cutpoints for the intervals.
  cuts_sum <- unique(sort(c(pw1, pw2)))
  haz_sum <- NULL
  # Get sum of hazards for all intervals.
  for (i in seq_along(cuts_sum)) {
    haz_sum[i] <- getPWCHazard(haz1, pw1, cuts_sum[i]) +
      getPWCHazard(haz2, pw2, cuts_sum[i])
  }
  list(
    hazards = haz_sum,
    intervals = cuts_sum
  )
}

#' Assertion for Positive Number
#'
#' @param x what to check.
#' @param zero_ok (`flag`)\cr whether `x` can be zero or not.
#'
#' @return Raises an error if `x` is not a single positive (or non-negative) number.
#' @export
#'
#' @examples
#' assert_positive_number(3.2)
#' assert_positive_number(0, zero_ok = TRUE)
assert_positive_number <- function(x,
                                   zero_ok = FALSE) {
  assert_number(x)
  assert_flag(zero_ok)
  if (zero_ok) {
    assert_true(x >= 0)
  } else {
    assert_true(x > 0)
  }
}

#' Assertion for vector describing intervals
#'
#' We define an intervals vector to always start with 0, and contain
#' unique ordered time points.
#'
#' @param x what to check.
#' @param y (`count`)\cr required length of `y`.
#'
#' @return Raises an error if `x` is not an intervals vector starting with 0.
#' @export
#'
#' @examples
#' assert_intervals(c(0, 5, 7), 3)
assert_intervals <- function(x, y) {
  assert_numeric(
    x,
    lower = 0,
    any.missing = FALSE,
    all.missing = FALSE,
    len = y,
    unique = TRUE,
    sorted = TRUE
  )
  assert_true(x[1] == 0)
}

1		#' Log-Rank Test for a Single Trial
2		#'
3		#' This function evaluates the significance of either PFS or OS endpoints in a trial,
4		#' based on a pre-specified critical value.
5		#'
6		#' @param data (`data.frame`)\cr data frame containing entry and exit times of an
7		#' illness-death model. See [getSimulatedData()] for details.
8		#' @param typeEvent (`string`)\cr endpoint to be evaluated, possible values are `PFS` and `OS`.
9		#' @param critical (positive `number`)\cr critical value of the log-rank test.
10		#'
11		#' @return Logical value indicating log-rank test significance.
12		#' @export
13		#'
14		#' @examples
15		#' transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
16		#' transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
17		#' simTrial <- getClinicalTrials(
18		#' nRep = 1, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
19		#' transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
20		#' accrual = list(param = "intensity", value = 7)
21		#' )[[1]]
22		#' logRankTest(data = simTrial, typeEvent = "OS", critical = 3.4)
23		logRankTest <- function(data, typeEvent = c("PFS", "OS"), critical) {
24		# Must be have the format of one row per patient (`datType` must be 1rowPatient`
25		# in getClinicalTrials()), i.e. have 10 (if censored) or 11 columns.
26	108x	assert_data_frame(data, min.cols = 10, max.cols = 11)
27	108x	typeEvent <- match.arg(typeEvent)
28	108x	assert_positive_number(critical)
29
30	108x	if (typeEvent == "OS") {
31	54x	time <- data$OStime
32	54x	event <- data$OSevent
33	54x	} else if (typeEvent == "PFS") {
34	54x	time <- data$PFStime
35	54x	event <- data$PFSevent
36		}
37
38	108x	logRank <- survival::survdiff(survival::Surv(time, event) ~ trt, data)
39	108x	sqrt(logRank$chisq) > critical
40		}
41
42
43		#' Helper function to conduct log-rank tests for either PFS or OS
44		#'
45		#' This function evaluates the significance of either PFS or OS endpoints for each trial
46		#' in a list of trials, based on a pre-specified critical value.
47		#'
48		#' @param simTrials (`list`)\cr simulated trial data sets, see [getClinicalTrials()].
49		#' @param typeEvent (`string`)\cr endpoint to be evaluated, possible values are `PFS` and `OS`.
50		#' @param eventNum (`integer`)\cr number of events required to trigger analysis.
51		#' @param critical (positive `number`)\cr critical value of the log-rank test.
52		#'
53		#' @return Logical vector indicating log-rank test significance for each trial.
54		#' @export
55		#'
56		#' @examples
57		#' transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
58		#' transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
59		#' simTrials <- getClinicalTrials(
60		#' nRep = 50, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
61		#' transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
62		#' accrual = list(param = "intensity", value = 7)
63		#' )
64		#' passedLogRank(simTrials = simTrials, typeEvent = "PFS", eventNum = 300, critical = 2.4)
65		#' @keywords internal
66		passedLogRank <- function(simTrials, typeEvent, eventNum, critical) {
67	4x	assert_list(simTrials, null.ok = FALSE)
68	4x	assert_positive_number(critical)
69	4x	assert_count(eventNum, positive = TRUE)
70
71		# Censor simulated trials at time-point of OS/PFS analysis.
72	4x	trialsAna <- lapply(
73	4x	X = simTrials,
74	4x	FUN = censoringByNumberEvents,
75	4x	eventNum = eventNum,
76	4x	typeEvent = typeEvent
77		)
78		# Compute log-rank test for all trials for OS/PFS.
79	4x	unlist(lapply(
80	4x	X = trialsAna,
81	4x	FUN = logRankTest,
82	4x	typeEvent = typeEvent,
83	4x	critical = critical
84		))
85		}
86
87		#' Empirical Significance for a List of Simulated Trials
88		#'
89		#' This function computes four types of empirical significance — PFS, OS, at-least (significant
90		#' in at least one of PFS/OS), and joint (significant in both PFS and OS) — using
91		#' the log-rank test. Empirical significance is calculated as the proportion of significant
92		#' results in simulated trials, each ending when a set number of PFS/OS events occur.
93		#' Critical values for PFS and OS test significance must be specified. If trials simulate equal
94		#' transition hazards across groups (H0), empirical significance estimates type I error;
95		#' if they simulate differing transition hazards (H1), it estimates power.
96		#'
97		#' @param simTrials (`list`)\cr simulated trial data sets, see [getClinicalTrials()].
98		#' @param criticalPFS (positive `number`)\cr critical value of the log-rank test for PFS.
99		#' @param criticalOS (positive `number`)\cr critical value of the log-rank test for OS.
100		#' @param eventNumPFS (`integer`)\cr number of PFS events required to trigger PFS analysis.
101		#' @param eventNumOS (`integer`)\cr number of OS events required to trigger OS analysis.
102		#'
103		#' @return This returns values of four measures of empirical significance.
104		#' @export
105		#'
106		#' @examples
107		#' transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
108		#' transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
109		#' simTrials <- getClinicalTrials(
110		#' nRep = 50, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
111		#' transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
112		#' accrual = list(param = "intensity", value = 7)
113		#' )
114		#' empSignificant(
115		#' simTrials = simTrials, criticalPFS = 2.4, criticalOS = 2.2,
116		#' eventNumPFS = 300, eventNumOS = 500
117		#' )
118		empSignificant <- function(simTrials, criticalPFS, criticalOS, eventNumPFS, eventNumOS) {
119	1x	assert_list(simTrials, null.ok = FALSE)
120	1x	assert_positive_number(criticalPFS)
121	1x	assert_positive_number(criticalOS)
122	1x	assert_count(eventNumPFS, positive = TRUE)
123	1x	assert_count(eventNumOS, positive = TRUE)
124
125	1x	nRep <- length(simTrials)
126	1x	simTrials <- lapply(
127	1x	X = simTrials,
128	1x	FUN = function(x) if (ncol(x) == 9) getDatasetWideFormat(x) else x
129		)
130
131		# Which trials passed the log-rank test for PFS/OS?
132	1x	passedLogRankPFS <- passedLogRank(
133	1x	simTrials = simTrials,
134	1x	typeEvent = "PFS",
135	1x	eventNum = eventNumPFS,
136	1x	critical = criticalPFS
137		)
138	1x	passedLogRankOS <- passedLogRank(
139	1x	simTrials = simTrials,
140	1x	typeEvent = "OS",
141	1x	eventNum = eventNumOS,
142	1x	critical = criticalOS
143		)
144
145		# Empirical significance is the fraction of trials with significant log-rank test:
146	1x	significantPFS <- sum(passedLogRankPFS) / nRep
147	1x	significantOS <- sum(passedLogRankOS) / nRep
148
149		# Derived measures of significance:
150	1x	sumPassed <- passedLogRankPFS + passedLogRankOS
151		# At least one of PFS/OS has significant log-rank test.
152	1x	significantAtLeastOne <- sum(sumPassed >= 1) / nRep
153		# Both PFS/OS have significant log-rank tests.
154	1x	significantBoth <- sum(sumPassed == 2) / nRep
155
156	1x	list(
157	1x	"significantPFS" = significantPFS,
158	1x	"significantOS" = significantOS,
159	1x	"significantAtLeastOne" = significantAtLeastOne,
160	1x	"significantBoth" = significantBoth
161		)
162		}

1		#' Transition Hazards for Exponential Event Times
2		#'
3		#' This creates a list with class `TransitionParameters` containing
4		#' hazards, time intervals and Weibull rates for exponential event times
5		#' in an illness-death model.
6		#'
7		#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition.
8		#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition.
9		#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition.
10		#'
11		#' @return List with elements `hazards`, `intervals`, `weibull_rates` and `family`
12		#' (exponential).
13		#' @export
14		#'
15		#' @examples
16		#' exponential_transition(1, 1.6, 0.3)
17		exponential_transition <- function(h01 = 1, h02 = 1, h12 = 1) {
18	137x	assert_positive_number(h01, zero_ok = TRUE)
19	136x	assert_positive_number(h02, zero_ok = TRUE)
20	136x	assert_positive_number(h12, zero_ok = TRUE)
21	136x	structure(
22	136x	list(
23	136x	hazards = list(h01 = h01, h02 = h02, h12 = h12),
24	136x	intervals = list(pw01 = 0, pw02 = 0, pw12 = 0),
25	136x	weibull_rates = list(p01 = 1, p02 = 1, p12 = 1),
26	136x	family = "exponential"
27		),
28	136x	class = c("ExponentialTransition", "TransitionParameters")
29		)
30		}
31
32		#' Transition Hazards for Piecewise Exponential Event Times
33		#'
34		#' This creates a list with class `TransitionParameters` containing
35		#' hazards, time intervals and Weibull rates for piecewise exponential event times
36		#' in an illness-death model.
37		#'
38		#' @param h01 (`numeric vector`)\cr constant transition hazards for 0 to 1 transition
39		#' @param h02 (`numeric vector`)\cr constant transition hazards for 0 to 2 transition
40		#' @param h12 (`numeric vector`)\cr constant transition hazards for 1 to 2 transition
41		#' @param pw01 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h01`
42		#' @param pw02 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h02`
43		#' @param pw12 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h12`
44		#'
45		#' @return List with elements `hazards`, `intervals`, `weibull_rates` and `family`
46		#' (piecewise exponential).
47		#' @export
48		#'
49		#' @examples
50		#' piecewise_exponential(
51		#' h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1),
52		#' pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9)
53		#' )
54		piecewise_exponential <- function(h01, h02, h12, pw01, pw02, pw12) {
55	15x	assert_numeric(h01, lower = 0, any.missing = FALSE, all.missing = FALSE)
56	14x	assert_numeric(h02, lower = 0, any.missing = FALSE, all.missing = FALSE)
57	14x	assert_numeric(h12, lower = 0, any.missing = FALSE, all.missing = FALSE)
58
59	14x	assert_intervals(pw01, length(h01))
60	13x	assert_intervals(pw02, length(h02))
61	12x	assert_intervals(pw12, length(h12))
62
63	12x	structure(
64	12x	list(
65	12x	hazards = list(h01 = h01, h02 = h02, h12 = h12),
66	12x	intervals = list(pw01 = pw01, pw02 = pw02, pw12 = pw12),
67	12x	weibull_rates = list(p01 = 1, p02 = 1, p12 = 1),
68	12x	family = "piecewise exponential"
69		),
70	12x	class = c("PWCTransition", "TransitionParameters")
71		)
72		}
73
74		#' Transition Hazards for Weibull Distributed Event Times
75		#'
76		#' This creates a list with class `TransitionParameters` containing
77		#' hazards, time intervals and Weibull rates for Weibull distributed event times
78		#' in an illness-death model.
79		#'
80		#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition
81		#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition
82		#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition
83		#' @param p01 (positive `number`)\cr rate parameter of Weibull distribution for `h01`
84		#' @param p02 (positive `number`)\cr rate parameter of Weibull distribution for `h02`
85		#' @param p12 (positive `number`)\cr rate parameter of Weibull distribution for `h12`
86		#'
87		#' @return List with elements `hazards`, `intervals`, `weibull_rates` and `family`
88		#' (Weibull).
89		#' @export
90		#'
91		#' @examples
92		#' weibull_transition(h01 = 1, h02 = 1.3, h12 = 0.5, p01 = 1.2, p02 = 1.3, p12 = 0.5)
93		weibull_transition <- function(h01 = 1, h02 = 1, h12 = 1, p01 = 1, p02 = 1, p12 = 1) {
94	236x	assert_positive_number(h01, zero_ok = TRUE)
95	236x	assert_positive_number(h02, zero_ok = TRUE)
96	235x	assert_positive_number(h12, zero_ok = TRUE)
97	235x	assert_positive_number(p01)
98	235x	assert_positive_number(p02)
99	235x	assert_positive_number(p12)
100	234x	structure(
101	234x	list(
102	234x	hazards = list(h01 = h01, h02 = h02, h12 = h12),
103	234x	weibull_rates = list(p01 = p01, p02 = p02, p12 = p12),
104	234x	intervals = list(pw01 = 0, pw02 = 0, pw12 = 0),
105	234x	family = "Weibull"
106		),
107	234x	class = c("WeibullTransition", "TransitionParameters")
108		)
109		}

1		#' Piecewise Exponentially Distributed Event Times
2		#'
3		#' This returns event times with a distribution resulting from piece-wise constant hazards
4		#' using the inversion method.
5		#'
6		#' @param U (`numeric`)\cr uniformly distributed random variables.
7		#' @param haz (`numeric`)\cr piecewise constant hazard.
8		#' @param pw (`numeric`)\cr time intervals for the piecewise constant hazard.
9		#' @param t_0 (`numeric`)\cr the starting times.
10		#'
11		#' @return This returns a vector with event times.
12		#' @export
13		#'
14		#' @examples
15		#' getPCWDistr(U = runif(3), haz = c(1.1, 0.5, 0.4), pw = c(0, 7, 10), t_0 = c(0, 1, 4.2))
16		getPCWDistr <- function(U, haz, pw, t_0) {
17	4013x	assert_numeric(U, lower = 0, upper = 1, any.missing = FALSE, all.missing = FALSE)
18	4013x	assert_numeric(haz, lower = 0, any.missing = FALSE, all.missing = FALSE)
19	4013x	assert_intervals(pw, length(haz))
20	4013x	assert_numeric(t_0, lower = 0, len = length(U), any.missing = FALSE, all.missing = FALSE)
21
22	4013x	N <- length(U)
23	4013x	t1 <- rep(NA, N) # Initialize event times.
24	4013x	n2 <- length(haz)
25	4013x	for (kk in seq_len(N)) {
26		# Shift if t_0 !=0.
27	556838x	cuts_temp <- c(t_0[kk], pw[pw > t_0[kk]])
28	556838x	cuts_temp <- cuts_temp - t_0[kk]
29		# Number of cut-points after shift.
30	556838x	n <- length(cuts_temp)
31		# Number of hazards to remove for shift (t_0 > times).
32	556838x	remov <- n2 - n
33	556838x	haz_temp <- haz[(remov + 1):n2]
34	556838x	LogU <- log(1 - U[kk])
35
36	556838x	t1[kk] <- if (n != 1) {
37	556833x	PCWInversionMethod(haz = haz_temp, pw = cuts_temp, LogU = LogU)
38	556838x	} else if (n == 1) { # I.e. exponential distribution.
39	5x	-LogU / haz_temp
40		}
41		}
42	4013x	t1
43		}
44
45		#' Single Piecewise Exponentially Distributed Event Time
46		#'
47		#' This returns an event time with a distribution resulting from piece-wise constant hazards
48		#' using the inversion method.
49		#'
50		#' @param pw (`numeric`)\cr time intervals for the piecewise constant hazard.
51		#' @param haz (`numeric`)\cr piecewise constant hazard.
52		#' @param LogU (`numeric`)\cr transformed uniformly distributed random variables (log(1-U)).
53		#'
54		#' @return This returns one single event time.
55		#' @export
56		#'
57		#' @examples
58		#' PCWInversionMethod(haz = c(1.1, 0.5, 0.4), pw = c(0, 7, 10), LogU = log(1 - runif(1)))
59		PCWInversionMethod <- function(haz, pw, LogU) {
60	556834x	assert_numeric(haz, lower = 0, any.missing = FALSE, all.missing = FALSE)
61	556834x	assert_intervals(pw, length(haz))
62	556834x	assert_number(LogU)
63
64	556834x	n <- length(pw)
65	556834x	t1 <- 0
66		# Determine sum of alpha*time-interval for all i.
67	556834x	dt <- pw[2:n] - pw[1:(n - 1)]
68
69		# Helping matrix.
70	556834x	tempMatrix <- matrix(0, nrow = n, ncol = n - 1)
71	556834x	tempMatrix[lower.tri(tempMatrix)] <- 1
72
73	556834x	sumA <- -as.vector(tempMatrix %% (haz[1:(n - 1)] dt))
74		# Find the appropriate time interval.
75	556834x	rev_index <- findInterval(LogU, rev(sumA), rightmost.closed = FALSE, left.open = TRUE, all.inside = FALSE)
76		# Use reverse index.
77	556834x	index <- n - rev_index
78		# Calculate event time.
79	556834x	pw[index] + (sumA[index] - LogU) / haz[index]
80		}

1		#' Event Times Distributed as Sum of Weibull
2		#'
3		#' This returns event times with a distribution resulting from the sum of two Weibull distributed random variables
4		#' using the inversion method.
5		#'
6		#' @param U (`numeric`)\cr uniformly distributed random variables.
7		#' @param haz1 (positive `number`)\cr first summand (constant hazard).
8		#' @param haz2 (positive `number`)\cr second summand (constant hazard).
9		#' @param p1 (positive `number`)\cr rate parameter of Weibull distribution for `haz1`.
10		#' @param p2 (positive `number`)\cr rate parameter of Weibull distribution for `haz2`.
11		#' @param entry (`numeric`)\cr the entry times in the current state.
12		#'
13		#' @return This returns a vector with event times.
14		#' @export
15		#'
16		#' @examples
17		#' getWaitTimeSum(U = c(0.4, 0.5), haz1 = 0.8, haz2 = 1, p1 = 1.1, p2 = 1.5, entry = c(0, 0))
18		getWaitTimeSum <- function(U, haz1, haz2, p1, p2, entry) {
19	6232x	assert_numeric(U, lower = 0, upper = 1, any.missing = FALSE, all.missing = FALSE)
20	6232x	assert_positive_number(haz1, zero_ok = TRUE)
21	6232x	assert_positive_number(haz2, zero_ok = TRUE)
22	6232x	assert_positive_number(p1)
23	6232x	assert_positive_number(p2)
24	6232x	assert_numeric(entry, lower = 0, len = length(U), any.missing = FALSE, all.missing = FALSE)
25
26	6232x	N <- length(U)
27		# Exponential distributed survival times.
28	6232x	if (p1 == 1 && p2 == 1) {
29	2187x	return(-log(1 - U) / (haz1 + haz2))
30		}
31		# Weibull distributed survival times.
32	4045x	temp <- function(x, y, t0) {
33	15888829x	return(haz1 * (x + t0)^p1 - haz1 * t0^p1 - haz2 * t0^p2 + haz2 * (x + t0)^p2 + y)
34		}
35	4045x	stime <- NULL
36	4045x	i <- 1
37	4045x	while (length(stime) < N) {
38	696412x	u <- U[i]
39	696412x	t0 <- entry[i]
40	696412x	if (temp(0, log(1 - u), t0) * temp(10000, log(1 - u), t0) < 0) {
41	696412x	res <- stats::uniroot(temp, c(0, 10000), tol = 10^-16, y = log(1 - u), t0 = t0)
42	696412x	stime[i] <- res$root
43	696412x	i <- i + 1
44	696412x	if (res$root == 0) {
45	!	warning("uniroot: accuracy (tol=10^-16) is not enough.")
46		}
47		} else {
48	!	warning("uniroot: Values at interval endpoints (interval=c(0,10000)) are not of opposite signs. \n")
49		}
50		}
51	4045x	stime
52		}

1		#' Time-point by which a specified number of events occurred.
2		#'
3		#' This returns the study time-point by which a specified number of events (PFS or OS) occurred.
4		#'
5		#' @param data (`data.frame`)\cr illness-death data set in `1rowPatient` format.
6		#' @param eventNum (`int`)\cr number of events.
7		#' @param typeEvent (`string`)\cr type of event. Possible values are `PFS` and `OS`.
8		#' @param byArm (`logical`)\cr if `TRUE` time-point per treatment arm, else joint evaluation
9		#' of treatment arms.
10		#'
11		#' @return This returns the time-point by which `eventNum` of `typeEvent`-events occurred.
12		#' @export
13		#'
14		#' @examples
15		#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
16		#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
17		#'
18		#' simStudy <- getOneClinicalTrial(
19		#' nPat = c(20, 20), transitionByArm = list(transition1, transition2),
20		#' dropout = list(rate = 0.3, time = 10),
21		#' accrual = list(param = "time", value = 0)
22		#' )
23		#' simStudyWide <- getDatasetWideFormat(simStudy)
24		#' getTimePoint(simStudyWide, eventNum = 10, typeEvent = "OS", byArm = FALSE)
25		getTimePoint <- function(data, eventNum, typeEvent, byArm = FALSE) {
26	112x	assert_data_frame(data, ncols = 11)
27	112x	assert_int(eventNum, lower = 1L)
28	112x	assert_choice(typeEvent, c("OS", "PFS"))
29	112x	assert_logical(byArm)
30
31	112x	if (!byArm) {
32	110x	data$trt <- "all"
33		}
34	112x	byVar <- unique(data$trt)
35
36	112x	if (typeEvent == "OS") {
37	56x	data$event <- data$OSevent
38	56x	data$time <- data$OStimeCal
39	56x	} else if (typeEvent == "PFS") {
40	56x	data$event <- data$PFSevent
41	56x	data$time <- data$PFStimeCal
42		}
43
44	112x	sapply(byVar, function(trt) {
45	114x	dataTemp <- data[data$trt == trt, ]
46	114x	sortedTimes <- sort(dataTemp$time[dataTemp$event == 1])
47	114x	if (eventNum > length(sortedTimes)) {
48	!	warning("Less events occur till EOS than specified \n")
49		}
50	114x	timePoint <- sortedTimes[eventNum]
51		})
52		}
53
54		#' Helper function for `censoringByNumberEvents`
55		#'
56		#' @param time (`numeric`) \cr event times.
57		#' @param event (`numeric`)\cr event indicator.
58		#' @param data (`data.frame`)\cr data frame including patient id `id`, recruiting time `recruitTime`
59		#' and individual censoring time `censTimeInd`.
60		#'
61		#' @return This function returns a data frame with columns:
62		#' event time, censoring indicator, event indicator and event time
63		#' in calendar time.
64		#' @export
65		#'
66		#' @examples
67		#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
68		#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
69		#'
70		#' simStudy <- getOneClinicalTrial(
71		#' nPat = c(20, 20), transitionByArm = list(transition1, transition2),
72		#' dropout = list(rate = 0.3, time = 10),
73		#' accrual = list(param = "time", value = 7)
74		#' )
75		#' simStudyWide <- getDatasetWideFormat(simStudy)
76		#' simStudyWide$censTimeInd <- 5 - simStudyWide$recruitTime
77		#' NotRecruited <- simStudyWide$id[simStudyWide$censTimeInd < 0]
78		#' censoredData <- simStudyWide[!(simStudyWide$id %in% NotRecruited), ]
79		#' getCensoredData(time = censoredData$OStime, event = censoredData$OSevent, data = censoredData)
80		getCensoredData <- function(time, event, data) {
81	217x	assert_numeric(time, lower = 0)
82	217x	assert_numeric(event, lower = 0, upper = 1)
83	217x	assert_data_frame(data)
84
85		# Event time censored?
86	217x	epsilon <- 1e-10
87	217x	Censored <- data$id[(time - epsilon) > data$censTimeInd]
88		# keep minimum of censoring and event time. epsilon)
89	217x	Censoredtime <- pmin(time, data$censTimeInd)
90		# adjust event and censoring indicators.
91	217x	CensoredEvent <- ifelse(data$id %in% Censored, 0, event)
92	217x	Censored <- abs(1 - CensoredEvent)
93		# calculate corresponding calendar times.
94	217x	timeCal <- Censoredtime + data$recruitTime
95
96	217x	data.frame(
97	217x	time = Censoredtime,
98	217x	Censored = Censored,
99	217x	event = CensoredEvent,
100	217x	timeCal = timeCal
101		)
102		}
103
104		#' Event-driven censoring.
105		#'
106		#' This function censors a study after a pre-specified number of events occurred.
107		#'
108		#' @param data (`data.frame`)\cr illness-death data set in `1rowPatient` format.
109		#' @param eventNum (`int`)\cr number of events.
110		#' @param typeEvent (`string`)\cr type of event. Possible values are `PFS` and `OS`.
111		#'
112		#' @return This function returns a data set that is censored after `eventNum` of
113		#' `typeEvent`-events occurred.
114		#' @export
115		#'
116		#' @examples
117		#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
118		#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
119		#'
120		#' simStudy <- getOneClinicalTrial(
121		#' nPat = c(20, 20), transitionByArm = list(transition1, transition2),
122		#' dropout = list(rate = 0.3, time = 10),
123		#' accrual = list(param = "time", value = 7)
124		#' )
125		#' simStudyWide <- getDatasetWideFormat(simStudy)
126		#' censoringByNumberEvents(data = simStudyWide, eventNum = 20, typeEvent = "PFS")
127		censoringByNumberEvents <- function(data, eventNum, typeEvent) {
128	108x	assert_data_frame(data, ncols = 11)
129	108x	assert_int(eventNum, lower = 1L)
130	108x	assert_choice(typeEvent, c("OS", "PFS"))
131
132		# first step get timePoint (calendar time) for censoring.
133	108x	censTime <- getTimePoint(data, eventNum, typeEvent, byArm = FALSE)
134
135		# censoring time at individual time scale.
136	108x	data$censTimeInd <- censTime - data$recruitTime
137		# patients that are not yet recruited at censoring time has to be deleted.
138	108x	NotRecruited <- data$id[data$censTimeInd < 0]
139	108x	censoredData <- data[!(data$id %in% NotRecruited), ]
140
141	108x	OSdata <- getCensoredData(time = censoredData$OStime, event = censoredData$OSevent, data = censoredData)
142	108x	PFSdata <- getCensoredData(time = censoredData$PFStime, event = censoredData$PFSevent, data = censoredData)
143
144	108x	data.frame(
145	108x	id = censoredData$id, trt = censoredData$trt,
146	108x	PFStime = PFSdata$time, PFSevent = PFSdata$event,
147	108x	OStime = OSdata$time, CensoredOS = OSdata$Censored, OSevent = OSdata$event,
148	108x	recruitTime = censoredData$recruitTime, OStimeCal = OSdata$timeCal, PFStimeCal = PFSdata$timeCal
149		)
150		}
151
152		#' Number of recruited/censored/ongoing Patients.
153		#'
154		#'
155		#' @param data (`data.frame`)\cr illness-death data set in `1rowPatient` format.
156		#' @param t (`numeric`)\cr study time-point.
157		#'
158		#' @return This function returns number of recruited patients,
159		#' number of censored and number of patients under observations.
160		#' @export
161		#'
162		#' @examples
163		#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
164		#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
165		#'
166		#' simStudy <- getOneClinicalTrial(
167		#' nPat = c(20, 20), transitionByArm = list(transition1, transition2),
168		#' dropout = list(rate = 0.6, time = 10),
169		#' accrual = list(param = "time", value = 0)
170		#' )
171		#' simStudyWide <- getDatasetWideFormat(simStudy)
172		#' getEventsAll(data = simStudyWide, t = 1.5)
173		getEventsAll <- function(data, t) {
174	3x	assert_data_frame(data, ncols = 11)
175	3x	assert_positive_number(t, zero_ok = FALSE)
176
177	3x	allRecruited <- data$id[data$recruitTime <= t]
178	3x	allCensored <- data$id[(data$OStime + data$recruitTime) <= t & data$OSevent == 0]
179	3x	allDeath <- data$id[(data$OStime + data$recruitTime) <= t & data$OSevent == 1]
180		## recruited, not death, not censored.
181	3x	allUnderObs <- data$id[data$id %in% allRecruited & !(data$id %in% allCensored) &
182	3x	!(data$id %in% allDeath)]
183	3x	numRecruited <- length(unique(allRecruited))
184	3x	numCensored <- length(unique(allCensored))
185	3x	numUnderObs <- length(unique(allUnderObs))
186
187	3x	c(
188	3x	Recruited = numRecruited,
189	3x	Censored = numCensored,
190	3x	UnderObs = numUnderObs
191		)
192		}
193
194		#' Helper Function for `trackEventsPerTrial`
195		#'
196		#' @param event (`numeric`)\cr event indicator.
197		#' @param time (`numeric`) \cr event times.
198		#' @param t (`numeric`)\cr study time-point.
199		#'
200		#' @return This function returns the number of events occurred until time t.
201		#' @export
202		#'
203		#' @examples
204		#' event <- c(0, 1, 1, 1, 0)
205		#' time <- c(3, 3.4, 5, 6, 5.5)
206		#' getNumberEvents(event = event, time = time, t = 5)
207		getNumberEvents <- function(event, time, t) {
208	5x	assert_numeric(event, lower = 0, upper = 1)
209	5x	assert_numeric(time, lower = 0)
210	5x	assert_number(t, lower = 0)
211
212	5x	sum(event[time <= t])
213		}
214
215
216		#' Event tracking in an oncology trial.
217		#'
218		#' @param data (`data.frame`)\cr illness-death data set in `1rowPatient` format.
219		#' @param timeP (`numeric`)\cr vector of study time-points.
220		#' @param byArm (`logical`)\cr if `TRUE` time-point per treatment arm, else joint evaluation of treatment arms.
221		#'
222		#' @return This function returns a data frame including number of PFS events, number of OS events,
223		#' number of recruited patients, number of censored patients and number of ongoing patients at `timeP`.
224		#'
225		#' @export
226		#'
227		#' @examples
228		#' transition1 <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 0.8, p02 = 0.9, p12 = 1)
229		#' transition2 <- weibull_transition(h01 = 1, h02 = 1.3, h12 = 1.7, p01 = 1.1, p02 = 0.9, p12 = 1.1)
230		#'
231		#' simStudy <- getOneClinicalTrial(
232		#' nPat = c(20, 20), transitionByArm = list(transition1, transition2),
233		#' dropout = list(rate = 0.3, time = 10),
234		#' accrual = list(param = "time", value = 0)
235		#' )
236		#' simStudyWide <- getDatasetWideFormat(simStudy)
237		#' trackEventsPerTrial(data = simStudyWide, timeP = 1.5, byArm = FALSE)
238		trackEventsPerTrial <- function(data, timeP, byArm = FALSE) {
239	1x	assert_data_frame(data, ncols = 11)
240	1x	assert_numeric(timeP, lower = 0)
241	1x	assert_logical(byArm)
242
243	1x	if (!byArm) {
244	!	data$trt <- "all"
245		}
246	1x	byVar <- unique(data$trt)
247	1x	allNumbers <- lapply(byVar, function(j) {
248	2x	datTemp <- data[data$trt == j, ]
249	2x	eventsPFS <- sapply(timeP, getNumberEvents, event = datTemp$PFSevent, time = datTemp$PFStimeCal)
250	2x	eventsOS <- sapply(timeP, getNumberEvents, event = datTemp$OSevent, time = datTemp$OStimeCal)
251	2x	eventsTrial <- sapply(timeP, getEventsAll, data = datTemp)
252
253	2x	allNumbers <- rbind(eventsPFS, eventsOS, eventsTrial)
254	2x	rownames(allNumbers) <- c("PFS", "OS", "Recruited", "Censored", "Ongoing")
255	2x	colnames(allNumbers) <- paste0("Timepoint: ", round(timeP, 2))
256	2x	allNumbers
257		})
258	1x	names(allNumbers) <- byVar
259	1x	allNumbers
260		}

1		# survPFS ----
2
3		#' PFS Survival Function for Different Transition Models
4		#'
5		#' @param transition (`TransitionParameters`)\cr
6		#' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
7		#' @param t (`numeric`)\cr time at which the value of the PFS survival function is to be computed.
8		#'
9		#' @return The value of the survival function for the specified transition and time.
10		#' @export
11		#'
12		#' @examples
13		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
14		#' survPFS(transition, 0.4)
15		survPFS <- function(transition, t) {
16	697x	UseMethod("survPFS")
17		}
18
19		#' @describeIn survPFS Survival Function for an exponential transition model.
20		#' @export
21		#'
22		#' @examples
23		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
24		#' survPFS(transition, 0.4)
25		survPFS.ExponentialTransition <- function(transition, t) {
26	533x	ExpSurvPFS(t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02)
27		}
28
29		#' @describeIn survPFS Survival Function for a Weibull transition model.
30		#' @export
31		#'
32		#' @examples
33		#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
34		#' survPFS(transition, 0.4)
35		survPFS.WeibullTransition <- function(transition, t) {
36	153x	WeibSurvPFS(
37	153x	t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02,
38	153x	p01 = transition$weibull_rates$p01, p02 = transition$weibull_rates$p02
39		)
40		}
41
42		#' @describeIn survPFS Survival Function for a piecewise constant transition model.
43		#' @export
44		#'
45		#' @examples
46		#' transition <- piecewise_exponential(
47		#' h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1),
48		#' pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9)
49		#' )
50		#' survPFS(transition, 0.4)
51		survPFS.PWCTransition <- function(transition, t) {
52	11x	PWCsurvPFS(t,
53	11x	h01 = transition$hazards$h01, h02 = transition$hazards$h02,
54	11x	pw01 = transition$intervals$pw01, pw02 = transition$intervals$pw02
55		)
56		}
57
58		# survOS ----
59
60		#' OS Survival Function for Different Transition Models
61		#'
62		#' @param transition (`TransitionParameters`)\cr
63		#' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
64		#' @param t (`numeric`)\cr time at which the value of the OS survival function is to be computed.
65		#'
66		#' @return The value of the survival function for the specified transition and time.
67		#' @export
68		#'
69		#' @examples
70		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
71		#' survOS(transition, 0.4)
72		survOS <- function(transition, t) {
73	25x	UseMethod("survOS")
74		}
75
76		#' @describeIn survOS Survival Function for an exponential transition model.
77		#' @export
78		#'
79		#' @examples
80		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
81		#' survOS(transition, 0.4)
82		survOS.ExponentialTransition <- function(transition, t) {
83	13x	ExpSurvOS(
84	13x	t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02,
85	13x	h12 = transition$hazards$h12
86		)
87		}
88
89		#' @describeIn survOS Survival Function for a Weibull transition model.
90		#' @export
91		#'
92		#' @examples
93		#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
94		#' survOS(transition, 0.4)
95		survOS.WeibullTransition <- function(transition, t) {
96	11x	WeibSurvOS(
97	11x	t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02,
98	11x	h12 = transition$hazards$h12, p01 = transition$weibull_rates$p01,
99	11x	p02 = transition$weibull_rates$p02, p12 = transition$weibull_rates$p12
100		)
101		}
102
103		#' @describeIn survOS Survival Function for a piecewise constant transition model.
104		#' @export
105		#'
106		#' @examples
107		#' transition <- piecewise_exponential(
108		#' h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1),
109		#' pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9)
110		#' )
111		#' survOS(transition, 0.4)
112		survOS.PWCTransition <- function(transition, t) {
113	1x	PWCsurvOS(
114	1x	t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02,
115	1x	h12 = transition$hazards$h12, pw01 = transition$intervals$pw01,
116	1x	pw02 = transition$intervals$pw02, pw12 = transition$intervals$pw12
117		)
118		}
119
120		# expval ----
121
122		#' Helper Function for Computing E(PFS^2)
123		#'
124		#' @param x (`numeric`)\cr variable of integration.
125		#' @param transition (`TransitionParameters`)\cr
126		#' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
127		#'
128		#' @return Numeric results of the integrand used to calculate E(PFS^2).
129		#' @export
130		#'
131		#' @examples
132		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
133		#' expvalPFSInteg(0.4, transition)
134		expvalPFSInteg <- function(x, transition) {
135	9x	x * survPFS(transition, x)
136		}
137
138		#' Helper Function for Computing E(OS^2)
139		#'
140		#' @param x (`numeric`)\cr variable of integration.
141		#' @param transition (`TransitionParameters`)\cr
142		#' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
143		#'
144		#' @return Numeric results of the integrand used to calculate E(OS^2).
145		#' @export
146		#'
147		#' @examples
148		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
149		#' expvalOSInteg(0.4, transition)
150		expvalOSInteg <- function(x, transition) {
151	11x	x * survOS(transition = transition, t = x)
152		}
153
154		# p11 ----
155
156		#' Helper Function for `log_p11()`
157		#'
158		#' @param x (`numeric`)\cr variable of integration.
159		#' @param transition (`TransitionParameters`)\cr
160		#' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
161		#'
162		#' @return Hazard rate at the specified time for the transition from progression to death.
163		#' @keywords internal
164		p11Integ <- function(x, transition) {
165	8974x	haz(transition = transition, t = x, trans = 3)
166		}
167
168		#' Probability of Remaining in Progression Between Two Time Points for Different Transition Models
169		#'
170		#' @param transition (`TransitionParameters`)\cr
171		#' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
172		#' @param s (`numeric`)\cr lower time points.
173		#' @param t (`numeric`)\cr higher time points.
174		#' @return This returns the natural logarithm of the probability of remaining in progression (state 1)
175		#' between two time points, conditional on being in state 1 at the lower time point.
176		#'
177		#' @export
178		#'
179		#' @examples
180		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
181		#' log_p11(transition, 1, 3)
182		log_p11 <- function(transition, s, t) {
183	431x	assert_numeric(s, finite = TRUE, any.missing = FALSE, lower = 0)
184	431x	assert_numeric(t, finite = TRUE, any.missing = FALSE, lower = 0)
185	431x	assert_true(identical(length(s), length(t)))
186	431x	assert_true(all(t > s))
187
188	431x	intval <- mapply(function(s, t) {
189	8973x	stats::integrate(p11Integ,
190	8973x	lower = s,
191	8973x	upper = t,
192	8973x	transition
193	8973x	)$value
194	431x	}, s, t)
195	431x	-intval
196		}
197
198		# PFSOS ----
199
200		#' Helper Function for `survPFSOS()`
201		#'
202		#' @param u (`numeric`)\cr variable of integration.
203		#' @param t (`numeric`)\cr time at which the value of the PFS*OS survival function is to be computed.
204		#' @param transition (`TransitionParameters`)\cr
205		#' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
206		#'
207		#' @note Not all vectors `u` and `t` work here due to assertions in [log_p11()].
208		#'
209		#' @return Numeric result of the integrand used to calculate the PFS*OS survival function.
210		#' @keywords internal
211		PFSOSInteg <- function(u, t, transition) {
212	429x	exp(log_p11(transition, u, t / u) + log(survPFS(transition, u)) + log(haz(transition, u, 1)))
213		}
214
215		#' Survival Function of the Product PFS*OS for Different Transition Models
216		#'
217		#' @param t (`numeric`)\cr time at which the value of the PFS*OS survival function is to be computed.
218		#' @param transition (`TransitionParameters`)\cr
219		#' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
220		#'
221		#' @return This returns the value of PFS*OS survival function at time t.
222		#' @export
223		#'
224		#' @examples
225		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
226		#' survPFSOS(0.4, transition)
227		survPFSOS <- function(t, transition) {
228	19x	sapply(t, function(x) {
229	249x	intval <- stats::integrate(PFSOSInteg, lower = 0, upper = sqrt(x), x, transition)$value
230	249x	survPFS(transition, sqrt(x)) + intval
231		})
232		}
233
234		# correlation ----
235
236		#' Correlation of PFS and OS event times for Different Transition Models
237		#'
238		#' @param transition (`TransitionParameters`)\cr
239		#' see [exponential_transition()], [weibull_transition()] or [piecewise_exponential()] for details.
240		#'
241		#' @return The correlation of PFS and OS.
242		#' @export
243		#'
244		#' @examples
245		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
246		#' corTrans(transition)
247		corTrans <- function(transition) {
248		# E(PFS) & E(OS).
249	2x	expvalPFS <- stats::integrate(survPFS,
250	2x	lower = 0, upper = Inf,
251	2x	transition = transition
252	2x	)$value
253
254	2x	expvalOS <- stats::integrate(survOS,
255	2x	lower = 0, upper = Inf,
256	2x	transition = transition
257	2x	)$value
258
259		# Var(PFS) & Var(OS).
260	2x	expvalPFS2 <- 2 * stats::integrate(expvalPFSInteg,
261	2x	lower = 0, upper = Inf,
262	2x	transition = transition
263	2x	)$value
264
265	2x	expvalOS2 <- 2 * stats::integrate(expvalOSInteg,
266	2x	lower = 0, upper = Inf,
267	2x	transition = transition
268	2x	)$value
269
270	2x	varPFS <- expvalPFS2 - expvalPFS^2
271
272	2x	varOS <- expvalOS2 - expvalOS^2
273
274		# E(PFS*OS).
275	2x	expvalPFSOS <- stats::integrate(survPFSOS,
276	2x	lower = 0, upper = Inf,
277	2x	transition
278	2x	)$value
279
280		# Cor(PFS, OS).
281	2x	(expvalPFSOS - expvalPFS * expvalOS) / sqrt(varPFS * varOS)
282		}
283
284		#' Correlation of PFS and OS event times for data from the IDM
285		#'
286		#' @param data (`data.frame`)\cr in the format produced by [getOneClinicalTrial()].
287		#' @param transition (`TransitionParameters` object)\cr specifying the assumed distribution of transition hazards.
288		#' Initial parameters for optimization can be specified here.
289		#' See [exponential_transition()] or [weibull_transition()] for details.
290		#' @param bootstrap (`flag`)\cr if `TRUE` computes confidence interval via bootstrap.
291		#' @param bootstrap_n (`count`)\cr number of bootstrap samples.
292		#' @param conf_level (`proportion`)\cr confidence level for the confidence interval.
293		#'
294		#' @return The correlation of PFS and OS.
295		#' @export
296		#'
297		#' @examples
298		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
299		#' data <- getClinicalTrials(
300		#' nRep = 1, nPat = c(100), seed = 1234, datType = "1rowTransition",
301		#' transitionByArm = list(transition), dropout = list(rate = 0.5, time = 12),
302		#' accrual = list(param = "intensity", value = 7)
303		#' )[[1]]
304		#' corPFSOS(data, transition = exponential_transition(), bootstrap = FALSE)
305		#' \dontrun{
306		#' corPFSOS(data, transition = exponential_transition(), bootstrap = TRUE)
307		#' }
308		corPFSOS <- function(data, transition, bootstrap = TRUE, bootstrap_n = 100, conf_level = 0.95) {
309	1x	assert_data_frame(data)
310	1x	assert_flag(bootstrap)
311	1x	assert_count(bootstrap_n)
312	1x	assert_number(conf_level, lower = 0.01, upper = 0.999)
313
314	1x	trans <- estimateParams(data, transition)
315	1x	res <- list("corPFSOS" = corTrans(trans))
316	1x	if (bootstrap) {
317	!	oplan <- future::plan(future::multisession, workers = parallelly::availableCores(omit = 1))
318	!	on.exit(future::plan(oplan), add = TRUE)
319	!	ids <- lapply(1:bootstrap_n, function(x) sample(seq_len(nrow(data)), nrow(data), replace = TRUE))
320	!	corBootstrap <- furrr::future_map_dbl(ids, ~ {
321	!	furrr::furrr_options(
322	!	globals = list(data = data, transition = transition),
323	!	packages = c("simIDM")
324		)
325	!	b_sample <- data[.x, , drop = FALSE]
326	!	b_transition <- estimateParams(b_sample, transition)
327	!	corTrans(b_transition)
328		})
329	!	lowerQuantile <- (1 - conf_level) / 2
330	!	upperQuantile <- lowerQuantile + conf_level
331	!	c(stats::quantile(corBootstrap, lowerQuantile),
332	!	"corPFSOS" = res,
333	!	stats::quantile(corBootstrap, upperQuantile)
334		)
335	!	res$lower <- stats::quantile(corBootstrap, lowerQuantile)
336	!	res$upper <- stats::quantile(corBootstrap, upperQuantile)
337		}
338	1x	res
339		}

1		#' Transitions from the Intermediate State to the Absorbing State
2		#'
3		#' This function creates transition entry and exit times from the intermediate state to the absorbing state
4		#' for an existing data frame containing the exit times out of the initial state.
5		#'
6		#' @param simDataOne (`data.frame`)\cr a data frame containing all patients with transitions
7		#' into the intermediate state. See [getSimulatedData()] for details.
8		#' @param transition (`TransitionParameters`)\cr transition parameters comprising
9		#' `hazards`, corresponding `intervals` and `weibull_rates`, see [exponential_transition()], [piecewise_exponential()]
10		#' and [weibull_transition()] for details.
11		#'
12		#' @return This returns a data frame with one row per patient for the second transition,
13		#' i.e. the transition out of the intermediate
14		#' state. This is a helper function of [getSimulatedData()].
15		#'
16		#' @export
17		#'
18		#' @examples
19		#' simDataOne <- data.frame(
20		#' id = c(1:3), to = c(1, 1, 1), from = c(0, 0, 0), entry = c(0, 0, 0),
21		#' exit = c(3, 5.6, 7.2), censTime = c(6.8, 5.9, 9.4)
22		#' )
23		#' transition <- exponential_transition(1, 1.6, 0.3)
24		#' getOneToTwoRows(simDataOne, transition)
25		getOneToTwoRows <- function(simDataOne, transition) {
26	6207x	assert_data_frame(simDataOne, ncols = 6)
27	6207x	assert_class(transition, "TransitionParameters")
28
29	6207x	id1 <- simDataOne$id
30	6207x	N1 <- nrow(simDataOne)
31	6207x	U1 <- stats::runif(N1)
32	6207x	to1 <- rep(2, N1)
33	6207x	from1 <- rep(1, N1)
34	6207x	entry1 <- simDataOne$exit
35	6207x	h12 <- transition$hazard$h12
36	6207x	p12 <- transition$weibull_rates$p12
37	6207x	pw12 <- transition$intervals$pw12
38
39		# Create waiting time in state 1.
40	6207x	wait_time1 <- if (transition$family == "exponential") {
41	2178x	(-log(1 - U1)) / (h12)
42	6207x	} else if (transition$family == "Weibull") {
43	2026x	getWaitTimeSum(U = U1, haz1 = h12, haz2 = 0, p1 = p12, p2 = 1, entry = entry1)
44	6207x	} else if (transition$family == "piecewise exponential") {
45	2003x	getPCWDistr(U = U1, haz = h12, pw = pw12, t_0 = entry1)
46		}
47	6207x	exit1 <- entry1 + wait_time1
48
49		# Add censoring.
50	6207x	censTime1 <- simDataOne$censTime
51	6207x	to1 <- ifelse(censTime1 < exit1, "cens", to1)
52	6207x	exit1 <- pmin(censTime1, exit1)
53
54	6207x	data.frame(
55	6207x	id = id1,
56	6207x	from = from1,
57	6207x	to = to1,
58	6207x	entry = entry1,
59	6207x	exit = exit1,
60	6207x	censTime = censTime1,
61	6207x	stringsAsFactors = FALSE
62		)
63		}
64
65		#' Staggered Study Entry
66		#'
67		#' This function adds staggered study entry times to a simulated data set with illness-death model structure.
68		#'
69		#' @param simData (`data.frame`)\cr simulated data frame containing entry and exit times
70		#' at individual study time scale. See [getSimulatedData()] for details.
71		#' @param N (`int`)\cr number of patients.
72		#' @param accrualParam (`string`)\cr possible values are 'time' or 'intensity'.
73		#' @param accrualValue (`number`)\cr specifies the accrual intensity. For `accrualParam` equal time,
74		#' it describes the length of the accrual period. For `accrualParam` equal intensity, it describes
75		#' the number of patients recruited per time unit. If `accrualValue` is equal to 0,
76		#' all patients start at calendar time 0
77		#' in the initial state.
78		#'
79		#' @return This returns a data set containing a single simulated study containing accrual times,
80		#' i.e. staggered study entry.
81		#' This is a helper function of [getSimulatedData()].
82		#'
83		#' @export
84		#'
85		#' @examples
86		#' simData <- data.frame(
87		#' id = c(1, 1, 2, 3), from = c(0, 1, 0, 0), to = c(1, 2, "cens", 2),
88		#' entry = c(0, 3, 0, 0),
89		#' exit = c(3, 5.3, 5.6, 7.2), censTime = c(6.8, 6.8, 5.6, 9.4)
90		#' )
91		#' addStaggeredEntry(simData, 3, accrualParam = "time", accrualValue = 5)
92		addStaggeredEntry <- function(simData, N, accrualParam, accrualValue) {
93	6214x	assert_choice(accrualParam, c("time", "intensity"))
94	6214x	assert_number(accrualValue, lower = 0)
95
96		# Get accrual times in calendar time per individual.
97		# If no staggered study entry is present, all individuals have the same entry time 0.
98	6214x	entry_act <- if (accrualValue != 0) {
99	6200x	if (accrualParam == "time") {
100	13x	stats::runif(N, 0, accrualValue)
101	6187x	} else if (accrualParam == "intensity") {
102	6187x	accrualTime <- N / accrualValue
103	6187x	stats::runif(N, 0, accrualTime)
104		}
105	6214x	} else if (accrualValue == 0) {
106	14x	rep(0, N)
107		}
108	6214x	entryAct <- cbind(id = 1:N, entry_act)
109		# Combine simulated data with actual entry time.
110		# Generate actual times (actual entry time + individual time).
111	6214x	simData <- merge(simData, entryAct)
112	6214x	simData$entryAct <- simData$entry + simData$entry_act
113	6214x	simData$exitAct <- simData$exit + simData$entry_act
114	6214x	simData$censAct <- simData$censTime + simData$entry_act
115		# Delete temporary helper columns.
116	6214x	simData$entry_act <- NULL
117	6214x	simData[, c("censTime")] <- list(NULL)
118	6214x	simData
119		}
120
121		#' Simulate Data Set from an Illness-Death Model
122		#'
123		#' This function creates a single simulated data set for a single treatment arm. It simulates data
124		#' from an illness-death model with one row per transition and subject.
125		#'
126		#' @param N (`int`)\cr number of patients.
127		#' @param transition (`TransitionParameters`)\cr transition parameters comprising
128		#' `hazards`, corresponding `intervals` and `weibull_rates`, see [exponential_transition()], [piecewise_exponential()]
129		#' and [weibull_transition()] for details.
130		#' @param dropout (`list`)\cr specifies drop-out probability.
131		#' Random censoring times are generated using exponential distribution. `dropout$rate` specifies
132		#' the drop-out probability per `dropout$time` time units.
133		#' If `dropout$rate` is equal to 0, then no censoring is applied.
134		#' @param accrual (`list`)\cr specifies accrual intensity. See [addStaggeredEntry()] for details.
135		#'
136		#' @return This returns a data frame with one row per transition per individual.
137		#' @details The output data set contains the following columns:
138		#' - id (`integer`): patient id.
139		#' - from (`numeric`): starting state of the transition.
140		#' - to (`character`): final state of the transition.
141		#' - entry (`numeric`): entry time of the transition on the individual time scale.
142		#' - exit (`numeric`): exit time of the transition on the individual time scale.
143		#' - entryAct (`numeric`): entry time of the transition on study time scale.
144		#' - exitAct (`numeric`): exit time of the transition on study time scale.
145		#' - censAct (`numeric`): censoring time of the individual on study time scale.
146		#' @export
147		#'
148		#' @examples
149		#' getSimulatedData(
150		#' N = 10,
151		#' transition = exponential_transition(h01 = 1, h02 = 1.5, h12 = 1),
152		#' dropout = list(rate = 0.3, time = 1),
153		#' accrual = list(param = "time", value = 5)
154		#' )
155		getSimulatedData <- function(N,
156		transition = exponential_transition(h01 = 1, h02 = 1, h12 = 1),
157		dropout = list(rate = 0, time = 12),
158		accrual = list(param = "time", value = 0)) {
159		# Check input parameters.
160	6210x	assert_int(N, lower = 1L)
161	6210x	assert_class(transition, "TransitionParameters")
162	6210x	assert_list(dropout)
163
164	6210x	assert(check_number(dropout$rate, lower = 0, upper = 1),
165	6210x	check_number(dropout$time),
166	6210x	check_true(dropout$time > 0, ),
167	6210x	combine = "and", .var.name = "dropout"
168		)
169	6210x	assert_list(accrual)
170		# Initialize transition hazards, Weibull rates and intervals.
171	6210x	h <- transition$hazard
172	6210x	p <- transition$weibull_rates
173	6210x	pw <- transition$intervals
174
175		# Get rate parameter for exponential distributed censoring times.
176		# if rate=0, no censoring is applied.
177	6210x	censRate <- if (dropout$rate > 0) {
178	198x	-log(1 - dropout$rate) / dropout$time
179		} else {
180	6012x	0
181		}
182		# Censoring times for all individuals, infinity if no censoring is applied.
183	6210x	censTime <- if (dropout$rate > 0) {
184	198x	stats::rexp(N, censRate)
185		} else {
186	6012x	rep(Inf, N)
187		}
188
189		# All individuals start in the initial state.
190	6210x	entry <- rep(0, N)
191	6210x	from <- rep(0, N)
192		# Waiting time in the initial state 0.
193	6210x	U <- stats::runif(N)
194		# Exponential or Weibull distributed survival times.
195	6210x	if (transition$family %in% c("exponential", "Weibull")) {
196		# For distribution of waiting time in the initial state the all-cause hazard (h01+h02) is needed.
197	4203x	wait_time <- getWaitTimeSum(U, h$h01, h$h02, p$p01, p$p02, entry)
198		# A binomial experiment decides on death or progression.
199	4203x	numerator <- p$p01 * h$h01 * wait_time^(p$p01 - 1) # nolint
200	4203x	denumerator <- numerator + p$p02 * h$h02 * wait_time^(p$p02 - 1) # nolint
201
202		# Piecewise exponential distributed survival times.
203	2007x	} else if (transition$family == "piecewise exponential") {
204	2007x	Sum_PCW <- getSumPCW(h$h01, h$h02, pw$pw01, pw$pw02)
205	2007x	wait_time <- getPCWDistr(U, Sum_PCW$hazards, Sum_PCW$intervals, entry)
206		# A binomial experiment decides on death or progression.
207	2007x	numerator <- getPWCHazard(h$h01, pw$pw01, wait_time)
208	2007x	denumerator <- getPWCHazard(h$h01, pw$pw01, wait_time) +
209	2007x	getPWCHazard(h$h02, pw$pw02, wait_time)
210		}
211
212	6210x	to_prob <- stats::rbinom(N, 1, numerator / denumerator)
213	6210x	to <- ifelse(to_prob == 0, 2, 1)
214	6210x	exit <- wait_time
215
216		# Add censoring.
217	6210x	to <- ifelse(censTime < wait_time, "cens", to)
218	6210x	exit <- pmin(censTime, wait_time)
219
220	6210x	simData <- data.frame(
221	6210x	id = 1:N, from = from, to = to, entry = entry, exit = exit,
222	6210x	censTime = censTime, stringsAsFactors = FALSE
223		)
224	6210x	is_intermediate_state <- simData$to == 1
225	6210x	if (any(is_intermediate_state)) {
226		# Add 1 -> 2 transition only for patients that are in the intermediate state 1.
227	6206x	simDataOne <- simData[is_intermediate_state, ]
228	6206x	newrows <- getOneToTwoRows(simDataOne, transition)
229	6206x	simData <- rbind(simData, newrows)
230	6206x	simData <- simData[order(simData$id), ]
231		}
232		# Add staggered study entry, i.e. study entry at calendar time.
233	6210x	simData <- addStaggeredEntry(
234	6210x	simData = simData,
235	6210x	N = N,
236	6210x	accrualParam = accrual$param,
237	6210x	accrualValue = accrual$value
238		)
239	6210x	simData$to <- as.character(simData$to)
240	6210x	simData
241		}

1		#' Preparation of a Data Set to Compute Log-likelihood
2		#'
3		#' @param data (`data.frame`)\cr containing entry and exit times of an illness-death model.
4		#' See [getOneClinicalTrial()] for details.
5		#'
6		#' @return This function returns a data set with one row per patient and transition, when the patient is at risk.
7		#' @export
8		#'
9		#' @details
10		#' The output data set contains the following columns:
11		#' - id (`integer`): patient id.
12		#' - from (`integer`): start event state.
13		#' - to (`integer`): end event state.
14		#' - trans (`integer`): transition (1, 2 or 3) identifier
15		#' - `1`: Transition from state 0 (stable) to 1 (progression).
16		#' - `2`: Transition from state 0 (stable) to 2 (death).
17		#' - `3`: Transition from state 1 (progression) to 2 (death).
18		#' - entry (`numeric`): time at which the patient begins to be at risk for the transition.
19		#' - exit (`numeric`): time at which the patient ends to be at risk for the transition.
20		#' - status (`logical`): event indicator for the transition.
21		#'
22		#' @examples
23		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
24		#' simData <- getOneClinicalTrial(
25		#' nPat = c(30), transitionByArm = list(transition),
26		#' dropout = list(rate = 0.8, time = 12),
27		#' accrual = list(param = "time", value = 1)
28		#' )
29		#' prepareData(simData)
30		prepareData <- function(data) {
31	8x	assert_data_frame(data, min.cols = 9, max.cols = 11)
32	8x	colNames <- c(
33	8x	"id", "trt", "PFStime", "CensoredPFS", "PFSevent", "OStime",
34	8x	"CensoredOS", "OSevent", "recruitTime",
35	8x	"OStimeCal", "PFStimeCal"
36		)
37	8x	if (!all(names(data) %in% colNames)) {
38	1x	data <- getDatasetWideFormat(data)
39		}
40
41		# Transform simIDM trial data to log-likelihood-compatible format.
42		# Suppress warning about how msprep handles 1 -> 3 transitions.
43	8x	dataNew <- suppressWarnings(mstate::msprep(
44	8x	time = c("recruitTime", "PFStime", "OStime"),
45	8x	status = c("trt", "PFSevent", "OSevent"),
46	8x	data = data,
47	8x	trans = mstate::trans.illdeath(),
48	8x	id = data$id
49		))
50	8x	cols <- which(names(dataNew) %in% c("Tstart", "Tstop"))
51	8x	names(dataNew)[cols] <- c("entry", "exit")
52		# Correct msprep results for uncensored PFS=OS events.
53	8x	ids <- data$id[data$PFStimeCal == data$OStimeCal & data$CensoredPFS == 0]
54	8x	dataNew <- dataNew[!(dataNew$id %in% ids & dataNew$trans == 3), ]
55	8x	dataNew$status[dataNew$id %in% ids] <- abs(dataNew$status[dataNew$id %in% ids] - 1)
56
57	8x	as.data.frame(dataNew[, -which(names(dataNew) == "time")], row.names = seq_len(nrow(dataNew)))
58		}
59
60		#' Compute the Negative Log-Likelihood for a Given Data Set and Transition Model
61		#'
62		#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr
63		#' see [exponential_transition()] or [weibull_transition()] for details.
64		#' @param data (`data.frame`)\cr in the format created by [prepareData()].
65		#'
66		#' @return The value of the negative log-likelihood.
67		#' @export
68		#'
69		#' @details
70		#' Calculates the negative log-likelihood for a given data set and transition model. It uses the hazard
71		#' and survival functions specific to the transition model.
72		#'
73		#' @examples
74		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
75		#' simData <- getOneClinicalTrial(
76		#' nPat = c(30), transitionByArm = list(transition),
77		#' dropout = list(rate = 0.8, time = 12),
78		#' accrual = list(param = "time", value = 1)
79		#' )
80		#' negLogLik(transition, prepareData(simData))
81		negLogLik <- function(transition, data) {
82	290x	with(
83	290x	data,
84	290x	-sum(log(
85	290x	haz(transition, exit, trans)^status * survTrans(transition, exit, trans) /
86	290x	survTrans(transition, entry, trans)
87		))
88		)
89		}
90
91		#' Hazard Function for Different Transition Models
92		#'
93		#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr
94		#' see [exponential_transition()] or [weibull_transition()] for details.
95		#' @param t (`numeric`)\cr time at which hazard is to be computed.
96		#' @param trans (`integer`)\cr index specifying the transition type.
97		#'
98		#' @details
99		#' The transition types are:
100		#' - `1`: Transition from state 0 (stable) to 1 (progression).
101		#' - `2`: Transition from state 0 (stable) to 2 (death).
102		#' - `3`: Transition from state 1 (progression) to 2 (death).
103		#'
104		#' @return The hazard rate for the specified transition and time.
105		#' @export
106		#'
107		#' @examples
108		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
109		#' haz(transition, 0.4, 2)
110		haz <- function(transition, t, trans) {
111	9698x	assert_class(transition, "TransitionParameters")
112	9698x	assert_numeric(t, lower = 0)
113	9698x	assert_subset(trans, c(1, 2, 3))
114	9698x	UseMethod("haz")
115		}
116
117		#' @describeIn haz for an exponential transition model.
118		#' @export
119		#'
120		#' @examples
121		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
122		#' haz(transition, 0.4, 2)
123		haz.ExponentialTransition <- function(transition, t, trans) {
124		# params (in this order): h01, h02, h12.
125	7233x	params <- unlist(transition$hazards, use.names = FALSE)
126	7233x	rep(params[trans], length(t) / length(trans))
127		}
128
129		#' @describeIn haz for the Weibull transition model.
130		#' @export
131		#'
132		#' @examples
133		#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
134		#' haz(transition, 0.4, 2)
135		haz.WeibullTransition <- function(transition, t, trans) {
136		# params (in this order): h01, h02, h12, p01, p02, p12.
137	2311x	params <- c(unlist(transition$hazards, use.names = FALSE), unlist(transition$weibull_rates, use.names = FALSE))
138	2311x	params[trans] * params[trans + 3] * t^(params[trans + 3] - 1)
139		}
140
141		#' @describeIn haz for the piecewise constant transition model.
142		#' @export
143		#'
144		#' @examples
145		#' transition <- piecewise_exponential(
146		#' h01 = c(1, 1, 1), h02 = c(1.5, 0.5, 1), h12 = c(1, 1, 1),
147		#' pw01 = c(0, 3, 8), pw02 = c(0, 6, 7), pw12 = c(0, 8, 9)
148		#' )
149		#' haz(transition, 6, 2)
150		haz.PWCTransition <- function(transition, t, trans) {
151	154x	getPWCHazard(unlist(transition$hazards[trans], use.names = FALSE),
152	154x	unlist(transition$intervals[trans], use.names = FALSE),
153	154x	x = t
154		)
155		}
156
157		#' Survival Function for Different Transition Models
158		#'
159		#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr
160		#' see [exponential_transition()] or [weibull_transition()] for details.
161		#' @param t (`numeric`)\cr time at which survival probability is to be computed.
162		#' @param trans (`integer`)\cr index specifying the transition type.
163		#'
164		#' @return The survival probability for the specified transition and time.
165		#' @export
166		#'
167		#' @examples
168		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
169		#' survTrans(transition, 0.4, 2)
170		survTrans <- function(transition, t, trans) {
171	582x	assert_class(transition, "TransitionParameters")
172	582x	assert_numeric(t, lower = 0)
173	582x	assert_subset(trans, c(1, 2, 3))
174	582x	UseMethod("survTrans")
175		}
176
177		#' @describeIn survTrans for the Exponential Transition Model
178		#' @export
179		#'
180		#' @examples
181		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
182		#' survTrans(transition, 0.4, 2)
183		survTrans.ExponentialTransition <- function(transition, t, trans) {
184		# params (in this order): h01, h02, h12.
185	187x	params <- unlist(transition$hazards, use.names = FALSE)
186	187x	exp(-params[trans] * t)
187		}
188
189		#' @describeIn survTrans for the Weibull Transition Model
190		#' @export
191		#'
192		#' @examples
193		#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
194		#' survTrans(transition, 0.4, 2)
195		survTrans.WeibullTransition <- function(transition, t, trans) {
196		# params (in this order): h01, h02, h12, p01, p02, p12.
197	395x	params <- c(unlist(transition$hazards, use.names = FALSE), unlist(transition$weibull_rates, use.names = FALSE))
198	395x	exp(-params[trans] * t^params[trans + 3])
199		}
200
201		#' Retrieve Initial Parameter Vectors for Likelihood Maximization
202		#'
203		#' @param transition (`ExponentialTransition` or `WeibullTransition`)\cr containing the initial parameters.
204		#' See [exponential_transition()] or [weibull_transition()] for details.
205		#'
206		#' @return The numeric vector of initial parameters for likelihood maximization.
207		#' @export
208		#'
209		#' @examples
210		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
211		#' getInit(transition)
212		getInit <- function(transition) {
213	5x	assert_class(transition, "TransitionParameters")
214	5x	UseMethod("getInit")
215		}
216
217		#' @describeIn getInit for the Exponential Transition Model
218		#' @export
219		#'
220		#' @examples
221		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
222		#' getInit(transition)
223		getInit.ExponentialTransition <- function(transition) {
224	3x	unlist(transition$hazards, use.names = FALSE)
225		}
226
227		#' @describeIn getInit for the Weibull Transition Model
228		#' @export
229		#'
230		#' @examples
231		#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
232		#' getInit(transition)
233		getInit.WeibullTransition <- function(transition) {
234	2x	c(unlist(transition$hazards, use.names = FALSE), unlist(transition$weibull_rates, use.names = FALSE))
235		}
236
237		#' Generate the Target Function for Optimization
238		#'
239		#' @param transition (`TransitionParameters`)\cr
240		#' specifying the distribution family. See [exponential_transition()] or [weibull_transition()] for details.
241		#'
242		#' @return Function that calculates the negative log-likelihood for the given parameters.
243		#' @export
244		#'
245		#' @details
246		#' This function creates a target function for optimization, computing the negative log-likelihood for given
247		#' parameters, data, and transition model type.
248		#'
249		#' @examples
250		#' transition <- exponential_transition(2, 1.3, 0.8)
251		#' simData <- getOneClinicalTrial(
252		#' nPat = c(30), transitionByArm = list(transition),
253		#' dropout = list(rate = 0.8, time = 12),
254		#' accrual = list(param = "time", value = 1)
255		#' )
256		#' params <- c(1.2, 1.5, 1.6) # For ExponentialTransition
257		#' data <- prepareData(simData)
258		#' transition <- exponential_transition()
259		#' fun <- getTarget(transition)
260		#' fun(params, data)
261		getTarget <- function(transition) {
262	5x	assert_class(transition, "TransitionParameters")
263	5x	UseMethod("getTarget", transition)
264		}
265
266		#' @describeIn getTarget for the Exponential Transition Model
267		#' @export
268		#'
269		#' @examples
270		#' transition <- exponential_transition(2, 1.3, 0.8)
271		#' simData <- getOneClinicalTrial(
272		#' nPat = c(30), transitionByArm = list(transition),
273		#' dropout = list(rate = 0.8, time = 12),
274		#' accrual = list(param = "time", value = 1)
275		#' )
276		#' params <- c(1.2, 1.5, 1.6)
277		#' data <- prepareData(simData)
278		#' transition <- exponential_transition()
279		#' target <- getTarget(transition)
280		#' target(params, data)
281		getTarget.ExponentialTransition <- function(transition) {
282	3x	function(params, data) {
283	92x	negLogLik(transition = exponential_transition(h01 = params[1], h02 = params[2], h12 = params[3]), data = data)
284		}
285		}
286
287		#' @describeIn getTarget for the Weibull Transition Model
288		#' @export
289		#'
290		#' @examples
291		#' transition <- weibull_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6, p01 = 2, p02 = 2.5, p12 = 3)
292		#' simData <- getOneClinicalTrial(
293		#' nPat = c(30), transitionByArm = list(transition),
294		#' dropout = list(rate = 0.8, time = 12),
295		#' accrual = list(param = "time", value = 1)
296		#' )
297		#' params <- c(1.2, 1.5, 1.6, 0.8, 1.3, 1.1)
298		#' data <- prepareData(simData)
299		#' transition <- weibull_transition()
300		#' target <- getTarget(transition)
301		#' target(params, data)
302		getTarget.WeibullTransition <- function(transition) {
303	2x	function(params, data) {
304	196x	negLogLik(transition = weibull_transition(
305	196x	h01 = params[1], h02 = params[2], h12 = params[3],
306	196x	p01 = params[4], p02 = params[5], p12 = params[6]
307	196x	), data = data)
308		}
309		}
310
311		#' Format Results of Parameter Estimation for Different Transition Models
312		#'
313		#' @param transition (`TransitionParameters`)\cr
314		#' see [exponential_transition()] or [weibull_transition()] for details.
315		#' @param res (`numeric` vector)\cr vector of parameter estimates from the likelihood maximization procedure.
316		#'
317		#' @return Returns a `TransitionParameters` object with parameter estimates.
318		#' @export
319		#'
320		#' @examples
321		#' results <- c(1.2, 1.5, 1.6)
322		#' getResults(exponential_transition(), results)
323		getResults <- function(transition, res) {
324	5x	UseMethod("getResults")
325		}
326
327		#' @describeIn getResults for the Exponential Transition Model
328		#' @export
329		#'
330		#' @examples
331		#' results <- c(1.2, 1.5, 1.6)
332		#' getResults(exponential_transition(), results)
333		getResults.ExponentialTransition <- function(transition, res) {
334	3x	exponential_transition(h01 = res[1], h02 = res[2], h12 = res[3])
335		}
336
337		#' @describeIn getResults for the Weibull Transition Model
338		#' @export
339		#'
340		#' @examples
341		#' results <- c(1.2, 1.5, 1.6, 2, 2.5, 1)
342		#' getResults(weibull_transition(), results)
343		getResults.WeibullTransition <- function(transition, res) {
344	2x	weibull_transition(
345	2x	h01 = res[1], h02 = res[2], h12 = res[3],
346	2x	p01 = res[4], p02 = res[5], p12 = res[6]
347		)
348		}
349
350		#' Estimate Parameters of the Multistate Model Using Clinical Trial Data
351		#'
352		#' @param data (`data.frame`)\cr in the format produced by [getOneClinicalTrial()].
353		#' @param transition (`TransitionParameters` object)\cr specifying the assumed distribution of transition hazards.
354		#' Initial parameters for optimization can be specified here.
355		#' See [exponential_transition()] or [weibull_transition()] for details.
356		#'
357		#' @return Returns a `TransitionParameters` object with the estimated parameters.
358		#' @export
359		#'
360		#' @details
361		#' This function estimates parameters for transition models using clinical trial data.
362		#' The `transition` object can be initialized with starting values for parameter estimation.
363		#' It uses [stats::optim()] to optimize the parameters.
364		#'
365		#' @examples
366		#' transition <- exponential_transition(h01 = 2, h02 = 1.4, h12 = 1.6)
367		#' simData <- getOneClinicalTrial(
368		#' nPat = c(30), transitionByArm = list(transition),
369		#' dropout = list(rate = 0.3, time = 12),
370		#' accrual = list(param = "time", value = 1)
371		#' )
372		#' # Initialize transition with desired starting values for optimization:
373		#' transition <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
374		#' estimate <- estimateParams(simData, transition)
375		estimateParams <- function(data, transition) {
376	3x	data <- prepareData(data)
377	3x	par <- getInit(transition)
378	3x	target <- getTarget(transition)
379
380	3x	res <- stats::optim(
381	3x	par = par,
382	3x	fn = target,
383	3x	method = "L-BFGS-B",
384	3x	lower = 1e-3,
385	3x	data = data
386	3x	)$par
387
388	3x	getResults(transition, res)
389		}

1		#' Simulation of a Single Oncology Clinical Trial
2		#'
3		#' This function creates a data set with a single simulated oncology clinical trial with one row per transition
4		#' based on an illness-death model. Studies with an arbitrary number of treatment arms are possible.
5		#'
6		#' @param nPat (`integer`)\cr numbers of patients per treatment arm.
7		#' @param transitionByArm (`list`) \cr transition parameters for each treatment group.
8		#' See [exponential_transition()], [piecewise_exponential()] and [weibull_transition()] for details.
9		#' @param dropout dropout (`list`)\cr specifies drop-out probability. See [getSimulatedData()] for details.
10		#' Can be specified either as one list that should be applied to all treatment groups or a separate list
11		#' for each treatment group.
12		#' @param accrual accrual (`list`)\cr specifies accrual intensity. See [addStaggeredEntry()] for details.
13		#' Can be specified either as one list that should be applied to all treatment groups or a separate list
14		#' for each treatment group.
15		#'
16		#' @return This returns a data frame with one simulated clinical trial and multiple treatment arms.
17		#' See [getSimulatedData()] for the explanation of the columns. The column `trt` contains the treatment indicator.
18		#' This is a helper function of [getClinicalTrials()].
19		#' @export
20		#'
21		#' @examples
22		#' transition1 <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
23		#' transition2 <- exponential_transition(h01 = 1, h02 = 1.3, h12 = 1.7)
24		#' transition3 <- exponential_transition(h01 = 1.1, h02 = 1, h12 = 1.5)
25		#' getOneClinicalTrial(
26		#' nPat = c(30, 20, 30), transitionByArm = list(transition1, transition2, transition3),
27		#' dropout = list(rate = 0, time = 12),
28		#' accrual = list(param = "time", value = 0)
29		#' )
30		getOneClinicalTrial <- function(nPat, transitionByArm,
31		dropout = list(rate = 0, time = 12),
32		accrual = list(param = "time", value = 0)) {
33	3105x	assert_list(transitionByArm)
34	3105x	nPat <- as.integer(nPat)
35	3105x	nArm <- length(transitionByArm)
36	3105x	assert_integer(nPat,
37	3105x	lower = 1,
38	3105x	any.missing = FALSE,
39	3105x	all.missing = FALSE,
40	3105x	len = nArm,
41		)
42	3105x	assert_list(dropout)
43	3105x	assert_list(accrual)
44
45		# Same accrual and dropout parameters for each group?
46	3105x	if (is.list(dropout[[1]])) {
47	11x	assert_list(dropout, len = nArm, types = "list")
48		} else {
49	3094x	dropout <- rep(list(dropout), nArm)
50		}
51
52	3105x	if (is.list(accrual[[1]])) {
53	11x	assert_list(accrual, len = nArm, types = "list")
54		} else {
55	3094x	accrual <- rep(list(accrual), nArm)
56		}
57
58		# Each loop simulates a single trial arm.
59		# Starting values for the loop.
60	3105x	simdata <- NULL
61	3105x	previousPts <- 0
62	3105x	for (i in seq_len(nArm)) {
63	6205x	group <- getSimulatedData(nPat[i], transitionByArm[[i]], dropout[[i]], accrual[[i]])
64	6205x	group$trt <- i
65	6205x	group$id <- group$id + previousPts
66	6205x	simdata <- rbind(simdata, group)
67	6205x	previousPts <- previousPts + nPat[i]
68		}
69	3105x	if (!any(simdata$from == 1)) {
70	1x	warning("no progression to death transitions included in the simulated data")
71		}
72	3105x	simdata
73		}
74
75		#' Conversion of a Data Set from One Row per Transition to One Row per Patient
76		#'
77		#' @param data (`data.frame`)\cr data frame containing entry and exit times of an illness-death model.
78		#' See [getSimulatedData()] for details.
79		#'
80		#' @return This function returns a data set with one row per patient and endpoints PFS and OS.
81		#' @export
82		#'
83		#' @details
84		#' The output data set contains the following columns:
85		#' - id (`integer`): patient id.
86		#' - trt `integer`): treatment id.
87		#' - PFStime (`numeric`): event time of PFS event.
88		#' - CensoredPFS (`logical`): censoring indicator for PFS event.
89		#' - PFSevent (`logical`): event indicator for PFS event.
90		#' - OStime (`numeric`): event time of OS event.
91		#' - CensoredOS (`logical`): censoring indicator for OS event.
92		#' - OSevent (`logical`): event indicator for OS event.
93		#' - recruitTime (`numeric`): time of recruitment.
94		#' - OStimeCal (`numeric`): OS event time at calendar time scale.
95		#' - PFStimeCal (`numeric`): PFS event time at calendar time scale.
96		#'
97		#' @examples
98		#' transition1 <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
99		#' transition2 <- exponential_transition(h01 = 1, h02 = 1.3, h12 = 1.7)
100		#' transition3 <- exponential_transition(h01 = 1.1, h02 = 1, h12 = 1.5)
101		#' simData <- getOneClinicalTrial(
102		#' nPat = c(30, 20, 30), transitionByArm = list(transition1, transition2, transition3),
103		#' dropout = list(rate = 0, time = 12),
104		#' accrual = list(param = "time", value = 0)
105		#' )
106		#' getDatasetWideFormat(simData)
107		getDatasetWideFormat <- function(data) {
108	90x	assert_data_frame(data, ncols = 9)
109	90x	assert_subset(c("id", "from", "to", "entry", "exit", "entryAct", "exitAct", "censAct", "trt"), names(data))
110
111		# Recruitment time is the actual entry time of the initial state.
112	90x	recruitTime <- subset(data[, c("id", "entryAct")], data$from == 0)
113	90x	names(recruitTime)[names(recruitTime) == "entryAct"] <- "recruitTime"
114
115		# The OS time is the entry time into state 2 or the censoring time whatever occurs first.
116	90x	OStime <- subset(data[, c("id", "exit")], data$to == 2 \| data$to == "cens")
117	90x	names(OStime)[names(OStime) == "exit"] <- "OStime"
118
119		# The PFS time is the entry time into state 1 or state 2 or the censoring time whatever occurs first.
120	90x	PFStime <- subset(
121	90x	data[, c("id", "exit")],
122	90x	(data$to == 2 & data$from == 0) \| (data$to == 1) \| (data$to == "cens" & data$from == 0)
123		)
124	90x	names(PFStime)[names(PFStime) == "exit"] <- "PFStime"
125
126		# Add censoring indicators.
127	90x	censoredIdsPFS <- data$id[data$to == "cens" & data$from == 0]
128	90x	censoredIdsOS <- data$id[data$to == "cens"]
129	90x	id <- unique(data$id)
130	90x	CensoredOS <- cbind(id = id, CensoredOS = as.integer(id %in% censoredIdsOS))
131	90x	CensoredPFS <- cbind(id = id, CensoredPFS = as.integer(id %in% censoredIdsPFS))
132
133		# Merge all data sets to one.
134	90x	newdata <- unique(data[, c("id", "trt")])
135	90x	newdata <- merge(x = newdata, y = PFStime, by = "id")
136	90x	newdata <- merge(x = newdata, y = CensoredPFS, by = "id")
137
138		# Do we have an observed PFS event?
139	90x	newdata$PFSevent <- abs(1 - newdata$CensoredPFS)
140	90x	newdata <- merge(x = newdata, y = OStime, by = "id")
141	90x	newdata <- merge(x = newdata, y = CensoredOS, by = "id")
142
143		# Do we have an observed OS event?
144	90x	newdata$OSevent <- abs(1 - newdata$CensoredOS)
145	90x	newdata <- merge(x = newdata, y = recruitTime, by = "id")
146
147		# Add variables with event times in calendar time.
148	90x	newdata$OStimeCal <- newdata$OStime + newdata$recruitTime
149	90x	newdata$PFStimeCal <- newdata$PFStime + newdata$recruitTime
150
151	90x	newdata
152		}
153
154		#' Helper Function for Adding Progress Bar to Trial Simulation
155		#'
156		#' @param x (`int`)\cr iteration index within lapply.
157		#' @param ... parameters transferred to [getOneClinicalTrial()], see [getOneClinicalTrial()] for details.
158		#'
159		#' @return This returns the same as [getOneClinicalTrial()], but updates the progress bar.
160		#' @keywords internal
161		runTrial <- function(x, pb, ...) {
162	3094x	utils::setTxtProgressBar(pb, x)
163	3094x	getOneClinicalTrial(...)
164		}
165
166		#' Simulation of a Large Number of Oncology Clinical Trials
167		#'
168		#' @param nRep (`int`)\cr number of simulated trials.
169		#' @param ... parameters transferred to [getOneClinicalTrial()], see [getOneClinicalTrial()] for details.
170		#' @param seed (`int`)\cr random seed used for this simulation.
171		#' @param datType (`string`)\cr possible values are `1rowTransition` and `1rowPatient`.
172		#'
173		#' @return This function returns a list with `nRep` simulated data sets in the format specified by `datType`.
174		#' See [getDatasetWideFormat()] [getOneClinicalTrial()] for details.
175		#' @export
176		#'
177		#' @examples
178		#' transition1 <- exponential_transition(h01 = 1.2, h02 = 1.5, h12 = 1.6)
179		#' transition2 <- exponential_transition(h01 = 1, h02 = 1.3, h12 = 1.7)
180		#' getClinicalTrials(
181		#' nRep = 10, nPat = c(20, 20), seed = 1234, datType = "1rowTransition",
182		#' transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
183		#' accrual = list(param = "intensity", value = 7)
184		#' )
185		getClinicalTrials <- function(nRep, ..., seed = 1234, datType = "1rowTransition") {
186	19x	assert_number(nRep, lower = 1)
187	19x	assert_choice(datType, c("1rowTransition", "1rowPatient"))
188
189	19x	set.seed(seed)
190	19x	cat("Simulating", nRep, "trials:\n")
191	19x	pb <- utils::txtProgressBar(min = 0, max = nRep, style = 3)
192		# getOneClinicalTrial generates a single clinical trial with multiple arms. Generate nRep simulated trials:
193	19x	simulatedTrials <- lapply(
194	19x	seq_len(nRep),
195	19x	FUN = function(x, ...) runTrial(x, pb, ...),
196		...
197		)
198	19x	close(pb)
199		# Final data set format: one row per patient or one row per transition?
200	19x	if (datType == "1rowPatient") {
201	14x	simulatedTrials <- lapply(simulatedTrials, getDatasetWideFormat)
202		}
203	19x	simulatedTrials
204		}

1		#' PFS Survival Function from Constant Transition Hazards
2		#'
3		#' @inheritParams WeibSurvOS
4		#' @return This returns the value of PFS survival function at time t.
5		#' @export
6		#'
7		#' @examples
8		#' ExpSurvPFS(c(1:5), 0.2, 0.4)
9		ExpSurvPFS <- function(t, h01, h02) {
10	800x	assert_numeric(t, lower = 0, any.missing = FALSE)
11	800x	assert_positive_number(h01, zero_ok = TRUE)
12	800x	assert_positive_number(h02, zero_ok = TRUE)
13
14	800x	exp(-(h01 + h02) * t)
15		}
16
17		#' OS Survival Function from Constant Transition Hazards
18		#'
19		#' @inheritParams WeibSurvOS
20		#'
21		#' @return This returns the value of OS survival function at time t.
22		#' @export
23		#'
24		#' @examples
25		#' ExpSurvOS(c(1:5), 0.2, 0.4, 0.1)
26		ExpSurvOS <- function(t, h01, h02, h12) {
27	130x	assert_numeric(t, lower = 0, any.missing = FALSE)
28	130x	assert_positive_number(h01, zero_ok = TRUE)
29	130x	assert_positive_number(h02, zero_ok = TRUE)
30	130x	assert_positive_number(h12, zero_ok = TRUE)
31
32	130x	h012 <- h12 - h01 - h02
33	130x	ExpSurvPFS(t, h01, h02) + h01 * h012^-1 * (ExpSurvPFS(t, h01, h02) - exp(-h12 * t))
34		}
35
36		#' PFS Survival Function from Weibull Transition Hazards
37		#'
38		#' @inheritParams WeibSurvOS
39		#' @return This returns the value of PFS survival function at time t.
40		#' @export
41		#' @export
42		#'
43		#' @examples
44		#' WeibSurvPFS(c(1:5), 0.2, 0.5, 1.2, 0.9)
45		WeibSurvPFS <- function(t, h01, h02, p01, p02) {
46	173x	assert_numeric(t, lower = 0, any.missing = FALSE)
47	173x	assert_positive_number(h01, zero_ok = TRUE)
48	173x	assert_positive_number(h02, zero_ok = TRUE)
49	173x	assert_positive_number(p01)
50	173x	assert_positive_number(p02)
51
52	173x	exp(-h01 * t^p01 - h02 * t^p02)
53		}
54
55		#' Helper for Efficient Integration
56		#'
57		#' @param integrand (`function`)\cr to be integrated.
58		#' @param upper (`numeric`)\cr upper limits of integration.
59		#' @param ... additional arguments to be passed to `integrand`.
60		#'
61		#' @return This function returns for each upper limit the estimates of the integral.
62		#' @keywords internal
63		integrateVector <- function(integrand, upper, ...) {
64	4175x	assert_true(all(upper >= 0))
65	4175x	boundaries <- sort(unique(upper))
66	4175x	nIntervals <- length(boundaries)
67	4175x	intervals <- mapply(
68	4175x	FUN = function(f, lower, upper, ...) stats::integrate(f, ..., lower, upper)$value,
69	4175x	MoreArgs = list(f = integrand, ...),
70	4175x	lower = c(0, boundaries[-nIntervals]),
71	4175x	upper = boundaries
72		)
73	4175x	cumIntervals <- cumsum(intervals)
74	4175x	cumIntervals[match(upper, boundaries)]
75		}
76
77		#' Helper Function for `WeibSurvOS()`
78		#'
79		#' @param x (`numeric`)\cr variable of integration.
80		#' @inheritParams WeibSurvOS
81		#'
82		#' @return Numeric results of the integrand used to calculate
83		#' the OS survival function for Weibull transition hazards, see `WeibSurvOS()`.
84		#'
85		#' @keywords internal
86		WeibOSInteg <- function(x, t, h01, h02, h12, p01, p02, p12) {
87	4212x	assert_numeric(x, finite = TRUE, any.missing = FALSE)
88	4212x	assert_numeric(t, finite = TRUE, any.missing = FALSE)
89	4212x	assert_true(test_scalar(x) \|\| identical(length(x), length(t)) \|\| test_scalar(t))
90	4212x	assert_positive_number(h01, zero_ok = TRUE)
91	4212x	assert_positive_number(h02, zero_ok = TRUE)
92	4212x	assert_positive_number(h12, zero_ok = TRUE)
93	4212x	assert_positive_number(p01)
94	4212x	assert_positive_number(p02)
95	4212x	assert_positive_number(p12)
96
97	4212x	x^(p01 - 1) * exp(-h01 * x^p01 - h02 * x^p02 - h12 * (t^p12 - x^p12))
98		}
99
100		#' OS Survival Function from Weibull Transition Hazards
101		#'
102		#' @param t (`numeric`)\cr study time-points.
103		#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition.
104		#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition.
105		#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition.
106		#' @param p01 (positive `number`)\cr rate parameter of Weibull distribution for `h01`.
107		#' @param p02 (positive `number`)\cr rate parameter of Weibull distribution for `h02`.
108		#' @param p12 (positive `number`)\cr rate parameter of Weibull distribution for `h12`.
109		#'
110		#' @return This returns the value of OS survival function at time t.
111		#' @export
112		#'
113		#' @examples WeibSurvOS(c(1:5), 0.2, 0.5, 2.1, 1.2, 0.9, 1)
114		WeibSurvOS <- function(t, h01, h02, h12, p01, p02, p12) {
115	16x	assert_numeric(t, lower = 0, any.missing = FALSE)
116	16x	assert_positive_number(h01, zero_ok = TRUE)
117	16x	assert_positive_number(p01)
118
119	16x	WeibSurvPFS(t, h01, h02, p01, p02) +
120	16x	h01 * p01 *
121	16x	sapply(t, function(t) {
122	4161x	integrateVector(WeibOSInteg,
123	4161x	upper = t,
124	4161x	t = t,
125	4161x	h01 = h01,
126	4161x	h02 = h02,
127	4161x	h12 = h12,
128	4161x	p01 = p01,
129	4161x	p02 = p02,
130	4161x	p12 = p12
131		)
132		})
133		}
134
135		#' Cumulative Hazard for Piecewise Constant Hazards
136		#'
137		#' @param t (`numeric`)\cr study time-points.
138		#' @param haz (`numeric vector`)\cr constant transition hazards.
139		#' @param pw (`numeric vector`)\cr time intervals for the piecewise constant hazards.
140		#'
141		#' @return This returns the value of cumulative hazard at time t.
142		#' @export
143		#'
144		#' @examples
145		#' pwA(1:5, c(0.5, 0.9), c(0, 4))
146		pwA <- function(t, haz, pw) {
147	12581x	assert_numeric(t, lower = 0, any.missing = FALSE)
148	12581x	assert_numeric(haz, lower = 0, any.missing = FALSE, all.missing = FALSE)
149	12581x	assert_intervals(pw, length(haz))
150
151		# Time intervals: time-points + cut points.
152	12581x	pw_new <- c(pw[pw < max(t)])
153	12581x	pw_time <- sort(unique(c(pw_new, t)))
154	12581x	haz_all <- getPWCHazard(haz, pw, pw_time)
155
156		# Cumulative hazard.
157	12581x	i <- 1:(length(pw_time) - 1)
158	12581x	dA_pw <- haz_all[i] * (pw_time[i + 1] - pw_time[i])
159	12581x	A_pw <- c(0, cumsum(dA_pw))
160
161		# Match result with input.
162	12581x	A_pw[match(t, pw_time)]
163		}
164
165		#' PFS Survival Function from Piecewise Constant Hazards
166		#'
167		#' @inheritParams PWCsurvOS
168		#'
169		#' @return This returns the value of PFS survival function at time t.
170		#' @export
171		#'
172		#' @examples
173		#' PWCsurvPFS(1:5, c(0.3, 0.5), c(0.5, 0.8), c(0, 4), c(0, 8))
174		PWCsurvPFS <- function(t, h01, h02, pw01, pw02) {
175	56x	assert_numeric(t, lower = 0, any.missing = FALSE)
176	56x	assert_numeric(h01, lower = 0, any.missing = FALSE, all.missing = FALSE)
177	56x	assert_numeric(h02, lower = 0, any.missing = FALSE, all.missing = FALSE)
178
179	56x	assert_intervals(pw01, length(h01))
180	56x	assert_intervals(pw02, length(h02))
181
182	56x	A01 <- pwA(t, h01, pw01)
183	56x	A02 <- pwA(t, h02, pw02)
184
185	56x	exp(-A01 - A02)
186		}
187
188		#' Helper Function of `PWCsurvOS()`
189		#'
190		#' @param x (`numeric`)\cr variable of integration.
191		#' @inheritParams PWCsurvOS
192		#'
193		#' @return Numeric results of the integrand used to calculate
194		#' the OS survival function for piecewise constant transition hazards, see `PWCsurvOS`.
195		#'
196		#' @keywords internal
197		PwcOSInt <- function(x, t, h01, h02, h12, pw01, pw02, pw12) {
198	30x	PWCsurvPFS(x, h01, h02, pw01, pw02) *
199	30x	getPWCHazard(h01, pw01, x) *
200	30x	exp(pwA(x, h12, pw12) - pwA(t, h12, pw12))
201		}
202
203		#' OS Survival Function from Piecewise Constant Hazards
204		#'
205		#' @param t (`numeric`)\cr study time-points.
206		#' @param h01 (`numeric vector`)\cr constant transition hazards for 0 to 1 transition.
207		#' @param h02 (`numeric vector`)\cr constant transition hazards for 0 to 2 transition.
208		#' @param h12 (`numeric vector`)\cr constant transition hazards for 1 to 2 transition.
209		#' @param pw01 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h01`.
210		#' @param pw02 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h02`.
211		#' @param pw12 (`numeric vector`)\cr time intervals for the piecewise constant hazards `h12`.
212		#'
213		#' @return This returns the value of OS survival function at time t.
214		#' @export
215		#'
216		#' @examples
217		#' PWCsurvOS(1:5, c(0.3, 0.5), c(0.5, 0.8), c(0.7, 1), c(0, 4), c(0, 8), c(0, 3))
218		PWCsurvOS <- function(t, h01, h02, h12, pw01, pw02, pw12) {
219	9x	assert_numeric(t, lower = 0, any.missing = FALSE)
220	9x	assert_numeric(h01, lower = 0, any.missing = FALSE, all.missing = FALSE)
221	9x	assert_numeric(h02, lower = 0, any.missing = FALSE, all.missing = FALSE)
222	9x	assert_numeric(h12, lower = 0, any.missing = FALSE, all.missing = FALSE)
223	9x	assert_intervals(pw01, length(h01))
224	9x	assert_intervals(pw02, length(h02))
225	9x	assert_intervals(pw12, length(h12))
226
227		# Assemble unique time points and corresponding hazard values once.
228	9x	unique_pw_times <- sort(unique(c(pw01, pw02, pw12)))
229	9x	h01_at_times <- getPWCHazard(h01, pw01, unique_pw_times)
230	9x	h02_at_times <- getPWCHazard(h02, pw02, unique_pw_times)
231	9x	h12_at_times <- getPWCHazard(h12, pw12, unique_pw_times)
232
233		# We work for the integral parts with sorted and unique time points.
234	9x	t_sorted <- sort(unique(t))
235	9x	t_sorted_zero <- c(0, t_sorted)
236	9x	n_t <- length(t_sorted)
237	9x	int_sums <- numeric(n_t)
238	9x	cum_haz_12 <- pwA(t_sorted, h12, pw12)
239	9x	for (i in seq_len(n_t)) {
240	4123x	t_start <- t_sorted_zero[i]
241	4123x	t_end <- t_sorted_zero[i + 1]
242		# Determine the indices of the time intervals we are working with here.
243	4123x	start_index <- findInterval(t_start, unique_pw_times)
244		# We use here left.open = TRUE, so that when t_end is identical to a change point,
245		# then we don't need an additional part.
246	4123x	end_index <- max(findInterval(t_end, unique_pw_times, left.open = TRUE), 1L)
247	4123x	index_bounds <- start_index:end_index
248		# Determine corresponding time intervals.
249	4123x	time_bounds <- if (length(index_bounds) == 1L) {
250	4116x	rbind(c(t_start, t_end))
251	4123x	} else if (length(index_bounds) == 2L) {
252	5x	rbind(
253	5x	c(t_start, unique_pw_times[end_index]),
254	5x	c(unique_pw_times[end_index], t_end)
255		)
256		} else {
257	2x	n_ind_bounds <- length(index_bounds)
258	2x	rbind(
259	2x	c(t_start, unique_pw_times[start_index + 1L]),
260	2x	cbind(
261	2x	unique_pw_times[index_bounds[-c(1, n_ind_bounds)]],
262	2x	unique_pw_times[index_bounds[-c(1, 2)]]
263		),
264	2x	c(unique_pw_times[end_index], t_end)
265		)
266		}
267	4123x	for (j in seq_along(index_bounds)) {
268	4133x	time_j <- time_bounds[j, 1]
269	4133x	time_jp1 <- time_bounds[j, 2]
270	4133x	intercept_a <- pwA(time_j, h12, pw12) - pwA(time_j, h01, pw01) - pwA(time_j, h02, pw02)
271	4133x	slope_b <- h12_at_times[index_bounds[j]] - h01_at_times[index_bounds[j]] - h02_at_times[index_bounds[j]]
272	4133x	h01_j <- h01_at_times[index_bounds[j]]
273	4133x	int_js <- if (slope_b != 0) {
274	4126x	h01_j / slope_b * exp(intercept_a) *
275	4126x	(exp(slope_b * (time_jp1 - time_j) - cum_haz_12[i:n_t]) - exp(-cum_haz_12[i:n_t]))
276		} else {
277	7x	h01_j * exp(intercept_a - cum_haz_12[i:n_t]) * (time_jp1 - time_j)
278		}
279	4133x	int_sums[i:n_t] <- int_sums[i:n_t] + int_js
280		}
281		}
282
283		# Match back to all time points,
284		# which might not be unique or ordered.
285	9x	int_sums <- int_sums[match(t, t_sorted)]
286	9x	result <- PWCsurvPFS(t, h01, h02, pw01, pw02) + int_sums
287
288		# Cap at 1 in a safe way - i.e. first check.
289	9x	assert_numeric(result, finite = TRUE, any.missing = FALSE)
290	9x	above_one <- result > 1
291	9x	if (any(above_one)) {
292	1x	assert_true(all((result[above_one] - 1) < sqrt(.Machine$double.eps)))
293	1x	result[above_one] <- 1
294		}
295	9x	result
296		}
297
298		#' Helper Function for Single Quantile for OS Survival Function
299		#'
300		#' @param q (`number`)\cr single quantile at which to compute event time.
301		#' @inheritParams ExpQuantOS
302		#'
303		#' @return Single time t such that the OS survival function at t equals q.
304		#' @keywords internal
305		singleExpQuantOS <- function(q, h01, h02, h12) {
306	5x	assert_number(q, finite = TRUE, lower = 0, upper = 1)
307	5x	toroot <- function(x) {
308	81x	return(q - ExpSurvOS(t = x, h01, h02, h12))
309		}
310	5x	stats::uniroot(
311	5x	toroot,
312	5x	interval = c(0, 10^3),
313	5x	extendInt = "yes"
314	5x	)$root
315		}
316
317		#' Quantile function for OS survival function induced by an illness-death model
318		#'
319		#' @param q (`numeric`)\cr quantile(s) at which to compute event time (q = 1 / 2 corresponds to median).
320		#' @param h01 (`numeric vector`)\cr constant transition hazards for 0 to 1 transition.
321		#' @param h02 (`numeric vector`)\cr constant transition hazards for 0 to 2 transition.
322		#' @param h12 (`numeric vector`)\cr constant transition hazards for 1 to 2 transition.
323		#'
324		#' @return This returns the time(s) t such that the OS survival function at t equals q.
325		#' @export
326		#'
327		#' @examples ExpQuantOS(1 / 2, 0.2, 0.5, 2.1)
328		ExpQuantOS <- function(q = 1 / 2, h01, h02, h12) {
329	2x	assert_numeric(q, lower = 0, upper = 1, any.missing = FALSE)
330	2x	assert_numeric(h01, lower = 0, any.missing = FALSE)
331	2x	assert_numeric(h02, lower = 0, any.missing = FALSE)
332	2x	assert_numeric(h12, lower = 0, any.missing = FALSE)
333
334	2x	vapply(
335	2x	q,
336	2x	FUN = singleExpQuantOS,
337	2x	FUN.VALUE = numeric(1),
338	2x	h01 = h01,
339	2x	h02 = h02,
340	2x	h12 = h12
341		)
342		}

1		#' OS Hazard Function from Constant Transition Hazards
2		#'
3		#' @param t (`numeric`)\cr study time-points.
4		#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition.
5		#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition.
6		#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition.
7		#'
8		#' @return This returns the value of the OS hazard function at time t.
9		#' @export
10		#'
11		#' @examples
12		#' ExpHazOS(c(1:5), 0.2, 1.1, 0.8)
13		ExpHazOS <- function(t, h01, h02, h12) {
14	17x	assert_numeric(t, lower = 0, any.missing = FALSE)
15	17x	assert_positive_number(h01)
16	17x	assert_positive_number(h02)
17	17x	assert_positive_number(h12)
18
19	17x	h012 <- h12 - h01 - h02
20	17x	((h12 - h02) * (h01 + h02) - h01 * h12 * exp(-h012 * t)) / ((h12 - h02) - h01 * exp(-h012 * t))
21		}
22
23		#' Helper Function for `avgHRExpOS()`
24		#'
25		#' It is an integrand of the form OS hazard function with intensities h01, h02, h12
26		#' at time point t multiplied with a weighted product of the two OS Survival functions
27		#' at t (one for intensities h0 and one for h1).
28		#'
29		#' @param x (`numeric`)\cr variable of integration.
30		#' @param h01 (positive `number`)\cr transition hazard for 0 to 1 transition.
31		#' @param h02 (positive `number`)\cr transition hazard for 0 to 2 transition.
32		#' @param h12 (positive `number`)\cr transition hazard for 1 to 2 transition.
33		#' @param h0 (`list`)\cr transition parameters for the first treatment group.
34		#' @param h1 (`list`)\cr transition parameters for the second treatment group.
35		#' @param alpha (`number`)\cr weight parameter, see [avgHRExpOS()].
36		#' @return This returns the value of the integrand used to calculate
37		#' the average hazard ratio for constant transition hazards, see [avgHRExpOS()].
38		#' @export
39		#'
40		#' @examples
41		#' h0 <- list(h01 = 0.18, h02 = 0.06, h12 = 0.17)
42		#' h1 <- list(h01 = 0.23, h02 = 0.07, h12 = 0.19)
43		#' avgHRIntegExpOS(x = 5, h01 = 0.2, h02 = 0.5, h12 = 0.7, h0 = h0, h1 = h1, alpha = 0.5)
44		avgHRIntegExpOS <- function(x, h01, h02, h12, h0, h1, alpha) {
45	14x	assert_positive_number(alpha)
46
47	14x	weightedSurv <- (ExpSurvOS(t = x, h01 = h0$h01, h02 = h0$h02, h12 = h0$h12) *
48	14x	ExpSurvOS(t = x, h01 = h1$h01, h02 = h1$h02, h12 = h1$h12))^alpha
49	14x	ExpHazOS(x, h01, h02, h12) * weightedSurv
50		}
51
52		#' Average OS Hazard Ratio from Constant Transition Hazards
53		#'
54		#' @param transitionByArm (`list`)\cr transition parameters for each treatment group.
55		#' See [exponential_transition()], [piecewise_exponential()] and [weibull_transition()] for details.
56		#' @param alpha (`number`)\cr assigns weights to time points, where values higher than 0.5 assign greater weight
57		#' to earlier times and values lower than 0.5 assign greater weight to later times.
58		#' @param upper (`number`)\cr upper (time) limit of integration.
59		#'
60		#' @return This returns the value of the average hazard ratio.
61		#' @export
62		#'
63		#' @examples
64		#' transitionTrt <- exponential_transition(h01 = 0.18, h02 = 0.06, h12 = 0.17)
65		#' transitionCtl <- exponential_transition(h01 = 0.23, h02 = 0.07, h12 = 0.19)
66		#' transitionList <- list(transitionCtl, transitionTrt)
67		#' avgHRExpOS(transitionByArm = transitionList, alpha = 0.5, upper = 100)
68		avgHRExpOS <- function(transitionByArm, alpha = 0.5, upper = Inf) {
69	2x	assert_list(transitionByArm, len = 2)
70	2x	assert_class(transitionByArm[[1]], "TransitionParameters")
71	2x	assert_class(transitionByArm[[2]], "TransitionParameters")
72	2x	assert_positive_number(alpha)
73	2x	assert_positive_number(upper)
74
75	2x	h0 <- transitionByArm[[1]]$hazard
76	2x	h1 <- transitionByArm[[2]]$hazard
77
78	2x	num <- stats::integrate(avgHRIntegExpOS,
79	2x	lower = 0, upper = upper, h01 = h1$h01, h02 = h1$h02, h12 = h1$h12,
80	2x	h0 = h0, h1 = h1, alpha = alpha
81	2x	)$value
82	2x	denom <- stats::integrate(avgHRIntegExpOS,
83	2x	lower = 0, upper = upper, h01 = h0$h01, h02 = h0$h02, h12 = h0$h12,
84	2x	h0 = h0, h1 = h1, alpha = alpha
85	2x	)$value
86	2x	num / denom
87		}

1		#' Piecewise Constant Hazard Values
2		#'
3		#' This returns piecewise constant hazard values at specified time points.
4		#'
5		#' @param haz (`numeric`)\cr piecewise constant input hazard.
6		#' @param pw (`numeric`)\cr time intervals for the piecewise constant hazard.
7		#' @param x (`numeric`)\cr time-points.
8		#'
9		#' @return Hazard values at input time-points.
10		#' @export
11		#'
12		#' @examples
13		#' getPWCHazard(c(1, 1.2, 1.4), c(0, 2, 3), c(1, 4, 6))
14		getPWCHazard <- function(haz, pw, x) { # nolint
15	38866x	sapply(x, function(jj) {
16	1283790x	y <- NULL
17		# Find interval and corresponding hazard value for time x[jj].
18	1283790x	for (ii in seq_along(haz)) {
19	3851341x	if (jj >= pw[ii]) {
20	1475722x	y <- haz[ii]
21		}
22		}
23	1283790x	y
24		})
25		}
26
27		#' Sum of Two Piecewise Constant Hazards
28		#'
29		#' This returns the sum of two piecewise constant hazards per interval.
30		#'
31		#' @param haz1 (`numeric`)\cr first summand (piecewise constant hazard).
32		#' @param haz2 (`numeric`)\cr second summand (piecewise constant hazard).
33		#' @param pw1 (`numeric`)\cr time intervals of first summand.
34		#' @param pw2 (`numeric`)\cr time intervals of second summand.
35		#'
36		#' @return List with elements `hazards` and `intervals` for the sum of two piecewise constant hazards.
37		#' @export
38		#'
39		#' @examples
40		#' getSumPCW(c(1.2, 0.3, 0.6), c(1.2, 0.7, 1), c(0, 8, 9), c(0, 1, 4))
41		getSumPCW <- function(haz1, haz2, pw1, pw2) {
42		# Get all cutpoints for the intervals.
43	2008x	cuts_sum <- unique(sort(c(pw1, pw2)))
44	2008x	haz_sum <- NULL
45		# Get sum of hazards for all intervals.
46	2008x	for (i in seq_along(cuts_sum)) {
47	10026x	haz_sum[i] <- getPWCHazard(haz1, pw1, cuts_sum[i]) +
48	10026x	getPWCHazard(haz2, pw2, cuts_sum[i])
49		}
50	2008x	list(
51	2008x	hazards = haz_sum,
52	2008x	intervals = cuts_sum
53		)
54		}

1		#' Assertion for Positive Number
2		#'
3		#' @param x what to check.
4		#' @param zero_ok (`flag`)\cr whether `x` can be zero or not.
5		#'
6		#' @return Raises an error if `x` is not a single positive (or non-negative) number.
7		#' @export
8		#'
9		#' @examples
10		#' assert_positive_number(3.2)
11		#' assert_positive_number(0, zero_ok = TRUE)
12		assert_positive_number <- function(x,
13		zero_ok = FALSE) {
14	54926x	assert_number(x)
15	54925x	assert_flag(zero_ok)
16	54925x	if (zero_ok) {
17	28569x	assert_true(x >= 0)
18		} else {
19	26356x	assert_true(x > 0)
20		}
21		}
22
23		#' Assertion for vector describing intervals
24		#'
25		#' We define an intervals vector to always start with 0, and contain
26		#' unique ordered time points.
27		#'
28		#' @param x what to check.
29		#' @param y (`count`)\cr required length of `y`.
30		#'
31		#' @return Raises an error if `x` is not an intervals vector starting with 0.
32		#' @export
33		#'
34		#' @examples
35		#' assert_intervals(c(0, 5, 7), 3)
36		assert_intervals <- function(x, y) {
37	573612x	assert_numeric(
38	573612x	x,
39	573612x	lower = 0,
40	573612x	any.missing = FALSE,
41	573612x	all.missing = FALSE,
42	573612x	len = y,
43	573612x	unique = TRUE,
44	573612x	sorted = TRUE
45		)
46	573607x	assert_true(x[1] == 0)
47		}